Simulation of decay photons through the D1S method (#3235)

D1S FTW!

CMakeLists.txt

@@ -342,6 +342,7 @@ list(APPEND libopenmc_SOURCES
   src/boundary_condition.cpp
   src/bremsstrahlung.cpp
   src/cell.cpp
+  src/chain.cpp
   src/cmfd_solver.cpp
   src/cross_sections.cpp
   src/dagmc.cpp
@@ -424,6 +425,7 @@ list(APPEND libopenmc_SOURCES
   src/tallies/filter_meshsurface.cpp
   src/tallies/filter_mu.cpp
   src/tallies/filter_musurface.cpp
+  src/tallies/filter_parent_nuclide.cpp
   src/tallies/filter_particle.cpp
   src/tallies/filter_polar.cpp
   src/tallies/filter_sph_harm.cpp

docs/source/conf.py

@@ -52,7 +52,7 @@ if not on_rtd:
 templates_path = ['_templates']
 
 # The suffix of source filenames.
-source_suffix = '.rst'
+source_suffix = {'.rst': 'restructuredtext'}
 
 # The encoding of source files.
 #source_encoding = 'utf-8'

docs/source/pythonapi/base.rst

@@ -145,6 +145,7 @@ Constructing Tallies
    openmc.TimeFilter
    openmc.ZernikeFilter
    openmc.ZernikeRadialFilter
+   openmc.ParentNuclideFilter
    openmc.ParticleFilter
    openmc.RegularMesh
    openmc.RectilinearMesh

docs/source/pythonapi/capi.rst

@@ -79,6 +79,7 @@ Classes
    MeshSurfaceFilter
    MuFilter
    Nuclide
+   ParentNuclideFilter
    ParticleFilter
    PolarFilter
    RectilinearMesh

docs/source/pythonapi/deplete.rst

@@ -287,3 +287,15 @@ the following abstract base classes:
    abc.Integrator
    abc.SIIntegrator
    abc.DepSystemSolver
+
+D1S Functions
+-------------
+
+.. autosummary::
+   :toctree: generated
+   :nosignatures:
+   :template: myfunction.rst
+
+   d1s.prepare_tallies
+   d1s.time_correction_factors
+   d1s.apply_time_correction

docs/source/usersguide/decay_sources.rst

@@ -0,0 +1,92 @@
+.. usersguide_decay_sources:
+
+=============
+Decay Sources
+=============
+
+Through the :ref:`depletion <usersguide_depletion>` capabilities in OpenMC, it
+is possible to simulate radiation emitted from the decay of activated materials.
+For fusion energy systems, this is commonly done using what is known as the
+`rigorous 2-step <https://doi.org/10.1016/S0920-3796(02)00144-8>`_ (R2S) method.
+In this method, a neutron transport calculation is used to determine the neutron
+flux and reaction rates over a cell- or mesh-based spatial discretization of the
+model. Then, the neutron flux in each discrete region is used to predict the
+activated material composition using a depletion solver. Finally, a photon
+transport calculation with a source based on the activity and energy spectrum of
+the activated materials is used to determine a desired physical response (e.g.,
+a dose rate) at one or more locations of interest.
+
+Once a depletion simulation has been completed in OpenMC, the intrinsic decay
+source can be determined as follows. First the activated material composition
+can be determined using the :class:`openmc.deplete.Results` object. Indexing an
+instance of this class with the timestep index returns a
+:class:`~openmc.deplete.StepResult` object, which itself has a
+:meth:`~openmc.deplete.StepResult.get_material` method. Once the activated
+:class:`~openmc.Material` has been obtained, the
+:meth:`~openmc.Material.get_decay_photon_energy` method will give the energy
+spectrum of the decay photon source. The integral of the spectrum also indicates
+the intensity of the source in units of [Bq]. Altogether, the workflow looks as
+follows::
+
+    results = openmc.deplete.Results("depletion_results.h5")
+
+    # Get results at last timestep
+    step = results[-1]
+
+    # Get activated material composition for ID=1
+    activated_mat = step.get_material('1')
+
+    # Determine photon source
+    photon_energy = activated_mat.get_decay_photon_energy()
+
+By default, the :meth:`~openmc.Material.get_decay_photon_energy` method will
+eliminate spectral lines with very low intensity, but this behavior can be
+configured with the ``clip_tolerance`` argument.
+
+Direct 1-Step (D1S) Calculations
+================================
+
+OpenMC also includes built-in capability for performing shutdown dose rate
+calculations using the `direct 1-step <https://10.1016/S0920-3796(01)00188-0>`_
+(D1S) method. In this method, a single coupled neutron--photon transport
+calculation is used where the prompt photon production is replaced with photons
+produced from the decay of radionuclides in an activated material. To obtain
+properly scaled results, it is also necessary to apply time correction factors.
+A normal neutron transport calculation can be extended to a D1S calculation with
+a few helper functions. First, import the ``d1s`` submodule, which is part of
+:mod:`openmc.deplete`::
+
+    from openmc.deplete import d1s
+
+First, you need to instruct OpenMC to use decay photon data instead of prompt
+photon data. This is done with an attribute on the :class:`~openmc.Settings`
+class::
+
+    model = openmc.Model()
+    ...
+    model.settings.use_decay_photons = True
+
+To prepare any tallies for use of the D1S method, you should call the
+:func:`~openmc.deplete.d1s.prepare_tallies` function, which adds a
+:class:`openmc.ParentNuclideFilter` (used later for assigning time correction
+factors) to any applicable tally and returns a list of possible radionuclides
+based on the :ref:`chain file <usersguide_data>`. Once the tallies are prepared,
+the model can be simulated::
+
+    output_path = model.run()
+
+Finally, the time correction factors need to be computed and applied to the
+relevant tallies. This can be done with the aid of the
+:func:`~openmc.deplete.d1s.time_correction_factors` and
+:func:`~openmc.deplete.d1s.apply_time_correction` functions::
+
+    # Compute time correction factors based on irradiation schedule
+    factors = d1s.time_correction_factors(nuclides, timesteps, source_rates)
+
+    # Get tally from statepoint
+    with openmc.StatePoint(output_path) as sp:
+        dose_tally = sp.get_tally(name='dose tally')
+
+    # Apply time correction factors
+    tally = d1s.apply_time_correction(tally, factors, time_index)
+

docs/source/usersguide/index.rst

@@ -21,6 +21,7 @@ essential aspects of using OpenMC to perform simulations.
     tallies
     plots
     depletion
+    decay_sources
     scripts
     processing
     parallel
@@ -28,4 +29,3 @@ essential aspects of using OpenMC to perform simulations.
     variance_reduction
     random_ray
     troubleshoot
-    

include/openmc/chain.h

@@ -0,0 +1,97 @@
+//! \file chain.h
+//! \brief Depletion chain and associated information
+
+#ifndef OPENMC_CHAIN_H
+#define OPENMC_CHAIN_H
+
+#include <cmath>
+#include <string>
+#include <unordered_map>
+
+#include "pugixml.hpp"
+
+#include "openmc/angle_energy.h" // for AngleEnergy
+#include "openmc/distribution.h" // for UPtrDist
+#include "openmc/memory.h"       // for unique_ptr
+#include "openmc/vector.h"
+
+namespace openmc {
+
+//==============================================================================
+// Data for a nuclide in the depletion chain
+//==============================================================================
+
+class ChainNuclide {
+public:
+  // Types
+  struct Product {
+    std::string name;       //!< Reaction product name
+    double branching_ratio; //!< Branching ratio
+  };
+
+  // Constructors, destructors
+  ChainNuclide(pugi::xml_node node);
+  ~ChainNuclide();
+
+  //! Compute the decay constant for the nuclide
+  //! \return Decay constant in [1/s]
+  double decay_constant() const { return std::log(2.0) / half_life_; }
+
+  const Distribution* photon_energy() const { return photon_energy_.get(); }
+  const std::unordered_map<int, vector<Product>>& reaction_products() const
+  {
+    return reaction_products_;
+  }
+
+private:
+  // Data members
+  std::string name_;          //!< Name of nuclide
+  double half_life_ {0.0};    //!< Half-life in [s]
+  double decay_energy_ {0.0}; //!< Decay energy in [eV]
+  std::unordered_map<int, vector<Product>>
+    reaction_products_;    //!< Map of MT to reaction products
+  UPtrDist photon_energy_; //!< Decay photon energy distribution
+};
+
+//==============================================================================
+// Angle-energy distribution for decay photon
+//==============================================================================
+
+class DecayPhotonAngleEnergy : public AngleEnergy {
+public:
+  explicit DecayPhotonAngleEnergy(const Distribution* dist)
+    : photon_energy_(dist)
+  {}
+
+  //! Sample distribution for an angle and energy
+  //! \param[in] E_in Incoming energy in [eV]
+  //! \param[out] E_out Outgoing energy in [eV]
+  //! \param[out] mu Outgoing cosine with respect to current direction
+  //! \param[inout] seed Pseudorandom seed pointer
+  void sample(
+    double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+
+private:
+  const Distribution* photon_energy_;
+};
+
+//==============================================================================
+// Global variables
+//==============================================================================
+
+namespace data {
+
+extern std::unordered_map<std::string, int> chain_nuclide_map;
+extern vector<unique_ptr<ChainNuclide>> chain_nuclides;
+
+} // namespace data
+
+//==============================================================================
+// Non-member functions
+//==============================================================================
+
+void read_chain_file_xml();
+
+} // namespace openmc
+
+#endif // OPENMC_CHAIN_H

include/openmc/endf.h

@@ -53,6 +53,10 @@ public:
   //! \param[in] dset Dataset containing coefficients
   explicit Polynomial(hid_t dset);
 
+  //! Construct polynomial from coefficients
+  //! \param[in] coef Polynomial coefficients
+  explicit Polynomial(vector<double> coef) : coef_(coef) {}
+
   //! Evaluate the polynomials
   //! \param[in] x independent variable
   //! \return Polynomial evaluated at x

include/openmc/particle.h

@@ -47,7 +47,8 @@ public:
   //! \param u Direction of the secondary particle
   //! \param E Energy of the secondary particle in [eV]
   //! \param type Particle type
-  void create_secondary(double wgt, Direction u, double E, ParticleType type);
+  //! \return Whether a secondary particle was created
+  bool create_secondary(double wgt, Direction u, double E, ParticleType type);
 
   //! split a particle
   //

include/openmc/particle_data.h

@@ -49,6 +49,9 @@ struct SourceSite {
   int delayed_group {0};
   int surf_id {0};
   ParticleType particle;
+
+  // Extra attributes that don't show up in source written to file
+  int parent_nuclide {-1};
   int64_t parent_id;
   int64_t progeny_id;
 };
@@ -441,6 +444,7 @@ private:
   int event_nuclide_;
   int event_mt_;
   int delayed_group_ {0};
+  int parent_nuclide_ {-1};
 
   int n_bank_ {0};
   double bank_second_E_ {0.0};
@@ -555,6 +559,8 @@ public:
   const int& event_nuclide() const { return event_nuclide_; }
   int& event_mt() { return event_mt_; }           // MT number of collision
   int& delayed_group() { return delayed_group_; } // delayed group
+  const int& parent_nuclide() const { return parent_nuclide_; }
+  int& parent_nuclide() { return parent_nuclide_; } // Parent nuclide
 
   // Post-collision data
   double& bank_second_E()

include/openmc/reaction.h

@@ -26,7 +26,9 @@ public:
   //! Construct reaction from HDF5 data
   //! \param[in] group HDF5 group containing reaction data
   //! \param[in] temperatures Desired temperatures for cross sections
-  explicit Reaction(hid_t group, const vector<int>& temperatures);
+  //! \param[in] name Name of the nuclide
+  explicit Reaction(
+    hid_t group, const vector<int>& temperatures, std::string name);
 
   //! Calculate cross section given temperautre/grid index, interpolation factor
   //

include/openmc/reaction_product.h

@@ -7,6 +7,7 @@
 #include "hdf5.h"
 
 #include "openmc/angle_energy.h"
+#include "openmc/chain.h"
 #include "openmc/endf.h"
 #include "openmc/memory.h" // for unique_ptr
 #include "openmc/particle.h"
@@ -37,6 +38,10 @@ public:
   //! \param[in] group HDF5 group containing data
   explicit ReactionProduct(hid_t group);
 
+  //! Construct reaction product for decay photon from chain nuclide product
+  //! \param[in] product Chain nuclide product
+  explicit ReactionProduct(const ChainNuclide::Product& product);
+
   //! Sample an outgoing angle and energy
   //! \param[in] E_in Incoming energy in [eV]
   //! \param[out] E_out Outgoing energy in [eV]
@@ -50,6 +55,7 @@ public:
   unique_ptr<Function1D> yield_;      //!< Yield as a function of energy
   vector<Tabulated1D> applicability_; //!< Applicability of distribution
   vector<Secondary> distribution_;    //!< Secondary angle-energy distribution
+  int parent_nuclide_ = -1;           //!< Index of chain nuclide that is parent
 };
 
 } // namespace openmc

include/openmc/settings.h

@@ -67,6 +67,7 @@ extern bool trigger_predict;         //!< predict batches for triggers?
 extern bool uniform_source_sampling; //!< sample sources uniformly?
 extern bool ufs_on;                  //!< uniform fission site method on?
 extern bool urr_ptables_on;          //!< use unresolved resonance prob. tables?
+extern bool use_decay_photons;       //!< use decay photons for D1S
 extern "C" bool weight_windows_on;   //!< are weight windows are enabled?
 extern bool weight_window_checkpoint_surface;   //!< enable weight window check
                                                 //!< upon surface crossing?

include/openmc/tallies/filter.h

@@ -36,6 +36,7 @@ enum class FilterType {
   MESH_SURFACE,
   MU,
   MUSURFACE,
+  PARENT_NUCLIDE,
   PARTICLE,
   POLAR,
   SPHERICAL_HARMONICS,

include/openmc/tallies/filter_parent_nuclide.h

@@ -0,0 +1,56 @@
+#ifndef OPENMC_TALLIES_FILTER_PARENT_NUCLIDE_H
+#define OPENMC_TALLIES_FILTER_PARENT_NUCLIDE_H
+
+#include <string>
+#include <unordered_map>
+
+#include "openmc/span.h"
+#include "openmc/tallies/filter.h"
+#include "openmc/vector.h"
+
+namespace openmc {
+
+//==============================================================================
+//! Bins events by parent nuclide (for decay photons)
+//==============================================================================
+
+class ParentNuclideFilter : public Filter {
+public:
+  //----------------------------------------------------------------------------
+  // Constructors, destructors
+
+  ~ParentNuclideFilter() = default;
+
+  //----------------------------------------------------------------------------
+  // Methods
+
+  std::string type_str() const override { return "parentnuclide"; }
+  FilterType type() const override { return FilterType::PARENT_NUCLIDE; }
+
+  void from_xml(pugi::xml_node node) override;
+
+  void get_all_bins(const Particle& p, TallyEstimator estimator,
+    FilterMatch& match) const override;
+
+  void to_statepoint(hid_t filter_group) const override;
+
+  std::string text_label(int bin) const override;
+
+  //----------------------------------------------------------------------------
+  // Accessors
+
+  const vector<int>& bins() const { return bins_; }
+  void set_bins(span<const int> bins);
+
+protected:
+  //----------------------------------------------------------------------------
+  // Data members
+
+  vector<int> bins_;
+  vector<std::string> nuclides_;
+
+  std::unordered_map<int, int> map_;
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLIES_FILTER_PARENT_NUCLIDE_H

openmc/deplete/abc.py

@@ -41,6 +41,63 @@ _SECONDS_PER_HOUR = 60*60
 _SECONDS_PER_DAY = 24*60*60
 _SECONDS_PER_JULIAN_YEAR = 365.25*24*60*60
 
+
+def _normalize_timesteps(
+        timesteps: Sequence[float] | Sequence[tuple[float, str]],
+        source_rates: float | Sequence[float],
+        timestep_units: str = 's',
+        operator: TransportOperator | None = None,
+):
+    if not isinstance(source_rates, Sequence):
+        # Ensure that rate is single value if that is the case
+        source_rates = [source_rates] * len(timesteps)
+
+    if len(source_rates) != len(timesteps):
+        raise ValueError(
+            "Number of time steps ({}) != number of powers ({})".format(
+                len(timesteps), len(source_rates)))
+
+    # Get list of times / units
+    if isinstance(timesteps[0], Sequence):
+        times, units = zip(*timesteps)
+    else:
+        times = timesteps
+        units = [timestep_units] * len(timesteps)
+
+    # Determine number of seconds for each timestep
+    seconds = []
+    for timestep, unit, rate in zip(times, units, source_rates):
+        # Make sure values passed make sense
+        check_type('timestep', timestep, Real)
+        check_greater_than('timestep', timestep, 0.0, False)
+        check_type('timestep units', unit, str)
+        check_type('source rate', rate, Real)
+        check_greater_than('source rate', rate, 0.0, True)
+
+        if unit in ('s', 'sec'):
+            seconds.append(timestep)
+        elif unit in ('min', 'minute'):
+            seconds.append(timestep*_SECONDS_PER_MINUTE)
+        elif unit in ('h', 'hr', 'hour'):
+            seconds.append(timestep*_SECONDS_PER_HOUR)
+        elif unit in ('d', 'day'):
+            seconds.append(timestep*_SECONDS_PER_DAY)
+        elif unit in ('a', 'year'):
+            seconds.append(timestep*_SECONDS_PER_JULIAN_YEAR)
+        elif unit.lower() == 'mwd/kg':
+            watt_days_per_kg = 1e6*timestep
+            kilograms = 1e-3*operator.heavy_metal
+            if rate == 0.0:
+                raise ValueError("Cannot specify a timestep in [MWd/kg] when"
+                                 " the power is zero.")
+            days = watt_days_per_kg * kilograms / rate
+            seconds.append(days*_SECONDS_PER_DAY)
+        else:
+            raise ValueError(f"Invalid timestep unit '{unit}'")
+
+    return (np.asarray(seconds), np.asarray(source_rates))
+
+
 OperatorResult = namedtuple('OperatorResult', ['k', 'rates'])
 OperatorResult.__doc__ = """\
 Result of applying transport operator
@@ -551,10 +608,10 @@ class Integrator(ABC):
     def __init__(
             self,
             operator: TransportOperator,
-            timesteps: Sequence[float],
+            timesteps: Sequence[float] | Sequence[tuple[float, str]],
             power: Optional[Union[float, Sequence[float]]] = None,
             power_density: Optional[Union[float, Sequence[float]]] = None,
-            source_rates: Optional[Sequence[float]] = None,
+            source_rates: Optional[Union[float, Sequence[float]]] = None,
             timestep_units: str = 's',
             solver: str = "cram48"
         ):
@@ -582,53 +639,9 @@ class Integrator(ABC):
         elif source_rates is None:
             raise ValueError("Either power, power_density, or source_rates must be set")
 
-        if not isinstance(source_rates, Iterable):
-            # Ensure that rate is single value if that is the case
-            source_rates = [source_rates] * len(timesteps)
-
-        if len(source_rates) != len(timesteps):
-            raise ValueError(
-                "Number of time steps ({}) != number of powers ({})".format(
-                    len(timesteps), len(source_rates)))
-
-        # Get list of times / units
-        if isinstance(timesteps[0], Iterable):
-            times, units = zip(*timesteps)
-        else:
-            times = timesteps
-            units = [timestep_units] * len(timesteps)
-
-        # Determine number of seconds for each timestep
-        seconds = []
-        for timestep, unit, rate in zip(times, units, source_rates):
-            # Make sure values passed make sense
-            check_type('timestep', timestep, Real)
-            check_greater_than('timestep', timestep, 0.0, False)
-            check_type('timestep units', unit, str)
-            check_type('source rate', rate, Real)
-            check_greater_than('source rate', rate, 0.0, True)
-
-            if unit in ('s', 'sec'):
-                seconds.append(timestep)
-            elif unit in ('min', 'minute'):
-                seconds.append(timestep*_SECONDS_PER_MINUTE)
-            elif unit in ('h', 'hr', 'hour'):
-                seconds.append(timestep*_SECONDS_PER_HOUR)
-            elif unit in ('d', 'day'):
-                seconds.append(timestep*_SECONDS_PER_DAY)
-            elif unit in ('a', 'year'):
-                seconds.append(timestep*_SECONDS_PER_JULIAN_YEAR)
-            elif unit.lower() == 'mwd/kg':
-                watt_days_per_kg = 1e6*timestep
-                kilograms = 1e-3*operator.heavy_metal
-                if rate == 0.0:
-                    raise ValueError("Cannot specify a timestep in [MWd/kg] when"
-                                     " the power is zero.")
-                days = watt_days_per_kg * kilograms / rate
-                seconds.append(days*_SECONDS_PER_DAY)
-            else:
-                raise ValueError(f"Invalid timestep unit '{unit}'")
-
+        # Normalize timesteps and source rates
+        seconds, source_rates = _normalize_timesteps(
+            timesteps, source_rates, timestep_units, operator)
         self.timesteps = np.asarray(seconds)
         self.source_rates = np.asarray(source_rates)
 

openmc/deplete/d1s.py

@@ -0,0 +1,246 @@
+"""D1S module
+
+This module contains functionality to support the direct 1-step (D1S) method for
+shutdown dose rate calculations.
+
+"""
+
+from copy import deepcopy
+from typing import Sequence
+from math import log, prod
+
+import numpy as np
+
+import openmc
+from openmc.data import half_life
+from .abc import _normalize_timesteps
+from .chain import Chain
+
+
+def get_radionuclides(model: openmc.Model, chain_file: str | None = None) -> list[str]:
+    """Determine all radionuclides that can be produced during D1S.
+
+    Parameters
+    ----------
+    model : openmc.Model
+        Model that should be used for determining what nuclides are present
+    chain_file : str, optional
+        Which chain file to use for inspecting decay data. If None is passed,
+        defaults to ``openmc.config['chain_file']``
+
+    Returns
+    -------
+    List of nuclide names
+
+    """
+    # Determine what nuclides appear in model
+    model_nuclides = set(model.geometry.get_all_nuclides())
+
+    # Load chain file
+    if chain_file is None:
+        chain_file = openmc.config['chain_file']
+    chain = Chain.from_xml(chain_file)
+
+    radionuclides = set()
+    for nuclide in chain.nuclides:
+        # Restrict to set of nuclides present in model
+        if nuclide.name not in model_nuclides:
+            continue
+
+        # Loop over reactions and add any targets that are unstable
+        for rx_tuple in nuclide.reactions:
+            target = rx_tuple.target
+            if target is None:
+                continue
+            target_nuclide = chain[target]
+            if target_nuclide.half_life is not None:
+                radionuclides.add(target_nuclide.name)
+
+    return list(radionuclides)
+
+
+def time_correction_factors(
+        nuclides: list[str],
+        timesteps: Sequence[float] | Sequence[tuple[float, str]],
+        source_rates: float | Sequence[float],
+        timestep_units: str = 's'
+) -> dict[str, np.ndarray]:
+    """Calculate time correction factors for the D1S method.
+
+    This function determines the time correction factor that should be applied
+    to photon tallies as part of the D1S method.
+
+    Parameters
+    ----------
+    nuclides : list of str
+        The name of the nuclide to find the time correction for, e.g., 'Ni65'
+    timesteps : iterable of float or iterable of tuple
+        Array of timesteps. Note that values are not cumulative. The units are
+        specified by the `timestep_units` argument when `timesteps` is an
+        iterable of float. Alternatively, units can be specified for each step
+        by passing a sequence of (value, unit) tuples.
+    source_rates : float or iterable of float
+        Source rate in [neutron/sec] for each interval in `timesteps`
+    timestep_units : {'s', 'min', 'h', 'd', 'a'}, optional
+        Units for values specified in the `timesteps` argument. 's' means
+        seconds, 'min' means minutes, 'h' means hours, and 'a' means Julian
+        years.
+
+    Returns
+    -------
+    dict
+        Dictionary mapping nuclide to an array of time correction factors for
+        each time.
+
+    """
+
+    # Determine normalized timesteps and source rates
+    timesteps, source_rates = _normalize_timesteps(
+        timesteps, source_rates, timestep_units)
+
+    # Calculate decay rate for each nuclide
+    decay_rate = np.array([log(2.0) / half_life(x) for x in nuclides])
+
+    n_timesteps = len(timesteps) + 1
+    n_nuclides = len(nuclides)
+
+    # Create a 2D array for the time correction factors
+    h = np.zeros((n_timesteps, n_nuclides))
+
+    for i, (dt, rate) in enumerate(zip(timesteps, source_rates)):
+        # Precompute the exponential terms. Since (1 - exp(-x)) is susceptible to
+        # roundoff error, use expm1 instead (which computes exp(x) - 1)
+        g = np.exp(-decay_rate*dt)
+        one_minus_g = -np.expm1(-decay_rate*dt)
+
+        # Eq. (4) in doi:10.1016/j.fusengdes.2019.111399
+        h[i + 1] = rate*one_minus_g + h[i]*g
+
+    return {nuclides[i]: h[:, i] for i in range(n_nuclides)}
+
+
+def apply_time_correction(
+        tally: openmc.Tally,
+        time_correction_factors: dict[str, np.ndarray],
+        index: int = -1,
+        sum_nuclides: bool = True
+) -> openmc.Tally:
+    """Apply time correction factors to a tally.
+
+    This function applies the time correction factors at the given index to a
+    tally that contains a :class:`~openmc.ParentNuclideFilter`. When
+    `sum_nuclides` is True, values over all parent nuclides will be summed,
+    leaving a single value for each filter combination.
+
+    Parameters
+    ----------
+    tally : openmc.Tally
+        Tally to apply the time correction factors to
+    time_correction_factors : dict
+        Time correction factors as returned by :func:`time_correction_factors`
+    index : int, optional
+        Index to use for the correction factors
+    sum_nuclides : bool
+        Whether to sum over the parent nuclides
+
+    Returns
+    -------
+    openmc.Tally
+        Derived tally with time correction factors applied
+
+    """
+    # Make sure the tally contains a ParentNuclideFilter
+    for i_filter, filter in enumerate(tally.filters):
+        if isinstance(filter, openmc.ParentNuclideFilter):
+            break
+    else:
+        raise ValueError('Tally must contain a ParentNuclideFilter')
+
+    # Get list of radionuclides based on tally filter
+    radionuclides = [str(x) for x in tally.filters[i_filter].bins]
+    tcf = np.array([time_correction_factors[x][index] for x in radionuclides])
+
+    # Create copy of tally
+    new_tally = deepcopy(tally)
+
+    # Determine number of bins in other filters
+    n_bins_before = prod([f.num_bins for f in tally.filters[:i_filter]])
+    n_bins_after = prod([f.num_bins for f in tally.filters[i_filter + 1:]])
+
+    # Reshape sum and sum_sq, apply TCF, and sum along that axis
+    _, n_nuclides, n_scores = new_tally.shape
+    n_radionuclides = len(radionuclides)
+    shape = (n_bins_before, n_radionuclides, n_bins_after, n_nuclides, n_scores)
+    tally_sum = new_tally.sum.reshape(shape)
+    tally_sum_sq = new_tally.sum_sq.reshape(shape)
+
+    # Apply TCF, broadcasting to the correct dimensions
+    tcf.shape = (1, -1, 1, 1, 1)
+    new_tally._sum = tally_sum * tcf
+    new_tally._sum_sq = tally_sum_sq * (tcf*tcf)
+    new_tally._mean = None
+    new_tally._std_dev = None
+
+    shape = (-1, n_nuclides, n_scores)
+
+    if sum_nuclides:
+        # Query the mean and standard deviation
+        mean = new_tally.mean
+        std_dev = new_tally.std_dev
+
+        # Sum over parent nuclides (note that when combining different bins for
+        # parent nuclide, we can't work directly on sum_sq)
+        new_tally._mean = mean.sum(axis=1).reshape(shape)
+        new_tally._std_dev = np.linalg.norm(std_dev, axis=1).reshape(shape)
+        new_tally._derived = True
+
+        # Remove ParentNuclideFilter
+        new_tally.filters.pop(i_filter)
+    else:
+        new_tally._sum.shape = shape
+        new_tally._sum_sq.shape = shape
+
+    return new_tally
+
+
+def prepare_tallies(
+        model: openmc.Model,
+        nuclides: list[str] | None = None,
+        chain_file: str | None = None
+) -> list[str]:
+    """Prepare tallies for the D1S method.
+
+    This function adds a :class:`~openmc.ParentNuclideFilter` to any tally that
+    has a particle filter with a single 'photon' bin.
+
+    Parameters
+    ----------
+    model : openmc.Model
+        Model to prepare tallies for
+    nuclides : list of str, optional
+        Nuclides to use for the parent nuclide filter. If None, radionuclides
+        are determined from :func:`get_radionuclides`.
+    chain_file : str, optional
+        Chain file to use for inspecting decay data. If None, defaults to
+        ``openmc.config['chain_file']``
+
+    Returns
+    -------
+    list of str
+        List of parent nuclides being filtered on
+
+    """
+    if nuclides is None:
+        nuclides = get_radionuclides(model, chain_file)
+    filter = openmc.ParentNuclideFilter(nuclides)
+
+    # Apply parent nuclide filter to any tally that has a particle filter with a
+    # single 'photon' bin
+    for tally in model.tallies:
+        for f in tally.filters:
+            if isinstance(f, openmc.ParticleFilter):
+                if list(f.bins) == ['photon']:
+                    tally.filters.append(filter)
+                    break
+
+    return nuclides

openmc/filter.py

@@ -25,7 +25,7 @@ _FILTER_TYPES = (
     'energyout', 'mu', 'musurface', 'polar', 'azimuthal', 'distribcell', 'delayedgroup',
     'energyfunction', 'cellfrom', 'materialfrom', 'legendre', 'spatiallegendre',
     'sphericalharmonics', 'zernike', 'zernikeradial', 'particle', 'cellinstance',
-    'collision', 'time'
+    'collision', 'time', 'parentnuclide'
 )
 
 _CURRENT_NAMES = (
@@ -732,7 +732,7 @@ class SurfaceFilter(WithIDFilter):
 
 
 class ParticleFilter(Filter):
-    """Bins tally events based on the Particle type.
+    """Bins tally events based on the particle type.
 
     Parameters
     ----------
@@ -788,6 +788,33 @@ class ParticleFilter(Filter):
         return cls(bins, filter_id=filter_id)
 
 
+class ParentNuclideFilter(ParticleFilter):
+    """Bins tally events based on the parent nuclide
+
+    Parameters
+    ----------
+    bins : str, or iterable of str
+        Names of nuclides (e.g., 'Ni65')
+    filter_id : int
+        Unique identifier for the filter
+
+    Attributes
+    ----------
+    bins : iterable of str
+        Names of nuclides
+    id : int
+        Unique identifier for the filter
+    num_bins : Integral
+        The number of filter bins
+
+    """
+    @Filter.bins.setter
+    def bins(self, bins):
+        bins = np.atleast_1d(bins)
+        cv.check_iterable_type('filter bins', bins, str)
+        self._bins = bins
+
+
 class MeshFilter(Filter):
     """Bins tally event locations by mesh elements.
 

openmc/lib/core.py

@@ -25,6 +25,7 @@ class _SourceSite(Structure):
                 ('delayed_group', c_int),
                 ('surf_id', c_int),
                 ('particle', c_int),
+                ('parent_nuclide', c_int),
                 ('parent_id', c_int64),
                 ('progeny_id', c_int64)]
 

openmc/lib/filter.py

@@ -21,9 +21,10 @@ __all__ = [
     'CellInstanceFilter', 'CollisionFilter', 'DistribcellFilter', 'DelayedGroupFilter',
     'EnergyFilter', 'EnergyoutFilter', 'EnergyFunctionFilter', 'LegendreFilter',
     'MaterialFilter', 'MaterialFromFilter', 'MeshFilter', 'MeshBornFilter',
-    'MeshSurfaceFilter', 'MuFilter', 'MuSurfaceFilter', 'ParticleFilter', 'PolarFilter',
-    'SphericalHarmonicsFilter', 'SpatialLegendreFilter', 'SurfaceFilter', 'UniverseFilter',
-    'ZernikeFilter', 'ZernikeRadialFilter', 'filters'
+    'MeshSurfaceFilter', 'MuFilter', 'MuSurfaceFilter', 'ParentNuclideFilter',
+    'ParticleFilter', 'PolarFilter', 'SphericalHarmonicsFilter',
+    'SpatialLegendreFilter', 'SurfaceFilter', 'UniverseFilter', 'ZernikeFilter',
+    'ZernikeRadialFilter', 'filters'
 ]
 
 # Tally functions
@@ -544,6 +545,10 @@ class MuSurfaceFilter(Filter):
     filter_type = 'musurface'
 
 
+class ParentNuclideFilter(Filter):
+    filter_type = 'parentnuclide'
+
+
 class ParticleFilter(Filter):
     filter_type = 'particle'
 
@@ -647,6 +652,7 @@ _FILTER_TYPE_MAP = {
     'meshsurface': MeshSurfaceFilter,
     'mu': MuFilter,
     'musurface': MuSurfaceFilter,
+    'parentnuclide': ParentNuclideFilter,
     'particle': ParticleFilter,
     'polar': PolarFilter,
     'sphericalharmonics': SphericalHarmonicsFilter,

openmc/material.py

@@ -308,8 +308,9 @@ class Material(IDManagerMixin):
 
         Returns
         -------
-        Decay photon energy distribution. The integral of this distribution is
-        the total intensity of the photon source in the requested units.
+        Univariate or None
+            Decay photon energy distribution. The integral of this distribution
+            is the total intensity of the photon source in the requested units.
 
         """
         cv.check_value('units', units, {'Bq', 'Bq/g', 'Bq/cm3'})

openmc/settings.py

@@ -275,6 +275,8 @@ class Settings:
     ufs_mesh : openmc.RegularMesh
         Mesh to be used for redistributing source sites via the uniform fission
         site (UFS) method.
+    use_decay_photons : bool
+        Produce decay photons from neutron reactions instead of prompt
     verbosity : int
         Verbosity during simulation between 1 and 10. Verbosity levels are
         described in :ref:`verbosity`.
@@ -393,6 +395,7 @@ class Settings:
         self._weight_window_checkpoints = {}
         self._max_history_splits = None
         self._max_tracks = None
+        self._use_decay_photons = None
 
         self._random_ray = {}
 
@@ -1143,6 +1146,15 @@ class Settings:
 
         self._random_ray = random_ray
 
+    @property
+    def use_decay_photons(self) -> bool:
+        return self._use_decay_photons
+
+    @use_decay_photons.setter
+    def use_decay_photons(self, value):
+        cv.check_type('use decay photons', value, bool)
+        self._use_decay_photons = value
+
     def _create_run_mode_subelement(self, root):
         elem = ET.SubElement(root, "run_mode")
         elem.text = self._run_mode.value
@@ -1431,6 +1443,11 @@ class Settings:
             root.append(self.ufs_mesh.to_xml_element())
             if mesh_memo is not None: mesh_memo.add(self.ufs_mesh.id)
 
+    def _create_use_decay_photons_subelement(self, root):
+        if self._use_decay_photons is not None:
+            element = ET.SubElement(root, "use_decay_photons")
+            element.text = str(self._use_decay_photons).lower()
+
     def _create_resonance_scattering_subelement(self, root):
         res = self.resonance_scattering
         if res:
@@ -1957,6 +1974,11 @@ class Settings:
                 elif child.tag == 'sample_method':
                     self.random_ray['sample_method'] = child.text
 
+    def _use_decay_photons_from_xml_element(self, root):
+        text = get_text(root, 'use_decay_photons')
+        if text is not None:
+            self.use_decay_photons = text in ('true', '1')
+
     def to_xml_element(self, mesh_memo=None):
         """Create a 'settings' element to be written to an XML file.
 
@@ -2021,6 +2043,7 @@ class Settings:
         self._create_max_history_splits_subelement(element)
         self._create_max_tracks_subelement(element)
         self._create_random_ray_subelement(element)
+        self._create_use_decay_photons_subelement(element)
 
         # Clean the indentation in the file to be user-readable
         clean_indentation(element)
@@ -2126,8 +2149,8 @@ class Settings:
         settings._max_history_splits_from_xml_element(elem)
         settings._max_tracks_from_xml_element(elem)
         settings._random_ray_from_xml_element(elem)
+        settings._use_decay_photons_from_xml_element(elem)
 
-        # TODO: Get volume calculations
         return settings
 
     @classmethod

src/chain.cpp

@@ -0,0 +1,116 @@
+//! \file chain.cpp
+//! \brief Depletion chain and associated information
+
+#include "openmc/chain.h"
+
+#include <cstdlib> // for getenv
+#include <memory>  // for make_unique
+#include <string>  // for stod
+
+#include <fmt/core.h>
+#include <pugixml.hpp>
+
+#include "openmc/distribution.h" // for distribution_from_xml
+#include "openmc/error.h"
+#include "openmc/reaction.h"
+#include "openmc/xml_interface.h" // for get_node_value
+
+namespace openmc {
+
+//==============================================================================
+// ChainNuclide implementation
+//==============================================================================
+
+ChainNuclide::ChainNuclide(pugi::xml_node node)
+{
+  name_ = get_node_value(node, "name");
+  if (check_for_node(node, "half_life")) {
+    half_life_ = std::stod(get_node_value(node, "half_life"));
+  }
+  if (check_for_node(node, "decay_energy")) {
+    decay_energy_ = std::stod(get_node_value(node, "decay_energy"));
+  }
+
+  // Read reactions to store MT -> product map
+  for (pugi::xml_node reaction_node : node.children("reaction")) {
+    std::string rx_name = get_node_value(reaction_node, "type");
+    if (!reaction_node.attribute("target"))
+      continue;
+    std::string rx_target = get_node_value(reaction_node, "target");
+    double branching_ratio = 1.0;
+    if (reaction_node.attribute("branching_ratio")) {
+      branching_ratio =
+        std::stod(get_node_value(reaction_node, "branching_ratio"));
+    }
+    int mt = reaction_type(rx_name);
+    reaction_products_[mt].push_back({rx_target, branching_ratio});
+  }
+
+  for (pugi::xml_node source_node : node.children("source")) {
+    auto particle = get_node_value(source_node, "particle");
+    if (particle == "photon") {
+      photon_energy_ = distribution_from_xml(source_node);
+      break;
+    }
+  }
+
+  // Set entry in mapping
+  data::chain_nuclide_map[name_] = data::chain_nuclides.size();
+}
+
+ChainNuclide::~ChainNuclide()
+{
+  data::chain_nuclide_map.erase(name_);
+}
+
+//==============================================================================
+// DecayPhotonAngleEnergy implementation
+//==============================================================================
+
+void DecayPhotonAngleEnergy::sample(
+  double E_in, double& E_out, double& mu, uint64_t* seed) const
+{
+  E_out = photon_energy_->sample(seed);
+  mu = Uniform(-1., 1.).sample(seed);
+}
+
+//==============================================================================
+// Global variables
+//==============================================================================
+
+namespace data {
+
+std::unordered_map<std::string, int> chain_nuclide_map;
+vector<unique_ptr<ChainNuclide>> chain_nuclides;
+
+} // namespace data
+
+//==============================================================================
+// Non-member functions
+//==============================================================================
+
+void read_chain_file_xml()
+{
+  char* chain_file_path = std::getenv("OPENMC_CHAIN_FILE");
+  if (!chain_file_path) {
+    return;
+  }
+
+  write_message(5, "Reading chain file: {}...", chain_file_path);
+
+  pugi::xml_document doc;
+  auto result = doc.load_file(chain_file_path);
+  if (!result) {
+    fatal_error(
+      fmt::format("Error processing chain file: {}", chain_file_path));
+  }
+
+  // Get root element
+  pugi::xml_node root = doc.document_element();
+
+  for (auto node : root.children("nuclide")) {
+    data::chain_nuclides.push_back(std::make_unique<ChainNuclide>(node));
+  }
+}
+
+} // namespace openmc

src/initialize.cpp

@@ -12,6 +12,7 @@
 #include <fmt/core.h>
 
 #include "openmc/capi.h"
+#include "openmc/chain.h"
 #include "openmc/constants.h"
 #include "openmc/cross_sections.h"
 #include "openmc/error.h"
@@ -164,16 +165,17 @@ void initialize_mpi(MPI_Comm intracomm)
   MPI_Get_address(&b.delayed_group, &disp[5]);
   MPI_Get_address(&b.surf_id, &disp[6]);
   MPI_Get_address(&b.particle, &disp[7]);
-  MPI_Get_address(&b.parent_id, &disp[8]);
-  MPI_Get_address(&b.progeny_id, &disp[9]);
-  for (int i = 9; i >= 0; --i) {
+  MPI_Get_address(&b.parent_nuclide, &disp[8]);
+  MPI_Get_address(&b.parent_id, &disp[9]);
+  MPI_Get_address(&b.progeny_id, &disp[10]);
+  for (int i = 10; i >= 0; --i) {
     disp[i] -= disp[0];
   }
 
-  int blocks[] {3, 3, 1, 1, 1, 1, 1, 1, 1, 1};
+  int blocks[] {3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1};
   MPI_Datatype types[] {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE,
-    MPI_DOUBLE, MPI_INT, MPI_INT, MPI_INT, MPI_LONG, MPI_LONG};
-  MPI_Type_create_struct(10, blocks, disp, types, &mpi::source_site);
+    MPI_DOUBLE, MPI_INT, MPI_INT, MPI_INT, MPI_INT, MPI_LONG, MPI_LONG};
+  MPI_Type_create_struct(11, blocks, disp, types, &mpi::source_site);
   MPI_Type_commit(&mpi::source_site);
 }
 #endif // OPENMC_MPI
@@ -371,6 +373,9 @@ bool read_model_xml()
     }
   }
 
+  // Read data from chain file
+  read_chain_file_xml();
+
   // Read materials and cross sections
   if (!check_for_node(root, "materials")) {
     fatal_error(fmt::format(
@@ -423,6 +428,10 @@ void read_separate_xml_files()
   if (settings::run_mode != RunMode::PLOTTING) {
     read_cross_sections_xml();
   }
+
+  // Read data from chain file
+  read_chain_file_xml();
+
   read_materials_xml();
   read_geometry_xml();
 

src/nuclide.cpp

@@ -248,7 +248,8 @@ Nuclide::Nuclide(hid_t group, const vector<double>& temperature)
   for (auto name : group_names(rxs_group)) {
     if (starts_with(name, "reaction_")) {
       hid_t rx_group = open_group(rxs_group, name.c_str());
-      reactions_.push_back(make_unique<Reaction>(rx_group, temps_to_read));
+      reactions_.push_back(
+        make_unique<Reaction>(rx_group, temps_to_read, name_));
 
       // Check for 0K elastic scattering
       const auto& rx = reactions_.back();
@@ -375,15 +376,15 @@ void Nuclide::create_derived(
       int j = rx->xs_[t].threshold;
       int n = rx->xs_[t].value.size();
       auto xs = xt::adapt(rx->xs_[t].value);
+      auto pprod = xt::view(xs_[t], xt::range(j, j + n), XS_PHOTON_PROD);
 
       for (const auto& p : rx->products_) {
         if (p.particle_ == ParticleType::photon) {
-          auto pprod = xt::view(xs_[t], xt::range(j, j + n), XS_PHOTON_PROD);
           for (int k = 0; k < n; ++k) {
             double E = grid_[t].energy[k + j];
 
-            // For fission, artificially increase the photon yield to account
-            // for delayed photons
+            // For fission, artificially increase the photon yield to
+            // account for delayed photons
             double f = 1.0;
             if (settings::delayed_photon_scaling) {
               if (is_fission(rx->mt_)) {
@@ -469,8 +470,8 @@ void Nuclide::create_derived(
         }
       }
     } else {
-      // Otherwise, assume that any that have 0 K elastic scattering data are
-      // resonant
+      // Otherwise, assume that any that have 0 K elastic scattering data
+      // are resonant
       resonant_ = !energy_0K_.empty();
     }
 
@@ -781,8 +782,8 @@ void Nuclide::calculate_xs(
       }
 
       for (int j = 0; j < DEPLETION_RX.size(); ++j) {
-        // If reaction is present and energy is greater than threshold, set the
-        // reaction xs appropriately
+        // If reaction is present and energy is greater than threshold, set
+        // the reaction xs appropriately
         int i_rx = reaction_index_[DEPLETION_RX[j]];
         if (i_rx >= 0) {
           const auto& rx = reactions_[i_rx];
@@ -819,9 +820,9 @@ void Nuclide::calculate_xs(
   // Initialize URR probability table treatment to false
   micro.use_ptable = false;
 
-  // If there is S(a,b) data for this nuclide, we need to set the sab_scatter
-  // and sab_elastic cross sections and correct the total and elastic cross
-  // sections.
+  // If there is S(a,b) data for this nuclide, we need to set the
+  // sab_scatter and sab_elastic cross sections and correct the total and
+  // elastic cross sections.
 
   if (i_sab >= 0)
     this->calculate_sab_xs(i_sab, sab_frac, p);
@@ -984,8 +985,8 @@ void Nuclide::calculate_urr_xs(int i_temp, Particle& p) const
   }
 
   // Set elastic, absorption, fission, total, and capture x/s. Note that the
-  // total x/s is calculated as a sum of partials instead of the table-provided
-  // value
+  // total x/s is calculated as a sum of partials instead of the
+  // table-provided value
   micro.elastic = elastic;
   micro.absorption = capture + fission;
   micro.fission = fission;

src/particle.cpp

@@ -75,13 +75,13 @@ double Particle::speed() const
   }
 }
 
-void Particle::create_secondary(
+bool Particle::create_secondary(
   double wgt, Direction u, double E, ParticleType type)
 {
   // If energy is below cutoff for this particle, don't create secondary
   // particle
   if (E < settings::energy_cutoff[static_cast<int>(type)]) {
-    return;
+    return false;
   }
 
   auto& bank = secondary_bank().emplace_back();
@@ -92,6 +92,7 @@ void Particle::create_secondary(
   bank.E = settings::run_CE ? E : g();
   bank.time = time();
   bank_second_E() += bank.E;
+  return true;
 }
 
 void Particle::split(double wgt)
@@ -137,6 +138,7 @@ void Particle::from_source(const SourceSite* src)
   E_last() = E();
   time() = src->time;
   time_last() = src->time;
+  parent_nuclide() = src->parent_nuclide;
 }
 
 void Particle::event_calculate_xs()

src/physics.cpp

@@ -2,6 +2,7 @@
 
 #include "openmc/bank.h"
 #include "openmc/bremsstrahlung.h"
+#include "openmc/chain.h"
 #include "openmc/constants.h"
 #include "openmc/distribution_multi.h"
 #include "openmc/eigenvalue.h"
@@ -1145,9 +1146,21 @@ void sample_secondary_photons(Particle& p, int i_nuclide)
   // Sample the number of photons produced
   double y_t =
     p.neutron_xs(i_nuclide).photon_prod / p.neutron_xs(i_nuclide).total;
-  int y = static_cast<int>(y_t);
-  if (prn(p.current_seed()) <= y_t - y)
-    ++y;
+  double photon_wgt = p.wgt();
+  int y = 1;
+
+  if (settings::use_decay_photons) {
+    // For decay photons, sample a single photon and modify the weight
+    if (y_t <= 0.0)
+      return;
+    photon_wgt *= y_t;
+  } else {
+    // For prompt photons, sample an integral number of photons with weight
+    // equal to the neutron's weight
+    y = static_cast<int>(y_t);
+    if (prn(p.current_seed()) <= y_t - y)
+      ++y;
+  }
 
   // Sample each secondary photon
   for (int i = 0; i < y; ++i) {
@@ -1170,15 +1183,19 @@ void sample_secondary_photons(Particle& p, int i_nuclide)
     // release and deposition. See D. P. Griesheimer, S. J. Douglass, and M. H.
     // Stedry, "Self-consistent energy normalization for quasistatic reactor
     // calculations", Proc. PHYSOR, Cambridge, UK, Mar 29-Apr 2, 2020.
-    double wgt;
+    double wgt = photon_wgt;
     if (settings::run_mode == RunMode::EIGENVALUE && !is_fission(rx->mt_)) {
-      wgt = simulation::keff * p.wgt();
-    } else {
-      wgt = p.wgt();
+      wgt *= simulation::keff;
     }
 
     // Create the secondary photon
-    p.create_secondary(wgt, u, E, ParticleType::photon);
+    bool created_photon = p.create_secondary(wgt, u, E, ParticleType::photon);
+
+    // Tag secondary particle with parent nuclide
+    if (created_photon && settings::use_decay_photons) {
+      p.secondary_bank().back().parent_nuclide =
+        rx->products_[i_product].parent_nuclide_;
+    }
   }
 }
 

src/reaction.cpp

@@ -1,17 +1,20 @@
 #include "openmc/reaction.h"
 
+#include <algorithm> // for remove_if
 #include <string>
 #include <unordered_map>
 #include <utility> // for move
 
 #include <fmt/core.h>
 
+#include "openmc/chain.h"
 #include "openmc/constants.h"
 #include "openmc/endf.h"
 #include "openmc/hdf5_interface.h"
 #include "openmc/random_lcg.h"
 #include "openmc/search.h"
 #include "openmc/secondary_uncorrelated.h"
+#include "openmc/settings.h"
 
 namespace openmc {
 
@@ -19,7 +22,8 @@ namespace openmc {
 // Reaction implementation
 //==============================================================================
 
-Reaction::Reaction(hid_t group, const vector<int>& temperatures)
+Reaction::Reaction(
+  hid_t group, const vector<int>& temperatures, std::string name)
 {
   read_attribute(group, "Q_value", q_value_);
   read_attribute(group, "mt", mt_);
@@ -63,6 +67,34 @@ Reaction::Reaction(hid_t group, const vector<int>& temperatures)
       close_group(pgroup);
     }
   }
+
+  if (settings::use_decay_photons) {
+    // Remove photon products for D1S method
+    products_.erase(
+      std::remove_if(products_.begin(), products_.end(),
+        [](const auto& p) { return p.particle_ == ParticleType::photon; }),
+      products_.end());
+
+    // Determine product for D1S method
+    auto nuclide_it = data::chain_nuclide_map.find(name);
+    if (nuclide_it != data::chain_nuclide_map.end()) {
+      const auto& chain_nuc = data::chain_nuclides[nuclide_it->second];
+      const auto& rx_products = chain_nuc->reaction_products();
+      auto product_it = rx_products.find(mt_);
+      if (product_it != rx_products.end()) {
+        auto decay_products = product_it->second;
+        for (const auto& decay_product : decay_products) {
+          auto product_it = data::chain_nuclide_map.find(decay_product.name);
+          if (product_it != data::chain_nuclide_map.end()) {
+            const auto& product_nuc = data::chain_nuclides[product_it->second];
+            if (product_nuc->photon_energy()) {
+              products_.emplace_back(decay_product);
+            }
+          }
+        }
+      }
+    }
+  }
 }
 
 double Reaction::xs(int64_t i_temp, int64_t i_grid, double interp_factor) const

src/reaction_product.cpp

@@ -83,6 +83,29 @@ ReactionProduct::ReactionProduct(hid_t group)
   }
 }
 
+ReactionProduct::ReactionProduct(const ChainNuclide::Product& product)
+{
+  particle_ = ParticleType::photon;
+  emission_mode_ = EmissionMode::delayed;
+
+  // Get chain nuclide object for radionuclide
+  parent_nuclide_ = data::chain_nuclide_map.at(product.name);
+  const auto& chain_nuc = data::chain_nuclides[parent_nuclide_].get();
+
+  // Determine decay constant in [s^-1]
+  decay_rate_ = chain_nuc->decay_constant();
+
+  // Determine number of photons per decay and set yield
+  double photon_per_sec = chain_nuc->photon_energy()->integral();
+  double photon_per_decay = photon_per_sec / decay_rate_;
+  vector<double> coef = {product.branching_ratio * photon_per_decay};
+  yield_ = make_unique<Polynomial>(coef);
+
+  // Set decay photon angle-energy distribution
+  distribution_.push_back(
+    make_unique<DecayPhotonAngleEnergy>(chain_nuc->photon_energy()));
+}
+
 void ReactionProduct::sample(
   double E_in, double& E_out, double& mu, uint64_t* seed) const
 {

src/settings.cpp

@@ -75,6 +75,7 @@ bool trigger_predict {false};
 bool uniform_source_sampling {false};
 bool ufs_on {false};
 bool urr_ptables_on {true};
+bool use_decay_photons {false};
 bool weight_windows_on {false};
 bool weight_window_checkpoint_surface {false};
 bool weight_window_checkpoint_collision {true};
@@ -1124,6 +1125,11 @@ void read_settings_xml(pugi::xml_node root)
         get_node_value_bool(ww_checkpoints, "surface");
     }
   }
+
+  if (check_for_node(root, "use_decay_photons")) {
+    settings::use_decay_photons =
+      get_node_value_bool(root, "use_decay_photons");
+  }
 }
 
 void free_memory_settings()

src/tallies/filter.cpp

@@ -28,6 +28,7 @@
 #include "openmc/tallies/filter_meshsurface.h"
 #include "openmc/tallies/filter_mu.h"
 #include "openmc/tallies/filter_musurface.h"
+#include "openmc/tallies/filter_parent_nuclide.h"
 #include "openmc/tallies/filter_particle.h"
 #include "openmc/tallies/filter_polar.h"
 #include "openmc/tallies/filter_sph_harm.h"
@@ -137,6 +138,8 @@ Filter* Filter::create(const std::string& type, int32_t id)
     return Filter::create<MuFilter>(id);
   } else if (type == "musurface") {
     return Filter::create<MuSurfaceFilter>(id);
+  } else if (type == "parentnuclide") {
+    return Filter::create<ParentNuclideFilter>(id);
   } else if (type == "particle") {
     return Filter::create<ParticleFilter>(id);
   } else if (type == "polar") {

src/tallies/filter_parent_nuclide.cpp

@@ -0,0 +1,79 @@
+#include "openmc/tallies/filter_parent_nuclide.h"
+
+#include <cstdint> // for int64_t
+
+#include <fmt/core.h>
+
+#include "openmc/capi.h"
+#include "openmc/chain.h"
+#include "openmc/search.h"
+#include "openmc/settings.h"
+#include "openmc/xml_interface.h"
+
+namespace openmc {
+
+//==============================================================================
+// ParentNuclideFilter implementation
+//==============================================================================
+
+void ParentNuclideFilter::from_xml(pugi::xml_node node)
+{
+  nuclides_ = get_node_array<std::string>(node, "bins");
+
+  // Convert nuclides to indices in data::chain_nuclides
+  std::vector<int> bins;
+  for (const auto& nuclide : nuclides_) {
+    auto it = data::chain_nuclide_map.find(nuclide);
+    if (it != data::chain_nuclide_map.end()) {
+      bins.push_back(it->second);
+    } else {
+      // The default value of parent_nuclide is -1, so to prevent a score to
+      // this bin assign the value -2.
+      bins.push_back(-2);
+    }
+  }
+  this->set_bins(bins);
+}
+
+void ParentNuclideFilter::set_bins(span<const int> bins)
+{
+  // Clear existing bins
+  bins_.clear();
+  bins_.reserve(bins.size());
+  map_.clear();
+
+  // Set bins based on chain nuclide indexing
+  for (int64_t i = 0; i < bins.size(); ++i) {
+    bins_.push_back(bins[i]);
+    map_[bins[i]] = i;
+  }
+
+  n_bins_ = bins_.size();
+}
+
+void ParentNuclideFilter::get_all_bins(
+  const Particle& p, TallyEstimator estimator, FilterMatch& match) const
+{
+  // Get the particle's parent nuclide
+  int parent_nuclide = p.parent_nuclide();
+
+  // Find bin matching parent nuclide
+  auto search = map_.find(parent_nuclide);
+  if (search != map_.end()) {
+    match.bins_.push_back(search->second);
+    match.weights_.push_back(1.0);
+  }
+}
+
+void ParentNuclideFilter::to_statepoint(hid_t filter_group) const
+{
+  Filter::to_statepoint(filter_group);
+  write_dataset(filter_group, "bins", nuclides_);
+}
+
+std::string ParentNuclideFilter::text_label(int bin) const
+{
+  return fmt::format("Parent Nuclide {}", nuclides_[bin]);
+}
+
+} // namespace openmc

tests/chain_ni.xml

@@ -0,0 +1,173 @@
+<depletion_chain>
+  <nuclide name="H1" reactions="1">
+    <reaction type="(n,gamma)" Q="2224631.0"/>
+  </nuclide>
+  <nuclide name="He4" reactions="0"/>
+  <nuclide name="Fe55" half_life="86594050.0" decay_modes="1" decay_energy="5842.102" reactions="0">
+    <decay type="ec/beta+" branching_ratio="1.0"/>
+    <source type="discrete" particle="photon">
+      <parameters>557.6039 640.4896 655.72 678.8113 5847.93 5858.68 6448.81 6450.12 6499.05 6499.18 126000.0 4.281895544855688e-11 8.113092795419058e-12 2.881383486075155e-13 8.900165070824795e-12 6.680021392530635e-10 1.3098276096378766e-09 7.960187308411459e-11 1.5666266690172673e-10 2.8061406316612144e-14 4.106367573258232e-14 1.0245835494664239e-17</parameters>
+    </source>
+    <source type="discrete" particle="positron">
+      <parameters>105210.0 231210.0 1.040592667426837e-17 8.004558980206438e-09</parameters>
+    </source>
+    <source type="discrete" particle="electron">
+      <parameters>36.67648 65.20168 562.7812 610.8647 5155.619 5764.491 6377.571 2.257651440763611e-08 2.654244519541021e-09 1.1943434358627441e-08 2.0711244046714524e-10 3.863926291321637e-09 9.499586470057554e-10 5.918783030777615e-11</parameters>
+    </source>
+  </nuclide>
+  <nuclide name="Fe57" reactions="4">
+    <reaction type="(n,2n)" Q="-7646000.0"/>
+    <reaction type="(n,gamma)" Q="10044540.0" target="Fe58"/>
+    <reaction type="(n,p)" Q="-1909000.0"/>
+    <reaction type="(n,a)" Q="2400000.0"/>
+  </nuclide>
+  <nuclide name="Fe58" reactions="4">
+    <reaction type="(n,2n)" Q="-10044000.0" target="Fe57"/>
+    <reaction type="(n,gamma)" Q="6580960.0" target="Fe59"/>
+    <reaction type="(n,p)" Q="-5159000.0"/>
+    <reaction type="(n,a)" Q="-1399000.0"/>
+  </nuclide>
+  <nuclide name="Fe59" half_life="3844368.0" decay_modes="1" decay_energy="1306009.3" reactions="0">
+    <decay type="beta-" branching_ratio="1.0"/>
+    <source type="discrete" particle="photon">
+      <parameters>780.0 6915.0 6930.0 7649.0 142651.0 189000.0 192343.0 334800.0 382000.0 1099245.0 1291590.0 1481700.0 1.0526870077721752e-12 1.0971242349369013e-11 2.1648606887732874e-11 3.9280951004469605e-12 1.8390802445841403e-09 1.622717862868359e-12 5.553301130705051e-09 4.868153588605077e-10 3.245435725736718e-11 1.018706213911803e-07 7.789045741768123e-08 1.0637817101025909e-10</parameters>
+    </source>
+    <source type="discrete" particle="electron">
+      <parameters>750.0 6070.0 83569.95 130946.9 134942.1 141725.4 184634.1 191417.4 273599.0 327091.1 333874.4 465943.0 1091536.0 1283881.0 1565200.0 1.4515384373262843e-10 5.994999523918275e-11 1.4063554811525778e-10 2.3619560003972782e-09 2.703447959538686e-11 2.7218387619845277e-12 4.520385317374064e-11 4.498173915871091e-12 8.16767990977074e-08 8.519268780058884e-13 8.275861100628631e-14 9.574035390923319e-08 1.6299297619569002e-11 8.567950315944936e-12 3.245435725736718e-10</parameters>
+    </source>
+  </nuclide>
+  <nuclide name="Fe61" half_life="358.8" decay_modes="1" decay_energy="2515414.0" reactions="0">
+    <decay type="beta-" target="Co61" branching_ratio="1.0"/>
+    <source type="discrete" particle="photon">
+      <parameters>120340.0 177610.0 297900.0 333000.0 349700.0 440500.0 542600.0 561400.0 603300.0 618400.0 657300.0 686000.0 696900.0 748100.0 769400.0 806300.0 925600.0 945400.0 978000.0 984100.0 989200.0 1027420.0 1097800.0 1205070.0 1275000.0 1285700.0 1381400.0 1403900.0 1538800.0 1618900.0 1645950.0 1659300.0 1837200.0 1879400.0 1889000.0 1899300.0 1972700.0 1999800.0 2011600.0 2177100.0 2230800.0 2484400.0 2754400.0 2920000.0 3191000.0 3204200.0 3239100.0 3364900.0 0.00010275887633317896 3.874515009283797e-05 0.00042956579450755135 4.2956579450755134e-06 3.116457724858706e-06 4.21142935791717e-06 1.263428807375151e-06 1.0949716330584643e-06 1.431885981691838e-06 1.802491765188549e-05 4.2956579450755134e-06 7.749030018567594e-06 2.1899432661169287e-06 1.558228862429353e-05 3.116457724858706e-06 3.7060578349671097e-06 6.569829798350785e-06 2.105714678958585e-06 1.431885981691838e-06 1.1792002202168078e-05 1.1792002202168078e-05 0.0008254401541517653 1.3476573945334946e-05 0.000842285871583434 1.1792002202168078e-05 7.159429908459189e-06 7.749030018567594e-06 2.2741718532752718e-06 5.306400990975635e-06 7.075201321300846e-06 0.00013476573945334945 1.4992688514185127e-05 2.6953147890669893e-06 5.053715229500604e-06 3.4533720734920796e-06 1.431885981691838e-06 1.1792002202168076e-06 2.526857614750302e-06 8.507087302992683e-05 4.042972183600483e-06 2.105714678958585e-06 2.3584004404336153e-06 1.4824231339868439e-05 1.3476573945334947e-06 1.6003431560085245e-06 8.422858715834341e-07 1.0949716330584643e-06 8.422858715834341e-07</parameters>
+    </source>
+    <source type="discrete" particle="electron">
+      <parameters>613000.0 738800.1 773300.1 786899.9 873600.1 977699.9 1057900.0 1113610.0 1223500.0 1493500.0 1546600.0 1675000.0 1747100.0 1966500.0 2024800.0 2089000.0 2332110.0 2359120.0 2652590.0 2692270.0 2772900.0 2950510.0 3978000.0 7.969724371709087e-07 1.032441566335041e-06 1.1773456458206606e-05 1.5214928345990077e-06 1.267910695499173e-06 2.535821390998346e-05 1.267910695499173e-06 2.1735611922842964e-05 1.3947017650490904e-05 2.227900222091404e-06 2.173561192284297e-06 2.8980815897123955e-05 1.0867805961421485e-06 9.056504967851236e-05 1.7931879836345447e-05 4.709382583282643e-06 0.00014490407948561978 2.4452563413198337e-06 0.0004890512682639667 3.622601987140494e-07 0.0006882943775566939 0.00028980815897123956 9.056504967851236e-05</parameters>
+    </source>
+  </nuclide>
+  <nuclide name="Co58" half_life="6122304.0" decay_modes="1" decay_energy="1008776.54" reactions="4">
+    <decay type="ec/beta+" target="Fe58" branching_ratio="1.0"/>
+    <source type="discrete" particle="photon">
+      <parameters>619.9594 707.9362 726.02 741.741 6349.85 6362.71 7015.36 7016.95 7070.49 7070.66 510998.9 810759.2 863951.0 1674725.0 5.276288639112756e-10 1.3715809425619083e-10 3.759818598537813e-12 1.1767757246637665e-10 8.777123989590544e-09 1.719716697150835e-08 1.057271626846154e-09 2.0783510382283326e-09 5.201844116739253e-13 7.595270495013994e-13 3.374000922877054e-08 1.1259402850084962e-07 7.768987966558623e-10 5.854889482044179e-10</parameters>
+    </source>
+    <source type="discrete" particle="positron">
+      <parameters>632864.1 1496834.0 2307600.0 1.3699223176071195e-09 1.1185811981783752e-07 8.491254034754873e-13</parameters>
+    </source>
+    <source type="discrete" particle="electron">
+      <parameters>40.23492 70.00108 628.9213 675.1537 5583.843 6258.984 6938.274 803647.2 809913.1 810759.2 856839.0 863104.9 863951.0 1667613.0 1673879.0 1674725.0 2.3068291731266982e-07 2.4528685530196402e-08 1.3309519902564132e-07 1.9030055102434398e-09 4.481805021434518e-08 1.1035282426757231e-08 6.934054278993223e-10 3.366561452175404e-11 3.231448617974384e-12 4.640000480603615e-13 1.6314869068937087e-13 1.553797593311725e-14 2.2436836794554425e-15 3.3782714575729324e-14 3.2026247731116077e-15 4.614238015862169e-16</parameters>
+    </source>
+    <reaction type="(n,2n)" Q="-8572000.0"/>
+    <reaction type="(n,gamma)" Q="10454000.0"/>
+    <reaction type="(n,p)" Q="0.0" target="Fe58"/>
+    <reaction type="(n,a)" Q="0.0"/>
+  </nuclide>
+  <nuclide name="Co60" half_life="166344200.0" decay_modes="1" decay_energy="2600613.1" reactions="0">
+    <decay type="beta-" target="Ni60" branching_ratio="1.0"/>
+    <source type="discrete" particle="photon">
+      <parameters>751.0021 852.3381 876.89 883.6421 7417.82 7435.78 8222.319 8224.59 8287.88 8288.141 347140.0 826100.0 1173228.0 1332492.0 2158570.0 2505692.0 7.0614888325101086e-15 3.3564685966506987e-15 5.761418570090495e-17 2.3600648888691433e-15 1.3290046971181013e-13 2.5954235427719607e-13 1.6253028424014347e-14 3.186921843254948e-14 1.3932264913028456e-17 2.024904561736029e-17 3.125208966828774e-13 3.1668784197198245e-13 4.160694871171375e-09 4.166220240624728e-09 5.000334346926039e-14 8.333890578210064e-17</parameters>
+    </source>
+    <source type="discrete" particle="electron">
+      <parameters>47.97734 78.78671 772.594 818.3802 6497.883 7313.626 8133.5 318151.8 1164895.0 1172220.0 1173228.0 1324159.0 1331484.0 1332492.0 1491392.0 2497359.0 2504684.0 2505692.0 2.430147074908485e-12 3.047352927658078e-13 1.6375715794161468e-12 1.653965175572544e-14 5.18885111875044e-13 1.2976030141139345e-13 8.558901456876446e-15 4.161944954758106e-09 6.241042306757062e-13 6.095415902793419e-14 8.94050113917305e-15 4.736991738545178e-13 4.616171159887578e-14 6.779692031857936e-15 5.000334346926038e-12 6.46709908869101e-21 6.3170890582832295e-22 9.678980517533168e-23</parameters>
+    </source>
+  </nuclide>
+  <nuclide name="Co61" half_life="5940.0" decay_modes="1" decay_energy="563911.34" reactions="0">
+    <decay type="beta-" target="Ni61" branching_ratio="1.0"/>
+    <source type="discrete" particle="photon">
+      <parameters>751.028 852.3382 876.89 883.6886 7417.82 7435.78 8222.319 8224.59 8287.88 8288.141 67415.0 841700.0 909200.0 7.866502846971466e-08 3.756939402646081e-08 6.436568219422988e-10 2.637348005677196e-08 1.4671031864629482e-06 2.8654120987125446e-06 1.794509040062992e-07 3.5187803327436886e-07 1.538314140492697e-10 2.235774236842793e-10 9.880497991981766e-05 9.262966867482906e-07 4.219796017408879e-06</parameters>
+    </source>
+    <source type="discrete" particle="electron">
+      <parameters>48.00067 78.83981 772.8435 818.5774 6498.158 7313.756 8133.501 59082.2 66406.9 67415.0 413499.9 833367.2 840691.9 841700.0 900867.2 908191.9 909200.0 1255285.0 2.6847340674754763e-05 3.3502043712197297e-06 1.814135374119958e-05 1.8444716489566865e-07 5.715307736044776e-06 1.430796977323619e-06 9.449890224831685e-08 1.2103610040177662e-05 1.246919266677299e-06 1.902588515965013e-07 5.134423559703298e-06 2.917834563257115e-10 2.862256762052218e-11 4.385088515066407e-12 9.241353278125446e-10 9.030363477255e-11 1.3801262380678854e-11 0.00011155702097900803</parameters>
+    </source>
+  </nuclide>
+  <nuclide name="Co62" half_life="90.0" decay_modes="1" decay_energy="3236869.0" reactions="0">
+    <decay type="beta-" target="Ni62" branching_ratio="1.0"/>
+    <source type="discrete" particle="photon">
+      <parameters>1128900.0 1172900.0 1886300.0 1985000.0 2083000.0 2097000.0 2301800.0 2345900.0 3158200.0 3271000.0 3369500.0 3519000.0 4063100.0 0.0008553051126342747 0.006430865508528381 3.2154327542641905e-05 0.00012218644466203924 2.5723462034113525e-05 7.07395205938122e-05 0.0011318323295009952 0.0001028938481364541 6.430865508528381e-05 1.1254014639924668e-05 2.8938894788377717e-05 6.430865508528381e-06 2.5723462034113525e-05</parameters>
+    </source>
+    <source type="discrete" particle="electron">
+      <parameters>1251800.0 1796100.0 1945400.0 2044800.0 2059000.0 2156900.0 2255700.0 3013160.0 4142070.0 5315000.0 2.5415396620531326e-05 0.00010782289475376927 2.9266214290308802e-05 9.318978760861486e-05 2.5415396620531326e-05 0.00018483924814931876 3.080654135821979e-05 0.0019870219176051766 0.005221708760218255 1.925408834888737e-05</parameters>
+    </source>
+  </nuclide>
+  <nuclide name="Co64" half_life="0.3" decay_modes="1" decay_energy="3469635.0" reactions="0">
+    <decay type="beta-" target="Ni64" branching_ratio="1.0"/>
+    <source type="discrete" particle="photon">
+      <parameters>931100.0 1346100.0 0.11552453009332422 0.23104906018664845</parameters>
+    </source>
+    <source type="discrete" particle="electron">
+      <parameters>5029800.0 5960900.0 7307000.0 0.11552453009332422 0.11552453009332422 2.079441541679836</parameters>
+    </source>
+  </nuclide>
+  <nuclide name="Ni57" half_life="128160.0" decay_modes="1" decay_energy="2092651.3" reactions="0">
+    <decay type="ec/beta+" branching_ratio="1.0"/>
+    <source type="discrete" particle="photon">
+      <parameters>683.2094 778.5527 799.73 813.244 6873.16 6888.41 7606.53 7608.44 7666.77 7666.98 127164.0 161860.0 304100.0 379940.0 510998.9 541900.0 673440.0 696000.0 755300.0 906980.0 1046680.0 1224000.0 1279990.0 1350520.0 1377630.0 1603280.0 1730440.0 1757550.0 1897420.0 1919520.0 2133040.0 2730910.0 2804200.0 3177280.0 1.8638833690911883e-08 6.653293496572297e-09 1.4203876186268797e-10 5.409592966721624e-09 3.0832339673226584e-07 6.028687625463321e-07 3.74301586798566e-08 7.349918311628966e-08 2.4690572137876036e-11 3.5966442331457405e-11 9.014158418349326e-07 1.2284000197554473e-09 1.0604892256881559e-10 3.6233381877678664e-09 4.684800934777033e-06 1.9884172981652923e-10 2.6556417693274235e-09 4.860575617737382e-11 2.916345370642429e-10 3.3140288302754875e-09 7.246676375535733e-09 3.402402932416167e-09 5.2140720263001e-10 1.0604892256881559e-10 4.418705107033983e-06 2.1209784513763118e-10 2.827971268501749e-09 3.110768395351924e-07 1.502359736391554e-09 6.628057660550975e-07 1.546546787461894e-09 1.0737453410092578e-09 5.30244612844078e-09 6.009438945566218e-10</parameters>
+    </source>
+    <source type="discrete" particle="positron">
+      <parameters>84790.04 154079.8 457929.9 531290.1 1129120.0 1342650.0 1504620.0 1757390.0 1884550.0 1.124957970946228e-09 3.2450710700371965e-09 1.5738594689680403e-08 1.0762819048956704e-09 1.8388736063544118e-09 6.652395693576253e-07 3.0611837094017556e-07 9.194368031772059e-07 3.488451400289987e-06</parameters>
+    </source>
+    <source type="discrete" particle="electron">
+      <parameters>43.92123 75.12743 700.3887 749.9867 6033.611 6776.5 7523.621 119455.1 126238.4 127164.0 154151.1 160934.4 161860.0 665731.1 672514.4 673440.0 747591.1 754374.4 755300.0 1216291.0 1223074.0 1224000.0 1369921.0 1376704.0 1377630.0 1749841.0 1756624.0 1757550.0 1889711.0 1896494.0 1897420.0 1911811.0 1918594.0 1919520.0 2125331.0 2132114.0 2133040.0 7.252641897852818e-06 9.885065183674125e-07 4.453862207006767e-06 5.575973168934215e-08 1.3531437967587472e-06 3.3801747364315905e-07 2.1887831296943824e-08 1.739732899248527e-08 1.7397323584033486e-09 2.636641229198142e-10 1.272622225762379e-11 1.2627950688719498e-12 1.9160584871075314e-13 9.58686478614484e-13 9.268187990163621e-14 1.4090839543995835e-14 8.632381215411233e-14 8.34074667834699e-15 1.2155339619769639e-15 3.6814000810433283e-13 3.538499049712814e-14 5.158725716665081e-15 4.21544467211042e-10 4.0607899933642303e-11 2.1687004665322787e-10 1.822910020070542e-11 1.7482517516525527e-12 6.153697369032366e-11 7.181277376570916e-14 6.8808054292926055e-15 3.195035376537882e-13 3.115187100458958e-11 2.9826259472479385e-12 1.482406852545487e-10 6.031532471101388e-14 5.768616813006974e-15 4.833646927113489e-13</parameters>
+    </source>
+  </nuclide>
+  <nuclide name="Ni58" reactions="4">
+    <reaction type="(n,2n)" Q="-12196000.0" target="Ni57"/>
+    <reaction type="(n,gamma)" Q="8999510.0" target="Ni59"/>
+    <reaction type="(n,p)" Q="402000.0" target="Co58"/>
+    <reaction type="(n,a)" Q="2901000.0" target="Fe55"/>
+  </nuclide>
+  <nuclide name="Ni59" half_life="2398380000000.0" decay_modes="1" decay_energy="6915.5689999999995" reactions="4">
+    <decay type="ec/beta+" branching_ratio="1.0"/>
+    <source type="discrete" particle="positron">
+      <parameters>1072500.0 2.890064045563861e-13</parameters>
+    </source>
+    <source type="discrete" particle="electron">
+      <parameters>44.06754 75.28693 701.5928 750.7433 6035.851 6777.615 7523.68 6.689278658453111e-13 8.913674252402247e-14 4.1657845563004787e-13 5.378516121164032e-15 1.2583954438026756e-13 3.16650156158609e-14 2.0701132819599697e-15</parameters>
+    </source>
+    <source type="discrete" particle="photon">
+      <parameters>683.3271 778.5529 799.73 813.5109 6873.16 6888.41 7606.53 7608.44 7666.77 7666.98 510998.9 1.7758787515451584e-15 6.581346326607483e-16 1.368589539770362e-17 5.248596182059754e-16 2.920762305855841e-14 5.711006429349971e-14 3.5457762265416177e-15 6.962612245690403e-15 2.338949640535961e-18 3.4071167337635547e-18 2.1386473937172572e-19</parameters>
+    </source>
+    <reaction type="(n,2n)" Q="-8997200.0" target="Ni58"/>
+    <reaction type="(n,gamma)" Q="11240000.0" target="Ni60"/>
+    <reaction type="(n,p)" Q="1584000.0"/>
+    <reaction type="(n,a)" Q="4878000.0"/>
+  </nuclide>
+  <nuclide name="Ni60" reactions="4">
+    <reaction type="(n,2n)" Q="-11389000.0" target="Ni59"/>
+    <reaction type="(n,gamma)" Q="7820050.0" target="Ni61"/>
+    <reaction type="(n,p)" Q="-2042000.0" target="Co60"/>
+    <reaction type="(n,a)" Q="1355000.0" target="Fe57"/>
+  </nuclide>
+  <nuclide name="Ni61" reactions="4">
+    <reaction type="(n,2n)" Q="-7820300.0" target="Ni60"/>
+    <reaction type="(n,gamma)" Q="10597330.0" target="Ni62"/>
+    <reaction type="(n,p)" Q="-539700.0" target="Co61"/>
+    <reaction type="(n,a)" Q="3579000.0" target="Fe58"/>
+  </nuclide>
+  <nuclide name="Ni62" reactions="4">
+    <reaction type="(n,2n)" Q="-10598000.0" target="Ni61"/>
+    <reaction type="(n,gamma)" Q="6837910.0" target="Ni63"/>
+    <reaction type="(n,p)" Q="-4459000.0" target="Co62"/>
+    <reaction type="(n,a)" Q="-437400.0" target="Fe59"/>
+  </nuclide>
+  <nuclide name="Ni63" half_life="3193630000.0" decay_modes="1" decay_energy="17425.0" reactions="0">
+    <decay type="beta-" branching_ratio="1.0"/>
+    <source type="discrete" particle="electron">
+      <parameters>66945.0 2.1704054025041888e-10</parameters>
+    </source>
+  </nuclide>
+  <nuclide name="Ni64" reactions="4">
+    <reaction type="(n,2n)" Q="-9657000.0" target="Ni63"/>
+    <reaction type="(n,gamma)" Q="6098061.0" target="Ni65"/>
+    <reaction type="(n,p)" Q="-6524000.0" target="Co64"/>
+    <reaction type="(n,a)" Q="-2442000.0" target="Fe61"/>
+  </nuclide>
+  <nuclide name="Ni65" half_life="9061.884" decay_modes="1" decay_energy="1186322.7999999998" reactions="0">
+    <decay type="beta-" branching_ratio="1.0"/>
+    <source type="discrete" particle="photon">
+      <parameters>820.0389 927.6652 951.42 959.5527 7984.67 8005.709 8862.7 8865.34 8933.11 8933.4 366270.0 507900.0 609500.0 770600.0 852700.0 954500.0 1115530.0 1481840.0 1623420.0 1724920.0 6.930967733828076e-11 4.4848940639261224e-11 6.225110448569713e-13 3.0091057755445835e-11 1.4120610720530457e-09 2.7533685191404717e-09 1.7112326402414573e-10 3.3482649978514024e-10 1.9318959861680173e-13 2.7929997279698636e-13 3.6755804457910027e-06 2.2374667416695354e-07 1.181887674027053e-07 7.92135402897521e-08 7.398075516810562e-08 1.3533064969775417e-09 1.1800832653644161e-05 1.804408662636722e-05 3.807302278163484e-07 3.0494506398560603e-07</parameters>
+    </source>
+    <source type="discrete" particle="electron">
+      <parameters>56.92597 90.91714 855.942 897.216 6991.356 7880.445 8773.007 357291.1 365173.4 366270.0 412919.9 498921.1 506803.4 507900.0 514460.0 600521.1 608403.4 609500.0 656069.9 761621.1 769503.4 770600.0 1022350.0 1106551.0 1114434.0 1115530.0 1472861.0 1480743.0 1481840.0 1614441.0 1622324.0 1623420.0 1715941.0 1723824.0 1724920.0 2137900.0 2.167548391930537e-08 2.8616659297729836e-09 1.5876524466163717e-08 5.919283450210211e-11 4.8238674145478966e-09 1.1910611683428657e-09 8.251524863747208e-11 6.6528004309537836e-09 6.652800430953784e-10 9.941708721258623e-11 4.2452174979371803e-07 1.968970885650002e-10 1.9600207433178643e-11 2.840463936760989e-12 6.807646077773136e-07 6.748578656939676e-11 6.701302996997783e-12 9.996408692926368e-13 2.172327512457944e-05 2.732867560693677e-11 2.7011821751739386e-12 4.0351373644973695e-13 7.786723266486576e-06 1.953037934211798e-09 1.9247158318161006e-10 3.950919133474779e-11 1.818843931937816e-09 1.7899733933356285e-10 1.4705028396157966e-09 3.076300301948419e-11 3.022997962967563e-12 4.574485034565062e-11 2.213900812679635e-11 2.1742581302894385e-12 5.0020380532292804e-11 4.589424322094249e-05</parameters>
+    </source>
+  </nuclide>
+</depletion_chain>

tests/regression_tests/surface_source/test.py

@@ -10,6 +10,17 @@ from tests.testing_harness import PyAPITestHarness
 from tests.regression_tests import config
 
 
+def assert_structured_arrays_close(arr1, arr2, rtol=1e-5, atol=1e-8):
+    assert arr1.dtype == arr2.dtype
+
+    for field in arr1.dtype.names:
+        data1, data2 = arr1[field], arr2[field]
+        if data1.dtype.names:
+            assert_structured_arrays_close(data1, data2, rtol=rtol, atol=atol)
+        else:
+            np.testing.assert_allclose(data1, data2, rtol=rtol, atol=atol)
+
+
 @pytest.fixture
 def model(request):
     openmc.reset_auto_ids()
@@ -76,16 +87,10 @@ class SurfaceSourceTestHarness(PyAPITestHarness):
         """Make sure the current surface_source.h5 agree with the reference."""
         if self._model.settings.surf_source_write:
             with h5py.File("surface_source_true.h5", 'r') as f:
-                source_true = f['source_bank'][()]
-                # Convert dtye from mixed to a float for comparison assertion
-                source_true.dtype = 'float64'
+                source_true = np.sort(f['source_bank'][()])
             with h5py.File("surface_source.h5", 'r') as f:
-                source_test = f['source_bank'][()]
-                # Convert dtye from mixed to a float for comparison assertion
-                source_test.dtype = 'float64'
-            np.testing.assert_allclose(np.sort(source_true),
-                                       np.sort(source_test),
-                                       atol=1e-07)
+                source_test = np.sort(f['source_bank'][()])
+            assert_structured_arrays_close(source_true, source_test, atol=1e-07)
 
     def execute_test(self):
         """Build input XMLs, run OpenMC, check output and results."""

tests/unit_tests/test_d1s.py

@@ -0,0 +1,132 @@
+from pathlib import Path
+from math import exp
+
+import numpy as np
+import pytest
+import openmc
+import openmc.deplete
+from openmc.deplete import d1s
+
+
+CHAIN_PATH = Path(__file__).parents[1] / "chain_ni.xml"
+
+
+@pytest.fixture
+def model():
+    """Simple model with natural Ni"""
+    mat = openmc.Material()
+    mat.add_element('Ni', 1.0)
+    geom = openmc.Geometry([openmc.Cell(fill=mat)])
+    return openmc.Model(geometry=geom)
+
+
+def test_get_radionuclides(model):
+    # Check that radionuclides are correct and are unstable
+    nuclides = d1s.get_radionuclides(model, CHAIN_PATH)
+    assert sorted(nuclides) == [
+        'Co58', 'Co60', 'Co61', 'Co62', 'Co64',
+        'Fe55', 'Fe59', 'Fe61', 'Ni57', 'Ni59', 'Ni63', 'Ni65'
+    ]
+    for nuc in nuclides:
+        assert openmc.data.half_life(nuc) is not None
+
+
+@pytest.mark.parametrize("nuclide", ['Co60', 'Ni63', 'H3', 'Na24', 'K40'])
+def test_time_correction_factors(nuclide):
+    # Irradiation schedule turning unit neutron source on and off
+    timesteps = [1.0, 1.0, 1.0]
+    source_rates = [1.0, 0.0, 1.0]
+
+    # Compute expected solution
+    decay_rate = openmc.data.decay_constant(nuclide)
+    g = exp(-decay_rate)
+    expected = [0.0, (1 - g), (1 - g)*g, (1 - g)*(1 + g*g)]
+
+    # Test against expected solution
+    tcf = d1s.time_correction_factors([nuclide], timesteps, source_rates)
+    assert tcf[nuclide] == pytest.approx(expected)
+
+    # Make sure all values at first timestep and onward are positive (K40 case
+    # has very small decay constant that stresses this)
+    assert np.all(tcf[nuclide][1:] > 0.0)
+
+    # Timesteps as a tuple
+    timesteps = [(1.0, 's'), (1.0, 's'), (1.0, 's')]
+    tcf = d1s.time_correction_factors([nuclide], timesteps, source_rates)
+    assert tcf[nuclide] == pytest.approx(expected)
+
+    # Test changing units
+    timesteps = [1.0/60.0, 1.0/60.0, 1.0/60.0]
+    tcf = d1s.time_correction_factors([nuclide], timesteps, source_rates,
+                                      timestep_units='min')
+    assert tcf[nuclide] == pytest.approx(expected)
+
+
+def test_prepare_tallies(model):
+    tally = openmc.Tally()
+    tally.filters = [openmc.ParticleFilter('photon')]
+    tally.scores = ['flux']
+    model.tallies = [tally]
+
+    # Check that prepare_tallies adds a ParentNuclideFilter
+    nuclides = ['Co58', 'Co60', 'Fe55']
+    d1s.prepare_tallies(model, nuclides, chain_file=CHAIN_PATH)
+    assert tally.contains_filter(openmc.ParentNuclideFilter)
+    assert list(tally.filters[-1].bins) == nuclides
+
+    # Get rid of parent nuclide filter
+    tally.filters.pop()
+
+    # With no nuclides specified, filter should use get_radionuclides
+    radionuclides = d1s.get_radionuclides(model, CHAIN_PATH)
+    d1s.prepare_tallies(model, chain_file=CHAIN_PATH)
+    assert tally.contains_filter(openmc.ParentNuclideFilter)
+    assert sorted(tally.filters[-1].bins) == sorted(radionuclides)
+
+
+def test_apply_time_correction(run_in_tmpdir):
+    # Make simple sphere model with elemental Ni
+    mat = openmc.Material()
+    mat.add_element('Ni', 1.0)
+    sphere = openmc.Sphere(r=10.0, boundary_type='vacuum')
+    cell = openmc.Cell(fill=mat, region=-sphere)
+    model = openmc.Model()
+    model.geometry = openmc.Geometry([cell])
+    model.settings.run_mode = 'fixed source'
+    model.settings.batches = 3
+    model.settings.particles = 10
+    model.settings.photon_transport = True
+    model.settings.use_decay_photons = True
+    particle_filter = openmc.ParticleFilter('photon')
+    tally = openmc.Tally()
+    tally.filters = [particle_filter]
+    tally.scores = ['flux']
+    model.tallies = [tally]
+
+    # Prepare tallies for D1S and compute time correction factors
+    nuclides = d1s.prepare_tallies(model, chain_file=CHAIN_PATH)
+    factors = d1s.time_correction_factors(nuclides, [1.0e10], [1.0])
+
+    # Run OpenMC and get tally result
+    with openmc.config.patch('chain_file', CHAIN_PATH):
+        output_path = model.run()
+    with openmc.StatePoint(output_path) as sp:
+        tally = sp.tallies[tally.id]
+        flux = tally.mean.flatten()
+
+    # Apply TCF and make sure results are consistent
+    result = d1s.apply_time_correction(tally, factors, sum_nuclides=False)
+    tcf = np.array([factors[nuc][-1] for nuc in nuclides])
+    assert result.mean.flatten() == pytest.approx(tcf * flux)
+
+    # Make sure summed results match a manual sum
+    result_summed = d1s.apply_time_correction(tally, factors)
+    assert result_summed.mean.flatten()[0] == pytest.approx(result.mean.sum())
+
+    # Make sure various tally methods work
+    result.get_values()
+    result_summed.get_values()
+    result.get_reshaped_data()
+    result_summed.get_reshaped_data()
+    result.get_pandas_dataframe()
+    result_summed.get_pandas_dataframe()