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 @@
 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,11 +21,11 @@ essential aspects of using OpenMC to perform simulations.
     tallies
     plots
     depletion
+    decay_sources
     scripts
     processing
     parallel
     volume
     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 @@ class Polynomial : public Function1D {
   //! \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 @@ class Particle : public ParticleData {
   //! \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 @@ class ParticleData : public GeometryState {
   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 @@ class ParticleData : public GeometryState {
   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 @@ class Reaction {
   //! 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 @@ class ReactionProduct {
   //! \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 @@ class ReactionProduct {
   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,