diff --git a/CMakeLists.txt b/CMakeLists.txt index 35694a14bc8..f76fcb247d9 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/docs/source/conf.py b/docs/source/conf.py index 726269f0b7c..60bd406fba5 100644 --- a/docs/source/conf.py +++ b/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' diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst index d6d04444f2a..bfd94272748 100644 --- a/docs/source/pythonapi/base.rst +++ b/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 diff --git a/docs/source/pythonapi/capi.rst b/docs/source/pythonapi/capi.rst index 9ceff83fde8..c8e0e874dbf 100644 --- a/docs/source/pythonapi/capi.rst +++ b/docs/source/pythonapi/capi.rst @@ -79,6 +79,7 @@ Classes MeshSurfaceFilter MuFilter Nuclide + ParentNuclideFilter ParticleFilter PolarFilter RectilinearMesh diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst index 8731a9a1383..d7a779b1298 100644 --- a/docs/source/pythonapi/deplete.rst +++ b/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 diff --git a/docs/source/usersguide/decay_sources.rst b/docs/source/usersguide/decay_sources.rst new file mode 100644 index 00000000000..e612206969c --- /dev/null +++ b/docs/source/usersguide/decay_sources.rst @@ -0,0 +1,92 @@ +.. usersguide_decay_sources: + +============= +Decay Sources +============= + +Through the :ref:`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 `_ (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 `_ +(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 `. 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) + diff --git a/docs/source/usersguide/index.rst b/docs/source/usersguide/index.rst index 74651c0114b..aef9b1a1c1a 100644 --- a/docs/source/usersguide/index.rst +++ b/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 - \ No newline at end of file diff --git a/include/openmc/chain.h b/include/openmc/chain.h new file mode 100644 index 00000000000..a3bc6f3a364 --- /dev/null +++ b/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 +#include +#include + +#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>& 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> + 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 chain_nuclide_map; +extern vector> chain_nuclides; + +} // namespace data + +//============================================================================== +// Non-member functions +//============================================================================== + +void read_chain_file_xml(); + +} // namespace openmc + +#endif // OPENMC_CHAIN_H diff --git a/include/openmc/endf.h b/include/openmc/endf.h index e580874b4b1..4a737eb8816 100644 --- a/include/openmc/endf.h +++ b/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 coef) : coef_(coef) {} + //! Evaluate the polynomials //! \param[in] x independent variable //! \return Polynomial evaluated at x diff --git a/include/openmc/particle.h b/include/openmc/particle.h index aaac864e11e..0f37719b94f 100644 --- a/include/openmc/particle.h +++ b/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 // diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index e2fac79e372..6b318385ea7 100644 --- a/include/openmc/particle_data.h +++ b/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() diff --git a/include/openmc/reaction.h b/include/openmc/reaction.h index d5f04d136d0..3314d18666f 100644 --- a/include/openmc/reaction.h +++ b/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& temperatures); + //! \param[in] name Name of the nuclide + explicit Reaction( + hid_t group, const vector& temperatures, std::string name); //! Calculate cross section given temperautre/grid index, interpolation factor // diff --git a/include/openmc/reaction_product.h b/include/openmc/reaction_product.h index ce4fa8fc766..4fbbc1b626a 100644 --- a/include/openmc/reaction_product.h +++ b/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 yield_; //!< Yield as a function of energy vector applicability_; //!< Applicability of distribution vector distribution_; //!< Secondary angle-energy distribution + int parent_nuclide_ = -1; //!< Index of chain nuclide that is parent }; } // namespace openmc diff --git a/include/openmc/settings.h b/include/openmc/settings.h index 9a4ce56ec8f..12835fdef9a 100644 --- a/include/openmc/settings.h +++ b/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? diff --git a/include/openmc/tallies/filter.h b/include/openmc/tallies/filter.h index 31dd609ed10..3e982d0cf45 100644 --- a/include/openmc/tallies/filter.h +++ b/include/openmc/tallies/filter.h @@ -36,6 +36,7 @@ enum class FilterType { MESH_SURFACE, MU, MUSURFACE, + PARENT_NUCLIDE, PARTICLE, POLAR, SPHERICAL_HARMONICS, diff --git a/include/openmc/tallies/filter_parent_nuclide.h b/include/openmc/tallies/filter_parent_nuclide.h new file mode 100644 index 00000000000..53f8a5fa41d --- /dev/null +++ b/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 +#include + +#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& bins() const { return bins_; } + void set_bins(span bins); + +protected: + //---------------------------------------------------------------------------- + // Data members + + vector bins_; + vector nuclides_; + + std::unordered_map map_; +}; + +} // namespace openmc +#endif // OPENMC_TALLIES_FILTER_PARENT_NUCLIDE_H diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index 784023f26ff..fc98239a561 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -41,6 +41,63 @@ _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 @@ def __init__( 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) diff --git a/openmc/deplete/d1s.py b/openmc/deplete/d1s.py new file mode 100644 index 00000000000..22fc6ec23da --- /dev/null +++ b/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 diff --git a/openmc/filter.py b/openmc/filter.py index b5926af5ad9..35a7e7c7e96 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -25,7 +25,7 @@ '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 @@ def from_xml_element(cls, elem, **kwargs): 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. diff --git a/openmc/lib/core.py b/openmc/lib/core.py index 8561602e670..577913dcb8f 100644 --- a/openmc/lib/core.py +++ b/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)] diff --git a/openmc/lib/filter.py b/openmc/lib/filter.py index b30f5e66282..b6086c56711 100644 --- a/openmc/lib/filter.py +++ b/openmc/lib/filter.py @@ -21,9 +21,10 @@ '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 @@ class ZernikeRadialFilter(ZernikeFilter): 'meshsurface': MeshSurfaceFilter, 'mu': MuFilter, 'musurface': MuSurfaceFilter, + 'parentnuclide': ParentNuclideFilter, 'particle': ParticleFilter, 'polar': PolarFilter, 'sphericalharmonics': SphericalHarmonicsFilter, diff --git a/openmc/material.py b/openmc/material.py index 343b2fff293..64458f57133 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -308,8 +308,9 @@ def get_decay_photon_energy( 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'}) diff --git a/openmc/settings.py b/openmc/settings.py index 110d57c19f3..0e2a1839907 100644 --- a/openmc/settings.py +++ b/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 @@ def __init__(self, **kwargs): 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 @@ def random_ray(self, random_ray: dict): 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 @@ def _create_ufs_mesh_subelement(self, root, mesh_memo=None): 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 @@ def _random_ray_from_xml_element(self, root): 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 @@ def to_xml_element(self, mesh_memo=None): 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 @@ def from_xml_element(cls, elem, meshes=None): 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 diff --git a/src/chain.cpp b/src/chain.cpp new file mode 100644 index 00000000000..e279d1f5916 --- /dev/null +++ b/src/chain.cpp @@ -0,0 +1,116 @@ +//! \file chain.cpp +//! \brief Depletion chain and associated information + +#include "openmc/chain.h" + +#include // for getenv +#include // for make_unique +#include // for stod + +#include +#include + +#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 chain_nuclide_map; +vector> 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(node)); + } +} + +} // namespace openmc diff --git a/src/initialize.cpp b/src/initialize.cpp index cc1eac9cf35..5f717b137f0 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -12,6 +12,7 @@ #include #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(); diff --git a/src/nuclide.cpp b/src/nuclide.cpp index d78b7a0101d..eaac2f756e7 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -248,7 +248,8 @@ Nuclide::Nuclide(hid_t group, const vector& 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(rx_group, temps_to_read)); + reactions_.push_back( + make_unique(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; diff --git a/src/particle.cpp b/src/particle.cpp index 5b46e2e639b..c63c846741c 100644 --- a/src/particle.cpp +++ b/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(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() diff --git a/src/physics.cpp b/src/physics.cpp index c324ce6b9cf..e04f54c0b13 100644 --- a/src/physics.cpp +++ b/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(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(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_; + } } } diff --git a/src/reaction.cpp b/src/reaction.cpp index 60a14b2b3cc..97147343831 100644 --- a/src/reaction.cpp +++ b/src/reaction.cpp @@ -1,17 +1,20 @@ #include "openmc/reaction.h" +#include // for remove_if #include #include #include // for move #include +#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& temperatures) +Reaction::Reaction( + hid_t group, const vector& 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& 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 diff --git a/src/reaction_product.cpp b/src/reaction_product.cpp index 4cef8d3a958..3ba2c0cfb05 100644 --- a/src/reaction_product.cpp +++ b/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 coef = {product.branching_ratio * photon_per_decay}; + yield_ = make_unique(coef); + + // Set decay photon angle-energy distribution + distribution_.push_back( + make_unique(chain_nuc->photon_energy())); +} + void ReactionProduct::sample( double E_in, double& E_out, double& mu, uint64_t* seed) const { diff --git a/src/settings.cpp b/src/settings.cpp index 60263c84653..e2f5a033b50 100644 --- a/src/settings.cpp +++ b/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() diff --git a/src/tallies/filter.cpp b/src/tallies/filter.cpp index 0ce1ee20b37..6430d9ae100 100644 --- a/src/tallies/filter.cpp +++ b/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(id); } else if (type == "musurface") { return Filter::create(id); + } else if (type == "parentnuclide") { + return Filter::create(id); } else if (type == "particle") { return Filter::create(id); } else if (type == "polar") { diff --git a/src/tallies/filter_parent_nuclide.cpp b/src/tallies/filter_parent_nuclide.cpp new file mode 100644 index 00000000000..d0494107529 --- /dev/null +++ b/src/tallies/filter_parent_nuclide.cpp @@ -0,0 +1,79 @@ +#include "openmc/tallies/filter_parent_nuclide.h" + +#include // for int64_t + +#include + +#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(node, "bins"); + + // Convert nuclides to indices in data::chain_nuclides + std::vector 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 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 diff --git a/tests/chain_ni.xml b/tests/chain_ni.xml new file mode 100644 index 00000000000..fc96d3d8b36 --- /dev/null +++ b/tests/chain_ni.xml @@ -0,0 +1,173 @@ + + + + + + + + + 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 + + + 105210.0 231210.0 1.040592667426837e-17 8.004558980206438e-09 + + + 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 + + + + + + + + + + + + + + + + + + 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 + + + 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 + + + + + + 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 + + + 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 + + + + + + 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 + + + 632864.1 1496834.0 2307600.0 1.3699223176071195e-09 1.1185811981783752e-07 8.491254034754873e-13 + + + 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 + + + + + + + + + + 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 + + + 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 + + + + + + 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 + + + 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 + + + + + + 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 + + + 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 + + + + + + 931100.0 1346100.0 0.11552453009332422 0.23104906018664845 + + + 5029800.0 5960900.0 7307000.0 0.11552453009332422 0.11552453009332422 2.079441541679836 + + + + + + 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 + + + 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 + + + 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 + + + + + + + + + + + + 1072500.0 2.890064045563861e-13 + + + 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 + + + 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 + + + + + + + + + + + + + + + + + + + + + + + + + + + + 66945.0 2.1704054025041888e-10 + + + + + + + + + + + + 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 + + + 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 + + + diff --git a/tests/regression_tests/surface_source/test.py b/tests/regression_tests/surface_source/test.py index 251e9e9ad4a..83b78e61deb 100644 --- a/tests/regression_tests/surface_source/test.py +++ b/tests/regression_tests/surface_source/test.py @@ -10,6 +10,17 @@ 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 @@ def _compare_output(self): """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.""" @@ -132,4 +137,4 @@ def test_surface_source_read(model): harness = SurfaceSourceTestHarness('statepoint.10.h5', model, 'inputs_true_read.dat') - harness.main() \ No newline at end of file + harness.main() diff --git a/tests/unit_tests/test_d1s.py b/tests/unit_tests/test_d1s.py new file mode 100644 index 00000000000..7c92dd6d9b8 --- /dev/null +++ b/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()