From 074354e29dab9b718cda50d73268dc4c85450059 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 4 Mar 2026 23:41:11 -0600 Subject: [PATCH 01/37] Add MPI parallelism to photon source generation --- openmc/deplete/r2s.py | 97 ++++++++++++++++++++++++++----------------- 1 file changed, 60 insertions(+), 37 deletions(-) diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index 10d7fdd502a..b33576a9c39 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -9,6 +9,7 @@ import openmc from . import IndependentOperator, PredictorIntegrator from .microxs import get_microxs_and_flux, write_microxs_hdf5, read_microxs_hdf5 +from .pool import _distribute from .results import Results from ..checkvalue import PathLike from ..mpi import comm @@ -516,29 +517,39 @@ def step3_photon_transport( if different_photon_model: photon_cells = self.photon_model.geometry.get_all_cells() + # For cell-based calculations, pre-compute and distribute the eligible + # (cell, original_mat) pairs across MPI ranks once, outside the time + # loop, since the eligibility check doesn't depend on time. + if self.method == 'cell-based': + _domain_pairs = [] + for cell, original_mat in zip( + self.domains, self.results['activation_materials']): + if different_photon_model: + if cell.id not in photon_cells or \ + cell.fill.id != photon_cells[cell.id].fill.id: + continue + _domain_pairs.append((cell, original_mat)) + _my_domain_pairs = _distribute(_domain_pairs) + for time_index in time_indices: # Create decay photon source if self.method == 'mesh-based': self.photon_model.settings.source = \ self.get_decay_photon_source_mesh(time_index) else: - sources = [] results = self.results['depletion_results'] - for cell, original_mat in zip(self.domains, self.results['activation_materials']): - # Skip if the cell is not in the photon model or the - # material has changed - if different_photon_model: - if cell.id not in photon_cells or \ - cell.fill.id != photon_cells[cell.id].fill.id: - continue + # Build sources for this rank's assigned cells + local_sources = [] + for cell, original_mat in _my_domain_pairs: # Get bounding box for the cell bounding_box = bounding_boxes[cell.id] # Get activated material composition - activated_mat = results[time_index].get_material(str(original_mat.id)) + activated_mat = results[time_index].get_material( + str(original_mat.id)) - # Create decay photon source source + # Create decay photon source space = openmc.stats.Box(*bounding_box) energy = activated_mat.get_decay_photon_energy() strength = energy.integral() if energy is not None else 0.0 @@ -549,8 +560,13 @@ def step3_photon_transport( strength=strength, constraints={'domains': [cell]} ) - sources.append(source) - self.photon_model.settings.source = sources + local_sources.append(source) + + # Gather sources from all ranks and flatten into a single list + all_sources = comm.allgather(local_sources) + self.photon_model.settings.source = [ + s for rank_sources in all_sources for s in rank_sources + ] # Convert time_index (which may be negative) to a normal index if time_index < 0: @@ -598,18 +614,18 @@ def get_decay_photon_source_mesh( """ mat_dict = self.neutron_model._get_all_materials() - # List to hold all sources - sources = [] - - # Index in the overall list of activated materials - index_mat = 0 - # Get various results from previous steps mmv_list = self.results['mesh_material_volumes'] materials = self.results['activation_materials'] results = self.results['depletion_results'] photon_mmv_list = self.results.get('mesh_material_volumes_photon') + # Phase 1: Enumerate all processable work items across all mesh + # elements and materials. This pass is cheap (no depletion data access) + # and records the stable index_mat needed to look up each activated + # material in the distributed phase. + work_items = [] # list of (index_mat, mat_id, bbox) + index_mat = 0 for mesh_idx, mat_vols in enumerate(mmv_list): photon_mat_vols = photon_mmv_list[mesh_idx] \ if photon_mmv_list is not None else None @@ -636,27 +652,34 @@ def get_decay_photon_source_mesh( index_mat += 1 continue - # Get activated material composition - original_mat = materials[index_mat] - activated_mat = results[time_index].get_material(str(original_mat.id)) - - # Create decay photon source - energy = activated_mat.get_decay_photon_energy() - if energy is not None: - strength = energy.integral() - space = openmc.stats.Box(*bbox) - sources.append(openmc.IndependentSource( - space=space, - energy=energy, - particle='photon', - strength=strength, - constraints={'domains': [mat_dict[mat_id]]} - )) - - # Increment index of activated material + work_items.append((index_mat, mat_id, bbox)) index_mat += 1 - return sources + # Phase 2: Distribute work items across MPI ranks + my_items = _distribute(work_items) + + # Phase 3: Each rank builds decay photon sources for its assigned items + local_sources = [] + for item_index_mat, mat_id, bbox in my_items: + original_mat = materials[item_index_mat] + activated_mat = results[time_index].get_material(str(original_mat.id)) + + # Create decay photon source + energy = activated_mat.get_decay_photon_energy() + if energy is not None: + strength = energy.integral() + space = openmc.stats.Box(*bbox) + local_sources.append(openmc.IndependentSource( + space=space, + energy=energy, + particle='photon', + strength=strength, + constraints={'domains': [mat_dict[mat_id]]} + )) + + # Phase 4: Gather sources from all ranks and return the complete list + all_sources = comm.allgather(local_sources) + return [s for rank_sources in all_sources for s in rank_sources] def load_results(self, path: PathLike): """Load results from a previous R2S calculation. From e89d74be19f0db607fef9b9ea353770461254cbe Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 6 Mar 2026 08:50:19 -0600 Subject: [PATCH 02/37] Create decay photon source in C++ for R2S --- CMakeLists.txt | 1 + include/openmc/capi.h | 20 +++ include/openmc/r2s_source.h | 40 ++++++ include/openmc/source.h | 13 ++ openmc/deplete/r2s.py | 200 +++++++++++++++++++++++------- openmc/lib/core.py | 77 ++++++++++++ src/r2s_source.cpp | 237 ++++++++++++++++++++++++++++++++++++ 7 files changed, 547 insertions(+), 41 deletions(-) create mode 100644 include/openmc/r2s_source.h create mode 100644 src/r2s_source.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 9fe133a22e3..62bffc8110d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -415,6 +415,7 @@ list(APPEND libopenmc_SOURCES src/ray.cpp src/reaction.cpp src/reaction_product.cpp + src/r2s_source.cpp src/scattdata.cpp src/secondary_correlated.cpp src/secondary_kalbach.cpp diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 911654d318f..e2339719fe2 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -320,6 +320,26 @@ int openmc_properties_export(const char* filename); // \return Error code int openmc_properties_import(const char* filename); +//! Create decay photon sources for R2S calculations +//! +//! Clears existing external sources and creates IndependentSource objects +//! for each spatial region based on the activated nuclide compositions. +//! +//! \param n_regions Number of spatial regions +//! \param domain_ids Material or cell IDs for constraints [n_regions] +//! \param domain_type "material" or "cell" +//! \param lower_left Flattened lower-left bbox corners [n_regions * 3] +//! \param upper_right Flattened upper-right bbox corners [n_regions * 3] +//! \param n_nuclides Number of nuclides +//! \param nuclide_names Nuclide name strings [n_nuclides] +//! \param atom_densities Atom densities [n_regions * n_nuclides] row-major +//! \param volumes Region volumes in cm^3 [n_regions] +//! \return Error code (0 on success) +int openmc_create_decay_photon_sources(int n_regions, const int32_t* domain_ids, + const char* domain_type, const double* lower_left, const double* upper_right, + int n_nuclides, const char** nuclide_names, const double* atom_densities, + const double* volumes); + // Error codes extern int OPENMC_E_UNASSIGNED; extern int OPENMC_E_ALLOCATE; diff --git a/include/openmc/r2s_source.h b/include/openmc/r2s_source.h new file mode 100644 index 00000000000..6da76222af3 --- /dev/null +++ b/include/openmc/r2s_source.h @@ -0,0 +1,40 @@ +//! \file r2s_source.h +//! \brief Decay photon source generation for R2S calculations + +#ifndef OPENMC_R2S_SOURCE_H +#define OPENMC_R2S_SOURCE_H + +#include // for size_t +#include // for int32_t + +#include "openmc/source.h" + +namespace openmc { + +//! Create decay photon IndependentSource objects from activated materials +//! +//! For each spatial region, this function combines the decay photon energy +//! spectra from all nuclides present (weighted by atom count and decay +//! constant) and creates an IndependentSource with a SpatialBox spatial +//! distribution and isotropic angular distribution. The resulting sources +//! are stored in model::external_sources with the probability mass function +//! rebuilt accordingly. +//! +//! \param n_regions Number of spatial regions +//! \param domain_ids Material or cell IDs for source constraints [n_regions] +//! \param domain_type MATERIAL or CELL +//! \param lower_left Flattened lower-left bbox corners [n_regions * 3] +//! \param upper_right Flattened upper-right bbox corners [n_regions * 3] +//! \param n_nuclides Number of nuclides in the atom density matrix +//! \param nuclide_names Array of nuclide name strings [n_nuclides] +//! \param atom_densities Row-major atom densities in [atom/b-cm], +//! shape [n_regions * n_nuclides] +//! \param volumes Volume of each region in [cm^3], shape [n_regions] +void create_decay_photon_sources(int n_regions, const int32_t* domain_ids, + Source::DomainType domain_type, const double* lower_left, + const double* upper_right, int n_nuclides, const char** nuclide_names, + const double* atom_densities, const double* volumes); + +} // namespace openmc + +#endif // OPENMC_R2S_SOURCE_H diff --git a/include/openmc/source.h b/include/openmc/source.h index 1ef7eba2b24..7f6c36ec0b9 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -73,6 +73,16 @@ class Source { // Methods that can be overridden virtual double strength() const { return strength_; } + //! Set the source strength + void set_strength(double strength) { strength_ = strength; } + + //! Set the domain type and IDs for source constraints + void set_domain(DomainType domain_type, std::unordered_set ids) + { + domain_type_ = domain_type; + domain_ids_ = std::move(ids); + } + //! Sample a source site and apply constraints // //! \param[inout] seed Pseudorandom seed pointer @@ -138,6 +148,9 @@ class IndependentSource : public Source { // Properties ParticleType particle_type() const { return particle_; } + //! Set the particle type + void set_particle(ParticleType p) { particle_ = p; } + // Make observing pointers available SpatialDistribution* space() const { return space_.get(); } UnitSphereDistribution* angle() const { return angle_.get(); } diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index b33576a9c39..0adb6d46561 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -7,6 +7,7 @@ import numpy as np import openmc +import openmc.lib from . import IndependentOperator, PredictorIntegrator from .microxs import get_microxs_and_flux, write_microxs_hdf5, read_microxs_hdf5 from .pool import _distribute @@ -517,10 +518,12 @@ def step3_photon_transport( if different_photon_model: photon_cells = self.photon_model.geometry.get_all_cells() - # For cell-based calculations, pre-compute and distribute the eligible - # (cell, original_mat) pairs across MPI ranks once, outside the time - # loop, since the eligibility check doesn't depend on time. - if self.method == 'cell-based': + # Determine eligible work items upfront. For cell-based calculations, + # pre-compute the eligible (cell, original_mat) pairs since the + # eligibility check doesn't depend on time. + if self.method == 'mesh-based': + work_items = self._get_mesh_work_items() + else: _domain_pairs = [] for cell, original_mat in zip( self.domains, self.results['activation_materials']): @@ -529,45 +532,11 @@ def step3_photon_transport( cell.fill.id != photon_cells[cell.id].fill.id: continue _domain_pairs.append((cell, original_mat)) - _my_domain_pairs = _distribute(_domain_pairs) - - for time_index in time_indices: - # Create decay photon source - if self.method == 'mesh-based': - self.photon_model.settings.source = \ - self.get_decay_photon_source_mesh(time_index) - else: - results = self.results['depletion_results'] - - # Build sources for this rank's assigned cells - local_sources = [] - for cell, original_mat in _my_domain_pairs: - # Get bounding box for the cell - bounding_box = bounding_boxes[cell.id] - - # Get activated material composition - activated_mat = results[time_index].get_material( - str(original_mat.id)) - - # Create decay photon source - space = openmc.stats.Box(*bounding_box) - energy = activated_mat.get_decay_photon_energy() - strength = energy.integral() if energy is not None else 0.0 - source = openmc.IndependentSource( - space=space, - energy=energy, - particle='photon', - strength=strength, - constraints={'domains': [cell]} - ) - local_sources.append(source) - # Gather sources from all ranks and flatten into a single list - all_sources = comm.allgather(local_sources) - self.photon_model.settings.source = [ - s for rank_sources in all_sources for s in rank_sources - ] + # Ensure photon transport is enabled in settings + self.photon_model.settings.photon_transport = True + for time_index in time_indices: # Convert time_index (which may be negative) to a normal index if time_index < 0: time_index = len(self.results['depletion_results']) + time_index @@ -575,6 +544,13 @@ def step3_photon_transport( # Run photon transport calculation photon_dir = Path(output_dir) / f'time_{time_index}' with TemporarySession(self.photon_model, cwd=photon_dir): + # Create decay photon sources in C++ memory + if self.method == 'mesh-based': + self._create_cpp_sources_mesh(time_index, work_items) + else: + self._create_cpp_sources_cells( + time_index, bounding_boxes, _domain_pairs) + statepoint_path = self.photon_model.run(**run_kwargs) # Store tally results @@ -681,6 +657,148 @@ def get_decay_photon_source_mesh( all_sources = comm.allgather(local_sources) return [s for rank_sources in all_sources for s in rank_sources] + def _get_mesh_work_items(self): + """Enumerate mesh-based work items across all meshes. + + Returns a list of (index_mat, mat_id, bbox) tuples for each eligible + mesh element--material combination. This is the same enumeration used + by :meth:`get_decay_photon_source_mesh` but separated out so it can be + called once outside the time loop. + + Returns + ------- + list of tuple + Each tuple is (index_mat, mat_id, bbox). + """ + mmv_list = self.results['mesh_material_volumes'] + photon_mmv_list = self.results.get('mesh_material_volumes_photon') + + work_items = [] + index_mat = 0 + for mesh_idx, mat_vols in enumerate(mmv_list): + photon_mat_vols = photon_mmv_list[mesh_idx] \ + if photon_mmv_list is not None else None + + n_elements = mat_vols.num_elements + for index_elem in range(n_elements): + if photon_mat_vols is not None: + photon_materials = { + mat_id + for mat_id, _ in photon_mat_vols.by_element(index_elem) + if mat_id is not None + } + + for mat_id, _, bbox in mat_vols.by_element( + index_elem, include_bboxes=True): + if mat_id is None: + continue + if photon_mat_vols is not None \ + and mat_id not in photon_materials: + index_mat += 1 + continue + work_items.append((index_mat, mat_id, bbox)) + index_mat += 1 + + return work_items + + def _get_region_data(self, time_index, work_items, domain_type): + """Extract region data as flat numpy arrays for C++ source creation. + + Parameters + ---------- + time_index : int + Index into depletion results for the desired time. + work_items : list of tuple + For mesh-based: list of (index_mat, mat_id, bbox). + For cell-based: list of (cell, original_mat, bbox). + domain_type : str + Either ``'material'`` or ``'cell'``. + + Returns + ------- + domain_ids : numpy.ndarray of int32 + lower_left : numpy.ndarray of float64, shape (n, 3) + upper_right : numpy.ndarray of float64, shape (n, 3) + nuclide_names : list of str + atom_densities : numpy.ndarray of float64, shape (n, n_nuclides) + volumes : numpy.ndarray of float64, shape (n,) + """ + step_result = self.results['depletion_results'][time_index] + materials = self.results['activation_materials'] + + nuclide_names = list(step_result.index_nuc.keys()) + n_nuclides = len(nuclide_names) + n_regions = len(work_items) + + domain_ids = np.empty(n_regions, dtype=np.int32) + lower_left = np.empty((n_regions, 3), dtype=np.float64) + upper_right = np.empty((n_regions, 3), dtype=np.float64) + atom_densities = np.empty((n_regions, n_nuclides), dtype=np.float64) + volumes_arr = np.empty(n_regions, dtype=np.float64) + + for i, item in enumerate(work_items): + if domain_type == 'material': + index_mat, mat_id, bbox = item + domain_ids[i] = mat_id + original_mat = materials[index_mat] + mat_key = str(original_mat.id) + vol = step_result.volume[mat_key] + else: + cell, original_mat, bbox = item + domain_ids[i] = cell.id + mat_key = str(original_mat.id) + vol = step_result.volume[mat_key] + + lower_left[i] = bbox[0] + upper_right[i] = bbox[1] + volumes_arr[i] = vol + + # Get atom counts for this material and convert to atom/b-cm + mat_idx = step_result.index_mat[mat_key] + atom_counts = step_result.data[mat_idx, :] + atom_densities[i, :] = atom_counts / vol * 1e-24 + + return (domain_ids, lower_left, upper_right, nuclide_names, + atom_densities, volumes_arr) + + def _create_cpp_sources_mesh(self, time_index, work_items): + """Create decay photon sources in C++ for mesh-based calculations. + + Parameters + ---------- + time_index : int + Index into depletion results. + work_items : list of tuple + List of (index_mat, mat_id, bbox) from :meth:`_get_mesh_work_items`. + """ + domain_ids, ll, ur, nuc_names, densities, vols = \ + self._get_region_data(time_index, work_items, 'material') + openmc.lib.create_decay_photon_sources( + domain_ids, 'material', ll, ur, nuc_names, densities, vols) + + def _create_cpp_sources_cells(self, time_index, bounding_boxes, + domain_pairs): + """Create decay photon sources in C++ for cell-based calculations. + + Parameters + ---------- + time_index : int + Index into depletion results. + bounding_boxes : dict + Mapping of cell ID to :class:`~openmc.BoundingBox`. + domain_pairs : list of tuple + List of (cell, original_mat) pairs. + """ + # Build work items with bounding boxes + work_items = [ + (cell, original_mat, bounding_boxes[cell.id]) + for cell, original_mat in domain_pairs + ] + domain_ids, ll, ur, nuc_names, densities, vols = \ + self._get_region_data(time_index, work_items, 'cell') + openmc.lib.create_decay_photon_sources( + domain_ids, 'cell', ll, ur, nuc_names, densities, vols) + def load_results(self, path: PathLike): """Load results from a previous R2S calculation. diff --git a/openmc/lib/core.py b/openmc/lib/core.py index 00876db9972..c25c324411d 100644 --- a/openmc/lib/core.py +++ b/openmc/lib/core.py @@ -112,6 +112,19 @@ class _SourceSite(Structure): _dll.openmc_sample_external_source.argtypes = [c_size_t, POINTER(c_uint64), POINTER(_SourceSite)] _dll.openmc_sample_external_source.restype = c_int _dll.openmc_sample_external_source.errcheck = _error_handler +_dll.openmc_create_decay_photon_sources.argtypes = [ + c_int, # n_regions + POINTER(c_int32), # domain_ids + c_char_p, # domain_type + POINTER(c_double), # lower_left + POINTER(c_double), # upper_right + c_int, # n_nuclides + POINTER(c_char_p), # nuclide_names + POINTER(c_double), # atom_densities + POINTER(c_double), # volumes +] +_dll.openmc_create_decay_photon_sources.restype = c_int +_dll.openmc_create_decay_photon_sources.errcheck = _error_handler def global_bounding_box(): """Calculate a global bounding box for the model""" @@ -491,6 +504,70 @@ def run_random_ray(output=True): with quiet_dll(output): _dll.openmc_run_random_ray() +def create_decay_photon_sources( + domain_ids: np.ndarray, + domain_type: str, + lower_left: np.ndarray, + upper_right: np.ndarray, + nuclide_names: list[str], + atom_densities: np.ndarray, + volumes: np.ndarray, +): + """Create decay photon sources from activated material compositions. + + This populates the internal external source bank with + :class:`~openmc.IndependentSource` objects for each activated region with + non-zero photon emission. Existing external sources are cleared before the + new sources are created. + + .. versionadded:: 0.15.1 + + Parameters + ---------- + domain_ids : numpy.ndarray of int32 + Material or cell IDs for source domain constraints, shape (n_regions,). + domain_type : str + Either "material" or "cell". + lower_left : numpy.ndarray of float64 + Lower-left bounding box corners, shape (n_regions, 3). + upper_right : numpy.ndarray of float64 + Upper-right bounding box corners, shape (n_regions, 3). + nuclide_names : list of str + Names of nuclides in the atom density matrix. + atom_densities : numpy.ndarray of float64 + Atom densities in [atom/b-cm], shape (n_regions, n_nuclides). + volumes : numpy.ndarray of float64 + Volume of each region in [cm^3], shape (n_regions,). + + """ + n_regions = len(domain_ids) + n_nuclides = len(nuclide_names) + + # Ensure arrays are contiguous with correct dtype + domain_ids = np.ascontiguousarray(domain_ids, dtype=np.int32) + lower_left = np.ascontiguousarray(lower_left.ravel(), dtype=np.float64) + upper_right = np.ascontiguousarray(upper_right.ravel(), dtype=np.float64) + atom_densities = np.ascontiguousarray(atom_densities.ravel(), dtype=np.float64) + volumes = np.ascontiguousarray(volumes, dtype=np.float64) + + # Convert nuclide names to array of c_char_p + names_arr = (c_char_p * n_nuclides)( + *(name.encode() for name in nuclide_names) + ) + + _dll.openmc_create_decay_photon_sources( + n_regions, + domain_ids.ctypes.data_as(POINTER(c_int32)), + domain_type.encode(), + lower_left.ctypes.data_as(POINTER(c_double)), + upper_right.ctypes.data_as(POINTER(c_double)), + n_nuclides, + names_arr, + atom_densities.ctypes.data_as(POINTER(c_double)), + volumes.ctypes.data_as(POINTER(c_double)), + ) + + def sample_external_source( n_samples: int = 1000, prn_seed: int | None = None, diff --git a/src/r2s_source.cpp b/src/r2s_source.cpp new file mode 100644 index 00000000000..ea7e7c08598 --- /dev/null +++ b/src/r2s_source.cpp @@ -0,0 +1,237 @@ +//! \file r2s_source.cpp +//! \brief Decay photon source generation for R2S calculations + +#include "openmc/r2s_source.h" + +#include // for sort +#include // for accumulate +#include + +#ifdef _OPENMP +#include +#endif + +#include "openmc/capi.h" +#include "openmc/chain.h" +#include "openmc/distribution.h" +#include "openmc/distribution_multi.h" +#include "openmc/distribution_spatial.h" +#include "openmc/error.h" +#include "openmc/memory.h" +#include "openmc/particle_type.h" +#include "openmc/settings.h" +#include "openmc/source.h" +#include "openmc/vector.h" + +namespace openmc { + +//============================================================================== +// DecayPhotonMixture — non-owning mixture of decay photon distributions +//============================================================================== + +//! Energy distribution formed by mixing multiple decay photon spectra. +//! +//! Unlike the general Mixture distribution, this class holds non-owning +//! pointers to the component distributions (which live in +//! data::chain_nuclides). Each component is weighted by the activity +//! (atoms * decay_constant) of the corresponding nuclide. + +class DecayPhotonMixture : public Distribution { +public: + //! Construct from non-owning distribution pointers and weights + //! + //! \param dists Non-owning pointers to component distributions + //! \param weights Activity-based weights for each component. The Mixture + //! probability for component i is weights[i] * dists[i]->integral() + //! (the integral encodes photons-per-decay). + DecayPhotonMixture(vector dists, vector weights) + : dists_(std::move(dists)) + { + vector probs; + probs.reserve(dists_.size()); + for (size_t i = 0; i < dists_.size(); ++i) { + probs.push_back(weights[i] * dists_[i]->integral()); + } + integral_ = std::accumulate(probs.begin(), probs.end(), 0.0); + di_.assign(probs); + } + + std::pair sample(uint64_t* seed) const override + { + size_t idx = di_.sample(seed); + return dists_[idx]->sample(seed); + } + + double integral() const override { return integral_; } + +protected: + double sample_unbiased(uint64_t* seed) const override + { + size_t idx = di_.sample(seed); + return dists_[idx]->sample(seed).first; + } + +private: + vector dists_; //!< Non-owning component distributions + DiscreteIndex di_; //!< Discrete index for component selection + double integral_; //!< Total photon emission rate +}; + +//============================================================================== +// create_decay_photon_sources implementation +//============================================================================== + +void create_decay_photon_sources(int n_regions, const int32_t* domain_ids, + Source::DomainType domain_type, const double* lower_left, + const double* upper_right, int n_nuclides, const char** nuclide_names, + const double* atom_densities, const double* volumes) +{ + // Pre-resolve nuclide names to chain nuclide indices. A value of -1 + // indicates the nuclide is not in the chain (or has no decay photon data). + vector chain_indices(n_nuclides, -1); + for (int j = 0; j < n_nuclides; ++j) { + auto it = data::chain_nuclide_map.find(nuclide_names[j]); + if (it != data::chain_nuclide_map.end()) { + chain_indices[j] = it->second; + } + } + + // Clear existing sources + model::external_sources.clear(); + + // Thread-local storage for sources built in parallel + vector>> thread_sources; + +#ifdef _OPENMP + int n_threads = omp_get_max_threads(); +#else + int n_threads = 1; +#endif + thread_sources.resize(n_threads); + +#pragma omp parallel + { +#ifdef _OPENMP + int tid = omp_get_thread_num(); +#else + int tid = 0; +#endif + +#pragma omp for schedule(dynamic, 64) + for (int i = 0; i < n_regions; ++i) { + // Collect distributions and weights for nuclides present in this region + vector dists; + vector weights; + + for (int j = 0; j < n_nuclides; ++j) { + double atom_dens = atom_densities[i * n_nuclides + j]; + if (atom_dens <= 0.0) + continue; + + int chain_idx = chain_indices[j]; + if (chain_idx < 0) + continue; + + const auto& chain_nuc = data::chain_nuclides[chain_idx]; + const Distribution* photon_dist = chain_nuc->photon_energy(); + if (!photon_dist) + continue; + + // Weight = number of atoms * decay_constant + // atom_dens is in [atom/b-cm], volumes in [cm^3] + // atoms = atom_dens * 1e24 * volume + double atoms = atom_dens * 1.0e24 * volumes[i]; + double activity = atoms * chain_nuc->decay_constant(); + + dists.push_back(photon_dist); + weights.push_back(activity); + } + + // Skip regions with no photon-emitting nuclides + if (dists.empty()) + continue; + + // Build combined energy distribution + auto energy = + make_unique(std::move(dists), std::move(weights)); + double strength = energy->integral(); + if (strength <= 0.0) + continue; + + // Build spatial distribution from bounding box + Position ll { + lower_left[i * 3], lower_left[i * 3 + 1], lower_left[i * 3 + 2]}; + Position ur { + upper_right[i * 3], upper_right[i * 3 + 1], upper_right[i * 3 + 2]}; + + // Build time distribution (delta at t=0) + double t0[] {0.0}; + double tp[] {1.0}; + + // Create IndependentSource + auto source = make_unique( + UPtrSpace {new SpatialBox(ll, ur)}, UPtrAngle {new Isotropic()}, + std::move(energy), UPtrDist {new Discrete(t0, tp, 1)}); + source->set_particle(ParticleType::photon()); + source->set_strength(strength); + source->set_domain(domain_type, {domain_ids[i]}); + + thread_sources[tid].push_back(std::move(source)); + } + } // end parallel + + // Merge thread-local sources into model::external_sources + for (auto& tvec : thread_sources) { + for (auto& src : tvec) { + model::external_sources.push_back(std::move(src)); + } + } + + // Enable photon transport since we are creating photon sources + settings::photon_transport = true; + + // Rebuild probability mass function for sampling + vector source_strengths; + source_strengths.reserve(model::external_sources.size()); + for (auto& s : model::external_sources) { + source_strengths.push_back(s->strength()); + } + model::external_sources_probability.assign(source_strengths); +} + +} // namespace openmc + +//============================================================================== +// C API +//============================================================================== + +extern "C" int openmc_create_decay_photon_sources(int n_regions, + const int32_t* domain_ids, const char* domain_type, const double* lower_left, + const double* upper_right, int n_nuclides, const char** nuclide_names, + const double* atom_densities, const double* volumes) +{ + using DT = openmc::Source::DomainType; + + // Validate domain type + std::string dt_str(domain_type); + DT dt; + if (dt_str == "material") { + dt = DT::MATERIAL; + } else if (dt_str == "cell") { + dt = DT::CELL; + } else { + openmc::set_errmsg( + "Invalid domain_type '" + dt_str + "'. Must be 'material' or 'cell'."); + return OPENMC_E_INVALID_ARGUMENT; + } + + try { + openmc::create_decay_photon_sources(n_regions, domain_ids, dt, lower_left, + upper_right, n_nuclides, nuclide_names, atom_densities, volumes); + } catch (const std::exception& e) { + openmc::set_errmsg(e.what()); + return OPENMC_E_DATA; + } + + return 0; +} From 71a38c8b8886bcacd5152084b436fefa31edc405 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 7 Mar 2026 11:10:37 -0600 Subject: [PATCH 03/37] Use functions from openmp_interface.h --- src/r2s_source.cpp | 17 +++-------------- 1 file changed, 3 insertions(+), 14 deletions(-) diff --git a/src/r2s_source.cpp b/src/r2s_source.cpp index ea7e7c08598..9ef938f8a0d 100644 --- a/src/r2s_source.cpp +++ b/src/r2s_source.cpp @@ -7,10 +7,6 @@ #include // for accumulate #include -#ifdef _OPENMP -#include -#endif - #include "openmc/capi.h" #include "openmc/chain.h" #include "openmc/distribution.h" @@ -18,6 +14,7 @@ #include "openmc/distribution_spatial.h" #include "openmc/error.h" #include "openmc/memory.h" +#include "openmc/openmp_interface.h" #include "openmc/particle_type.h" #include "openmc/settings.h" #include "openmc/source.h" @@ -102,20 +99,12 @@ void create_decay_photon_sources(int n_regions, const int32_t* domain_ids, // Thread-local storage for sources built in parallel vector>> thread_sources; -#ifdef _OPENMP - int n_threads = omp_get_max_threads(); -#else - int n_threads = 1; -#endif + int n_threads = num_threads(); thread_sources.resize(n_threads); #pragma omp parallel { -#ifdef _OPENMP - int tid = omp_get_thread_num(); -#else - int tid = 0; -#endif + int tid = thread_num(); #pragma omp for schedule(dynamic, 64) for (int i = 0; i < n_regions; ++i) { From 4a4d1787ca24eb3dc59f29bba337134a43f74525 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 7 Mar 2026 11:14:25 -0600 Subject: [PATCH 04/37] Remove unused get_decay_photon_source_mesh function --- openmc/deplete/r2s.py | 105 ++---------------------------------------- 1 file changed, 3 insertions(+), 102 deletions(-) diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index 0adb6d46561..2da1629460c 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -10,7 +10,6 @@ import openmc.lib from . import IndependentOperator, PredictorIntegrator from .microxs import get_microxs_and_flux, write_microxs_hdf5, read_microxs_hdf5 -from .pool import _distribute from .results import Results from ..checkvalue import PathLike from ..mpi import comm @@ -559,111 +558,13 @@ def step3_photon_transport( sp.tallies[tally.id] for tally in self.photon_model.tallies ] - def get_decay_photon_source_mesh( - self, - time_index: int = -1 - ) -> list[openmc.IndependentSource]: - """Create decay photon source for a mesh-based calculation. - - For each mesh element-material combination across all meshes, an - :class:`~openmc.IndependentSource` is created with a - :class:`~openmc.stats.Box` spatial distribution based on the bounding - box of the material within the mesh element. A material constraint is - also applied so that sampled source sites are limited to the correct - region. - - When the photon transport model is different from the neutron model, the - photon MeshMaterialVolumes is used to determine whether an (element, - material) combination exists in the photon model. - - Parameters - ---------- - time_index : int, optional - Time index for the decay photon source. Default is -1 (last time). - - Returns - ------- - list of openmc.IndependentSource - A list of IndependentSource objects for the decay photons, one for - each mesh element-material combination with non-zero source strength. - - """ - mat_dict = self.neutron_model._get_all_materials() - - # Get various results from previous steps - mmv_list = self.results['mesh_material_volumes'] - materials = self.results['activation_materials'] - results = self.results['depletion_results'] - photon_mmv_list = self.results.get('mesh_material_volumes_photon') - - # Phase 1: Enumerate all processable work items across all mesh - # elements and materials. This pass is cheap (no depletion data access) - # and records the stable index_mat needed to look up each activated - # material in the distributed phase. - work_items = [] # list of (index_mat, mat_id, bbox) - index_mat = 0 - for mesh_idx, mat_vols in enumerate(mmv_list): - photon_mat_vols = photon_mmv_list[mesh_idx] \ - if photon_mmv_list is not None else None - - # Total number of mesh elements for this mesh - n_elements = mat_vols.num_elements - - for index_elem in range(n_elements): - # Determine which materials exist in the photon model for this element - if photon_mat_vols is not None: - photon_materials = { - mat_id - for mat_id, _ in photon_mat_vols.by_element(index_elem) - if mat_id is not None - } - - for mat_id, _, bbox in mat_vols.by_element(index_elem, include_bboxes=True): - # Skip void volume - if mat_id is None: - continue - - # Skip if this material doesn't exist in photon model - if photon_mat_vols is not None and mat_id not in photon_materials: - index_mat += 1 - continue - - work_items.append((index_mat, mat_id, bbox)) - index_mat += 1 - - # Phase 2: Distribute work items across MPI ranks - my_items = _distribute(work_items) - - # Phase 3: Each rank builds decay photon sources for its assigned items - local_sources = [] - for item_index_mat, mat_id, bbox in my_items: - original_mat = materials[item_index_mat] - activated_mat = results[time_index].get_material(str(original_mat.id)) - - # Create decay photon source - energy = activated_mat.get_decay_photon_energy() - if energy is not None: - strength = energy.integral() - space = openmc.stats.Box(*bbox) - local_sources.append(openmc.IndependentSource( - space=space, - energy=energy, - particle='photon', - strength=strength, - constraints={'domains': [mat_dict[mat_id]]} - )) - - # Phase 4: Gather sources from all ranks and return the complete list - all_sources = comm.allgather(local_sources) - return [s for rank_sources in all_sources for s in rank_sources] - def _get_mesh_work_items(self): """Enumerate mesh-based work items across all meshes. Returns a list of (index_mat, mat_id, bbox) tuples for each eligible - mesh element--material combination. This is the same enumeration used - by :meth:`get_decay_photon_source_mesh` but separated out so it can be - called once outside the time loop. + mesh element--material combination, where index_mat is the index into + the activation materials list, mat_id is the material ID, and bbox is + the bounding box for that mesh element--material combination. Returns ------- From 85d21dfc28e2a8587e996ca1b4ec9fd6f889bb81 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 7 Mar 2026 11:37:57 -0600 Subject: [PATCH 05/37] Single _create_cpp_sources method --- openmc/deplete/r2s.py | 45 ++++++++++++++----------------------------- 1 file changed, 14 insertions(+), 31 deletions(-) diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index 2da1629460c..282266a30a3 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -545,10 +545,13 @@ def step3_photon_transport( with TemporarySession(self.photon_model, cwd=photon_dir): # Create decay photon sources in C++ memory if self.method == 'mesh-based': - self._create_cpp_sources_mesh(time_index, work_items) + self._create_cpp_sources(time_index, work_items, 'material') else: - self._create_cpp_sources_cells( - time_index, bounding_boxes, _domain_pairs) + cell_work_items = [ + (cell, original_mat, bounding_boxes[cell.id]) + for cell, original_mat in _domain_pairs + ] + self._create_cpp_sources(time_index, cell_work_items, 'cell') statepoint_path = self.photon_model.run(**run_kwargs) @@ -662,43 +665,23 @@ def _get_region_data(self, time_index, work_items, domain_type): return (domain_ids, lower_left, upper_right, nuclide_names, atom_densities, volumes_arr) - def _create_cpp_sources_mesh(self, time_index, work_items): - """Create decay photon sources in C++ for mesh-based calculations. + def _create_cpp_sources(self, time_index, work_items, domain_type): + """Create decay photon sources in C++ for a set of regions. Parameters ---------- time_index : int Index into depletion results. work_items : list of tuple - List of (index_mat, mat_id, bbox) from :meth:`_get_mesh_work_items`. - """ - domain_ids, ll, ur, nuc_names, densities, vols = \ - self._get_region_data(time_index, work_items, 'material') - openmc.lib.create_decay_photon_sources( - domain_ids, 'material', ll, ur, nuc_names, densities, vols) - - def _create_cpp_sources_cells(self, time_index, bounding_boxes, - domain_pairs): - """Create decay photon sources in C++ for cell-based calculations. - - Parameters - ---------- - time_index : int - Index into depletion results. - bounding_boxes : dict - Mapping of cell ID to :class:`~openmc.BoundingBox`. - domain_pairs : list of tuple - List of (cell, original_mat) pairs. + For mesh-based: list of (index_mat, mat_id, bbox). + For cell-based: list of (cell, original_mat, bbox). + domain_type : {'material', 'cell'} + The type of domain constraint to apply to each source. """ - # Build work items with bounding boxes - work_items = [ - (cell, original_mat, bounding_boxes[cell.id]) - for cell, original_mat in domain_pairs - ] domain_ids, ll, ur, nuc_names, densities, vols = \ - self._get_region_data(time_index, work_items, 'cell') + self._get_region_data(time_index, work_items, domain_type) openmc.lib.create_decay_photon_sources( - domain_ids, 'cell', ll, ur, nuc_names, densities, vols) + domain_ids, domain_type, ll, ur, nuc_names, densities, vols) def load_results(self, path: PathLike): """Load results from a previous R2S calculation. From bbe1df151a5c580ff0f7aa3a4ddde3aaeaf1f964 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 7 Mar 2026 11:40:50 -0600 Subject: [PATCH 06/37] Unify work items logic --- openmc/deplete/r2s.py | 19 ++++++------------- 1 file changed, 6 insertions(+), 13 deletions(-) diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index 282266a30a3..c59a34a7037 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -517,20 +517,20 @@ def step3_photon_transport( if different_photon_model: photon_cells = self.photon_model.geometry.get_all_cells() - # Determine eligible work items upfront. For cell-based calculations, - # pre-compute the eligible (cell, original_mat) pairs since the - # eligibility check doesn't depend on time. + # Determine eligible work items upfront (independent of time index). if self.method == 'mesh-based': work_items = self._get_mesh_work_items() + domain_type = 'material' else: - _domain_pairs = [] + work_items = [] for cell, original_mat in zip( self.domains, self.results['activation_materials']): if different_photon_model: if cell.id not in photon_cells or \ cell.fill.id != photon_cells[cell.id].fill.id: continue - _domain_pairs.append((cell, original_mat)) + work_items.append((cell, original_mat, bounding_boxes[cell.id])) + domain_type = 'cell' # Ensure photon transport is enabled in settings self.photon_model.settings.photon_transport = True @@ -544,14 +544,7 @@ def step3_photon_transport( photon_dir = Path(output_dir) / f'time_{time_index}' with TemporarySession(self.photon_model, cwd=photon_dir): # Create decay photon sources in C++ memory - if self.method == 'mesh-based': - self._create_cpp_sources(time_index, work_items, 'material') - else: - cell_work_items = [ - (cell, original_mat, bounding_boxes[cell.id]) - for cell, original_mat in _domain_pairs - ] - self._create_cpp_sources(time_index, cell_work_items, 'cell') + self._create_cpp_sources(time_index, work_items, domain_type) statepoint_path = self.photon_model.run(**run_kwargs) From d0a2117610ecd6034d797b76575ea69d2141dfad Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 7 Mar 2026 11:44:06 -0600 Subject: [PATCH 07/37] Simplify use of domain type --- openmc/deplete/r2s.py | 18 +++++++----------- 1 file changed, 7 insertions(+), 11 deletions(-) diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index c59a34a7037..da8270c0e9f 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -520,7 +520,6 @@ def step3_photon_transport( # Determine eligible work items upfront (independent of time index). if self.method == 'mesh-based': work_items = self._get_mesh_work_items() - domain_type = 'material' else: work_items = [] for cell, original_mat in zip( @@ -530,7 +529,6 @@ def step3_photon_transport( cell.fill.id != photon_cells[cell.id].fill.id: continue work_items.append((cell, original_mat, bounding_boxes[cell.id])) - domain_type = 'cell' # Ensure photon transport is enabled in settings self.photon_model.settings.photon_transport = True @@ -544,7 +542,7 @@ def step3_photon_transport( photon_dir = Path(output_dir) / f'time_{time_index}' with TemporarySession(self.photon_model, cwd=photon_dir): # Create decay photon sources in C++ memory - self._create_cpp_sources(time_index, work_items, domain_type) + self._create_cpp_sources(time_index, work_items) statepoint_path = self.photon_model.run(**run_kwargs) @@ -598,7 +596,7 @@ def _get_mesh_work_items(self): return work_items - def _get_region_data(self, time_index, work_items, domain_type): + def _get_region_data(self, time_index, work_items): """Extract region data as flat numpy arrays for C++ source creation. Parameters @@ -608,8 +606,6 @@ def _get_region_data(self, time_index, work_items, domain_type): work_items : list of tuple For mesh-based: list of (index_mat, mat_id, bbox). For cell-based: list of (cell, original_mat, bbox). - domain_type : str - Either ``'material'`` or ``'cell'``. Returns ------- @@ -622,6 +618,7 @@ def _get_region_data(self, time_index, work_items, domain_type): """ step_result = self.results['depletion_results'][time_index] materials = self.results['activation_materials'] + mesh_based = self.method == 'mesh-based' nuclide_names = list(step_result.index_nuc.keys()) n_nuclides = len(nuclide_names) @@ -634,7 +631,7 @@ def _get_region_data(self, time_index, work_items, domain_type): volumes_arr = np.empty(n_regions, dtype=np.float64) for i, item in enumerate(work_items): - if domain_type == 'material': + if mesh_based: index_mat, mat_id, bbox = item domain_ids[i] = mat_id original_mat = materials[index_mat] @@ -658,7 +655,7 @@ def _get_region_data(self, time_index, work_items, domain_type): return (domain_ids, lower_left, upper_right, nuclide_names, atom_densities, volumes_arr) - def _create_cpp_sources(self, time_index, work_items, domain_type): + def _create_cpp_sources(self, time_index, work_items): """Create decay photon sources in C++ for a set of regions. Parameters @@ -668,11 +665,10 @@ def _create_cpp_sources(self, time_index, work_items, domain_type): work_items : list of tuple For mesh-based: list of (index_mat, mat_id, bbox). For cell-based: list of (cell, original_mat, bbox). - domain_type : {'material', 'cell'} - The type of domain constraint to apply to each source. """ domain_ids, ll, ur, nuc_names, densities, vols = \ - self._get_region_data(time_index, work_items, domain_type) + self._get_region_data(time_index, work_items) + domain_type = 'material' if self.method == 'mesh-based' else 'cell' openmc.lib.create_decay_photon_sources( domain_ids, domain_type, ll, ur, nuc_names, densities, vols) From 68cc9d466c735be754fcb94413641584b2c7f3cc Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 9 Mar 2026 09:30:54 -0500 Subject: [PATCH 08/37] Eliminate _get_region_data --- openmc/deplete/r2s.py | 49 +++++++++++-------------------------------- 1 file changed, 12 insertions(+), 37 deletions(-) diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index da8270c0e9f..b46770bd27a 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -596,25 +596,16 @@ def _get_mesh_work_items(self): return work_items - def _get_region_data(self, time_index, work_items): - """Extract region data as flat numpy arrays for C++ source creation. + def _create_cpp_sources(self, time_index, work_items): + """Create decay photon sources in C++ for a set of regions. Parameters ---------- time_index : int - Index into depletion results for the desired time. + Index into depletion results. work_items : list of tuple For mesh-based: list of (index_mat, mat_id, bbox). For cell-based: list of (cell, original_mat, bbox). - - Returns - ------- - domain_ids : numpy.ndarray of int32 - lower_left : numpy.ndarray of float64, shape (n, 3) - upper_right : numpy.ndarray of float64, shape (n, 3) - nuclide_names : list of str - atom_densities : numpy.ndarray of float64, shape (n, n_nuclides) - volumes : numpy.ndarray of float64, shape (n,) """ step_result = self.results['depletion_results'][time_index] materials = self.results['activation_materials'] @@ -628,49 +619,33 @@ def _get_region_data(self, time_index, work_items): lower_left = np.empty((n_regions, 3), dtype=np.float64) upper_right = np.empty((n_regions, 3), dtype=np.float64) atom_densities = np.empty((n_regions, n_nuclides), dtype=np.float64) - volumes_arr = np.empty(n_regions, dtype=np.float64) + volumes = np.empty(n_regions, dtype=np.float64) for i, item in enumerate(work_items): if mesh_based: index_mat, mat_id, bbox = item domain_ids[i] = mat_id original_mat = materials[index_mat] - mat_key = str(original_mat.id) - vol = step_result.volume[mat_key] else: cell, original_mat, bbox = item domain_ids[i] = cell.id - mat_key = str(original_mat.id) - vol = step_result.volume[mat_key] - lower_left[i] = bbox[0] - upper_right[i] = bbox[1] - volumes_arr[i] = vol + mat_key = str(original_mat.id) + vol = step_result.volume[mat_key] + + lower_left[i] = bbox.lower_left + upper_right[i] = bbox.upper_right + volumes[i] = vol # Get atom counts for this material and convert to atom/b-cm mat_idx = step_result.index_mat[mat_key] atom_counts = step_result.data[mat_idx, :] atom_densities[i, :] = atom_counts / vol * 1e-24 - return (domain_ids, lower_left, upper_right, nuclide_names, - atom_densities, volumes_arr) - - def _create_cpp_sources(self, time_index, work_items): - """Create decay photon sources in C++ for a set of regions. - - Parameters - ---------- - time_index : int - Index into depletion results. - work_items : list of tuple - For mesh-based: list of (index_mat, mat_id, bbox). - For cell-based: list of (cell, original_mat, bbox). - """ - domain_ids, ll, ur, nuc_names, densities, vols = \ - self._get_region_data(time_index, work_items) domain_type = 'material' if self.method == 'mesh-based' else 'cell' openmc.lib.create_decay_photon_sources( - domain_ids, domain_type, ll, ur, nuc_names, densities, vols) + domain_ids, domain_type, lower_left, upper_right, nuclide_names, + atom_densities, volumes) def load_results(self, path: PathLike): """Load results from a previous R2S calculation. From 6b1b24f67eb443dc5c2b3e067ef4f77d33be1e1f Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 13 Mar 2026 09:35:07 -0500 Subject: [PATCH 09/37] DecayPhoton univariate type instead of C interface --- include/openmc/capi.h | 20 --- include/openmc/r2s_source.h | 62 ++++++--- openmc/deplete/r2s.py | 73 ++++++---- openmc/lib/core.py | 77 ----------- openmc/stats/univariate.py | 149 +++++++++++++++++++++ src/distribution.cpp | 3 + src/r2s_source.cpp | 258 +++++++++--------------------------- 7 files changed, 299 insertions(+), 343 deletions(-) diff --git a/include/openmc/capi.h b/include/openmc/capi.h index e2339719fe2..911654d318f 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -320,26 +320,6 @@ int openmc_properties_export(const char* filename); // \return Error code int openmc_properties_import(const char* filename); -//! Create decay photon sources for R2S calculations -//! -//! Clears existing external sources and creates IndependentSource objects -//! for each spatial region based on the activated nuclide compositions. -//! -//! \param n_regions Number of spatial regions -//! \param domain_ids Material or cell IDs for constraints [n_regions] -//! \param domain_type "material" or "cell" -//! \param lower_left Flattened lower-left bbox corners [n_regions * 3] -//! \param upper_right Flattened upper-right bbox corners [n_regions * 3] -//! \param n_nuclides Number of nuclides -//! \param nuclide_names Nuclide name strings [n_nuclides] -//! \param atom_densities Atom densities [n_regions * n_nuclides] row-major -//! \param volumes Region volumes in cm^3 [n_regions] -//! \return Error code (0 on success) -int openmc_create_decay_photon_sources(int n_regions, const int32_t* domain_ids, - const char* domain_type, const double* lower_left, const double* upper_right, - int n_nuclides, const char** nuclide_names, const double* atom_densities, - const double* volumes); - // Error codes extern int OPENMC_E_UNASSIGNED; extern int OPENMC_E_ALLOCATE; diff --git a/include/openmc/r2s_source.h b/include/openmc/r2s_source.h index 6da76222af3..cb0d336a9d3 100644 --- a/include/openmc/r2s_source.h +++ b/include/openmc/r2s_source.h @@ -7,33 +7,51 @@ #include // for size_t #include // for int32_t +#include "openmc/distribution.h" #include "openmc/source.h" +#include "openmc/vector.h" namespace openmc { -//! Create decay photon IndependentSource objects from activated materials -//! -//! For each spatial region, this function combines the decay photon energy -//! spectra from all nuclides present (weighted by atom count and decay -//! constant) and creates an IndependentSource with a SpatialBox spatial -//! distribution and isotropic angular distribution. The resulting sources -//! are stored in model::external_sources with the probability mass function -//! rebuilt accordingly. +//============================================================================== +// DecayPhotonMixture — non-owning mixture of decay photon distributions +//============================================================================== + +//! Energy distribution formed by mixing multiple decay photon spectra. //! -//! \param n_regions Number of spatial regions -//! \param domain_ids Material or cell IDs for source constraints [n_regions] -//! \param domain_type MATERIAL or CELL -//! \param lower_left Flattened lower-left bbox corners [n_regions * 3] -//! \param upper_right Flattened upper-right bbox corners [n_regions * 3] -//! \param n_nuclides Number of nuclides in the atom density matrix -//! \param nuclide_names Array of nuclide name strings [n_nuclides] -//! \param atom_densities Row-major atom densities in [atom/b-cm], -//! shape [n_regions * n_nuclides] -//! \param volumes Volume of each region in [cm^3], shape [n_regions] -void create_decay_photon_sources(int n_regions, const int32_t* domain_ids, - Source::DomainType domain_type, const double* lower_left, - const double* upper_right, int n_nuclides, const char** nuclide_names, - const double* atom_densities, const double* volumes); +//! Unlike the general Mixture distribution, this class holds non-owning +//! pointers to the component distributions (which live in +//! data::chain_nuclides). Each component is weighted by the activity +//! (atoms * decay_constant) of the corresponding nuclide. + +class DecayPhotonMixture : public Distribution { +public: + //! Construct from non-owning distribution pointers and weights + //! + //! \param dists Non-owning pointers to component distributions + //! \param weights Activity-based weights for each component. The Mixture + //! probability for component i is weights[i] * dists[i]->integral() + //! (the integral encodes photons-per-decay). + DecayPhotonMixture(vector dists, vector weights); + + //! Construct from an XML node containing nuclide names and atom densities. + //! + //! Reads child ```` elements with ``name`` and ``density`` + //! attributes, resolves them against the loaded depletion chain, and + //! constructs the mixed distribution. + explicit DecayPhotonMixture(pugi::xml_node node); + + std::pair sample(uint64_t* seed) const override; + double integral() const override; + +protected: + double sample_unbiased(uint64_t* seed) const override; + +private: + vector dists_; //!< Non-owning component distributions + DiscreteIndex di_; //!< Discrete index for component selection + double integral_; //!< Total photon emission rate +}; } // namespace openmc diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index b46770bd27a..bb25ad3356d 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -7,7 +7,6 @@ import numpy as np import openmc -import openmc.lib from . import IndependentOperator, PredictorIntegrator from .microxs import get_microxs_and_flux, write_microxs_hdf5, read_microxs_hdf5 from .results import Results @@ -538,12 +537,13 @@ def step3_photon_transport( if time_index < 0: time_index = len(self.results['depletion_results']) + time_index + # Build decay photon sources and assign to the photon model + sources = self._create_photon_sources(time_index, work_items) + self.photon_model.settings.source = sources + # Run photon transport calculation photon_dir = Path(output_dir) / f'time_{time_index}' with TemporarySession(self.photon_model, cwd=photon_dir): - # Create decay photon sources in C++ memory - self._create_cpp_sources(time_index, work_items) - statepoint_path = self.photon_model.run(**run_kwargs) # Store tally results @@ -596,8 +596,13 @@ def _get_mesh_work_items(self): return work_items - def _create_cpp_sources(self, time_index, work_items): - """Create decay photon sources in C++ for a set of regions. + def _create_photon_sources(self, time_index, work_items): + """Create decay photon sources for a set of regions. + + Builds :class:`openmc.IndependentSource` objects with + :class:`openmc.stats.DecayPhoton` energy distributions that will be + serialized to XML and resolved against the depletion chain by the C++ + solver. Parameters ---------- @@ -606,46 +611,58 @@ def _create_cpp_sources(self, time_index, work_items): work_items : list of tuple For mesh-based: list of (index_mat, mat_id, bbox). For cell-based: list of (cell, original_mat, bbox). + + Returns + ------- + list of openmc.IndependentSource + Photon sources for each activated region. """ step_result = self.results['depletion_results'][time_index] materials = self.results['activation_materials'] mesh_based = self.method == 'mesh-based' nuclide_names = list(step_result.index_nuc.keys()) - n_nuclides = len(nuclide_names) - n_regions = len(work_items) - - domain_ids = np.empty(n_regions, dtype=np.int32) - lower_left = np.empty((n_regions, 3), dtype=np.float64) - upper_right = np.empty((n_regions, 3), dtype=np.float64) - atom_densities = np.empty((n_regions, n_nuclides), dtype=np.float64) - volumes = np.empty(n_regions, dtype=np.float64) - for i, item in enumerate(work_items): + sources = [] + for item in work_items: if mesh_based: - index_mat, mat_id, bbox = item - domain_ids[i] = mat_id + index_mat, domain_id, bbox = item original_mat = materials[index_mat] else: cell, original_mat, bbox = item - domain_ids[i] = cell.id + domain_id = cell.id mat_key = str(original_mat.id) vol = step_result.volume[mat_key] - lower_left[i] = bbox.lower_left - upper_right[i] = bbox.upper_right - volumes[i] = vol - - # Get atom counts for this material and convert to atom/b-cm + # Get atom counts and convert to atom/b-cm mat_idx = step_result.index_mat[mat_key] atom_counts = step_result.data[mat_idx, :] - atom_densities[i, :] = atom_counts / vol * 1e-24 + atom_densities = atom_counts / vol * 1e-24 + + # Build nuclide dict with only non-zero densities + nuclides = { + name: dens for name, dens in zip(nuclide_names, atom_densities) + if dens > 0.0 + } + if not nuclides: + continue + + energy = openmc.stats.DecayPhoton(nuclides) + + if mesh_based: + domains = [openmc.Material(material_id=domain_id)] + else: + domains = [openmc.Cell(cell_id=domain_id)] + source = openmc.IndependentSource( + space=openmc.stats.Box(bbox.lower_left, bbox.upper_right), + energy=energy, + particle='photon', + constraints={'domains': domains}, + ) + sources.append(source) - domain_type = 'material' if self.method == 'mesh-based' else 'cell' - openmc.lib.create_decay_photon_sources( - domain_ids, domain_type, lower_left, upper_right, nuclide_names, - atom_densities, volumes) + return sources def load_results(self, path: PathLike): """Load results from a previous R2S calculation. diff --git a/openmc/lib/core.py b/openmc/lib/core.py index c25c324411d..00876db9972 100644 --- a/openmc/lib/core.py +++ b/openmc/lib/core.py @@ -112,19 +112,6 @@ class _SourceSite(Structure): _dll.openmc_sample_external_source.argtypes = [c_size_t, POINTER(c_uint64), POINTER(_SourceSite)] _dll.openmc_sample_external_source.restype = c_int _dll.openmc_sample_external_source.errcheck = _error_handler -_dll.openmc_create_decay_photon_sources.argtypes = [ - c_int, # n_regions - POINTER(c_int32), # domain_ids - c_char_p, # domain_type - POINTER(c_double), # lower_left - POINTER(c_double), # upper_right - c_int, # n_nuclides - POINTER(c_char_p), # nuclide_names - POINTER(c_double), # atom_densities - POINTER(c_double), # volumes -] -_dll.openmc_create_decay_photon_sources.restype = c_int -_dll.openmc_create_decay_photon_sources.errcheck = _error_handler def global_bounding_box(): """Calculate a global bounding box for the model""" @@ -504,70 +491,6 @@ def run_random_ray(output=True): with quiet_dll(output): _dll.openmc_run_random_ray() -def create_decay_photon_sources( - domain_ids: np.ndarray, - domain_type: str, - lower_left: np.ndarray, - upper_right: np.ndarray, - nuclide_names: list[str], - atom_densities: np.ndarray, - volumes: np.ndarray, -): - """Create decay photon sources from activated material compositions. - - This populates the internal external source bank with - :class:`~openmc.IndependentSource` objects for each activated region with - non-zero photon emission. Existing external sources are cleared before the - new sources are created. - - .. versionadded:: 0.15.1 - - Parameters - ---------- - domain_ids : numpy.ndarray of int32 - Material or cell IDs for source domain constraints, shape (n_regions,). - domain_type : str - Either "material" or "cell". - lower_left : numpy.ndarray of float64 - Lower-left bounding box corners, shape (n_regions, 3). - upper_right : numpy.ndarray of float64 - Upper-right bounding box corners, shape (n_regions, 3). - nuclide_names : list of str - Names of nuclides in the atom density matrix. - atom_densities : numpy.ndarray of float64 - Atom densities in [atom/b-cm], shape (n_regions, n_nuclides). - volumes : numpy.ndarray of float64 - Volume of each region in [cm^3], shape (n_regions,). - - """ - n_regions = len(domain_ids) - n_nuclides = len(nuclide_names) - - # Ensure arrays are contiguous with correct dtype - domain_ids = np.ascontiguousarray(domain_ids, dtype=np.int32) - lower_left = np.ascontiguousarray(lower_left.ravel(), dtype=np.float64) - upper_right = np.ascontiguousarray(upper_right.ravel(), dtype=np.float64) - atom_densities = np.ascontiguousarray(atom_densities.ravel(), dtype=np.float64) - volumes = np.ascontiguousarray(volumes, dtype=np.float64) - - # Convert nuclide names to array of c_char_p - names_arr = (c_char_p * n_nuclides)( - *(name.encode() for name in nuclide_names) - ) - - _dll.openmc_create_decay_photon_sources( - n_regions, - domain_ids.ctypes.data_as(POINTER(c_int32)), - domain_type.encode(), - lower_left.ctypes.data_as(POINTER(c_double)), - upper_right.ctypes.data_as(POINTER(c_double)), - n_nuclides, - names_arr, - atom_densities.ctypes.data_as(POINTER(c_double)), - volumes.ctypes.data_as(POINTER(c_double)), - ) - - def sample_external_source( n_samples: int = 1000, prn_seed: int | None = None, diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py index 1cf9a1ad507..4a86f461a69 100644 --- a/openmc/stats/univariate.py +++ b/openmc/stats/univariate.py @@ -124,6 +124,8 @@ def from_xml_element(cls, elem): return Legendre.from_xml_element(elem) elif distribution == 'mixture': return Mixture.from_xml_element(elem) + elif distribution == 'decay_photon': + return DecayPhoton.from_xml_element(elem) @abstractmethod def _sample_unbiased(self, n_samples: int = 1, seed: int | None = None): @@ -2196,6 +2198,153 @@ def clip(self, tolerance: float = 1e-6, inplace: bool = False) -> Mixture: return new_dist +class DecayPhoton(Univariate): + """Energy distribution from decay photon spectra of a mixture of nuclides. + + This distribution stores nuclide names and their atom densities. When + written to XML and read by the C++ solver, the nuclide names are resolved + against the depletion chain to obtain the decay photon energy spectra + and decay constants. The resulting distribution is a mixture of per-nuclide + photon spectra weighted by activity (atom_density * decay_constant * + photons_per_decay). + + .. versionadded:: 0.15.1 + + Parameters + ---------- + nuclides : dict + Dictionary mapping nuclide name (str) to atom density (float) in + units of atom/b-cm. + + Attributes + ---------- + nuclides : dict + Dictionary mapping nuclide name to atom density in atom/b-cm. + + """ + + def __init__(self, nuclides: dict[str, float]): + super().__init__(bias=None) + self.nuclides = nuclides + + def __len__(self): + return len(self.nuclides) + + @property + def nuclides(self): + return self._nuclides + + @nuclides.setter + def nuclides(self, nuclides): + cv.check_type('nuclides', nuclides, dict) + for name, density in nuclides.items(): + cv.check_type('nuclide name', name, str) + cv.check_type(f'atom density for {name}', density, Real) + cv.check_greater_than(f'atom density for {name}', density, 0.0) + self._nuclides = dict(nuclides) + + def to_xml_element(self, element_name: str): + """Return XML representation of the decay photon distribution + + Parameters + ---------- + element_name : str + XML element name + + Returns + ------- + element : lxml.etree._Element + XML element containing decay photon distribution data + + """ + element = ET.Element(element_name) + element.set("type", "decay_photon") + for name, density in self.nuclides.items(): + nuclide_elem = ET.SubElement(element, "nuclide") + nuclide_elem.set("name", name) + nuclide_elem.set("density", str(density)) + return element + + @classmethod + def from_xml_element(cls, elem: ET.Element): + """Generate decay photon distribution from an XML element + + Parameters + ---------- + elem : lxml.etree._Element + XML element + + Returns + ------- + openmc.stats.DecayPhoton + Decay photon distribution generated from XML element + + """ + nuclides = {} + for nuclide_elem in elem.findall('nuclide'): + name = nuclide_elem.get('name') + density = float(nuclide_elem.get('density')) + nuclides[name] = density + return cls(nuclides) + + def _sample_unbiased(self, n_samples=1, seed=None): + raise NotImplementedError( + "DecayPhoton distributions can only be sampled by the C++ solver") + + def integral(self): + """Return integral of the distribution + + Returns + ------- + float + Integral of the distribution (requires chain data, returns 0.0 + when chain data is not available) + """ + return 0.0 + + def clip(self, tolerance: float = 1e-6, inplace: bool = False): + """Remove nuclides with negligible contribution to photon emission. + + Nuclides are sorted by their atom density and those contributing the + least are removed until the cumulative removed fraction exceeds the + tolerance. + + Parameters + ---------- + tolerance : float + Maximum fraction of total atom density that can be discarded. + inplace : bool + Whether to modify the current object in-place or return a new one. + + Returns + ------- + openmc.stats.DecayPhoton + Distribution with low-density nuclides removed + + """ + densities = np.array(list(self.nuclides.values())) + names = list(self.nuclides.keys()) + indices = _intensity_clip(densities, tolerance=tolerance) + new_nuclides = {names[i]: densities[i] for i in indices} + + if inplace: + self._nuclides = new_nuclides + return self + return type(self)(new_nuclides) + + @property + def support(self): + return (0.0, np.inf) + + def evaluate(self, x): + raise NotImplementedError( + "evaluate() is undefined for DecayPhoton distributions") + + def mean(self): + raise NotImplementedError( + "mean() is undefined for DecayPhoton distributions") + + def combine_distributions( dists: Sequence[Discrete | Tabular | Mixture], probs: Sequence[float] diff --git a/src/distribution.cpp b/src/distribution.cpp index c722d0366b6..c3ec71268df 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -11,6 +11,7 @@ #include "openmc/constants.h" #include "openmc/error.h" #include "openmc/math_functions.h" +#include "openmc/r2s_source.h" #include "openmc/random_dist.h" #include "openmc/random_lcg.h" #include "openmc/xml_interface.h" @@ -758,6 +759,8 @@ UPtrDist distribution_from_xml(pugi::xml_node node) dist = UPtrDist {new Tabular(node)}; } else if (type == "mixture") { dist = UPtrDist {new Mixture(node)}; + } else if (type == "decay_photon") { + dist = UPtrDist {new DecayPhotonMixture(node)}; } else if (type == "muir") { openmc::fatal_error( "'muir' distributions are now specified using the openmc.stats.muir() " diff --git a/src/r2s_source.cpp b/src/r2s_source.cpp index 9ef938f8a0d..7947a7d5e7c 100644 --- a/src/r2s_source.cpp +++ b/src/r2s_source.cpp @@ -3,224 +3,90 @@ #include "openmc/r2s_source.h" -#include // for sort -#include // for accumulate +#include // for accumulate #include -#include "openmc/capi.h" #include "openmc/chain.h" -#include "openmc/distribution.h" -#include "openmc/distribution_multi.h" -#include "openmc/distribution_spatial.h" #include "openmc/error.h" #include "openmc/memory.h" -#include "openmc/openmp_interface.h" -#include "openmc/particle_type.h" -#include "openmc/settings.h" -#include "openmc/source.h" #include "openmc/vector.h" +#include "openmc/xml_interface.h" namespace openmc { //============================================================================== -// DecayPhotonMixture — non-owning mixture of decay photon distributions +// DecayPhotonMixture implementation //============================================================================== -//! Energy distribution formed by mixing multiple decay photon spectra. -//! -//! Unlike the general Mixture distribution, this class holds non-owning -//! pointers to the component distributions (which live in -//! data::chain_nuclides). Each component is weighted by the activity -//! (atoms * decay_constant) of the corresponding nuclide. - -class DecayPhotonMixture : public Distribution { -public: - //! Construct from non-owning distribution pointers and weights - //! - //! \param dists Non-owning pointers to component distributions - //! \param weights Activity-based weights for each component. The Mixture - //! probability for component i is weights[i] * dists[i]->integral() - //! (the integral encodes photons-per-decay). - DecayPhotonMixture(vector dists, vector weights) - : dists_(std::move(dists)) - { - vector probs; - probs.reserve(dists_.size()); - for (size_t i = 0; i < dists_.size(); ++i) { - probs.push_back(weights[i] * dists_[i]->integral()); - } - integral_ = std::accumulate(probs.begin(), probs.end(), 0.0); - di_.assign(probs); - } - - std::pair sample(uint64_t* seed) const override - { - size_t idx = di_.sample(seed); - return dists_[idx]->sample(seed); - } - - double integral() const override { return integral_; } - -protected: - double sample_unbiased(uint64_t* seed) const override - { - size_t idx = di_.sample(seed); - return dists_[idx]->sample(seed).first; - } - -private: - vector dists_; //!< Non-owning component distributions - DiscreteIndex di_; //!< Discrete index for component selection - double integral_; //!< Total photon emission rate -}; - -//============================================================================== -// create_decay_photon_sources implementation -//============================================================================== - -void create_decay_photon_sources(int n_regions, const int32_t* domain_ids, - Source::DomainType domain_type, const double* lower_left, - const double* upper_right, int n_nuclides, const char** nuclide_names, - const double* atom_densities, const double* volumes) +DecayPhotonMixture::DecayPhotonMixture( + vector dists, vector weights) + : dists_(std::move(dists)) { - // Pre-resolve nuclide names to chain nuclide indices. A value of -1 - // indicates the nuclide is not in the chain (or has no decay photon data). - vector chain_indices(n_nuclides, -1); - for (int j = 0; j < n_nuclides; ++j) { - auto it = data::chain_nuclide_map.find(nuclide_names[j]); - if (it != data::chain_nuclide_map.end()) { - chain_indices[j] = it->second; - } + vector probs; + probs.reserve(dists_.size()); + for (size_t i = 0; i < dists_.size(); ++i) { + probs.push_back(weights[i] * dists_[i]->integral()); } + integral_ = std::accumulate(probs.begin(), probs.end(), 0.0); + di_.assign(probs); +} - // Clear existing sources - model::external_sources.clear(); - - // Thread-local storage for sources built in parallel - vector>> thread_sources; - - int n_threads = num_threads(); - thread_sources.resize(n_threads); - -#pragma omp parallel - { - int tid = thread_num(); - -#pragma omp for schedule(dynamic, 64) - for (int i = 0; i < n_regions; ++i) { - // Collect distributions and weights for nuclides present in this region - vector dists; - vector weights; - - for (int j = 0; j < n_nuclides; ++j) { - double atom_dens = atom_densities[i * n_nuclides + j]; - if (atom_dens <= 0.0) - continue; - - int chain_idx = chain_indices[j]; - if (chain_idx < 0) - continue; - - const auto& chain_nuc = data::chain_nuclides[chain_idx]; - const Distribution* photon_dist = chain_nuc->photon_energy(); - if (!photon_dist) - continue; - - // Weight = number of atoms * decay_constant - // atom_dens is in [atom/b-cm], volumes in [cm^3] - // atoms = atom_dens * 1e24 * volume - double atoms = atom_dens * 1.0e24 * volumes[i]; - double activity = atoms * chain_nuc->decay_constant(); - - dists.push_back(photon_dist); - weights.push_back(activity); - } - - // Skip regions with no photon-emitting nuclides - if (dists.empty()) - continue; - - // Build combined energy distribution - auto energy = - make_unique(std::move(dists), std::move(weights)); - double strength = energy->integral(); - if (strength <= 0.0) - continue; - - // Build spatial distribution from bounding box - Position ll { - lower_left[i * 3], lower_left[i * 3 + 1], lower_left[i * 3 + 2]}; - Position ur { - upper_right[i * 3], upper_right[i * 3 + 1], upper_right[i * 3 + 2]}; - - // Build time distribution (delta at t=0) - double t0[] {0.0}; - double tp[] {1.0}; - - // Create IndependentSource - auto source = make_unique( - UPtrSpace {new SpatialBox(ll, ur)}, UPtrAngle {new Isotropic()}, - std::move(energy), UPtrDist {new Discrete(t0, tp, 1)}); - source->set_particle(ParticleType::photon()); - source->set_strength(strength); - source->set_domain(domain_type, {domain_ids[i]}); - - thread_sources[tid].push_back(std::move(source)); - } - } // end parallel - - // Merge thread-local sources into model::external_sources - for (auto& tvec : thread_sources) { - for (auto& src : tvec) { - model::external_sources.push_back(std::move(src)); - } +DecayPhotonMixture::DecayPhotonMixture(pugi::xml_node node) +{ + // Read nuclide names and atom densities from XML + vector dists; + vector weights; + + for (auto nuclide_node : node.children("nuclide")) { + std::string name = get_node_value(nuclide_node, "name"); + double density = std::stod(get_node_value(nuclide_node, "density")); + + // Look up nuclide in the depletion chain + auto it = data::chain_nuclide_map.find(name); + if (it == data::chain_nuclide_map.end()) + continue; + + const auto& chain_nuc = data::chain_nuclides[it->second]; + const Distribution* photon_dist = chain_nuc->photon_energy(); + if (!photon_dist) + continue; + + // Weight = atom_density * decay_constant + double weight = density * chain_nuc->decay_constant(); + if (weight <= 0.0) + continue; + + dists.push_back(photon_dist); + weights.push_back(weight); } - // Enable photon transport since we are creating photon sources - settings::photon_transport = true; + dists_ = std::move(dists); - // Rebuild probability mass function for sampling - vector source_strengths; - source_strengths.reserve(model::external_sources.size()); - for (auto& s : model::external_sources) { - source_strengths.push_back(s->strength()); + // Build alias table from weighted probabilities + vector probs; + probs.reserve(dists_.size()); + for (size_t i = 0; i < dists_.size(); ++i) { + probs.push_back(weights[i] * dists_[i]->integral()); } - model::external_sources_probability.assign(source_strengths); + integral_ = std::accumulate(probs.begin(), probs.end(), 0.0); + di_.assign(probs); } -} // namespace openmc - -//============================================================================== -// C API -//============================================================================== - -extern "C" int openmc_create_decay_photon_sources(int n_regions, - const int32_t* domain_ids, const char* domain_type, const double* lower_left, - const double* upper_right, int n_nuclides, const char** nuclide_names, - const double* atom_densities, const double* volumes) +std::pair DecayPhotonMixture::sample(uint64_t* seed) const { - using DT = openmc::Source::DomainType; - - // Validate domain type - std::string dt_str(domain_type); - DT dt; - if (dt_str == "material") { - dt = DT::MATERIAL; - } else if (dt_str == "cell") { - dt = DT::CELL; - } else { - openmc::set_errmsg( - "Invalid domain_type '" + dt_str + "'. Must be 'material' or 'cell'."); - return OPENMC_E_INVALID_ARGUMENT; - } + size_t idx = di_.sample(seed); + return dists_[idx]->sample(seed); +} - try { - openmc::create_decay_photon_sources(n_regions, domain_ids, dt, lower_left, - upper_right, n_nuclides, nuclide_names, atom_densities, volumes); - } catch (const std::exception& e) { - openmc::set_errmsg(e.what()); - return OPENMC_E_DATA; - } +double DecayPhotonMixture::integral() const +{ + return integral_; +} - return 0; +double DecayPhotonMixture::sample_unbiased(uint64_t* seed) const +{ + size_t idx = di_.sample(seed); + return dists_[idx]->sample(seed).first; } + +} // namespace openmc From 8d098f506777b3ebc91b362cf5641a1cfb782109 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 15 Mar 2026 16:50:15 -0500 Subject: [PATCH 10/37] Pass volume through DecayPhoton --- include/openmc/r2s_source.h | 1 - openmc/deplete/r2s.py | 2 +- openmc/stats/univariate.py | 40 +++++++++++++++++++++++++++---------- src/r2s_source.cpp | 11 ++++++++-- src/source.cpp | 8 ++++++++ 5 files changed, 48 insertions(+), 14 deletions(-) diff --git a/include/openmc/r2s_source.h b/include/openmc/r2s_source.h index cb0d336a9d3..2bb08ed7f4c 100644 --- a/include/openmc/r2s_source.h +++ b/include/openmc/r2s_source.h @@ -8,7 +8,6 @@ #include // for int32_t #include "openmc/distribution.h" -#include "openmc/source.h" #include "openmc/vector.h" namespace openmc { diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index bb25ad3356d..39b9407b1d7 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -648,7 +648,7 @@ def _create_photon_sources(self, time_index, work_items): if not nuclides: continue - energy = openmc.stats.DecayPhoton(nuclides) + energy = openmc.stats.DecayPhoton(nuclides, vol) if mesh_based: domains = [openmc.Material(material_id=domain_id)] diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py index 4a86f461a69..b670cd055fa 100644 --- a/openmc/stats/univariate.py +++ b/openmc/stats/univariate.py @@ -2201,31 +2201,39 @@ def clip(self, tolerance: float = 1e-6, inplace: bool = False) -> Mixture: class DecayPhoton(Univariate): """Energy distribution from decay photon spectra of a mixture of nuclides. - This distribution stores nuclide names and their atom densities. When - written to XML and read by the C++ solver, the nuclide names are resolved - against the depletion chain to obtain the decay photon energy spectra - and decay constants. The resulting distribution is a mixture of per-nuclide - photon spectra weighted by activity (atom_density * decay_constant * - photons_per_decay). + This distribution stores nuclide names, their atom densities, and the + volume of the region. When written to XML and read by the C++ solver, the + nuclide names are resolved against the depletion chain to obtain the decay + photon energy spectra and decay constants. The resulting distribution is a + mixture of per-nuclide photon spectra weighted by absolute activity + (atom_density * 1e24 * volume * decay_constant * photons_per_decay). + The volume is necessary so that the C++ solver can compute the total photon + emission rate in [photons/s], which is used as the source strength. - .. versionadded:: 0.15.1 + .. versionadded:: 0.15.4 Parameters ---------- nuclides : dict Dictionary mapping nuclide name (str) to atom density (float) in units of atom/b-cm. + volume : float + Volume of the source region in cm\ :sup:`3`. Used together with atom + densities to compute the absolute photon emission rate. Attributes ---------- nuclides : dict Dictionary mapping nuclide name to atom density in atom/b-cm. + volume : float + Volume of the source region in cm\ :sup:`3`. """ - def __init__(self, nuclides: dict[str, float]): + def __init__(self, nuclides: dict[str, float], volume: float): super().__init__(bias=None) self.nuclides = nuclides + self.volume = volume def __len__(self): return len(self.nuclides) @@ -2243,6 +2251,16 @@ def nuclides(self, nuclides): cv.check_greater_than(f'atom density for {name}', density, 0.0) self._nuclides = dict(nuclides) + @property + def volume(self): + return self._volume + + @volume.setter + def volume(self, volume): + cv.check_type('volume', volume, Real) + cv.check_greater_than('volume', volume, 0.0) + self._volume = float(volume) + def to_xml_element(self, element_name: str): """Return XML representation of the decay photon distribution @@ -2259,6 +2277,7 @@ def to_xml_element(self, element_name: str): """ element = ET.Element(element_name) element.set("type", "decay_photon") + element.set("volume", str(self.volume)) for name, density in self.nuclides.items(): nuclide_elem = ET.SubElement(element, "nuclide") nuclide_elem.set("name", name) @@ -2280,12 +2299,13 @@ def from_xml_element(cls, elem: ET.Element): Decay photon distribution generated from XML element """ + volume = float(elem.get('volume')) nuclides = {} for nuclide_elem in elem.findall('nuclide'): name = nuclide_elem.get('name') density = float(nuclide_elem.get('density')) nuclides[name] = density - return cls(nuclides) + return cls(nuclides, volume) def _sample_unbiased(self, n_samples=1, seed=None): raise NotImplementedError( @@ -2330,7 +2350,7 @@ def clip(self, tolerance: float = 1e-6, inplace: bool = False): if inplace: self._nuclides = new_nuclides return self - return type(self)(new_nuclides) + return type(self)(new_nuclides, self.volume) @property def support(self): diff --git a/src/r2s_source.cpp b/src/r2s_source.cpp index 7947a7d5e7c..96afd0635ac 100644 --- a/src/r2s_source.cpp +++ b/src/r2s_source.cpp @@ -33,6 +33,11 @@ DecayPhotonMixture::DecayPhotonMixture( DecayPhotonMixture::DecayPhotonMixture(pugi::xml_node node) { + // Read the region volume [cm^3] needed for absolute emission rate + if (!check_for_node(node, "volume")) + fatal_error("DecayPhotonMixture: 'volume' attribute is required."); + double volume = std::stod(get_node_value(node, "volume")); + // Read nuclide names and atom densities from XML vector dists; vector weights; @@ -51,8 +56,10 @@ DecayPhotonMixture::DecayPhotonMixture(pugi::xml_node node) if (!photon_dist) continue; - // Weight = atom_density * decay_constant - double weight = density * chain_nuc->decay_constant(); + // Weight = atom_count * decay_constant + // atom_count = atom_density [atom/b-cm] * 1e24 [b-cm/cm^3] * volume [cm^3] + double atom_count = density * 1.0e24 * volume; + double weight = atom_count * chain_nuc->decay_constant(); if (weight <= 0.0) continue; diff --git a/src/source.cpp b/src/source.cpp index f0b8f48edda..22492083efd 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -27,6 +27,7 @@ #include "openmc/message_passing.h" #include "openmc/mgxs_interface.h" #include "openmc/nuclide.h" +#include "openmc/r2s_source.h" #include "openmc/random_lcg.h" #include "openmc/search.h" #include "openmc/settings.h" @@ -352,6 +353,13 @@ IndependentSource::IndependentSource(pugi::xml_node node) : Source(node) if (check_for_node(node, "energy")) { pugi::xml_node node_dist = node.child("energy"); energy_ = distribution_from_xml(node_dist); + + // For decay photon sources the C++ computes the absolute photon emission + // rate [photons/s] as the distribution integral. Use it as the source + // strength so that sources in more active regions are sampled more often. + if (dynamic_cast(energy_.get())) { + strength_ = energy_->integral(); + } } else { // Default to a Watt spectrum with parameters 0.988 MeV and 2.249 MeV^-1 energy_ = UPtrDist {new Watt(0.988e6, 2.249e-6)}; From f22ecba0a07fe30021ef277f09be258eeeaa0f57 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 15 Mar 2026 17:23:05 -0500 Subject: [PATCH 11/37] Rename to DecaySpectrum --- docs/source/pythonapi/stats.rst | 1 + include/openmc/r2s_source.h | 8 ++++---- openmc/deplete/r2s.py | 4 ++-- openmc/stats/univariate.py | 14 +++++++------- src/distribution.cpp | 2 +- src/r2s_source.cpp | 14 +++++++------- src/source.cpp | 2 +- 7 files changed, 23 insertions(+), 22 deletions(-) diff --git a/docs/source/pythonapi/stats.rst b/docs/source/pythonapi/stats.rst index 203d4fea4a4..2f2e2835a29 100644 --- a/docs/source/pythonapi/stats.rst +++ b/docs/source/pythonapi/stats.rst @@ -22,6 +22,7 @@ Univariate Probability Distributions openmc.stats.Legendre openmc.stats.Mixture openmc.stats.Normal + openmc.stats.DecaySpectrum .. autosummary:: :toctree: generated diff --git a/include/openmc/r2s_source.h b/include/openmc/r2s_source.h index 2bb08ed7f4c..7a4ba0410b4 100644 --- a/include/openmc/r2s_source.h +++ b/include/openmc/r2s_source.h @@ -13,7 +13,7 @@ namespace openmc { //============================================================================== -// DecayPhotonMixture — non-owning mixture of decay photon distributions +// DecaySpectrum — non-owning mixture of decay photon distributions //============================================================================== //! Energy distribution formed by mixing multiple decay photon spectra. @@ -23,7 +23,7 @@ namespace openmc { //! data::chain_nuclides). Each component is weighted by the activity //! (atoms * decay_constant) of the corresponding nuclide. -class DecayPhotonMixture : public Distribution { +class DecaySpectrum : public Distribution { public: //! Construct from non-owning distribution pointers and weights //! @@ -31,14 +31,14 @@ class DecayPhotonMixture : public Distribution { //! \param weights Activity-based weights for each component. The Mixture //! probability for component i is weights[i] * dists[i]->integral() //! (the integral encodes photons-per-decay). - DecayPhotonMixture(vector dists, vector weights); + DecaySpectrum(vector dists, vector weights); //! Construct from an XML node containing nuclide names and atom densities. //! //! Reads child ```` elements with ``name`` and ``density`` //! attributes, resolves them against the loaded depletion chain, and //! constructs the mixed distribution. - explicit DecayPhotonMixture(pugi::xml_node node); + explicit DecaySpectrum(pugi::xml_node node); std::pair sample(uint64_t* seed) const override; double integral() const override; diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index 39b9407b1d7..2680533a86a 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -600,7 +600,7 @@ def _create_photon_sources(self, time_index, work_items): """Create decay photon sources for a set of regions. Builds :class:`openmc.IndependentSource` objects with - :class:`openmc.stats.DecayPhoton` energy distributions that will be + :class:`openmc.stats.DecaySpectrum` energy distributions that will be serialized to XML and resolved against the depletion chain by the C++ solver. @@ -648,7 +648,7 @@ def _create_photon_sources(self, time_index, work_items): if not nuclides: continue - energy = openmc.stats.DecayPhoton(nuclides, vol) + energy = openmc.stats.DecaySpectrum(nuclides, vol) if mesh_based: domains = [openmc.Material(material_id=domain_id)] diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py index b670cd055fa..926381a4851 100644 --- a/openmc/stats/univariate.py +++ b/openmc/stats/univariate.py @@ -125,7 +125,7 @@ def from_xml_element(cls, elem): elif distribution == 'mixture': return Mixture.from_xml_element(elem) elif distribution == 'decay_photon': - return DecayPhoton.from_xml_element(elem) + return DecaySpectrum.from_xml_element(elem) @abstractmethod def _sample_unbiased(self, n_samples: int = 1, seed: int | None = None): @@ -2198,7 +2198,7 @@ def clip(self, tolerance: float = 1e-6, inplace: bool = False) -> Mixture: return new_dist -class DecayPhoton(Univariate): +class DecaySpectrum(Univariate): """Energy distribution from decay photon spectra of a mixture of nuclides. This distribution stores nuclide names, their atom densities, and the @@ -2295,7 +2295,7 @@ def from_xml_element(cls, elem: ET.Element): Returns ------- - openmc.stats.DecayPhoton + openmc.stats.DecaySpectrum Decay photon distribution generated from XML element """ @@ -2309,7 +2309,7 @@ def from_xml_element(cls, elem: ET.Element): def _sample_unbiased(self, n_samples=1, seed=None): raise NotImplementedError( - "DecayPhoton distributions can only be sampled by the C++ solver") + "DecaySpectrum distributions can only be sampled by the C++ solver") def integral(self): """Return integral of the distribution @@ -2338,7 +2338,7 @@ def clip(self, tolerance: float = 1e-6, inplace: bool = False): Returns ------- - openmc.stats.DecayPhoton + openmc.stats.DecaySpectrum Distribution with low-density nuclides removed """ @@ -2358,11 +2358,11 @@ def support(self): def evaluate(self, x): raise NotImplementedError( - "evaluate() is undefined for DecayPhoton distributions") + "evaluate() is undefined for DecaySpectrum distributions") def mean(self): raise NotImplementedError( - "mean() is undefined for DecayPhoton distributions") + "mean() is undefined for DecaySpectrum distributions") def combine_distributions( diff --git a/src/distribution.cpp b/src/distribution.cpp index c3ec71268df..e90907afe35 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -760,7 +760,7 @@ UPtrDist distribution_from_xml(pugi::xml_node node) } else if (type == "mixture") { dist = UPtrDist {new Mixture(node)}; } else if (type == "decay_photon") { - dist = UPtrDist {new DecayPhotonMixture(node)}; + dist = UPtrDist {new DecaySpectrum(node)}; } else if (type == "muir") { openmc::fatal_error( "'muir' distributions are now specified using the openmc.stats.muir() " diff --git a/src/r2s_source.cpp b/src/r2s_source.cpp index 96afd0635ac..ea7fc3c6a1d 100644 --- a/src/r2s_source.cpp +++ b/src/r2s_source.cpp @@ -15,10 +15,10 @@ namespace openmc { //============================================================================== -// DecayPhotonMixture implementation +// DecaySpectrum implementation //============================================================================== -DecayPhotonMixture::DecayPhotonMixture( +DecaySpectrum::DecaySpectrum( vector dists, vector weights) : dists_(std::move(dists)) { @@ -31,11 +31,11 @@ DecayPhotonMixture::DecayPhotonMixture( di_.assign(probs); } -DecayPhotonMixture::DecayPhotonMixture(pugi::xml_node node) +DecaySpectrum::DecaySpectrum(pugi::xml_node node) { // Read the region volume [cm^3] needed for absolute emission rate if (!check_for_node(node, "volume")) - fatal_error("DecayPhotonMixture: 'volume' attribute is required."); + fatal_error("DecaySpectrum: 'volume' attribute is required."); double volume = std::stod(get_node_value(node, "volume")); // Read nuclide names and atom densities from XML @@ -79,18 +79,18 @@ DecayPhotonMixture::DecayPhotonMixture(pugi::xml_node node) di_.assign(probs); } -std::pair DecayPhotonMixture::sample(uint64_t* seed) const +std::pair DecaySpectrum::sample(uint64_t* seed) const { size_t idx = di_.sample(seed); return dists_[idx]->sample(seed); } -double DecayPhotonMixture::integral() const +double DecaySpectrum::integral() const { return integral_; } -double DecayPhotonMixture::sample_unbiased(uint64_t* seed) const +double DecaySpectrum::sample_unbiased(uint64_t* seed) const { size_t idx = di_.sample(seed); return dists_[idx]->sample(seed).first; diff --git a/src/source.cpp b/src/source.cpp index 22492083efd..9511b4675e5 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -357,7 +357,7 @@ IndependentSource::IndependentSource(pugi::xml_node node) : Source(node) // For decay photon sources the C++ computes the absolute photon emission // rate [photons/s] as the distribution integral. Use it as the source // strength so that sources in more active regions are sampled more often. - if (dynamic_cast(energy_.get())) { + if (dynamic_cast(energy_.get())) { strength_ = energy_->integral(); } } else { From 78cf20f5264b63320e27ee3ddd7bd6b5b2431de1 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 16 Mar 2026 00:23:03 -0500 Subject: [PATCH 12/37] Update docstring --- openmc/stats/univariate.py | 27 +++++++++++++-------------- 1 file changed, 13 insertions(+), 14 deletions(-) diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py index 926381a4851..46ed94d6e15 100644 --- a/openmc/stats/univariate.py +++ b/openmc/stats/univariate.py @@ -2201,32 +2201,31 @@ def clip(self, tolerance: float = 1e-6, inplace: bool = False) -> Mixture: class DecaySpectrum(Univariate): """Energy distribution from decay photon spectra of a mixture of nuclides. - This distribution stores nuclide names, their atom densities, and the - volume of the region. When written to XML and read by the C++ solver, the - nuclide names are resolved against the depletion chain to obtain the decay - photon energy spectra and decay constants. The resulting distribution is a - mixture of per-nuclide photon spectra weighted by absolute activity - (atom_density * 1e24 * volume * decay_constant * photons_per_decay). - The volume is necessary so that the C++ solver can compute the total photon - emission rate in [photons/s], which is used as the source strength. + This distribution stores nuclide names, their atom densities, and the volume + of the region. When written to XML and read by the C++ solver, the nuclide + names are resolved against the depletion chain to obtain the decay photon + energy spectra and decay constants. The resulting distribution is a mixture + of per-nuclide photon spectra weighted by absolute activity. The volume is + necessary so that the C++ solver can compute the total photon emission rate + in [photons/s], which is used as the source strength. .. versionadded:: 0.15.4 Parameters ---------- nuclides : dict - Dictionary mapping nuclide name (str) to atom density (float) in - units of atom/b-cm. + Dictionary mapping nuclide name (str) to atom density (float) in units + of [atom/b-cm]. volume : float - Volume of the source region in cm\ :sup:`3`. Used together with atom - densities to compute the absolute photon emission rate. + Volume of the source region in [cm³]. Used together with atom densities + to compute the absolute photon emission rate. Attributes ---------- nuclides : dict - Dictionary mapping nuclide name to atom density in atom/b-cm. + Dictionary mapping nuclide name to atom density in [atom/b-cm]. volume : float - Volume of the source region in cm\ :sup:`3`. + Volume of the source region in [cm³]. """ From 3411fa1068adb32bd4b4e6e202256e864d562d73 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 19 Mar 2026 07:23:49 -0500 Subject: [PATCH 13/37] Complete unimplemented methods by converting to mixture --- openmc/stats/univariate.py | 109 ++++++++++++++++++++++++++++++++++--- 1 file changed, 100 insertions(+), 9 deletions(-) diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py index 46ed94d6e15..5b44fbadc8b 100644 --- a/openmc/stats/univariate.py +++ b/openmc/stats/univariate.py @@ -2231,6 +2231,7 @@ class DecaySpectrum(Univariate): def __init__(self, nuclides: dict[str, float], volume: float): super().__init__(bias=None) + self._dist_cache = None self.nuclides = nuclides self.volume = volume @@ -2249,6 +2250,7 @@ def nuclides(self, nuclides): cv.check_type(f'atom density for {name}', density, Real) cv.check_greater_than(f'atom density for {name}', density, 0.0) self._nuclides = dict(nuclides) + self._dist_cache = None @property def volume(self): @@ -2259,6 +2261,43 @@ def volume(self, volume): cv.check_type('volume', volume, Real) cv.check_greater_than('volume', volume, 0.0) self._volume = float(volume) + self._dist_cache = None + + def to_distribution(self): + """Convert to a concrete distribution using decay chain data. + + Builds a combined photon energy distribution by looking up each nuclide + in the depletion chain via :func:`openmc.data.decay_photon_energy` and + weighting by absolute atom count (``density * 1e24 * volume``). The + result is cached on the object; the cache is invalidated automatically + when :attr:`nuclides` or :attr:`volume` are reassigned. + + Requires ``openmc.config['chain_file']`` to be set. + + Returns + ------- + openmc.stats.Univariate or None + Combined photon energy distribution, or ``None`` if no nuclide in + :attr:`nuclides` has a photon source in the chain. + + """ + if self._dist_cache is not None: + return self._dist_cache + + import openmc.data + dists = [] + weights = [] + for name, density in self.nuclides.items(): + dist = openmc.data.decay_photon_energy(name) + if dist is not None: + dists.append(dist) + weights.append(density * 1e24 * self.volume) + + if not dists: + return None + + self._dist_cache = combine_distributions(dists, weights) + return self._dist_cache def to_xml_element(self, element_name: str): """Return XML representation of the decay photon distribution @@ -2307,19 +2346,35 @@ def from_xml_element(cls, elem: ET.Element): return cls(nuclides, volume) def _sample_unbiased(self, n_samples=1, seed=None): - raise NotImplementedError( - "DecaySpectrum distributions can only be sampled by the C++ solver") + dist = self.to_distribution() + if dist is None: + raise RuntimeError( + "DecaySpectrum._sample_unbiased requires chain data but none " + "was found. Ensure openmc.config['chain_file'] is set and the " + "chain contains photon sources for the nuclides present." + ) + return dist.sample(n_samples, seed)[0] def integral(self): """Return integral of the distribution + Returns the total photon emission rate in [photons/s] by delegating to + :meth:`to_distribution`. Returns ``0.0`` when no chain data is + available (e.g., ``openmc.config['chain_file']`` is not set). + Returns ------- float - Integral of the distribution (requires chain data, returns 0.0 - when chain data is not available) + Total photon emission rate in [photons/s], or ``0.0`` if chain + data is unavailable. """ - return 0.0 + try: + dist = self.to_distribution() + except Exception: + return 0.0 + if dist is None: + return 0.0 + return dist.integral() def clip(self, tolerance: float = 1e-6, inplace: bool = False): """Remove nuclides with negligible contribution to photon emission. @@ -2356,12 +2411,48 @@ def support(self): return (0.0, np.inf) def evaluate(self, x): - raise NotImplementedError( - "evaluate() is undefined for DecaySpectrum distributions") + """Evaluate the probability density at a given value. + + Delegates to the combined distribution built from chain data. Raises + ``NotImplementedError`` if the combined distribution is a + :class:`~openmc.stats.Mixture` (which does not support + ``evaluate()``). + + Parameters + ---------- + x : float + Value at which to evaluate the PDF. + + Returns + ------- + float + Probability density at *x*. + """ + dist = self.to_distribution() + if dist is None: + raise RuntimeError( + "DecaySpectrum.evaluate requires chain data. Ensure " + "openmc.config['chain_file'] is set." + ) + return dist.evaluate(x) def mean(self): - raise NotImplementedError( - "mean() is undefined for DecaySpectrum distributions") + """Return the mean of the distribution. + + Delegates to the combined distribution built from chain data. + + Returns + ------- + float + Mean photon energy in [eV]. + """ + dist = self.to_distribution() + if dist is None: + raise RuntimeError( + "DecaySpectrum.mean requires chain data. Ensure " + "openmc.config['chain_file'] is set." + ) + return dist.mean() def combine_distributions( From 96d986120168560a2efbf80b96b74ccb5304bf24 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 19 Mar 2026 07:42:38 -0500 Subject: [PATCH 14/37] Clip based on photon emission rate --- openmc/deplete/r2s.py | 3 +++ openmc/stats/univariate.py | 48 +++++++++++++++++++++++++++----------- 2 files changed, 38 insertions(+), 13 deletions(-) diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index 2680533a86a..33e12dfbea5 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -649,6 +649,9 @@ def _create_photon_sources(self, time_index, work_items): continue energy = openmc.stats.DecaySpectrum(nuclides, vol) + energy.clip(inplace=True) + if not energy.nuclides: + continue if mesh_based: domains = [openmc.Material(material_id=domain_id)] diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py index 5b44fbadc8b..c44399207b6 100644 --- a/openmc/stats/univariate.py +++ b/openmc/stats/univariate.py @@ -13,7 +13,7 @@ import scipy import openmc.checkvalue as cv -from openmc.data import atomic_mass, NEUTRON_MASS +from openmc.data import atomic_mass, NEUTRON_MASS, decay_photon_energy from .._xml import get_elem_list, get_text from ..mixin import EqualityMixin @@ -2284,11 +2284,10 @@ def to_distribution(self): if self._dist_cache is not None: return self._dist_cache - import openmc.data dists = [] weights = [] for name, density in self.nuclides.items(): - dist = openmc.data.decay_photon_energy(name) + dist = decay_photon_energy(name) if dist is not None: dists.append(dist) weights.append(density * 1e24 * self.volume) @@ -2376,33 +2375,56 @@ def integral(self): return 0.0 return dist.integral() - def clip(self, tolerance: float = 1e-6, inplace: bool = False): + def clip(self, tolerance: float = 1e-9, inplace: bool = False): """Remove nuclides with negligible contribution to photon emission. - Nuclides are sorted by their atom density and those contributing the - least are removed until the cumulative removed fraction exceeds the - tolerance. + Nuclides that are stable or have no photon source in the depletion + chain are removed unconditionally. The remaining nuclides are ranked + by their photon emission rate (proportional to + ``atom_density * decay_constant * photon_yield``) and the least + important are discarded until the cumulative discarded fraction of the + total emission rate exceeds *tolerance*. + + Requires ``openmc.config['chain_file']`` to be set. Parameters ---------- tolerance : float - Maximum fraction of total atom density that can be discarded. + Maximum fraction of total photon emission rate that may be + discarded. inplace : bool Whether to modify the current object in-place or return a new one. Returns ------- openmc.stats.DecaySpectrum - Distribution with low-density nuclides removed + Distribution with negligible nuclides removed. """ - densities = np.array(list(self.nuclides.values())) - names = list(self.nuclides.keys()) - indices = _intensity_clip(densities, tolerance=tolerance) - new_nuclides = {names[i]: densities[i] for i in indices} + # Compute per-nuclide emission rate; drop non-emitters + emitting_names = [] + emitting_densities = [] + rates = [] + for name, density in self.nuclides.items(): + dist = decay_photon_energy(name) + if dist is None: + continue + rate = density * self.volume * dist.integral() + emitting_names.append(name) + emitting_densities.append(density) + rates.append(rate) + + if not emitting_names: + new_nuclides = {} + else: + indices = _intensity_clip(rates, tolerance=tolerance) + new_nuclides = { + emitting_names[i]: emitting_densities[i] for i in indices + } if inplace: self._nuclides = new_nuclides + self._dist_cache = None return self return type(self)(new_nuclides, self.volume) From 27303e947fbca77a94e81de4d7647eecf84798fd Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 19 Mar 2026 07:45:37 -0500 Subject: [PATCH 15/37] Use 'decay_spectrum' in XML --- openmc/stats/univariate.py | 4 ++-- src/distribution.cpp | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py index c44399207b6..63fb518919c 100644 --- a/openmc/stats/univariate.py +++ b/openmc/stats/univariate.py @@ -124,7 +124,7 @@ def from_xml_element(cls, elem): return Legendre.from_xml_element(elem) elif distribution == 'mixture': return Mixture.from_xml_element(elem) - elif distribution == 'decay_photon': + elif distribution == 'decay_spectrum': return DecaySpectrum.from_xml_element(elem) @abstractmethod @@ -2313,7 +2313,7 @@ def to_xml_element(self, element_name: str): """ element = ET.Element(element_name) - element.set("type", "decay_photon") + element.set("type", "decay_spectrum") element.set("volume", str(self.volume)) for name, density in self.nuclides.items(): nuclide_elem = ET.SubElement(element, "nuclide") diff --git a/src/distribution.cpp b/src/distribution.cpp index e90907afe35..04b851cbeec 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -759,7 +759,7 @@ UPtrDist distribution_from_xml(pugi::xml_node node) dist = UPtrDist {new Tabular(node)}; } else if (type == "mixture") { dist = UPtrDist {new Mixture(node)}; - } else if (type == "decay_photon") { + } else if (type == "decay_spectrum") { dist = UPtrDist {new DecaySpectrum(node)}; } else if (type == "muir") { openmc::fatal_error( From fca5d5e550ef3ab2ff0ea17f7fd98cdbdd7cbe2c Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 19 Mar 2026 07:53:09 -0500 Subject: [PATCH 16/37] Move DecaySpectrum to distribution.h/cpp --- CMakeLists.txt | 1 - include/openmc/distribution.h | 40 ++++++++++++++ include/openmc/r2s_source.h | 57 -------------------- src/distribution.cpp | 84 ++++++++++++++++++++++++++++- src/r2s_source.cpp | 99 ----------------------------------- src/source.cpp | 1 - 6 files changed, 123 insertions(+), 159 deletions(-) delete mode 100644 include/openmc/r2s_source.h delete mode 100644 src/r2s_source.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 62bffc8110d..9fe133a22e3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -415,7 +415,6 @@ list(APPEND libopenmc_SOURCES src/ray.cpp src/reaction.cpp src/reaction_product.cpp - src/r2s_source.cpp src/scattdata.cpp src/secondary_correlated.cpp src/secondary_kalbach.cpp diff --git a/include/openmc/distribution.h b/include/openmc/distribution.h index cd81d8801d6..62375e2e0b9 100644 --- a/include/openmc/distribution.h +++ b/include/openmc/distribution.h @@ -407,6 +407,46 @@ class Mixture : public Distribution { double integral_; //!< Integral of distribution }; +//============================================================================== +// DecaySpectrum — non-owning mixture of decay photon distributions +//============================================================================== + +//! Energy distribution formed by mixing multiple decay photon spectra. +//! +//! Unlike the general Mixture distribution, this class holds non-owning +//! pointers to the component distributions (which live in +//! data::chain_nuclides). Each component is weighted by the activity +//! (atoms * decay_constant) of the corresponding nuclide. + +class DecaySpectrum : public Distribution { +public: + //! Construct from non-owning distribution pointers and weights + //! + //! \param dists Non-owning pointers to component distributions + //! \param weights Activity-based weights for each component. The Mixture + //! probability for component i is weights[i] * dists[i]->integral() + //! (the integral encodes photons-per-decay). + DecaySpectrum(vector dists, vector weights); + + //! Construct from an XML node containing nuclide names and atom densities. + //! + //! Reads child ```` elements with ``name`` and ``density`` + //! attributes, resolves them against the loaded depletion chain, and + //! constructs the mixed distribution. + explicit DecaySpectrum(pugi::xml_node node); + + std::pair sample(uint64_t* seed) const override; + double integral() const override; + +protected: + double sample_unbiased(uint64_t* seed) const override; + +private: + vector dists_; //!< Non-owning component distributions + DiscreteIndex di_; //!< Discrete index for component selection + double integral_; //!< Total photon emission rate +}; + } // namespace openmc #endif // OPENMC_DISTRIBUTION_H diff --git a/include/openmc/r2s_source.h b/include/openmc/r2s_source.h deleted file mode 100644 index 7a4ba0410b4..00000000000 --- a/include/openmc/r2s_source.h +++ /dev/null @@ -1,57 +0,0 @@ -//! \file r2s_source.h -//! \brief Decay photon source generation for R2S calculations - -#ifndef OPENMC_R2S_SOURCE_H -#define OPENMC_R2S_SOURCE_H - -#include // for size_t -#include // for int32_t - -#include "openmc/distribution.h" -#include "openmc/vector.h" - -namespace openmc { - -//============================================================================== -// DecaySpectrum — non-owning mixture of decay photon distributions -//============================================================================== - -//! Energy distribution formed by mixing multiple decay photon spectra. -//! -//! Unlike the general Mixture distribution, this class holds non-owning -//! pointers to the component distributions (which live in -//! data::chain_nuclides). Each component is weighted by the activity -//! (atoms * decay_constant) of the corresponding nuclide. - -class DecaySpectrum : public Distribution { -public: - //! Construct from non-owning distribution pointers and weights - //! - //! \param dists Non-owning pointers to component distributions - //! \param weights Activity-based weights for each component. The Mixture - //! probability for component i is weights[i] * dists[i]->integral() - //! (the integral encodes photons-per-decay). - DecaySpectrum(vector dists, vector weights); - - //! Construct from an XML node containing nuclide names and atom densities. - //! - //! Reads child ```` elements with ``name`` and ``density`` - //! attributes, resolves them against the loaded depletion chain, and - //! constructs the mixed distribution. - explicit DecaySpectrum(pugi::xml_node node); - - std::pair sample(uint64_t* seed) const override; - double integral() const override; - -protected: - double sample_unbiased(uint64_t* seed) const override; - -private: - vector dists_; //!< Non-owning component distributions - DiscreteIndex di_; //!< Discrete index for component selection - double integral_; //!< Total photon emission rate -}; - -} // namespace openmc - -#endif // OPENMC_R2S_SOURCE_H diff --git a/src/distribution.cpp b/src/distribution.cpp index 04b851cbeec..bee7240faaf 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -8,10 +8,10 @@ #include // for runtime_error #include // for string, stod +#include "openmc/chain.h" #include "openmc/constants.h" #include "openmc/error.h" #include "openmc/math_functions.h" -#include "openmc/r2s_source.h" #include "openmc/random_dist.h" #include "openmc/random_lcg.h" #include "openmc/xml_interface.h" @@ -771,4 +771,86 @@ UPtrDist distribution_from_xml(pugi::xml_node node) return dist; } +//============================================================================== +// DecaySpectrum implementation +//============================================================================== + +DecaySpectrum::DecaySpectrum( + vector dists, vector weights) + : dists_(std::move(dists)) +{ + vector probs; + probs.reserve(dists_.size()); + for (size_t i = 0; i < dists_.size(); ++i) { + probs.push_back(weights[i] * dists_[i]->integral()); + } + integral_ = std::accumulate(probs.begin(), probs.end(), 0.0); + di_.assign(probs); +} + +DecaySpectrum::DecaySpectrum(pugi::xml_node node) +{ + // Read the region volume [cm^3] needed for absolute emission rate + if (!check_for_node(node, "volume")) + fatal_error("DecaySpectrum: 'volume' attribute is required."); + double volume = std::stod(get_node_value(node, "volume")); + + // Read nuclide names and atom densities from XML + vector dists; + vector weights; + + for (auto nuclide_node : node.children("nuclide")) { + std::string name = get_node_value(nuclide_node, "name"); + double density = std::stod(get_node_value(nuclide_node, "density")); + + // Look up nuclide in the depletion chain + auto it = data::chain_nuclide_map.find(name); + if (it == data::chain_nuclide_map.end()) + continue; + + const auto& chain_nuc = data::chain_nuclides[it->second]; + const Distribution* photon_dist = chain_nuc->photon_energy(); + if (!photon_dist) + continue; + + // Weight = atom_count * decay_constant + // atom_count = atom_density [atom/b-cm] * 1e24 [b-cm/cm^3] * volume [cm^3] + double atom_count = density * 1.0e24 * volume; + double weight = atom_count * chain_nuc->decay_constant(); + if (weight <= 0.0) + continue; + + dists.push_back(photon_dist); + weights.push_back(weight); + } + + dists_ = std::move(dists); + + // Build alias table from weighted probabilities + vector probs; + probs.reserve(dists_.size()); + for (size_t i = 0; i < dists_.size(); ++i) { + probs.push_back(weights[i] * dists_[i]->integral()); + } + integral_ = std::accumulate(probs.begin(), probs.end(), 0.0); + di_.assign(probs); +} + +std::pair DecaySpectrum::sample(uint64_t* seed) const +{ + size_t idx = di_.sample(seed); + return dists_[idx]->sample(seed); +} + +double DecaySpectrum::integral() const +{ + return integral_; +} + +double DecaySpectrum::sample_unbiased(uint64_t* seed) const +{ + size_t idx = di_.sample(seed); + return dists_[idx]->sample(seed).first; +} + } // namespace openmc diff --git a/src/r2s_source.cpp b/src/r2s_source.cpp deleted file mode 100644 index ea7fc3c6a1d..00000000000 --- a/src/r2s_source.cpp +++ /dev/null @@ -1,99 +0,0 @@ -//! \file r2s_source.cpp -//! \brief Decay photon source generation for R2S calculations - -#include "openmc/r2s_source.h" - -#include // for accumulate -#include - -#include "openmc/chain.h" -#include "openmc/error.h" -#include "openmc/memory.h" -#include "openmc/vector.h" -#include "openmc/xml_interface.h" - -namespace openmc { - -//============================================================================== -// DecaySpectrum implementation -//============================================================================== - -DecaySpectrum::DecaySpectrum( - vector dists, vector weights) - : dists_(std::move(dists)) -{ - vector probs; - probs.reserve(dists_.size()); - for (size_t i = 0; i < dists_.size(); ++i) { - probs.push_back(weights[i] * dists_[i]->integral()); - } - integral_ = std::accumulate(probs.begin(), probs.end(), 0.0); - di_.assign(probs); -} - -DecaySpectrum::DecaySpectrum(pugi::xml_node node) -{ - // Read the region volume [cm^3] needed for absolute emission rate - if (!check_for_node(node, "volume")) - fatal_error("DecaySpectrum: 'volume' attribute is required."); - double volume = std::stod(get_node_value(node, "volume")); - - // Read nuclide names and atom densities from XML - vector dists; - vector weights; - - for (auto nuclide_node : node.children("nuclide")) { - std::string name = get_node_value(nuclide_node, "name"); - double density = std::stod(get_node_value(nuclide_node, "density")); - - // Look up nuclide in the depletion chain - auto it = data::chain_nuclide_map.find(name); - if (it == data::chain_nuclide_map.end()) - continue; - - const auto& chain_nuc = data::chain_nuclides[it->second]; - const Distribution* photon_dist = chain_nuc->photon_energy(); - if (!photon_dist) - continue; - - // Weight = atom_count * decay_constant - // atom_count = atom_density [atom/b-cm] * 1e24 [b-cm/cm^3] * volume [cm^3] - double atom_count = density * 1.0e24 * volume; - double weight = atom_count * chain_nuc->decay_constant(); - if (weight <= 0.0) - continue; - - dists.push_back(photon_dist); - weights.push_back(weight); - } - - dists_ = std::move(dists); - - // Build alias table from weighted probabilities - vector probs; - probs.reserve(dists_.size()); - for (size_t i = 0; i < dists_.size(); ++i) { - probs.push_back(weights[i] * dists_[i]->integral()); - } - integral_ = std::accumulate(probs.begin(), probs.end(), 0.0); - di_.assign(probs); -} - -std::pair DecaySpectrum::sample(uint64_t* seed) const -{ - size_t idx = di_.sample(seed); - return dists_[idx]->sample(seed); -} - -double DecaySpectrum::integral() const -{ - return integral_; -} - -double DecaySpectrum::sample_unbiased(uint64_t* seed) const -{ - size_t idx = di_.sample(seed); - return dists_[idx]->sample(seed).first; -} - -} // namespace openmc diff --git a/src/source.cpp b/src/source.cpp index 9511b4675e5..f2b1015e324 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -27,7 +27,6 @@ #include "openmc/message_passing.h" #include "openmc/mgxs_interface.h" #include "openmc/nuclide.h" -#include "openmc/r2s_source.h" #include "openmc/random_lcg.h" #include "openmc/search.h" #include "openmc/settings.h" From 5c6da0b0beb32ea671df03079b5fedc801f5c55f Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 19 Mar 2026 07:54:43 -0500 Subject: [PATCH 17/37] Remove unused code --- include/openmc/source.h | 13 ------------- 1 file changed, 13 deletions(-) diff --git a/include/openmc/source.h b/include/openmc/source.h index 7f6c36ec0b9..1ef7eba2b24 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -73,16 +73,6 @@ class Source { // Methods that can be overridden virtual double strength() const { return strength_; } - //! Set the source strength - void set_strength(double strength) { strength_ = strength; } - - //! Set the domain type and IDs for source constraints - void set_domain(DomainType domain_type, std::unordered_set ids) - { - domain_type_ = domain_type; - domain_ids_ = std::move(ids); - } - //! Sample a source site and apply constraints // //! \param[inout] seed Pseudorandom seed pointer @@ -148,9 +138,6 @@ class IndependentSource : public Source { // Properties ParticleType particle_type() const { return particle_; } - //! Set the particle type - void set_particle(ParticleType p) { particle_ = p; } - // Make observing pointers available SpatialDistribution* space() const { return space_.get(); } UnitSphereDistribution* angle() const { return angle_.get(); } From 662f7f771a87590a56919fe2a2fed881b5709e9e Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 19 Mar 2026 08:13:51 -0500 Subject: [PATCH 18/37] Simplification in _create_photon_sources --- openmc/deplete/r2s.py | 33 ++++++++------------------------- 1 file changed, 8 insertions(+), 25 deletions(-) diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index 33e12dfbea5..6ef6eecdd71 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -621,49 +621,32 @@ def _create_photon_sources(self, time_index, work_items): materials = self.results['activation_materials'] mesh_based = self.method == 'mesh-based' - nuclide_names = list(step_result.index_nuc.keys()) - sources = [] for item in work_items: if mesh_based: index_mat, domain_id, bbox = item original_mat = materials[index_mat] + domain = openmc.Material(material_id=domain_id) else: cell, original_mat, bbox = item - domain_id = cell.id - - mat_key = str(original_mat.id) - vol = step_result.volume[mat_key] - - # Get atom counts and convert to atom/b-cm - mat_idx = step_result.index_mat[mat_key] - atom_counts = step_result.data[mat_idx, :] - atom_densities = atom_counts / vol * 1e-24 + domain = cell - # Build nuclide dict with only non-zero densities - nuclides = { - name: dens for name, dens in zip(nuclide_names, atom_densities) - if dens > 0.0 - } + activated_mat = step_result.get_material(str(original_mat.id)) + nuclides = activated_mat.get_nuclide_atom_densities() if not nuclides: continue - energy = openmc.stats.DecaySpectrum(nuclides, vol) + energy = openmc.stats.DecaySpectrum(nuclides, activated_mat.volume) energy.clip(inplace=True) if not energy.nuclides: continue - if mesh_based: - domains = [openmc.Material(material_id=domain_id)] - else: - domains = [openmc.Cell(cell_id=domain_id)] - source = openmc.IndependentSource( + sources.append(openmc.IndependentSource( space=openmc.stats.Box(bbox.lower_left, bbox.upper_right), energy=energy, particle='photon', - constraints={'domains': domains}, - ) - sources.append(source) + constraints={'domains': [domain]}, + )) return sources From ab915de11c424d2f3fdd98d8bb8168d487b5b1e6 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 19 Mar 2026 08:36:36 -0500 Subject: [PATCH 19/37] Cache photon spectrum integral --- openmc/stats/univariate.py | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py index 63fb518919c..a4c6ba8abb7 100644 --- a/openmc/stats/univariate.py +++ b/openmc/stats/univariate.py @@ -2,6 +2,7 @@ from abc import ABC, abstractmethod from collections import defaultdict from collections.abc import Iterable, Sequence +from functools import cache from math import sqrt, pi, exp, log from numbers import Real from warnings import warn @@ -2375,6 +2376,13 @@ def integral(self): return 0.0 return dist.integral() + @staticmethod + @cache + def _photon_integral(nuclide: str) -> float | None: + """Return the per-atom photon emission integral for a nuclide""" + dist = decay_photon_energy(nuclide) + return dist.integral() if dist is not None else None + def clip(self, tolerance: float = 1e-9, inplace: bool = False): """Remove nuclides with negligible contribution to photon emission. @@ -2406,13 +2414,12 @@ def clip(self, tolerance: float = 1e-9, inplace: bool = False): emitting_densities = [] rates = [] for name, density in self.nuclides.items(): - dist = decay_photon_energy(name) - if dist is None: + integral = DecaySpectrum._photon_integral(name) + if integral is None: continue - rate = density * self.volume * dist.integral() emitting_names.append(name) emitting_densities.append(density) - rates.append(rate) + rates.append(density * self.volume * integral) if not emitting_names: new_nuclides = {} From a22e84f9cda18ed4c7115c96c49d5cc4335a77ea Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 19 Mar 2026 08:52:12 -0500 Subject: [PATCH 20/37] Add tests for DecaySpectrum --- tests/unit_tests/test_stats.py | 132 ++++++++++++++++++++++++++++++++- 1 file changed, 131 insertions(+), 1 deletion(-) diff --git a/tests/unit_tests/test_stats.py b/tests/unit_tests/test_stats.py index ca961c8b0f1..88195b159d5 100644 --- a/tests/unit_tests/test_stats.py +++ b/tests/unit_tests/test_stats.py @@ -1,10 +1,11 @@ from math import pi +from pathlib import Path import numpy as np import pytest import openmc import openmc.stats -from openmc.stats.univariate import _INTERPOLATION_SCHEMES +from openmc.stats.univariate import _INTERPOLATION_SCHEMES, DecaySpectrum from scipy.integrate import trapezoid from tests.unit_tests import assert_sample_mean @@ -1013,3 +1014,132 @@ def test_fusion_spectrum_invalid(): # Temperature above 100 keV should raise an error with pytest.raises(ValueError): openmc.stats.fusion_neutron_spectrum(101e3, 'DT') + + +@pytest.fixture(autouse=False) +def decay_spectrum_chain(): + """Set chain_file for the duration of a test and clear the _photon_integral + cache so results from a different chain don't bleed across tests.""" + CHAIN_FILE = (Path(__file__).parents[1] / 'chain_simple.xml').resolve() + DecaySpectrum._photon_integral.cache_clear() + with openmc.config.patch('chain_file', CHAIN_FILE): + yield + DecaySpectrum._photon_integral.cache_clear() + + +def test_decay_spectrum_construction(): + nuclides = {'I135': 1.5e-3, 'Xe135': 8.2e-4} + d = openmc.stats.DecaySpectrum(nuclides, volume=100.0) + assert d.nuclides == nuclides + assert d.volume == pytest.approx(100.0) + assert len(d) == 2 + + +def test_decay_spectrum_validation(): + # nuclides must be a dict + with pytest.raises(TypeError): + openmc.stats.DecaySpectrum(['I135'], volume=1.0) + + # densities must be > 0 + with pytest.raises(ValueError): + openmc.stats.DecaySpectrum({'I135': -1.0}, volume=1.0) + + # volume must be > 0 + with pytest.raises(ValueError): + openmc.stats.DecaySpectrum({'I135': 1e-3}, volume=-1.0) + + with pytest.raises(ValueError): + openmc.stats.DecaySpectrum({'I135': 1e-3}, volume=0.0) + + +def test_decay_spectrum_xml_roundtrip(): + nuclides = {'I135': 1.5e-3, 'Xe135': 8.2e-4} + d = openmc.stats.DecaySpectrum(nuclides, volume=100.0) + + elem = d.to_xml_element('energy') + assert elem.get('type') == 'decay_spectrum' + assert float(elem.get('volume')) == pytest.approx(100.0) + + # Round-trip via DecaySpectrum.from_xml_element + d2 = openmc.stats.DecaySpectrum.from_xml_element(elem) + assert d2.nuclides == nuclides + assert d2.volume == pytest.approx(100.0) + + # Round-trip via the Univariate dispatcher + d3 = openmc.stats.Univariate.from_xml_element(elem) + assert isinstance(d3, openmc.stats.DecaySpectrum) + assert d3 == d + + +def test_decay_spectrum_to_distribution(decay_spectrum_chain): + # Single emitting nuclide -> concrete distribution, not None + d = openmc.stats.DecaySpectrum({'I135': 1e-3}, volume=10.0) + dist = d.to_distribution() + assert dist is not None + + # Result is cached on second call + dist2 = d.to_distribution() + assert dist2 is dist + + # Nuclide with no photon source -> None + d_stable = openmc.stats.DecaySpectrum({'Xe136': 1e-3}, volume=10.0) + assert d_stable.to_distribution() is None + + # Mixture of emitters -> non-None combined distribution + d_mix = openmc.stats.DecaySpectrum( + {'I135': 1e-3, 'Xe135': 5e-4}, volume=10.0 + ) + dist_mix = d_mix.to_distribution() + assert dist_mix is not None + + +def test_decay_spectrum_integral(decay_spectrum_chain): + # For an emitting nuclide, integral should be > 0 + d = openmc.stats.DecaySpectrum({'I135': 1e-3}, volume=10.0) + assert d.integral() > 0.0 + + # Proportional to density: doubling density doubles integral + d2 = openmc.stats.DecaySpectrum({'I135': 2e-3}, volume=10.0) + assert d2.integral() == pytest.approx(2.0 * d.integral()) + + # Proportional to volume + d3 = openmc.stats.DecaySpectrum({'I135': 1e-3}, volume=20.0) + assert d3.integral() == pytest.approx(2.0 * d.integral()) + + # Pure non-emitter -> 0.0 + d_stable = openmc.stats.DecaySpectrum({'Xe136': 1e-3}, volume=10.0) + assert d_stable.integral() == pytest.approx(0.0) + + +def test_decay_spectrum_clip(decay_spectrum_chain): + # Stable / non-emitting nuclides are removed unconditionally + d = openmc.stats.DecaySpectrum( + {'I135': 1e-3, 'Xe135': 5e-4, 'Xe136': 1.0, 'Cs135': 1.0}, + volume=10.0, + ) + d_clip = d.clip() + assert 'Xe136' not in d_clip.nuclides + assert 'Cs135' not in d_clip.nuclides + assert 'I135' in d_clip.nuclides + assert 'Xe135' in d_clip.nuclides + # Original is unchanged + assert 'Xe136' in d.nuclides + + # inplace=True modifies and returns the same object + d_same = d.clip(inplace=True) + assert d_same is d + assert 'Xe136' not in d.nuclides + + # A nuclide with negligible emission rate is removed by tolerance clipping. + # U235 has a very small integral (~4e-17 Bq/atom) compared with I135 (~4e-5) + d_tight = openmc.stats.DecaySpectrum( + {'I135': 1e-3, 'U235': 1e-3}, volume=10.0 + ) + d_tight_clip = d_tight.clip(tolerance=1e-9) + assert 'U235' not in d_tight_clip.nuclides + assert 'I135' in d_tight_clip.nuclides + + # All non-emitters -> empty nuclides dict + d_empty = openmc.stats.DecaySpectrum({'Xe136': 1e-3}, volume=10.0) + d_empty.clip(inplace=True) + assert d_empty.nuclides == {} From 95cf3cbbb651aa396c0013b8559a8a149e41109f Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 23 Apr 2026 00:34:05 -0500 Subject: [PATCH 21/37] Allow parent_nuclide to be set when using DecaySpectrum --- include/openmc/chain.h | 2 ++ include/openmc/distribution.h | 45 +++++++++++++++++++++----- src/chain.cpp | 8 +++++ src/distribution.cpp | 60 ++++++++++++++++++++++++----------- src/finalize.cpp | 1 + src/initialize.cpp | 14 ++++---- src/source.cpp | 17 +++++++--- 7 files changed, 110 insertions(+), 37 deletions(-) diff --git a/include/openmc/chain.h b/include/openmc/chain.h index 6f58303585a..e01f2a01c18 100644 --- a/include/openmc/chain.h +++ b/include/openmc/chain.h @@ -101,6 +101,8 @@ extern vector> chain_nuclides; void read_chain_file_xml(); +void free_memory_chain(); + } // namespace openmc #endif // OPENMC_CHAIN_H diff --git a/include/openmc/distribution.h b/include/openmc/distribution.h index 62375e2e0b9..442dbb04d27 100644 --- a/include/openmc/distribution.h +++ b/include/openmc/distribution.h @@ -420,13 +420,27 @@ class Mixture : public Distribution { class DecaySpectrum : public Distribution { public: - //! Construct from non-owning distribution pointers and weights + //============================================================================ + // Types, aliases + + struct Sample { + double energy; + double weight; + int parent_nuclide; + }; + + //============================================================================ + // Constructors + + //! Construct from depletion chain nuclide indices and weights //! - //! \param dists Non-owning pointers to component distributions + //! \param nuclide_indices Indices of decay photon emitters in + //! data::chain_nuclides //! \param weights Activity-based weights for each component. The Mixture - //! probability for component i is weights[i] * dists[i]->integral() - //! (the integral encodes photons-per-decay). - DecaySpectrum(vector dists, vector weights); + //! probability for component i is weights[i] * + //! chain_nuclides[nuclide_indices[i]]->photon_energy()->integral() (the + //! integral encodes photons-per-decay). + DecaySpectrum(vector nuclide_indices, vector weights); //! Construct from an XML node containing nuclide names and atom densities. //! @@ -435,16 +449,31 @@ class DecaySpectrum : public Distribution { //! constructs the mixed distribution. explicit DecaySpectrum(pugi::xml_node node); + //============================================================================ + // Methods + + //! Sample a value from the distribution and return the parent nuclide index + //! \param seed Pseudorandom number seed pointer + //! \return (Sampled energy, sample weight, chain nuclide index) + Sample sample_with_parent(uint64_t* seed) const; + + //! Sample a value from the distribution + //! \param seed Pseudorandom number seed pointer + //! \return (sampled value, sample weight) std::pair sample(uint64_t* seed) const override; + double integral() const override; protected: + //! Sample a value (unbiased) from the distribution + //! \param seed Pseudorandom number seed pointer + //! \return Sampled value double sample_unbiased(uint64_t* seed) const override; private: - vector dists_; //!< Non-owning component distributions - DiscreteIndex di_; //!< Discrete index for component selection - double integral_; //!< Total photon emission rate + vector nuclide_indices_; //!< Indices of emitting nuclides in the chain + DiscreteIndex di_; //!< Discrete index for component selection + double integral_; //!< Total photon emission rate }; } // namespace openmc diff --git a/src/chain.cpp b/src/chain.cpp index e4d0324d3b5..3915c016c29 100644 --- a/src/chain.cpp +++ b/src/chain.cpp @@ -98,6 +98,8 @@ vector> chain_nuclides; void read_chain_file_xml() { + free_memory_chain(); + char* chain_file_path = std::getenv("OPENMC_CHAIN_FILE"); if (!chain_file_path) { return; @@ -120,4 +122,10 @@ void read_chain_file_xml() } } +void free_memory_chain() +{ + data::chain_nuclides.clear(); + data::chain_nuclide_map.clear(); +} + } // namespace openmc diff --git a/src/distribution.cpp b/src/distribution.cpp index bee7240faaf..301896286cb 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -775,15 +775,31 @@ UPtrDist distribution_from_xml(pugi::xml_node node) // DecaySpectrum implementation //============================================================================== -DecaySpectrum::DecaySpectrum( - vector dists, vector weights) - : dists_(std::move(dists)) +// Helper function to compute weighted probabilities for decay spectrum sampling +vector decay_spectrum_probabilities( + const vector& nuclide_indices, const vector& weights) { + if (nuclide_indices.size() != weights.size()) { + fatal_error("DecaySpectrum nuclide index and weight arrays must have the " + "same length."); + } + vector probs; - probs.reserve(dists_.size()); - for (size_t i = 0; i < dists_.size(); ++i) { - probs.push_back(weights[i] * dists_[i]->integral()); + probs.reserve(nuclide_indices.size()); + for (size_t i = 0; i < nuclide_indices.size(); ++i) { + const auto* dist = + data::chain_nuclides[nuclide_indices[i]]->photon_energy(); + probs.push_back(weights[i] * dist->integral()); } + return probs; +} + +DecaySpectrum::DecaySpectrum( + vector nuclide_indices, vector weights) + : nuclide_indices_(std::move(nuclide_indices)) +{ + vector probs = + decay_spectrum_probabilities(nuclide_indices_, weights); integral_ = std::accumulate(probs.begin(), probs.end(), 0.0); di_.assign(probs); } @@ -796,7 +812,7 @@ DecaySpectrum::DecaySpectrum(pugi::xml_node node) double volume = std::stod(get_node_value(node, "volume")); // Read nuclide names and atom densities from XML - vector dists; + vector nuclide_indices; vector weights; for (auto nuclide_node : node.children("nuclide")) { @@ -808,7 +824,8 @@ DecaySpectrum::DecaySpectrum(pugi::xml_node node) if (it == data::chain_nuclide_map.end()) continue; - const auto& chain_nuc = data::chain_nuclides[it->second]; + int nuclide_index = it->second; + const auto& chain_nuc = data::chain_nuclides[nuclide_index]; const Distribution* photon_dist = chain_nuc->photon_energy(); if (!photon_dist) continue; @@ -820,26 +837,32 @@ DecaySpectrum::DecaySpectrum(pugi::xml_node node) if (weight <= 0.0) continue; - dists.push_back(photon_dist); + nuclide_indices.push_back(nuclide_index); weights.push_back(weight); } - dists_ = std::move(dists); + nuclide_indices_ = std::move(nuclide_indices); // Build alias table from weighted probabilities - vector probs; - probs.reserve(dists_.size()); - for (size_t i = 0; i < dists_.size(); ++i) { - probs.push_back(weights[i] * dists_[i]->integral()); - } + vector probs = + decay_spectrum_probabilities(nuclide_indices_, weights); integral_ = std::accumulate(probs.begin(), probs.end(), 0.0); di_.assign(probs); } -std::pair DecaySpectrum::sample(uint64_t* seed) const +DecaySpectrum::Sample DecaySpectrum::sample_with_parent(uint64_t* seed) const { size_t idx = di_.sample(seed); - return dists_[idx]->sample(seed); + int parent_nuclide = nuclide_indices_[idx]; + const auto* dist = data::chain_nuclides[parent_nuclide]->photon_energy(); + auto [energy, weight] = dist->sample(seed); + return {energy, weight, parent_nuclide}; +} + +std::pair DecaySpectrum::sample(uint64_t* seed) const +{ + auto sample = sample_with_parent(seed); + return {sample.energy, sample.weight}; } double DecaySpectrum::integral() const @@ -849,8 +872,7 @@ double DecaySpectrum::integral() const double DecaySpectrum::sample_unbiased(uint64_t* seed) const { - size_t idx = di_.sample(seed); - return dists_[idx]->sample(seed).first; + return sample_with_parent(seed).energy; } } // namespace openmc diff --git a/src/finalize.cpp b/src/finalize.cpp index 98f1125347d..1c4b3df6811 100644 --- a/src/finalize.cpp +++ b/src/finalize.cpp @@ -44,6 +44,7 @@ void free_memory() free_memory_photon(); free_memory_settings(); free_memory_thermal(); + free_memory_chain(); library_clear(); nuclides_clear(); free_memory_source(); diff --git a/src/initialize.cpp b/src/initialize.cpp index efb462f5c01..991877b2429 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -405,6 +405,10 @@ bool read_model_xml() write_message( fmt::format("Reading model XML file '{}' ...", model_filename), 5); + // Read chain data before settings so DecaySpectrum source distributions can + // resolve nuclides while sources are constructed. + read_chain_file_xml(); + read_settings_xml(settings_root); // If other XML files are present, display warning @@ -420,9 +424,6 @@ 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( @@ -475,14 +476,15 @@ bool read_model_xml() void read_separate_xml_files() { + // Read chain data before settings so DecaySpectrum source distributions can + // resolve nuclides while sources are constructed. + read_chain_file_xml(); + read_settings_xml(); 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/source.cpp b/src/source.cpp index f2b1015e324..d81205e2085 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -419,6 +419,7 @@ SourceSite IndependentSource::sample(uint64_t* seed) const // Check for monoenergetic source above maximum particle energy auto p = particle_.transport_index(); auto energy_ptr = dynamic_cast(energy_.get()); + auto decay_spectrum = dynamic_cast(energy_.get()); if (energy_ptr) { auto energies = tensor::Tensor(energy_ptr->x().data(), energy_ptr->x().size()); @@ -429,10 +430,18 @@ SourceSite IndependentSource::sample(uint64_t* seed) const } while (true) { - // Sample energy spectrum - auto [E, E_wgt_temp] = energy_->sample(seed); - site.E = E; - E_wgt = E_wgt_temp; + // Sample energy spectrum. For decay photon sources, also get the parent + // nuclide index to store in the source site for tallying purposes. + if (decay_spectrum) { + auto sample = decay_spectrum->sample_with_parent(seed); + site.E = sample.energy; + E_wgt = sample.weight; + site.parent_nuclide = sample.parent_nuclide; + } else { + auto [E, E_wgt_temp] = energy_->sample(seed); + site.E = E; + E_wgt = E_wgt_temp; + } // Resample if energy falls above maximum particle energy if (site.E < data::energy_max[p] && From 01441f8141e1099c270d14d2ac1bd2812110e0fa Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 23 Apr 2026 12:35:08 -0500 Subject: [PATCH 22/37] Add test for DecaySpectrum + ParentNuclideFilter --- tests/unit_tests/test_source.py | 69 +++++++++++++++++++++++++++++++++ 1 file changed, 69 insertions(+) diff --git a/tests/unit_tests/test_source.py b/tests/unit_tests/test_source.py index 6eb03f38810..394c09e739d 100644 --- a/tests/unit_tests/test_source.py +++ b/tests/unit_tests/test_source.py @@ -1,5 +1,6 @@ from collections import Counter from math import pi +from pathlib import Path import openmc import openmc.lib @@ -91,6 +92,74 @@ def test_point_cloud_strengths(run_in_tmpdir, sphere_box_model): assert sampled_strength == expected_strength, f'Strength incorrect for {positions[i]}' +def test_decay_spectrum_parent_nuclide(run_in_tmpdir): + chain_file = Path('chain_decay_spectrum_parent.xml') + chain_file.write_text(""" + + + + 1000000.0 1.0 + + + + + 2000000.0 1.0 + + + +""") + + inner_sphere = openmc.Sphere(r=10.0) + outer_sphere = openmc.Sphere(r=20.0, boundary_type='vacuum') + + shell_mat = openmc.Material() + shell_mat.add_nuclide('H1', 1.0) + shell_mat.set_density('atom/b-cm', 1.0e-12) + + void_cell = openmc.Cell(region=-inner_sphere) + shell_cell = openmc.Cell(fill=shell_mat, region=+inner_sphere & -outer_sphere) + + model = openmc.Model() + model.geometry = openmc.Geometry([void_cell, shell_cell]) + model.materials = [shell_mat] + model.settings.run_mode = 'fixed source' + model.settings.photon_transport = True + model.settings.particles = 1000 + model.settings.batches = 5 + model.settings.source = openmc.IndependentSource( + particle='photon', + space=openmc.stats.Point((0.0, 0.0, 0.0)), + energy=openmc.stats.DecaySpectrum( + {'ParentA': 1.0, 'ParentB': 1.0}, + volume=1.0 + ) + ) + + tally = openmc.Tally() + tally.filters = [ + openmc.CellFilter([void_cell]), + openmc.ParticleFilter(['photon']), + openmc.EnergyFilter([0.0, 1.5e6, 2.5e6]), + openmc.ParentNuclideFilter(['ParentA', 'ParentB']) + ] + tally.scores = ['flux'] + model.tallies = [tally] + + with openmc.config.patch('chain_file', chain_file): + sp_filename = model.run() + + with openmc.StatePoint(sp_filename) as sp: + tally_out = sp.tallies[tally.id] + mean = tally_out.get_reshaped_data('mean').squeeze() + + assert mean.shape == (2, 2) + assert mean[0, 0] > 0.0 + assert mean[1, 1] > 0.0 + assert mean[0, 1] == 0.0 + assert mean[1, 0] == 0.0 + assert np.count_nonzero(mean) == 2 + + def test_source_file(): filename = 'source.h5' src = openmc.FileSource(path=filename) From 5bc8f5e9a402b85a44a06e2660350030c76a469a Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 23 Apr 2026 16:10:31 -0500 Subject: [PATCH 23/37] Update comment --- src/source.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/source.cpp b/src/source.cpp index d81205e2085..50a6640c2ca 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -353,9 +353,8 @@ IndependentSource::IndependentSource(pugi::xml_node node) : Source(node) pugi::xml_node node_dist = node.child("energy"); energy_ = distribution_from_xml(node_dist); - // For decay photon sources the C++ computes the absolute photon emission - // rate [photons/s] as the distribution integral. Use it as the source - // strength so that sources in more active regions are sampled more often. + // For decay photon sources, use the absolute photon emission rate in + // [photons/s] as the source strength if (dynamic_cast(energy_.get())) { strength_ = energy_->integral(); } From 9b3b4e3408278ed798285aed35edf563aac9b172 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 23 Apr 2026 16:49:57 -0500 Subject: [PATCH 24/37] Handle zero density nuclides for R2S DecaySpectrum --- openmc/deplete/r2s.py | 4 ++++ openmc/stats/univariate.py | 2 +- 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index 6ef6eecdd71..d9c539f5419 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -636,6 +636,10 @@ def _create_photon_sources(self, time_index, work_items): if not nuclides: continue + # Eliminate nuclides with zero density + nuclides = {nuclide: density for nuclide, density in nuclides.items() + if density > 0} + energy = openmc.stats.DecaySpectrum(nuclides, activated_mat.volume) energy.clip(inplace=True) if not energy.nuclides: diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py index a4c6ba8abb7..b0203fbbc49 100644 --- a/openmc/stats/univariate.py +++ b/openmc/stats/univariate.py @@ -2249,7 +2249,7 @@ def nuclides(self, nuclides): for name, density in nuclides.items(): cv.check_type('nuclide name', name, str) cv.check_type(f'atom density for {name}', density, Real) - cv.check_greater_than(f'atom density for {name}', density, 0.0) + cv.check_greater_than(f'atom density for {name}', density, 0.0, True) self._nuclides = dict(nuclides) self._dist_cache = None From 315953c6730000f87eae4255a6f32d2e14fdcc4d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 23 Apr 2026 18:55:27 -0500 Subject: [PATCH 25/37] Fix bug in DecaySpectrum integral --- include/openmc/distribution.h | 6 +++--- src/distribution.cpp | 8 +++++--- 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/include/openmc/distribution.h b/include/openmc/distribution.h index 442dbb04d27..cf01ebcadf9 100644 --- a/include/openmc/distribution.h +++ b/include/openmc/distribution.h @@ -436,10 +436,10 @@ class DecaySpectrum : public Distribution { //! //! \param nuclide_indices Indices of decay photon emitters in //! data::chain_nuclides - //! \param weights Activity-based weights for each component. The Mixture - //! probability for component i is weights[i] * + //! \param weights Atom counts for each component. The probability for + //! component i is weights[i] * //! chain_nuclides[nuclide_indices[i]]->photon_energy()->integral() (the - //! integral encodes photons-per-decay). + //! integral encodes photons-per-second-per-atom, i.e., Bq/atom). DecaySpectrum(vector nuclide_indices, vector weights); //! Construct from an XML node containing nuclide names and atom densities. diff --git a/src/distribution.cpp b/src/distribution.cpp index 301896286cb..bd7ebafaf7e 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -830,10 +830,12 @@ DecaySpectrum::DecaySpectrum(pugi::xml_node node) if (!photon_dist) continue; - // Weight = atom_count * decay_constant - // atom_count = atom_density [atom/b-cm] * 1e24 [b-cm/cm^3] * volume [cm^3] + // Weight = atom_count: photon_energy()->integral() is in [photons/s/atom] + // (Bq/atom), so multiplying by atom count gives total emission rate + // [photons/s]. atom_count = atom_density [atom/b-cm] * 1e24 [b-cm/cm^3] * + // volume [cm^3] double atom_count = density * 1.0e24 * volume; - double weight = atom_count * chain_nuc->decay_constant(); + double weight = atom_count; if (weight <= 0.0) continue; From c4e650904f31006156e7ede90cd7097701d51a65 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 23 Apr 2026 20:49:56 -0500 Subject: [PATCH 26/37] Improve variable names --- include/openmc/distribution.h | 7 ++----- src/distribution.cpp | 35 +++++++++++++++-------------------- 2 files changed, 17 insertions(+), 25 deletions(-) diff --git a/include/openmc/distribution.h b/include/openmc/distribution.h index cf01ebcadf9..bfb3ca96378 100644 --- a/include/openmc/distribution.h +++ b/include/openmc/distribution.h @@ -436,11 +436,8 @@ class DecaySpectrum : public Distribution { //! //! \param nuclide_indices Indices of decay photon emitters in //! data::chain_nuclides - //! \param weights Atom counts for each component. The probability for - //! component i is weights[i] * - //! chain_nuclides[nuclide_indices[i]]->photon_energy()->integral() (the - //! integral encodes photons-per-second-per-atom, i.e., Bq/atom). - DecaySpectrum(vector nuclide_indices, vector weights); + //! \param atoms Number of atoms for each component. + DecaySpectrum(vector nuclide_indices, vector atoms); //! Construct from an XML node containing nuclide names and atom densities. //! diff --git a/src/distribution.cpp b/src/distribution.cpp index bd7ebafaf7e..0154df1b499 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -777,29 +777,29 @@ UPtrDist distribution_from_xml(pugi::xml_node node) // Helper function to compute weighted probabilities for decay spectrum sampling vector decay_spectrum_probabilities( - const vector& nuclide_indices, const vector& weights) + const vector& nuclide_indices, const vector& atoms) { - if (nuclide_indices.size() != weights.size()) { - fatal_error("DecaySpectrum nuclide index and weight arrays must have the " - "same length."); + if (nuclide_indices.size() != atoms.size()) { + fatal_error("DecaySpectrum nuclide index and atoms arrays must have " + "the same length."); } vector probs; probs.reserve(nuclide_indices.size()); for (size_t i = 0; i < nuclide_indices.size(); ++i) { + // Distribution integral is in [photons/s/atom]; multiplying by atoms gives + // the total emission rate [photons/s] for this nuclide. const auto* dist = data::chain_nuclides[nuclide_indices[i]]->photon_energy(); - probs.push_back(weights[i] * dist->integral()); + probs.push_back(atoms[i] * dist->integral()); } return probs; } -DecaySpectrum::DecaySpectrum( - vector nuclide_indices, vector weights) +DecaySpectrum::DecaySpectrum(vector nuclide_indices, vector atoms) : nuclide_indices_(std::move(nuclide_indices)) { - vector probs = - decay_spectrum_probabilities(nuclide_indices_, weights); + vector probs = decay_spectrum_probabilities(nuclide_indices_, atoms); integral_ = std::accumulate(probs.begin(), probs.end(), 0.0); di_.assign(probs); } @@ -813,7 +813,7 @@ DecaySpectrum::DecaySpectrum(pugi::xml_node node) // Read nuclide names and atom densities from XML vector nuclide_indices; - vector weights; + vector atoms; for (auto nuclide_node : node.children("nuclide")) { std::string name = get_node_value(nuclide_node, "name"); @@ -830,24 +830,19 @@ DecaySpectrum::DecaySpectrum(pugi::xml_node node) if (!photon_dist) continue; - // Weight = atom_count: photon_energy()->integral() is in [photons/s/atom] - // (Bq/atom), so multiplying by atom count gives total emission rate - // [photons/s]. atom_count = atom_density [atom/b-cm] * 1e24 [b-cm/cm^3] * - // volume [cm^3] - double atom_count = density * 1.0e24 * volume; - double weight = atom_count; - if (weight <= 0.0) + // atoms = density [atom/b-cm] * 1e24 [b/cm^2] * volume [cm^3] + double atoms_i = density * 1.0e24 * volume; + if (atoms_i <= 0.0) continue; nuclide_indices.push_back(nuclide_index); - weights.push_back(weight); + atoms.push_back(atoms_i); } nuclide_indices_ = std::move(nuclide_indices); // Build alias table from weighted probabilities - vector probs = - decay_spectrum_probabilities(nuclide_indices_, weights); + vector probs = decay_spectrum_probabilities(nuclide_indices_, atoms); integral_ = std::accumulate(probs.begin(), probs.end(), 0.0); di_.assign(probs); } From 5772de2e2056adfd7650e51c1fcc8e2e97573af1 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 24 Apr 2026 09:23:31 -0500 Subject: [PATCH 27/37] Avoid IDWarning when creating photon sources Co-authored-by: Copilot --- openmc/deplete/r2s.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index d9c539f5419..4ab007d5e41 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -620,13 +620,15 @@ def _create_photon_sources(self, time_index, work_items): step_result = self.results['depletion_results'][time_index] materials = self.results['activation_materials'] mesh_based = self.method == 'mesh-based' + if mesh_based: + mat_dict = self.neutron_model._get_all_materials() sources = [] for item in work_items: if mesh_based: index_mat, domain_id, bbox = item original_mat = materials[index_mat] - domain = openmc.Material(material_id=domain_id) + domain = mat_dict[domain_id] else: cell, original_mat, bbox = item domain = cell From c5060a5215de4603d918efa806c4987c365b685e Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 28 Apr 2026 14:58:20 -0500 Subject: [PATCH 28/37] Fix import, missing header Co-authored-by: Copilot --- openmc/stats/univariate.py | 7 ++++--- src/finalize.cpp | 1 + 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py index b0203fbbc49..834eacd91f6 100644 --- a/openmc/stats/univariate.py +++ b/openmc/stats/univariate.py @@ -14,7 +14,8 @@ import scipy import openmc.checkvalue as cv -from openmc.data import atomic_mass, NEUTRON_MASS, decay_photon_energy +from openmc.data import atomic_mass, NEUTRON_MASS +import openmc.data from .._xml import get_elem_list, get_text from ..mixin import EqualityMixin @@ -2288,7 +2289,7 @@ def to_distribution(self): dists = [] weights = [] for name, density in self.nuclides.items(): - dist = decay_photon_energy(name) + dist = openmc.data.decay_photon_energy(name) if dist is not None: dists.append(dist) weights.append(density * 1e24 * self.volume) @@ -2380,7 +2381,7 @@ def integral(self): @cache def _photon_integral(nuclide: str) -> float | None: """Return the per-atom photon emission integral for a nuclide""" - dist = decay_photon_energy(nuclide) + dist = openmc.data.decay_photon_energy(nuclide) return dist.integral() if dist is not None else None def clip(self, tolerance: float = 1e-9, inplace: bool = False): diff --git a/src/finalize.cpp b/src/finalize.cpp index 1c4b3df6811..59711dcc31d 100644 --- a/src/finalize.cpp +++ b/src/finalize.cpp @@ -2,6 +2,7 @@ #include "openmc/bank.h" #include "openmc/capi.h" +#include "openmc/chain.h" #include "openmc/cmfd_solver.h" #include "openmc/collision_track.h" #include "openmc/constants.h" From a9aa40fd177b9e189a0a7c6047dcfa230fcbcda0 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 28 Apr 2026 15:07:51 -0500 Subject: [PATCH 29/37] Avoid empty DecaySpectrum distribution --- src/distribution.cpp | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/distribution.cpp b/src/distribution.cpp index 0154df1b499..30589a4a1e9 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -801,6 +801,10 @@ DecaySpectrum::DecaySpectrum(vector nuclide_indices, vector atoms) { vector probs = decay_spectrum_probabilities(nuclide_indices_, atoms); integral_ = std::accumulate(probs.begin(), probs.end(), 0.0); + if (nuclide_indices_.empty() || integral_ <= 0.0) { + fatal_error("DecaySpectrum must contain at least one decay photon emitter " + "with a positive emission rate."); + } di_.assign(probs); } @@ -844,6 +848,12 @@ DecaySpectrum::DecaySpectrum(pugi::xml_node node) // Build alias table from weighted probabilities vector probs = decay_spectrum_probabilities(nuclide_indices_, atoms); integral_ = std::accumulate(probs.begin(), probs.end(), 0.0); + if (nuclide_indices_.empty() || integral_ <= 0.0) { + fatal_error("DecaySpectrum source did not resolve any nuclides with decay " + "photon spectra and positive atom densities. Ensure " + "OPENMC_CHAIN_FILE is set and matches the nuclides in the " + "source definition."); + } di_.assign(probs); } From 9d869e6cb3ade8cc243f73afc12ec597b1e1f23a Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 28 Apr 2026 17:49:19 -0500 Subject: [PATCH 30/37] Make the distribution cache in DecaySpectrum chain-specific --- openmc/stats/univariate.py | 28 +++++++++++++++++++++++++--- 1 file changed, 25 insertions(+), 3 deletions(-) diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py index 834eacd91f6..2a7d62bbd5a 100644 --- a/openmc/stats/univariate.py +++ b/openmc/stats/univariate.py @@ -5,6 +5,7 @@ from functools import cache from math import sqrt, pi, exp, log from numbers import Real +from pathlib import Path from warnings import warn import lxml.etree as ET @@ -2234,6 +2235,7 @@ class DecaySpectrum(Univariate): def __init__(self, nuclides: dict[str, float], volume: float): super().__init__(bias=None) self._dist_cache = None + self._dist_cache_key = None self.nuclides = nuclides self.volume = volume @@ -2253,6 +2255,7 @@ def nuclides(self, nuclides): cv.check_greater_than(f'atom density for {name}', density, 0.0, True) self._nuclides = dict(nuclides) self._dist_cache = None + self._dist_cache_key = None @property def volume(self): @@ -2264,6 +2267,21 @@ def volume(self, volume): cv.check_greater_than('volume', volume, 0.0) self._volume = float(volume) self._dist_cache = None + self._dist_cache_key = None + + @staticmethod + def _chain_file_cache_key(): + """Return a hashable key for the active depletion chain.""" + chain_file = openmc.config.get('chain_file') + if chain_file is None: + return None + + path = Path(chain_file).resolve() + try: + stat = path.stat() + except OSError: + return (path, None, None) + return (path, stat.st_mtime, stat.st_size) def to_distribution(self): """Convert to a concrete distribution using decay chain data. @@ -2283,7 +2301,8 @@ def to_distribution(self): :attr:`nuclides` has a photon source in the chain. """ - if self._dist_cache is not None: + chain_key = self._chain_file_cache_key() + if self._dist_cache is not None and self._dist_cache_key == chain_key: return self._dist_cache dists = [] @@ -2298,6 +2317,7 @@ def to_distribution(self): return None self._dist_cache = combine_distributions(dists, weights) + self._dist_cache_key = chain_key return self._dist_cache def to_xml_element(self, element_name: str): @@ -2379,7 +2399,7 @@ def integral(self): @staticmethod @cache - def _photon_integral(nuclide: str) -> float | None: + def _photon_integral(nuclide: str, chain_key) -> float | None: """Return the per-atom photon emission integral for a nuclide""" dist = openmc.data.decay_photon_energy(nuclide) return dist.integral() if dist is not None else None @@ -2414,8 +2434,9 @@ def clip(self, tolerance: float = 1e-9, inplace: bool = False): emitting_names = [] emitting_densities = [] rates = [] + chain_key = self._chain_file_cache_key() for name, density in self.nuclides.items(): - integral = DecaySpectrum._photon_integral(name) + integral = DecaySpectrum._photon_integral(name, chain_key) if integral is None: continue emitting_names.append(name) @@ -2433,6 +2454,7 @@ def clip(self, tolerance: float = 1e-9, inplace: bool = False): if inplace: self._nuclides = new_nuclides self._dist_cache = None + self._dist_cache_key = None return self return type(self)(new_nuclides, self.volume) From 4d085c26af032134cd3ef7063643139513879cd0 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 28 Apr 2026 23:05:22 -0500 Subject: [PATCH 31/37] Use chain_file consistently in R2SManager.run Co-authored-by: Copilot --- openmc/deplete/r2s.py | 56 ++++++++++++++++++++++++++++--------------- 1 file changed, 37 insertions(+), 19 deletions(-) diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index 4ab007d5e41..6dbb3adb2c0 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -1,5 +1,6 @@ from __future__ import annotations from collections.abc import Sequence +from contextlib import nullcontext import copy from datetime import datetime import json @@ -8,6 +9,7 @@ import numpy as np import openmc from . import IndependentOperator, PredictorIntegrator +from .chain import Chain from .microxs import get_microxs_and_flux, write_microxs_hdf5, read_microxs_hdf5 from .results import Results from ..checkvalue import PathLike @@ -149,7 +151,7 @@ def run( photon_time_indices: Sequence[int] | None = None, output_dir: PathLike | None = None, bounding_boxes: dict[int, openmc.BoundingBox] | None = None, - chain_file: PathLike | None = None, + chain_file: PathLike | Chain | None = None, micro_kwargs: dict | None = None, mat_vol_kwargs: dict | None = None, run_kwargs: dict | None = None, @@ -187,9 +189,10 @@ def run( Dictionary mapping cell IDs to bounding boxes used for spatial source sampling in cell-based R2S calculations. Required if method is 'cell-based'. - chain_file : PathLike, optional - Path to the depletion chain XML file to use during activation. If - not provided, the default configured chain file will be used. + chain_file : PathLike or openmc.deplete.Chain, optional + Path to the depletion chain XML file or depletion chain object to + use during activation. If not provided, the default configured + chain file will be used. micro_kwargs : dict, optional Additional keyword arguments passed to :func:`openmc.deplete.get_microxs_and_flux` during the neutron @@ -216,6 +219,8 @@ def run( # consistency (different ranks may have slightly different times) stamp = datetime.now().strftime('%Y-%m-%dT%H-%M-%S') output_dir = Path(comm.bcast(f'r2s_{stamp}')) + else: + output_dir = Path(output_dir) # Set run_kwargs for the neutron transport step if micro_kwargs is None: @@ -226,22 +231,35 @@ def run( operator_kwargs = {} run_kwargs.setdefault('output', False) micro_kwargs.setdefault('run_kwargs', run_kwargs) - # If a chain file is provided, prefer it for steps 1 and 2 - if chain_file is not None: - micro_kwargs.setdefault('chain_file', chain_file) - operator_kwargs.setdefault('chain_file', chain_file) - self.step1_neutron_transport( - output_dir / 'neutron_transport', mat_vol_kwargs, micro_kwargs - ) - self.step2_activation( - timesteps, source_rates, timestep_units, output_dir / 'activation', - operator_kwargs=operator_kwargs - ) - self.step3_photon_transport( - photon_time_indices, bounding_boxes, output_dir / 'photon_transport', - mat_vol_kwargs=mat_vol_kwargs, run_kwargs=run_kwargs + # DecaySpectrum distributions are resolved in the C++ solver using + # OPENMC_CHAIN_FILE. If a Chain object was passed, write an XML + # representation alongside the R2S outputs. + if isinstance(chain_file, Chain): + output_dir.mkdir(parents=True, exist_ok=True) + chain_path = output_dir / 'chain.xml' + if comm.rank == 0: + chain_file.export_to_xml(chain_path) + comm.barrier() + else: + chain_path = chain_file + + chain_context = ( + openmc.config.patch('chain_file', chain_path) + if chain_path is not None else nullcontext() ) + with chain_context: + self.step1_neutron_transport( + output_dir / 'neutron_transport', mat_vol_kwargs, micro_kwargs + ) + self.step2_activation( + timesteps, source_rates, timestep_units, + output_dir / 'activation', operator_kwargs=operator_kwargs + ) + self.step3_photon_transport( + photon_time_indices, bounding_boxes, output_dir / 'photon_transport', + mat_vol_kwargs=mat_vol_kwargs, run_kwargs=run_kwargs + ) return output_dir @@ -535,7 +553,7 @@ def step3_photon_transport( for time_index in time_indices: # Convert time_index (which may be negative) to a normal index if time_index < 0: - time_index = len(self.results['depletion_results']) + time_index + time_index += len(self.results['depletion_results']) # Build decay photon sources and assign to the photon model sources = self._create_photon_sources(time_index, work_items) From 65a502f9a69e8e3608b0ce106d0ddd51795aabdb Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 1 May 2026 17:38:59 -0500 Subject: [PATCH 32/37] Add warning about decay spectrum nuclide missing in chain --- src/distribution.cpp | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/src/distribution.cpp b/src/distribution.cpp index 30589a4a1e9..f937ade9b94 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -7,6 +7,7 @@ #include // for accumulate #include // for runtime_error #include // for string, stod +#include #include "openmc/chain.h" #include "openmc/constants.h" @@ -16,6 +17,10 @@ #include "openmc/random_lcg.h" #include "openmc/xml_interface.h" +namespace { +std::unordered_set decay_spectrum_missing_chain_nuclides; +} + namespace openmc { //============================================================================== @@ -825,8 +830,14 @@ DecaySpectrum::DecaySpectrum(pugi::xml_node node) // Look up nuclide in the depletion chain auto it = data::chain_nuclide_map.find(name); - if (it == data::chain_nuclide_map.end()) + if (it == data::chain_nuclide_map.end()) { + if (decay_spectrum_missing_chain_nuclides.insert(name).second) { + warning("Nuclide '" + name + + "' appears in a DecaySpectrum source but is not present in " + "the depletion chain; it will be ignored."); + } continue; + } int nuclide_index = it->second; const auto& chain_nuc = data::chain_nuclides[nuclide_index]; From 7031432defa0bb3a58505018aefa1c0df2be0d5b Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 1 May 2026 17:51:25 -0500 Subject: [PATCH 33/37] Simplify initialization --- include/openmc/distribution.h | 13 ++++---- src/distribution.cpp | 57 +++++++++++++---------------------- 2 files changed, 27 insertions(+), 43 deletions(-) diff --git a/include/openmc/distribution.h b/include/openmc/distribution.h index bfb3ca96378..f319dd19df4 100644 --- a/include/openmc/distribution.h +++ b/include/openmc/distribution.h @@ -432,13 +432,6 @@ class DecaySpectrum : public Distribution { //============================================================================ // Constructors - //! Construct from depletion chain nuclide indices and weights - //! - //! \param nuclide_indices Indices of decay photon emitters in - //! data::chain_nuclides - //! \param atoms Number of atoms for each component. - DecaySpectrum(vector nuclide_indices, vector atoms); - //! Construct from an XML node containing nuclide names and atom densities. //! //! Reads child ```` elements with ``name`` and ``density`` @@ -468,6 +461,12 @@ class DecaySpectrum : public Distribution { double sample_unbiased(uint64_t* seed) const override; private: + //! Initialize decay spectrum sampling data + //! \param nuclide_indices Indices of decay photon emitters in + //! data::chain_nuclides + //! \param atoms Number of atoms for each component. + void init(vector nuclide_indices, const vector& atoms); + vector nuclide_indices_; //!< Indices of emitting nuclides in the chain DiscreteIndex di_; //!< Discrete index for component selection double integral_; //!< Total photon emission rate diff --git a/src/distribution.cpp b/src/distribution.cpp index f937ade9b94..aa96c694015 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -780,39 +780,6 @@ UPtrDist distribution_from_xml(pugi::xml_node node) // DecaySpectrum implementation //============================================================================== -// Helper function to compute weighted probabilities for decay spectrum sampling -vector decay_spectrum_probabilities( - const vector& nuclide_indices, const vector& atoms) -{ - if (nuclide_indices.size() != atoms.size()) { - fatal_error("DecaySpectrum nuclide index and atoms arrays must have " - "the same length."); - } - - vector probs; - probs.reserve(nuclide_indices.size()); - for (size_t i = 0; i < nuclide_indices.size(); ++i) { - // Distribution integral is in [photons/s/atom]; multiplying by atoms gives - // the total emission rate [photons/s] for this nuclide. - const auto* dist = - data::chain_nuclides[nuclide_indices[i]]->photon_energy(); - probs.push_back(atoms[i] * dist->integral()); - } - return probs; -} - -DecaySpectrum::DecaySpectrum(vector nuclide_indices, vector atoms) - : nuclide_indices_(std::move(nuclide_indices)) -{ - vector probs = decay_spectrum_probabilities(nuclide_indices_, atoms); - integral_ = std::accumulate(probs.begin(), probs.end(), 0.0); - if (nuclide_indices_.empty() || integral_ <= 0.0) { - fatal_error("DecaySpectrum must contain at least one decay photon emitter " - "with a positive emission rate."); - } - di_.assign(probs); -} - DecaySpectrum::DecaySpectrum(pugi::xml_node node) { // Read the region volume [cm^3] needed for absolute emission rate @@ -854,10 +821,28 @@ DecaySpectrum::DecaySpectrum(pugi::xml_node node) atoms.push_back(atoms_i); } - nuclide_indices_ = std::move(nuclide_indices); + init(std::move(nuclide_indices), atoms); +} - // Build alias table from weighted probabilities - vector probs = decay_spectrum_probabilities(nuclide_indices_, atoms); +void DecaySpectrum::init( + vector nuclide_indices, const vector& atoms) +{ + if (nuclide_indices.size() != atoms.size()) { + fatal_error("DecaySpectrum nuclide index and atoms arrays must have " + "the same length."); + } + + vector probs; + probs.reserve(nuclide_indices.size()); + for (size_t i = 0; i < nuclide_indices.size(); ++i) { + // Distribution integral is in [photons/s/atom]; multiplying by atoms gives + // the total emission rate [photons/s] for this nuclide. + const auto* dist = + data::chain_nuclides[nuclide_indices[i]]->photon_energy(); + probs.push_back(atoms[i] * dist->integral()); + } + + nuclide_indices_ = std::move(nuclide_indices); integral_ = std::accumulate(probs.begin(), probs.end(), 0.0); if (nuclide_indices_.empty() || integral_ <= 0.0) { fatal_error("DecaySpectrum source did not resolve any nuclides with decay " From 8dc10bbe51bcb21af33719a3932cffcf325ddb60 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 1 May 2026 17:57:26 -0500 Subject: [PATCH 34/37] Show warning for negative densities --- src/distribution.cpp | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/src/distribution.cpp b/src/distribution.cpp index aa96c694015..2a80ad09b19 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -812,10 +812,18 @@ DecaySpectrum::DecaySpectrum(pugi::xml_node node) if (!photon_dist) continue; + // Skip non-positive densities and warn if negative + if (density <= 0.0) { + if (density < 0.0) { + warning("Nuclide '" + name + + "' has a negative density in a DecaySpectrum source; it will " + "be ignored."); + } + continue; + } + // atoms = density [atom/b-cm] * 1e24 [b/cm^2] * volume [cm^3] double atoms_i = density * 1.0e24 * volume; - if (atoms_i <= 0.0) - continue; nuclide_indices.push_back(nuclide_index); atoms.push_back(atoms_i); From f375e3c1b76ca81666df41d6e31afbf662840af3 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 1 May 2026 18:09:40 -0500 Subject: [PATCH 35/37] Warn about non-unity source strength getting overridden --- src/source.cpp | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/source.cpp b/src/source.cpp index 50a6640c2ca..38492056f92 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -356,6 +356,13 @@ IndependentSource::IndependentSource(pugi::xml_node node) : Source(node) // For decay photon sources, use the absolute photon emission rate in // [photons/s] as the source strength if (dynamic_cast(energy_.get())) { + if (strength_ != 1.0) { + warning(fmt::format( + "Source strength of {} is ignored because the source uses a " + "DecaySpectrum energy distribution. The source strength will be " + "set from the DecaySpectrum emission rate.", + strength_)); + } strength_ = energy_->integral(); } } else { From 163224d9afd379428a8f1d758805e86e2d9bc7cd Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 1 May 2026 18:18:42 -0500 Subject: [PATCH 36/37] Update documentation --- docs/source/io_formats/settings.rst | 34 +++++++++++++++++++++-------- 1 file changed, 25 insertions(+), 9 deletions(-) diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index a50922b041e..1bd0dd25f8d 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -1029,17 +1029,19 @@ variable and whose sub-elements/attributes are as follows: :type: The type of the distribution. Valid options are "uniform", "discrete", - "tabular", "maxwell", "watt", and "mixture". The "uniform" option produces - variates sampled from a uniform distribution over a finite interval. The - "discrete" option produces random variates that can assume a finite number - of values (i.e., a distribution characterized by a probability mass function). - The "tabular" option produces random variates sampled from a tabulated - distribution where the density function is either a histogram or + "tabular", "maxwell", "watt", "mixture", and "decay_spectrum". The "uniform" + option produces variates sampled from a uniform distribution over a finite + interval. The "discrete" option produces random variates that can assume a + finite number of values (i.e., a distribution characterized by a probability + mass function). The "tabular" option produces random variates sampled from a + tabulated distribution where the density function is either a histogram or linearly-interpolated between tabulated points. The "watt" option produces random variates is sampled from a Watt fission spectrum (only used for energies). The "maxwell" option produce variates sampled from a Maxwell - fission spectrum (only used for energies). The "mixture" option produces samples - from univariate sub-distributions with given probabilities. + fission spectrum (only used for energies). The "mixture" option produces + samples from univariate sub-distributions with given probabilities. The + "decay_spectrum" option produces photon energies sampled from decay photon + spectra in a depletion chain (only used for energies). *Default*: None @@ -1086,6 +1088,20 @@ variable and whose sub-elements/attributes are as follows: This sub-element of a ``pair`` element provides information on the corresponding univariate distribution. +:volume: + For a "decay_spectrum" distribution, this attribute specifies the source + region volume in cm\ :sup:`3`. It is used together with atom densities to + determine the absolute photon emission rate. When a source uses a + "decay_spectrum" energy distribution, the source strength is set from this + emission rate. + +:nuclide: + For a "decay_spectrum" distribution, each ``nuclide`` sub-element specifies a + nuclide contributing to the decay photon source. The nuclide name is given by + the ``name`` attribute, and the atom density in [atom/b-cm] is given by the + ``density`` attribute. Nuclides are resolved against the depletion chain, and + nuclides without decay photon spectra do not contribute to the distribution. + :bias: This optional element specifies a biased distribution for importance sampling. For continuous distributions, the ``bias`` element should contain another @@ -1271,7 +1287,7 @@ The ```` element specifies the surface flux cosine cutof ```` Element ----------------------------------- -The ```` element specifies the surface flux cosine +The ```` element specifies the surface flux cosine substitution ratio. *Default*: 0.5 From a33f4907ea7136009d56add7f235b2589760a197 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 2 May 2026 09:55:12 -0500 Subject: [PATCH 37/37] Make XML representation for DecaySpectrum more compact --- docs/source/io_formats/settings.rst | 17 +++++++++++------ openmc/stats/univariate.py | 16 +++++++--------- src/distribution.cpp | 12 +++++++++--- tests/unit_tests/test_stats.py | 3 +++ 4 files changed, 30 insertions(+), 18 deletions(-) diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 1bd0dd25f8d..1aeffc23441 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -1059,6 +1059,10 @@ variable and whose sub-elements/attributes are as follows: :math:`(x,p)` pairs defining the discrete/tabular distribution. All :math:`x` points are given first followed by corresponding :math:`p` points. + For a "decay_spectrum" distribution, ``parameters`` gives the atom densities + in [atom/b-cm] for the nuclides listed in the ``nuclides`` element, in the + same order. + For a "watt" distribution, ``parameters`` should be given as two real numbers :math:`a` and :math:`b` that parameterize the distribution :math:`p(x) dx = c e^{-x/a} \sinh \sqrt{b \, x} dx`. @@ -1095,12 +1099,13 @@ variable and whose sub-elements/attributes are as follows: "decay_spectrum" energy distribution, the source strength is set from this emission rate. -:nuclide: - For a "decay_spectrum" distribution, each ``nuclide`` sub-element specifies a - nuclide contributing to the decay photon source. The nuclide name is given by - the ``name`` attribute, and the atom density in [atom/b-cm] is given by the - ``density`` attribute. Nuclides are resolved against the depletion chain, and - nuclides without decay photon spectra do not contribute to the distribution. +:nuclides: + For a "decay_spectrum" distribution, this element specifies a + whitespace-separated list of nuclide names contributing to the decay photon + source. The atom densities for these nuclides are given by the ``parameters`` + element in the same order. Nuclides are resolved against the depletion chain, + and nuclides without decay photon spectra do not contribute to the + distribution. :bias: This optional element specifies a biased distribution for importance sampling. diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py index 2a7d62bbd5a..1b18bcb1745 100644 --- a/openmc/stats/univariate.py +++ b/openmc/stats/univariate.py @@ -2337,10 +2337,10 @@ def to_xml_element(self, element_name: str): element = ET.Element(element_name) element.set("type", "decay_spectrum") element.set("volume", str(self.volume)) - for name, density in self.nuclides.items(): - nuclide_elem = ET.SubElement(element, "nuclide") - nuclide_elem.set("name", name) - nuclide_elem.set("density", str(density)) + nuclides = ET.SubElement(element, "nuclides") + nuclides.text = ' '.join(self.nuclides) + parameters = ET.SubElement(element, "parameters") + parameters.text = ' '.join(str(density) for density in self.nuclides.values()) return element @classmethod @@ -2359,11 +2359,9 @@ def from_xml_element(cls, elem: ET.Element): """ volume = float(elem.get('volume')) - nuclides = {} - for nuclide_elem in elem.findall('nuclide'): - name = nuclide_elem.get('name') - density = float(nuclide_elem.get('density')) - nuclides[name] = density + names = get_elem_list(elem, 'nuclides', str) + densities = get_elem_list(elem, 'parameters', float) + nuclides = dict(zip(names, densities)) return cls(nuclides, volume) def _sample_unbiased(self, n_samples=1, seed=None): diff --git a/src/distribution.cpp b/src/distribution.cpp index 2a80ad09b19..7f5b498add9 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -790,10 +790,16 @@ DecaySpectrum::DecaySpectrum(pugi::xml_node node) // Read nuclide names and atom densities from XML vector nuclide_indices; vector atoms; + auto names = get_node_array(node, "nuclides"); + auto densities = get_node_array(node, "parameters"); + if (names.size() != densities.size()) { + fatal_error("DecaySpectrum nuclides and parameters must have the same " + "length."); + } - for (auto nuclide_node : node.children("nuclide")) { - std::string name = get_node_value(nuclide_node, "name"); - double density = std::stod(get_node_value(nuclide_node, "density")); + for (size_t i = 0; i < names.size(); ++i) { + const auto& name = names[i]; + double density = densities[i]; // Look up nuclide in the depletion chain auto it = data::chain_nuclide_map.find(name); diff --git a/tests/unit_tests/test_stats.py b/tests/unit_tests/test_stats.py index 88195b159d5..2754cf5d09a 100644 --- a/tests/unit_tests/test_stats.py +++ b/tests/unit_tests/test_stats.py @@ -1059,6 +1059,9 @@ def test_decay_spectrum_xml_roundtrip(): elem = d.to_xml_element('energy') assert elem.get('type') == 'decay_spectrum' assert float(elem.get('volume')) == pytest.approx(100.0) + assert elem.findtext('nuclides').split() == list(nuclides) + assert [float(x) for x in elem.findtext('parameters').split()] == pytest.approx( + list(nuclides.values())) # Round-trip via DecaySpectrum.from_xml_element d2 = openmc.stats.DecaySpectrum.from_xml_element(elem)