From e542b2f035e397a647cc3ad0ae39946b1ad81e7b Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Wed, 6 May 2026 16:07:35 +0100 Subject: [PATCH 01/10] Allow Mesh.volumes property for 1D and 2D RegularMesh (#3914) --- openmc/mesh.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 3c3c0a1ac..4129913d8 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -268,10 +268,11 @@ class MeshBase(IDManagerMixin, ABC): return string def _volume_dim_check(self): - if self.n_dimension != 3 or \ - any([d == 0 for d in self.dimension]): - raise RuntimeError(f'Mesh {self.id} is not 3D. ' - 'Volumes cannot be provided.') + if any(d == 0 for d in self.dimension): + raise RuntimeError( + f'Mesh {self.id} has a zero-size dimension. ' + 'Volumes cannot be provided.' + ) @classmethod def from_hdf5(cls, group: h5py.Group): From f3e1066d467a3c12a9c964b624608910bc2121ad Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 7 May 2026 18:50:23 -0500 Subject: [PATCH 02/10] Include mass_attenuation.h5 in package data (#3936) --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index bf4011344..098487c06 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -69,7 +69,7 @@ include = ['openmc*'] exclude = ['tests*'] [tool.setuptools.package-data] -"openmc.data.dose" = ["**/*.txt"] +"openmc.data.dose" = ["**/*.txt", "*.h5"] "openmc.data" = ["*.txt", "*.DAT", "*.json", "*.h5"] "openmc.lib" = ["libopenmc.dylib", "libopenmc.so"] From d56cda254408971c3f0a79736c46b080ca37f1b7 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 8 May 2026 20:53:12 -0500 Subject: [PATCH 03/10] Implement DecaySpectrum distribution type and utilize in R2S (#3930) Co-authored-by: Copilot --- docs/source/io_formats/settings.rst | 39 +++- docs/source/pythonapi/stats.rst | 1 + include/openmc/chain.h | 2 + include/openmc/distribution.h | 65 ++++++ openmc/deplete/r2s.py | 244 ++++++++++++---------- openmc/stats/univariate.py | 309 ++++++++++++++++++++++++++++ src/chain.cpp | 8 + src/distribution.cpp | 124 +++++++++++ src/finalize.cpp | 2 + src/initialize.cpp | 14 +- src/source.cpp | 30 ++- tests/unit_tests/test_source.py | 69 +++++++ tests/unit_tests/test_stats.py | 135 +++++++++++- 13 files changed, 908 insertions(+), 134 deletions(-) diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index d63376e2b..f8b064a37 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -1058,17 +1058,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,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`. @@ -1115,6 +1121,21 @@ 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. + +: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. For continuous distributions, the ``bias`` element should contain another @@ -1300,7 +1321,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 diff --git a/docs/source/pythonapi/stats.rst b/docs/source/pythonapi/stats.rst index 203d4fea4..2f2e2835a 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/chain.h b/include/openmc/chain.h index 6f5830358..e01f2a01c 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 cd81d8801..f319dd19d 100644 --- a/include/openmc/distribution.h +++ b/include/openmc/distribution.h @@ -407,6 +407,71 @@ private: 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: + //============================================================================ + // Types, aliases + + struct Sample { + double energy; + double weight; + int parent_nuclide; + }; + + //============================================================================ + // Constructors + + //! 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); + + //============================================================================ + // 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: + //! 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 +}; + } // namespace openmc #endif // OPENMC_DISTRIBUTION_H diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index 10d7fdd50..6dbb3adb2 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 @@ from pathlib import Path 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 @@ class R2SManager: 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 @@ class R2SManager: 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 @@ class R2SManager: # 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 @@ class R2SManager: 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 @@ -516,45 +534,30 @@ class R2SManager: if different_photon_model: photon_cells = self.photon_model.geometry.get_all_cells() - 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 \ + # Determine eligible work items upfront (independent of time index). + if self.method == 'mesh-based': + work_items = self._get_mesh_work_items() + else: + 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 + continue + work_items.append((cell, original_mat, bounding_boxes[cell.id])) - # 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 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]} - ) - sources.append(source) - self.photon_model.settings.source = 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 + 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) + self.photon_model.settings.source = sources # Run photon transport calculation photon_dir = Path(output_dir) / f'time_{time_index}' @@ -567,58 +570,30 @@ class R2SManager: 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. + def _get_mesh_work_items(self): + """Enumerate mesh-based work items across all meshes. - 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 a list of (index_mat, mat_id, bbox) tuples for each eligible + 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 ------- - 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. - + list of tuple + Each tuple is (index_mat, mat_id, bbox). """ - 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') + 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 - # 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 @@ -626,36 +601,77 @@ class R2SManager: if mat_id is not None } - for mat_id, _, bbox in mat_vols.by_element(index_elem, include_bboxes=True): - # Skip void volume + for mat_id, _, bbox in mat_vols.by_element( + index_elem, include_bboxes=True): 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: + if photon_mat_vols is not None \ + and mat_id not in photon_materials: 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 work_items + + 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.DecaySpectrum` energy distributions that will be + serialized to XML and resolved against the depletion chain by the C++ + solver. + + 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). + + 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' + 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 = mat_dict[domain_id] + else: + cell, original_mat, bbox = item + domain = cell + + activated_mat = step_result.get_material(str(original_mat.id)) + nuclides = activated_mat.get_nuclide_atom_densities() + 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: + continue + + sources.append(openmc.IndependentSource( + space=openmc.stats.Box(bbox.lower_left, bbox.upper_right), + energy=energy, + particle='photon', + constraints={'domains': [domain]}, + )) + return sources def load_results(self, path: PathLike): diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py index 1cf9a1ad5..1b18bcb17 100644 --- a/openmc/stats/univariate.py +++ b/openmc/stats/univariate.py @@ -2,8 +2,10 @@ from __future__ import annotations 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 pathlib import Path from warnings import warn import lxml.etree as ET @@ -14,6 +16,7 @@ import scipy import openmc.checkvalue as cv from openmc.data import atomic_mass, NEUTRON_MASS +import openmc.data from .._xml import get_elem_list, get_text from ..mixin import EqualityMixin @@ -124,6 +127,8 @@ class Univariate(EqualityMixin, ABC): return Legendre.from_xml_element(elem) elif distribution == 'mixture': return Mixture.from_xml_element(elem) + elif distribution == 'decay_spectrum': + return DecaySpectrum.from_xml_element(elem) @abstractmethod def _sample_unbiased(self, n_samples: int = 1, seed: int | None = None): @@ -2196,6 +2201,310 @@ class Mixture(Univariate): return new_dist +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. 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]. + volume : float + 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]. + volume : float + Volume of the source region in [cm³]. + + """ + + 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 + + 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, True) + self._nuclides = dict(nuclides) + self._dist_cache = None + self._dist_cache_key = None + + @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) + 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. + + 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. + + """ + 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 = [] + 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) + self._dist_cache_key = chain_key + return self._dist_cache + + 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_spectrum") + element.set("volume", str(self.volume)) + 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 + 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.DecaySpectrum + Decay photon distribution generated from XML element + + """ + volume = float(elem.get('volume')) + 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): + 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 + Total photon emission rate in [photons/s], or ``0.0`` if chain + data is unavailable. + """ + try: + dist = self.to_distribution() + except Exception: + return 0.0 + if dist is None: + return 0.0 + return dist.integral() + + @staticmethod + @cache + 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 + + def clip(self, tolerance: float = 1e-9, inplace: bool = False): + """Remove nuclides with negligible contribution to photon emission. + + 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 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 negligible nuclides removed. + + """ + # Compute per-nuclide emission rate; drop non-emitters + emitting_names = [] + emitting_densities = [] + rates = [] + chain_key = self._chain_file_cache_key() + for name, density in self.nuclides.items(): + integral = DecaySpectrum._photon_integral(name, chain_key) + if integral is None: + continue + emitting_names.append(name) + emitting_densities.append(density) + rates.append(density * self.volume * integral) + + 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 + self._dist_cache_key = None + return self + return type(self)(new_nuclides, self.volume) + + @property + def support(self): + return (0.0, np.inf) + + def evaluate(self, x): + """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): + """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( dists: Sequence[Discrete | Tabular | Mixture], probs: Sequence[float] diff --git a/src/chain.cpp b/src/chain.cpp index e4d0324d3..3915c016c 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 c722d0366..7f5b498ad 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -7,7 +7,9 @@ #include // for accumulate #include // for runtime_error #include // for string, stod +#include +#include "openmc/chain.h" #include "openmc/constants.h" #include "openmc/error.h" #include "openmc/math_functions.h" @@ -15,6 +17,10 @@ #include "openmc/random_lcg.h" #include "openmc/xml_interface.h" +namespace { +std::unordered_set decay_spectrum_missing_chain_nuclides; +} + namespace openmc { //============================================================================== @@ -758,6 +764,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_spectrum") { + dist = UPtrDist {new DecaySpectrum(node)}; } else if (type == "muir") { openmc::fatal_error( "'muir' distributions are now specified using the openmc.stats.muir() " @@ -768,4 +776,120 @@ UPtrDist distribution_from_xml(pugi::xml_node node) return dist; } +//============================================================================== +// DecaySpectrum implementation +//============================================================================== + +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 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 (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); + 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]; + const Distribution* photon_dist = chain_nuc->photon_energy(); + 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; + + nuclide_indices.push_back(nuclide_index); + atoms.push_back(atoms_i); + } + + init(std::move(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 " + "photon spectra and positive atom densities. Ensure " + "OPENMC_CHAIN_FILE is set and matches the nuclides in the " + "source definition."); + } + di_.assign(probs); +} + +DecaySpectrum::Sample DecaySpectrum::sample_with_parent(uint64_t* seed) const +{ + size_t idx = di_.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 +{ + return integral_; +} + +double DecaySpectrum::sample_unbiased(uint64_t* seed) const +{ + return sample_with_parent(seed).energy; +} + } // namespace openmc diff --git a/src/finalize.cpp b/src/finalize.cpp index 98f112534..59711dcc3 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" @@ -44,6 +45,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 efb462f5c..991877b24 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 524a72bf0..bee20d5c3 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -354,6 +354,19 @@ 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, 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 { // Default to a Watt spectrum with parameters 0.988 MeV and 2.249 MeV^-1 energy_ = UPtrDist {new Watt(0.988e6, 2.249e-6)}; @@ -414,6 +427,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()); @@ -424,10 +438,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] && diff --git a/tests/unit_tests/test_source.py b/tests/unit_tests/test_source.py index 6eb03f388..394c09e73 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) diff --git a/tests/unit_tests/test_stats.py b/tests/unit_tests/test_stats.py index ca961c8b0..2754cf5d0 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,135 @@ 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) + 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) + 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 195dff6d1d22cc8dd5ed840d846073ff22fed550 Mon Sep 17 00:00:00 2001 From: Hridoy Kabiraj Date: Mon, 18 May 2026 20:20:59 +0600 Subject: [PATCH 04/10] Fix missing volume checks in material activity/decay heat (#3939) --- openmc/material.py | 5 +++++ tests/unit_tests/test_material.py | 6 ++++++ 2 files changed, 11 insertions(+) diff --git a/openmc/material.py b/openmc/material.py index 2dbe691f5..86cc3d793 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -1418,6 +1418,9 @@ class Material(IDManagerMixin): if volume is None: volume = self.volume + if units in {'Bq', 'Ci'} and volume is None: + raise ValueError(f"Volume must be set in order to compute activity in '{units}'.") + if units == 'Bq': multiplier = volume elif units == 'Bq/cm3': @@ -1474,6 +1477,8 @@ class Material(IDManagerMixin): if units == 'W': multiplier = volume if volume is not None else self.volume + if multiplier is None: + raise ValueError("Volume must be set in order to compute total decay heat.") elif units == 'W/cm3': multiplier = 1 elif units == 'W/m3': diff --git a/tests/unit_tests/test_material.py b/tests/unit_tests/test_material.py index 911b8867f..89dfc03dd 100644 --- a/tests/unit_tests/test_material.py +++ b/tests/unit_tests/test_material.py @@ -561,6 +561,10 @@ def test_get_activity(): m1.add_element("Fe", 0.7) m1.add_element("Li", 0.3) m1.set_density('g/cm3', 1.5) + with pytest.raises(ValueError, match="Volume must be set"): + m1.get_activity(units='Bq') + with pytest.raises(ValueError, match="Volume must be set"): + m1.get_activity(units='Ci') # activity in Bq/cc and Bq/g should not require volume setting assert m1.get_activity(units='Bq/cm3') == 0 assert m1.get_activity(units='Bq/g') == 0 @@ -619,6 +623,8 @@ def test_get_decay_heat(): m1.add_nuclide("U235", 0.2) m1.add_nuclide("U238", 0.8) m1.set_density('g/cm3', 10.5) + with pytest.raises(ValueError, match="Volume must be set"): + m1.get_decay_heat(units='W') # decay heat in W/cc and W/g should not require volume setting assert m1.get_decay_heat(units='W/cm3') == 0 assert m1.get_decay_heat(units='W/g') == 0 From beed56e6ee44037a1543e7ff3cc8815376634c18 Mon Sep 17 00:00:00 2001 From: charliesheh <59520443+charliesheh@users.noreply.github.com> Date: Mon, 18 May 2026 11:35:56 -0700 Subject: [PATCH 05/10] Check for missing thread count argument for -s/--threads (#3940) --- src/initialize.cpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/initialize.cpp b/src/initialize.cpp index 991877b24..c04a2fcbd 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -300,6 +300,11 @@ int parse_command_line(int argc, char* argv[]) settings::run_mode = RunMode::VOLUME; } else if (arg == "-s" || arg == "--threads") { // Read number of threads + if (i + 1 >= argc) { + std::string msg {"Number of threads not specified."}; + strcpy(openmc_err_msg, msg.c_str()); + return OPENMC_E_INVALID_ARGUMENT; + } i += 1; #ifdef _OPENMP From 0169fd9226786863e499a367d71e5eaae6adabb8 Mon Sep 17 00:00:00 2001 From: John Tramm Date: Tue, 19 May 2026 23:23:10 -0500 Subject: [PATCH 06/10] Shared Secondary Particle Bank (#3863) Co-authored-by: John Tramm Co-authored-by: Claude Opus 4.6 Co-authored-by: Paul Romano Co-authored-by: Copilot --- docs/source/io_formats/settings.rst | 59 +-- include/openmc/bank.h | 8 +- include/openmc/event.h | 13 + include/openmc/particle.h | 3 +- include/openmc/particle_data.h | 34 +- include/openmc/settings.h | 3 +- include/openmc/shared_array.h | 30 +- include/openmc/simulation.h | 33 +- openmc/lib/core.py | 5 +- openmc/settings.py | 43 +- src/bank.cpp | 202 +++++++++- src/event.cpp | 51 ++- src/finalize.cpp | 3 + src/ifp.cpp | 4 +- src/initialize.cpp | 59 ++- src/output.cpp | 8 + src/particle.cpp | 178 ++++---- src/particle_restart.cpp | 23 +- src/physics.cpp | 14 +- src/physics_mg.cpp | 9 +- src/settings.cpp | 22 + src/simulation.cpp | 379 +++++++++++++++--- src/tallies/filter_particle_production.cpp | 2 +- src/tallies/tally_scoring.cpp | 9 +- .../__init__.py | 0 .../inputs_true.dat | 65 +++ .../results_true.dat | 11 + .../test.py | 91 +++++ .../particle_production_fission/__init__.py | 0 .../local/inputs_true.dat | 38 ++ .../local/results_true.dat | 3 + .../shared/inputs_true.dat | 38 ++ .../shared/results_true.dat | 3 + .../particle_production_fission/test.py | 49 +++ .../__init__.py | 0 .../geometry.xml | 9 + .../materials.xml | 9 + .../results_true.dat | 16 + .../settings.xml | 15 + .../test.py | 6 + .../pulse_height/{ => local}/inputs_true.dat | 4 +- .../pulse_height/{ => local}/results_true.dat | 0 .../pulse_height/shared/inputs_true.dat | 40 ++ .../pulse_height/shared/results_true.dat | 201 ++++++++++ tests/regression_tests/pulse_height/test.py | 80 ++-- .../weightwindows/local/inputs_true.dat | 94 +++++ .../weightwindows/local/results_true.dat | 1 + .../weightwindows/results_true.dat | 1 - .../{ => shared}/inputs_true.dat | 1 + .../weightwindows/shared/results_true.dat | 1 + .../survival_biasing/local/inputs_true.dat | 69 ++++ .../{ => local}/results_true.dat | 0 .../{ => shared}/inputs_true.dat | 1 + .../survival_biasing/shared/results_true.dat | 1 + .../weightwindows/survival_biasing/test.py | 28 +- tests/regression_tests/weightwindows/test.py | 28 +- .../weightwindows_pulse_height/__init__.py | 0 .../local/inputs_true.dat | 60 +++ .../local/results_true.dat | 201 ++++++++++ .../shared/inputs_true.dat | 60 +++ .../shared/results_true.dat | 201 ++++++++++ .../weightwindows_pulse_height/test.py | 73 ++++ tests/unit_tests/weightwindows/test.py | 16 +- 63 files changed, 2409 insertions(+), 299 deletions(-) create mode 100644 tests/regression_tests/mg_fixed_source_ww_fission_shared_secondary/__init__.py create mode 100644 tests/regression_tests/mg_fixed_source_ww_fission_shared_secondary/inputs_true.dat create mode 100644 tests/regression_tests/mg_fixed_source_ww_fission_shared_secondary/results_true.dat create mode 100644 tests/regression_tests/mg_fixed_source_ww_fission_shared_secondary/test.py create mode 100644 tests/regression_tests/particle_production_fission/__init__.py create mode 100644 tests/regression_tests/particle_production_fission/local/inputs_true.dat create mode 100644 tests/regression_tests/particle_production_fission/local/results_true.dat create mode 100644 tests/regression_tests/particle_production_fission/shared/inputs_true.dat create mode 100644 tests/regression_tests/particle_production_fission/shared/results_true.dat create mode 100644 tests/regression_tests/particle_production_fission/test.py create mode 100644 tests/regression_tests/particle_restart_fixed_shared_secondary/__init__.py create mode 100644 tests/regression_tests/particle_restart_fixed_shared_secondary/geometry.xml create mode 100644 tests/regression_tests/particle_restart_fixed_shared_secondary/materials.xml create mode 100644 tests/regression_tests/particle_restart_fixed_shared_secondary/results_true.dat create mode 100644 tests/regression_tests/particle_restart_fixed_shared_secondary/settings.xml create mode 100644 tests/regression_tests/particle_restart_fixed_shared_secondary/test.py rename tests/regression_tests/pulse_height/{ => local}/inputs_true.dat (95%) rename tests/regression_tests/pulse_height/{ => local}/results_true.dat (100%) create mode 100644 tests/regression_tests/pulse_height/shared/inputs_true.dat create mode 100644 tests/regression_tests/pulse_height/shared/results_true.dat create mode 100644 tests/regression_tests/weightwindows/local/inputs_true.dat create mode 100644 tests/regression_tests/weightwindows/local/results_true.dat delete mode 100644 tests/regression_tests/weightwindows/results_true.dat rename tests/regression_tests/weightwindows/{ => shared}/inputs_true.dat (99%) create mode 100644 tests/regression_tests/weightwindows/shared/results_true.dat create mode 100644 tests/regression_tests/weightwindows/survival_biasing/local/inputs_true.dat rename tests/regression_tests/weightwindows/survival_biasing/{ => local}/results_true.dat (100%) rename tests/regression_tests/weightwindows/survival_biasing/{ => shared}/inputs_true.dat (99%) create mode 100644 tests/regression_tests/weightwindows/survival_biasing/shared/results_true.dat create mode 100644 tests/regression_tests/weightwindows_pulse_height/__init__.py create mode 100644 tests/regression_tests/weightwindows_pulse_height/local/inputs_true.dat create mode 100644 tests/regression_tests/weightwindows_pulse_height/local/results_true.dat create mode 100644 tests/regression_tests/weightwindows_pulse_height/shared/inputs_true.dat create mode 100644 tests/regression_tests/weightwindows_pulse_height/shared/results_true.dat create mode 100644 tests/regression_tests/weightwindows_pulse_height/test.py diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index f8b064a37..1a436a75e 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -741,14 +741,16 @@ pseudo-random number generator. *Default*: 1 --------------------- -```` Element --------------------- +----------------------------------- +```` Element +----------------------------------- -The ``stride`` element is used to specify how many random numbers are allocated -for each source particle history. - - *Default*: 152,917 + The ``shared_secondary_bank`` element indicates whether to use a shared + secondary particle bank. When enabled, secondary particles are collected into + a global bank, sorted for reproducibility, and load-balanced across MPI ranks + between generations. If not specified, the shared secondary bank is enabled + automatically for fixed-source simulations with weight windows active, and + disabled otherwise. .. _source_element: @@ -1156,23 +1158,6 @@ based on constraints. *Default*: 0.05 -------------------------- -```` Element -------------------------- - -The ```` element indicates at what batches a state point file -should be written. A state point file can be used to restart a run or to get -tally results at any batch. The default behavior when using this tag is to -write out the source bank in the state_point file. This behavior can be -customized by using the ```` element. This element has the -following attributes/sub-elements: - - :batches: - A list of integers separated by spaces indicating at what batches a state - point file should be written. - - *Default*: Last batch only - -------------------------- ```` Element -------------------------- @@ -1222,6 +1207,32 @@ attributes/sub-elements: *Default*: false +------------------------- +```` Element +------------------------- + +The ```` element indicates at what batches a state point file +should be written. A state point file can be used to restart a run or to get +tally results at any batch. The default behavior when using this tag is to +write out the source bank in the state_point file. This behavior can be +customized by using the ```` element. This element has the +following attributes/sub-elements: + + :batches: + A list of integers separated by spaces indicating at what batches a state + point file should be written. + + *Default*: Last batch only + +-------------------- +```` Element +-------------------- + +The ``stride`` element is used to specify how many random numbers are allocated +for each source particle history. + + *Default*: 152,917 + ------------------------------ ```` Element ------------------------------ diff --git a/include/openmc/bank.h b/include/openmc/bank.h index c4e940bc8..6abcdd7f1 100644 --- a/include/openmc/bank.h +++ b/include/openmc/bank.h @@ -34,18 +34,24 @@ extern vector> ifp_fission_lifetime_bank; extern vector progeny_per_particle; +extern SharedArray shared_secondary_bank_read; +extern SharedArray shared_secondary_bank_write; + } // namespace simulation //============================================================================== // Non-member functions //============================================================================== -void sort_fission_bank(); +void sort_bank(SharedArray& bank, bool is_fission_bank); void free_memory_bank(); void init_fission_bank(int64_t max); +int64_t synchronize_global_secondary_bank( + SharedArray& shared_secondary_bank); + } // namespace openmc #endif // OPENMC_BANK_H diff --git a/include/openmc/event.h b/include/openmc/event.h index 2d215a10e..48a009f2b 100644 --- a/include/openmc/event.h +++ b/include/openmc/event.h @@ -112,6 +112,19 @@ void process_collision_events(); //! \param n_particles The number of particles in the particle buffer void process_death_events(int64_t n_particles); +//! Process event queues until all are empty. Each iteration processes the +//! longest queue first to maximize vectorization efficiency. +void process_transport_events(); + +//! Initialize secondary particles from a shared secondary bank for +//! event-based transport +// +//! \param n_particles The number of particles to initialize +//! \param offset The offset index in the shared secondary bank +//! \param shared_secondary_bank The shared secondary bank to read from +void process_init_secondary_events(int64_t n_particles, int64_t offset, + const SharedArray& shared_secondary_bank); + } // namespace openmc #endif // OPENMC_EVENT_H diff --git a/include/openmc/particle.h b/include/openmc/particle.h index e75f3785a..8db9721ba 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -71,7 +71,8 @@ public: void event_advance(); void event_cross_surface(); void event_collide(); - void event_revive_from_secondary(); + void event_revive_from_secondary(const SourceSite& site); + void event_check_limit_and_revive(); void event_death(); //! pulse-height recording diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index 5b632bbce..f72948f6e 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -50,8 +50,11 @@ struct SourceSite { // Extra attributes that don't show up in source written to file int parent_nuclide {-1}; - int64_t parent_id; - int64_t progeny_id; + int64_t parent_id {0}; + int64_t progeny_id {0}; + double wgt_born {1.0}; + double wgt_ww_born {-1.0}; + int64_t n_split {0}; }; struct CollisionTrackSite { @@ -533,14 +536,14 @@ private: uint64_t seeds_[N_STREAMS]; int stream_; - vector secondary_bank_; + vector local_secondary_bank_; // Keep track of how many secondary particles were created in the collision // and what the starting index is in the secondary bank for this particle int n_secondaries_ {0}; int secondary_bank_index_ {0}; - int64_t current_work_; + int64_t current_work_ {0}; vector flux_derivs_; @@ -563,7 +566,9 @@ private: int n_event_ {0}; - int n_split_ {0}; + int64_t n_tracks_ {0}; //!< number of tracks in this particle history + + int64_t n_split_ {0}; double ww_factor_ {0.0}; int64_t n_progeny_ {0}; @@ -693,12 +698,14 @@ public: int& stream() { return stream_; } // secondary particle bank - SourceSite& secondary_bank(int i) { return secondary_bank_[i]; } - const SourceSite& secondary_bank(int i) const { return secondary_bank_[i]; } - decltype(secondary_bank_)& secondary_bank() { return secondary_bank_; } - decltype(secondary_bank_) const& secondary_bank() const + SourceSite& local_secondary_bank(int i) { return local_secondary_bank_[i]; } + const SourceSite& local_secondary_bank(int i) const { - return secondary_bank_; + return local_secondary_bank_[i]; + } + decltype(local_secondary_bank_)& local_secondary_bank() + { + return local_secondary_bank_; } // Number of secondaries created in a collision @@ -747,13 +754,16 @@ public: int& n_event() { return n_event_; } // Number of times variance reduction has caused a particle split - int n_split() const { return n_split_; } - int& n_split() { return n_split_; } + int64_t n_split() const { return n_split_; } + int64_t& n_split() { return n_split_; } // Particle-specific factor for on-the-fly weight window adjustment double ww_factor() const { return ww_factor_; } double& ww_factor() { return ww_factor_; } + // Number of tracks in this particle history + int64_t& n_tracks() { return n_tracks_; } + // Number of progeny produced by this particle int64_t& n_progeny() { return n_progeny_; } diff --git a/include/openmc/settings.h b/include/openmc/settings.h index 0914a0958..3bba040f5 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -98,7 +98,8 @@ extern bool uniform_source_sampling; //!< sample sources uniformly? extern bool ufs_on; //!< uniform fission site method on? extern bool urr_ptables_on; //!< use unresolved resonance prob. tables? extern bool use_decay_photons; //!< use decay photons for D1S -extern "C" bool weight_windows_on; //!< are weight windows are enabled? +extern bool use_shared_secondary_bank; //!< Use shared bank for secondaries +extern "C" bool weight_windows_on; //!< are weight windows are enabled? extern bool weight_window_checkpoint_surface; //!< enable weight window check //!< upon surface crossing? extern bool weight_window_checkpoint_collision; //!< enable weight window check diff --git a/include/openmc/shared_array.h b/include/openmc/shared_array.h index 7e9ef28c5..b309ca3f1 100644 --- a/include/openmc/shared_array.h +++ b/include/openmc/shared_array.h @@ -4,6 +4,8 @@ //! \file shared_array.h //! \brief Shared array data structure +#include // for copy_n + #include "openmc/memory.h" namespace openmc { @@ -30,14 +32,12 @@ public: //! Default constructor. SharedArray() = default; - //! Construct a zero size container with space to hold capacity number of - //! elements. + //! Construct a container with `size` elements and capacity equal to `size`. // - //! \param capacity The number of elements for the container to allocate - //! space for - SharedArray(int64_t capacity) : capacity_(capacity) + //! \param size The number of elements to allocate and initialize + SharedArray(int64_t size) : size_(size), capacity_(size) { - data_ = make_unique(capacity); + data_ = make_unique(size); } //========================================================================== @@ -97,8 +97,26 @@ public: capacity_ = 0; } + //! Push back an element to the array, with capacity and reallocation behavior + //! as if this were a vector. This does not perform any thread safety checks. + //! If the size exceeds the capacity, then the capacity will double just as + //! with a vector. Data will be reallocated and moved to a new pointer and + //! copied in before the new item is appended. Old data will be freed. + void thread_unsafe_append(const T& value) + { + if (size_ == capacity_) { + int64_t new_capacity = capacity_ == 0 ? 8 : 2 * capacity_; + unique_ptr new_data = make_unique(new_capacity); + std::copy_n(data_.get(), size_, new_data.get()); + data_ = std::move(new_data); + capacity_ = new_capacity; + } + data_[size_++] = value; + } + //! Return the number of elements in the container int64_t size() { return size_; } + int64_t size() const { return size_; } //! Resize the container to contain a specified number of elements. This is //! useful in cases where the container is written to in a non-thread safe diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index 9a6cf1b21..454752cd2 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -49,6 +49,9 @@ extern const RegularMesh* ufs_mesh; extern vector k_generation; extern vector work_index; +extern int64_t + simulation_tracks_completed; //!< Number of tracks completed on this rank + } // namespace simulation //============================================================================== @@ -59,7 +62,7 @@ extern vector work_index; void allocate_banks(); //! Determine number of particles to transport per process -void calculate_work(); +void calculate_work(int64_t n_particles); //! Initialize nuclear data before a simulation void initialize_data(); @@ -70,8 +73,9 @@ void initialize_batch(); //! Initialize a fission generation void initialize_generation(); -//! Full initialization of a particle history -void initialize_history(Particle& p, int64_t index_source); +//! Full initialization of a particle track +void initialize_particle_track( + Particle& p, int64_t index_source, bool is_secondary); //! Finalize a batch //! @@ -92,16 +96,35 @@ void broadcast_results(); void free_memory_simulation(); -//! Simulate a single particle history (and all generated secondary particles, -//! if enabled), from birth to death +//! Compute unique particle ID from a 1-based source index +//! \param index_source 1-based source index within this rank's work +//! \return globally unique particle ID +int64_t compute_particle_id(int64_t index_source); + +//! Compute the transport RNG seed from a particle ID +//! \param particle_id the particle's globally unique ID +//! \return seed value passed to init_particle_seeds() +int64_t compute_transport_seed(int64_t particle_id); + +//! Simulate a single particle history from birth to death, inclusive of any +//! secondary particles. In shared secondary mode, only a single track is +//! transported and secondaries are deposited into a shared bank instead. void transport_history_based_single_particle(Particle& p); //! Simulate all particle histories using history-based parallelism void transport_history_based(); +//! Simulate all particles using history-based parallelism, with a shared +//! secondary bank +void transport_history_based_shared_secondary(); + //! Simulate all particle histories using event-based parallelism void transport_event_based(); +//! Simulate all particles using event-based parallelism, with a shared +//! secondary bank +void transport_event_based_shared_secondary(); + } // namespace openmc #endif // OPENMC_SIMULATION_H diff --git a/openmc/lib/core.py b/openmc/lib/core.py index 00876db99..0ee8c3ff0 100644 --- a/openmc/lib/core.py +++ b/openmc/lib/core.py @@ -31,7 +31,10 @@ class _SourceSite(Structure): ('particle', c_int32), ('parent_nuclide', c_int), ('parent_id', c_int64), - ('progeny_id', c_int64)] + ('progeny_id', c_int64), + ('wgt_born', c_double), + ('wgt_ww_born', c_double), + ('n_split', c_int64)] # Define input type for numpy arrays that will be passed into C++ functions # Must be an int or double array, with single dimension that is contiguous diff --git a/openmc/settings.py b/openmc/settings.py index 1090babda..8120eb073 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -67,6 +67,10 @@ class Settings: (ex: ["(n,fission)", 2, "(n,2n)"] ). (list of str or int) :deposited_E_threshold: Number to define the minimum deposited energy during per collision to trigger banking. (float) + create_delayed_neutrons : bool + Whether delayed neutrons are created in fission. + + .. versionadded:: 0.13.3 create_fission_neutrons : bool Indicate whether fission neutrons should be created or not. cutoff : dict @@ -256,8 +260,15 @@ class Settings: The type of calculation to perform (default is 'eigenvalue') seed : int Seed for the linear congruential pseudorandom number generator - stride : int - Number of random numbers allocated for each source particle history + shared_secondary_bank : bool + Whether to use a shared secondary particle bank. When enabled, + secondary particles are collected into a global bank, sorted for + reproducibility, and load-balanced across MPI ranks between + generations. If not specified, the shared secondary bank is + enabled automatically for fixed-source simulations with weight + windows active, and disabled otherwise. + + .. versionadded:: 0.15.4 source : Iterable of openmc.SourceBase Distribution of source sites in space, angle, and energy source_rejection_fraction : float @@ -277,6 +288,8 @@ class Settings: Options for writing state points. Acceptable keys are: :batches: list of batches at which to write statepoint files + stride : int + Number of random numbers allocated for each source particle history surf_source_read : dict Options for reading surface source points. Acceptable keys are: @@ -373,10 +386,6 @@ class Settings: .. versionadded:: 0.14.0 - create_delayed_neutrons : bool - Whether delayed neutrons are created in fission. - - .. versionadded:: 0.13.3 weight_windows_on : bool Whether weight windows are enabled @@ -482,6 +491,7 @@ class Settings: self._weight_window_generators = cv.CheckedList( WeightWindowGenerator, 'weight window generators') self._weight_windows_on = None + self._shared_secondary_bank = None self._weight_windows_file = None self._weight_window_checkpoints = {} self._max_history_splits = None @@ -1306,6 +1316,15 @@ class Settings: cv.check_type('weight windows on', value, bool) self._weight_windows_on = value + @property + def shared_secondary_bank(self) -> bool: + return self._shared_secondary_bank + + @shared_secondary_bank.setter + def shared_secondary_bank(self, value: bool): + cv.check_type('shared secondary bank', value, bool) + self._shared_secondary_bank = value + @property def weight_window_checkpoints(self) -> dict: return self._weight_window_checkpoints @@ -1924,6 +1943,11 @@ class Settings: elem = ET.SubElement(root, "weight_windows_on") elem.text = str(self._weight_windows_on).lower() + def _create_shared_secondary_bank_subelement(self, root): + if self._shared_secondary_bank is not None: + elem = ET.SubElement(root, "shared_secondary_bank") + elem.text = str(self._shared_secondary_bank).lower() + def _create_weight_window_generators_subelement(self, root, mesh_memo=None): if not self.weight_window_generators: return @@ -2430,6 +2454,11 @@ class Settings: if text is not None: self.weight_windows_on = text in ('true', '1') + def _shared_secondary_bank_from_xml_element(self, root): + text = get_text(root, 'shared_secondary_bank') + if text is not None: + self.shared_secondary_bank = text in ('true', '1') + def _weight_windows_file_from_xml_element(self, root): text = get_text(root, 'weight_windows_file') if text is not None: @@ -2597,6 +2626,7 @@ class Settings: self._create_write_initial_source_subelement(element) self._create_weight_windows_subelement(element, mesh_memo) self._create_weight_windows_on_subelement(element) + self._create_shared_secondary_bank_subelement(element) self._create_weight_window_generators_subelement(element, mesh_memo) self._create_weight_windows_file_element(element) self._create_weight_window_checkpoints_subelement(element) @@ -2714,6 +2744,7 @@ class Settings: settings._write_initial_source_from_xml_element(elem) settings._weight_windows_from_xml_element(elem, meshes) settings._weight_windows_on_from_xml_element(elem) + settings._shared_secondary_bank_from_xml_element(elem) settings._weight_windows_file_from_xml_element(elem) settings._weight_window_generators_from_xml_element(elem, meshes) settings._weight_window_checkpoints_from_xml_element(elem) diff --git a/src/bank.cpp b/src/bank.cpp index 5afdc2eca..5b12b48fd 100644 --- a/src/bank.cpp +++ b/src/bank.cpp @@ -7,6 +7,7 @@ #include "openmc/vector.h" #include +#include #include namespace openmc { @@ -43,6 +44,13 @@ vector> ifp_fission_lifetime_bank; // used to efficiently sort the fission bank after each iteration. vector progeny_per_particle; +// When shared secondary bank mode is enabled, secondaries produced during +// transport are collected in the write bank. When a secondary generation is +// complete, write is moved to read for transport, and a new empty write bank +// is created. This repeats until no secondaries remain. +SharedArray shared_secondary_bank_read; +SharedArray shared_secondary_bank_write; + } // namespace simulation //============================================================================== @@ -60,6 +68,8 @@ void free_memory_bank() simulation::ifp_source_lifetime_bank.clear(); simulation::ifp_fission_delayed_group_bank.clear(); simulation::ifp_fission_lifetime_bank.clear(); + simulation::shared_secondary_bank_read.clear(); + simulation::shared_secondary_bank_write.clear(); } void init_fission_bank(int64_t max) @@ -68,13 +78,13 @@ void init_fission_bank(int64_t max) simulation::progeny_per_particle.resize(simulation::work_per_rank); } -// Performs an O(n) sort on the fission bank, by leveraging +// Performs an O(n) sort on a fission or secondary bank, by leveraging // the parent_id and progeny_id fields of banked particles. See the following // paper for more details: // "Reproducibility and Monte Carlo Eigenvalue Calculations," F.B. Brown and // T.M. Sutton, 1992 ANS Annual Meeting, Transactions of the American Nuclear // Society, Volume 65, Page 235. -void sort_fission_bank() +void sort_bank(SharedArray& bank, bool is_fission_bank) { // Ensure we don't read off the end of the array if we ran with 0 particles if (simulation::progeny_per_particle.size() == 0) { @@ -96,44 +106,200 @@ void sort_fission_bank() vector> sorted_ifp_lifetime_bank; // If there is not enough space, allocate a temporary vector and point to it - if (simulation::fission_bank.size() > - simulation::fission_bank.capacity() / 2) { - sorted_bank_holder.resize(simulation::fission_bank.size()); + if (bank.size() > bank.capacity() / 2) { + sorted_bank_holder.resize(bank.size()); sorted_bank = sorted_bank_holder.data(); } else { // otherwise, point sorted_bank to unused portion of the fission bank - sorted_bank = &simulation::fission_bank[simulation::fission_bank.size()]; + sorted_bank = bank.data() + bank.size(); } - if (settings::ifp_on) { + if (settings::ifp_on && is_fission_bank) { allocate_temporary_vector_ifp( sorted_ifp_delayed_group_bank, sorted_ifp_lifetime_bank); } - // Use parent and progeny indices to sort fission bank - for (int64_t i = 0; i < simulation::fission_bank.size(); i++) { - const auto& site = simulation::fission_bank[i]; - int64_t offset = site.parent_id - 1 - simulation::work_index[mpi::rank]; - int64_t idx = simulation::progeny_per_particle[offset] + site.progeny_id; - if (idx >= simulation::fission_bank.size()) { + // Use parent and progeny indices to sort bank + for (int64_t i = 0; i < bank.size(); i++) { + const auto& site = bank[i]; + if (site.parent_id < 0 || + site.parent_id >= + static_cast(simulation::progeny_per_particle.size())) { + fatal_error(fmt::format("Invalid parent_id {} for banked site (expected " + "range [0, {})).", + site.parent_id, simulation::progeny_per_particle.size())); + } + int64_t idx = + simulation::progeny_per_particle[site.parent_id] + site.progeny_id; + if (idx < 0 || idx >= bank.size()) { fatal_error("Mismatch detected between sum of all particle progeny and " - "shared fission bank size."); + "bank size during sorting."); } sorted_bank[idx] = site; - if (settings::ifp_on) { + if (settings::ifp_on && is_fission_bank) { copy_ifp_data_from_fission_banks( i, sorted_ifp_delayed_group_bank[idx], sorted_ifp_lifetime_bank[idx]); } } // Copy sorted bank into the fission bank - std::copy(sorted_bank, sorted_bank + simulation::fission_bank.size(), - simulation::fission_bank.data()); - if (settings::ifp_on) { + std::copy(sorted_bank, sorted_bank + bank.size(), bank.data()); + if (settings::ifp_on && is_fission_bank) { copy_ifp_data_to_fission_banks( sorted_ifp_delayed_group_bank.data(), sorted_ifp_lifetime_bank.data()); } } +// This function redistributes SourceSite particles across MPI ranks to +// achieve load balancing while preserving the global ordering of particles. +// +// GUARANTEES: +// ----------- +// 1. Global Order Preservation: After redistribution, each rank holds a +// contiguous slice of the original global ordering. For example, if the +// input across 3 ranks was: +// - Rank 0: IDs 0-4 +// - Rank 1: IDs 5-6 +// - Rank 2: IDs 7-200 +// Then after redistribution (assuming ~67 particles per rank): +// - Rank 0: IDs 0-66 (contiguous) +// - Rank 1: IDs 67-133 (contiguous) +// - Rank 2: IDs 134-200 (contiguous) +// The global ordering is always preserved - no rank will ever hold +// non-contiguous ID ranges like "0-4 and 100-200". +// +// 2. Even Load Balancing: Particles are distributed as evenly as possible. +// If total % n_procs != 0, the first 'remainder' ranks each get one extra +// particle (i.e., floor division with remainder distributed to lower +// ranks). This follows the same logic as calculate_work(). +// +// HOW IT WORKS: +// ------------- +// The algorithm uses overlap-based redistribution: +// 1. Each rank's current data occupies a range [cumulative_before[rank], +// cumulative_before[rank+1]) in the global index space. +// 2. Each rank's target data should occupy [cumulative_target[rank], +// cumulative_target[rank+1]) in the same global index space. +// 3. For each pair of (source_rank, dest_rank), we calculate the overlap +// between what source_rank currently has and what dest_rank needs. +// 4. MPI_Alltoallv transfers exactly these overlapping regions, with +// displacements ensuring data lands at the correct position in the +// receiving buffer. +// +// EDGE CASES HANDLED: +// ------------------- +// - Single rank (n_procs == 1): Returns immediately with local size, no MPI. +// - Empty total (all ranks have 0 particles): Returns 0 immediately. +// - Imbalanced input (e.g., one rank has all particles): Works correctly; +// that rank will send portions to all other ranks based on target ranges. +// - Non-divisible totals: First 'remainder' ranks get one extra particle. +int64_t synchronize_global_secondary_bank( + SharedArray& shared_secondary_bank) +{ + // Get current size of local bank + int64_t local_size = shared_secondary_bank.size(); + + if (mpi::n_procs == 1) { + return local_size; + } + +#ifdef OPENMC_MPI + // Gather all sizes to all ranks + vector all_sizes(mpi::n_procs); + MPI_Allgather(&local_size, 1, MPI_INT64_T, all_sizes.data(), 1, MPI_INT64_T, + mpi::intracomm); + + // Calculate total and check for empty case + int64_t total = 0; + for (int64_t size : all_sizes) { + total += size; + } + + // If we don't have any items to distribute, return + if (total == 0) { + return total; + } + + int64_t base_count = total / mpi::n_procs; + int64_t remainder = total % mpi::n_procs; + + // Calculate target size for each rank + // First 'remainder' ranks get base_count + 1, rest get base_count + vector target_sizes(mpi::n_procs); + for (int i = 0; i < mpi::n_procs; ++i) { + target_sizes[i] = base_count + (i < remainder ? 1 : 0); + } + + // Calculate send and receive counts in terms of SourceSite objects + // (not bytes) + vector send_counts(mpi::n_procs, 0); + vector recv_counts(mpi::n_procs, 0); + vector send_displs(mpi::n_procs, 0); + vector recv_displs(mpi::n_procs, 0); + + // Calculate cumulative positions (starting index for each rank in the + // global array) + vector cumulative_before(mpi::n_procs + 1, 0); + vector cumulative_target(mpi::n_procs + 1, 0); + for (int i = 0; i < mpi::n_procs; ++i) { + cumulative_before[i + 1] = cumulative_before[i] + all_sizes[i]; + cumulative_target[i + 1] = cumulative_target[i] + target_sizes[i]; + } + + // Determine send and receive amounts for each rank + int64_t my_start = cumulative_before[mpi::rank]; + int64_t my_end = cumulative_before[mpi::rank + 1]; + int64_t my_target_start = cumulative_target[mpi::rank]; + int64_t my_target_end = cumulative_target[mpi::rank + 1]; + + for (int r = 0; r < mpi::n_procs; ++r) { + // Send: overlap between my current range and rank r's target range + int64_t send_overlap_start = std::max(my_start, cumulative_target[r]); + int64_t send_overlap_end = std::min(my_end, cumulative_target[r + 1]); + if (send_overlap_start < send_overlap_end) { + int64_t count = send_overlap_end - send_overlap_start; + int64_t displ = send_overlap_start - my_start; + if (count > std::numeric_limits::max() || + displ > std::numeric_limits::max()) { + fatal_error("Secondary bank size exceeds MPI_Alltoallv int limit."); + } + send_counts[r] = static_cast(count); + send_displs[r] = static_cast(displ); + } + + // Recv: overlap between rank r's current range and my target range + int64_t recv_overlap_start = + std::max(cumulative_before[r], my_target_start); + int64_t recv_overlap_end = + std::min(cumulative_before[r + 1], my_target_end); + if (recv_overlap_start < recv_overlap_end) { + int64_t count = recv_overlap_end - recv_overlap_start; + int64_t displ = recv_overlap_start - my_target_start; + if (count > std::numeric_limits::max() || + displ > std::numeric_limits::max()) { + fatal_error("Secondary bank size exceeds MPI_Alltoallv int limit."); + } + recv_counts[r] = static_cast(count); + recv_displs[r] = static_cast(displ); + } + } + + // Prepare receive buffer with target size + SharedArray new_bank(target_sizes[mpi::rank]); + + // Perform all-to-all redistribution using the custom MPI type + MPI_Alltoallv(shared_secondary_bank.data(), send_counts.data(), + send_displs.data(), mpi::source_site, new_bank.data(), recv_counts.data(), + recv_displs.data(), mpi::source_site, mpi::intracomm); + + // Replace old bank with redistributed data + shared_secondary_bank = std::move(new_bank); + + return total; +#else + return local_size; +#endif +} + //============================================================================== // C API //============================================================================== diff --git a/src/event.cpp b/src/event.cpp index f33e132d0..2d436bb9d 100644 --- a/src/event.cpp +++ b/src/event.cpp @@ -1,6 +1,9 @@ #include "openmc/event.h" +#include "openmc/bank.h" +#include "openmc/error.h" #include "openmc/material.h" +#include "openmc/settings.h" #include "openmc/simulation.h" #include "openmc/timer.h" @@ -64,7 +67,8 @@ void process_init_events(int64_t n_particles, int64_t source_offset) simulation::time_event_init.start(); #pragma omp parallel for schedule(runtime) for (int64_t i = 0; i < n_particles; i++) { - initialize_history(simulation::particles[i], source_offset + i + 1); + initialize_particle_track( + simulation::particles[i], source_offset + i + 1, false); dispatch_xs_event(i); } simulation::time_event_init.stop(); @@ -136,7 +140,7 @@ void process_surface_crossing_events() int64_t buffer_idx = simulation::surface_crossing_queue[i].idx; Particle& p = simulation::particles[buffer_idx]; p.event_cross_surface(); - p.event_revive_from_secondary(); + p.event_check_limit_and_revive(); if (p.alive()) dispatch_xs_event(buffer_idx); } @@ -155,7 +159,7 @@ void process_collision_events() int64_t buffer_idx = simulation::collision_queue[i].idx; Particle& p = simulation::particles[buffer_idx]; p.event_collide(); - p.event_revive_from_secondary(); + p.event_check_limit_and_revive(); if (p.alive()) dispatch_xs_event(buffer_idx); } @@ -176,4 +180,45 @@ void process_death_events(int64_t n_particles) simulation::time_event_death.stop(); } +void process_transport_events() +{ + while (true) { + int64_t max = std::max({simulation::calculate_fuel_xs_queue.size(), + simulation::calculate_nonfuel_xs_queue.size(), + simulation::advance_particle_queue.size(), + simulation::surface_crossing_queue.size(), + simulation::collision_queue.size()}); + + if (max == 0) { + break; + } else if (max == simulation::calculate_fuel_xs_queue.size()) { + process_calculate_xs_events(simulation::calculate_fuel_xs_queue); + } else if (max == simulation::calculate_nonfuel_xs_queue.size()) { + process_calculate_xs_events(simulation::calculate_nonfuel_xs_queue); + } else if (max == simulation::advance_particle_queue.size()) { + process_advance_particle_events(); + } else if (max == simulation::surface_crossing_queue.size()) { + process_surface_crossing_events(); + } else if (max == simulation::collision_queue.size()) { + process_collision_events(); + } + } +} + +void process_init_secondary_events(int64_t n_particles, int64_t offset, + const SharedArray& shared_secondary_bank) +{ + simulation::time_event_init.start(); +#pragma omp parallel for schedule(runtime) + for (int64_t i = 0; i < n_particles; i++) { + initialize_particle_track(simulation::particles[i], offset + i + 1, true); + const SourceSite& site = shared_secondary_bank[offset + i]; + simulation::particles[i].event_revive_from_secondary(site); + if (simulation::particles[i].alive()) { + dispatch_xs_event(i); + } + } + simulation::time_event_init.stop(); +} + } // namespace openmc diff --git a/src/finalize.cpp b/src/finalize.cpp index 59711dcc3..fd891d9dd 100644 --- a/src/finalize.cpp +++ b/src/finalize.cpp @@ -149,6 +149,8 @@ int openmc_finalize() settings::uniform_source_sampling = false; settings::ufs_on = false; settings::urr_ptables_on = true; + settings::use_decay_photons = false; + settings::use_shared_secondary_bank = false; settings::verbosity = -1; settings::weight_cutoff = 0.25; settings::weight_survive = 1.0; @@ -219,6 +221,7 @@ int openmc_reset() settings::cmfd_run = false; simulation::n_lost_particles = 0; + simulation::simulation_tracks_completed = 0; return 0; } diff --git a/src/ifp.cpp b/src/ifp.cpp index a5157440d..f0f98b8ef 100644 --- a/src/ifp.cpp +++ b/src/ifp.cpp @@ -32,13 +32,13 @@ void ifp(const Particle& p, int64_t idx) { if (is_beta_effective_or_both()) { const auto& delayed_groups = - simulation::ifp_source_delayed_group_bank[p.current_work() - 1]; + simulation::ifp_source_delayed_group_bank[p.current_work()]; simulation::ifp_fission_delayed_group_bank[idx] = _ifp(p.delayed_group(), delayed_groups); } if (is_generation_time_or_both()) { const auto& lifetimes = - simulation::ifp_source_lifetime_bank[p.current_work() - 1]; + simulation::ifp_source_lifetime_bank[p.current_work()]; simulation::ifp_fission_lifetime_bank[idx] = _ifp(p.lifetime(), lifetimes); } } diff --git a/src/initialize.cpp b/src/initialize.cpp index c04a2fcbd..78b414f78 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -161,7 +161,7 @@ void initialize_mpi(MPI_Comm intracomm) // Create bank datatype SourceSite b; - MPI_Aint disp[11]; + MPI_Aint disp[14]; MPI_Get_address(&b.r, &disp[0]); MPI_Get_address(&b.u, &disp[1]); MPI_Get_address(&b.E, &disp[2]); @@ -173,14 +173,35 @@ void initialize_mpi(MPI_Comm intracomm) MPI_Get_address(&b.parent_nuclide, &disp[8]); MPI_Get_address(&b.parent_id, &disp[9]); MPI_Get_address(&b.progeny_id, &disp[10]); - for (int i = 10; i >= 0; --i) { + MPI_Get_address(&b.wgt_born, &disp[11]); + MPI_Get_address(&b.wgt_ww_born, &disp[12]); + MPI_Get_address(&b.n_split, &disp[13]); + for (int i = 13; i >= 0; --i) { disp[i] -= disp[0]; } - int blocks[] {3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1}; - MPI_Datatype types[] {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, - MPI_DOUBLE, MPI_INT, MPI_INT, MPI_INT, MPI_INT, MPI_LONG, MPI_LONG}; - MPI_Type_create_struct(11, blocks, disp, types, &mpi::source_site); + // Block counts for each field + int blocks[] = {3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}; + + // Types for each field + MPI_Datatype types[] = { + MPI_DOUBLE, // r (3 doubles) + MPI_DOUBLE, // u (3 doubles) + MPI_DOUBLE, // E + MPI_DOUBLE, // time + MPI_DOUBLE, // wgt + MPI_INT, // delayed_group + MPI_INT, // surf_id + MPI_INT, // particle (enum) + MPI_INT, // parent_nuclide + MPI_INT64_T, // parent_id + MPI_INT64_T, // progeny_id + MPI_DOUBLE, // wgt_born + MPI_DOUBLE, // wgt_ww_born + MPI_INT64_T // n_split + }; + + MPI_Type_create_struct(14, blocks, disp, types, &mpi::source_site); MPI_Type_commit(&mpi::source_site); CollisionTrackSite bc; @@ -366,6 +387,28 @@ int parse_command_line(int argc, char* argv[]) return 0; } +// TODO: Pulse-height tallies require per-history scoring across the full +// particle tree (parent + all descendants). The shared secondary bank +// transports each secondary as an independent Particle, breaking this +// assumption. A proper fix would defer pulse-height scoring: save +// (root_source_id, cell, pht_storage) per particle, then aggregate by +// root_source_id after all secondary generations complete before scoring +// into the histogram. For now, disable shared secondary when pulse-height +// tallies are present. +static void check_pulse_height_compatibility() +{ + if (settings::use_shared_secondary_bank) { + for (const auto& t : model::tallies) { + if (t->type_ == TallyType::PULSE_HEIGHT) { + settings::use_shared_secondary_bank = false; + warning("Pulse-height tallies are not yet compatible with the shared " + "secondary bank. Disabling shared secondary bank."); + break; + } + } + } +} + bool read_model_xml() { std::string model_filename = settings::path_input; @@ -460,6 +503,8 @@ bool read_model_xml() if (check_for_node(root, "tallies")) read_tallies_xml(root.child("tallies")); + check_pulse_height_compatibility(); + // Initialize distribcell_filters prepare_distribcell(); @@ -505,6 +550,8 @@ void read_separate_xml_files() read_tallies_xml(); + check_pulse_height_compatibility(); + // Initialize distribcell_filters prepare_distribcell(); diff --git a/src/output.cpp b/src/output.cpp index f8c0f2a97..725e8d069 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -501,6 +501,14 @@ void print_runtime() show_rate("Calculation Rate (inactive)", speed_inactive); } show_rate("Calculation Rate (active)", speed_active); + + // Display track rate when weight windows are enabled + if (settings::weight_windows_on) { + double speed_tracks = + simulation::simulation_tracks_completed / time_active.elapsed(); + fmt::print( + " {:<33} = {:.6} tracks/second\n", "Track Rate (active)", speed_tracks); + } } //============================================================================== diff --git a/src/particle.cpp b/src/particle.cpp index 46df63cd1..779ae18da 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -94,7 +94,7 @@ bool Particle::create_secondary( // Increment number of secondaries created (for ParticleProductionFilter) n_secondaries()++; - auto& bank = secondary_bank().emplace_back(); + SourceSite bank; bank.particle = type; bank.wgt = wgt; bank.r = r(); @@ -102,12 +102,21 @@ bool Particle::create_secondary( bank.E = settings::run_CE ? E : g(); bank.time = time(); bank_second_E() += bank.E; + bank.parent_id = current_work(); + if (settings::use_shared_secondary_bank) { + bank.progeny_id = n_progeny()++; + } + bank.wgt_born = wgt_born(); + bank.wgt_ww_born = wgt_ww_born(); + bank.n_split = n_split(); + + local_secondary_bank().emplace_back(bank); return true; } void Particle::split(double wgt) { - auto& bank = secondary_bank().emplace_back(); + SourceSite bank; bank.particle = type(); bank.wgt = wgt; bank.r = r(); @@ -122,6 +131,16 @@ void Particle::split(double wgt) int surf_id = model::surfaces[surface_index()]->id_; bank.surf_id = (surface() > 0) ? surf_id : -surf_id; } + + bank.wgt_born = wgt_born(); + bank.wgt_ww_born = wgt_ww_born(); + bank.n_split = n_split(); + bank.parent_id = current_work(); + if (settings::use_shared_secondary_bank) { + bank.progeny_id = n_progeny()++; + } + + local_secondary_bank().emplace_back(bank); } void Particle::from_source(const SourceSite* src) @@ -168,6 +187,10 @@ void Particle::from_source(const SourceSite* src) int index_plus_one = model::surface_map[std::abs(src->surf_id)] + 1; surface() = (src->surf_id > 0) ? index_plus_one : -index_plus_one; } + + wgt_born() = src->wgt_born; + wgt_ww_born() = src->wgt_ww_born; + n_split() = src->n_split; } void Particle::event_calculate_xs() @@ -449,61 +472,72 @@ void Particle::event_collide() #endif } -void Particle::event_revive_from_secondary() +void Particle::event_revive_from_secondary(const SourceSite& site) +{ + // Write final position for the previous track (skip if this is a freshly + // constructed particle with no prior track, e.g., Phase 2 of shared + // secondary transport) + if (write_track() && n_event() > 0) { + write_particle_track(*this); + } + + from_source(&site); + + n_event() = 0; + if (!settings::use_shared_secondary_bank) { + n_tracks()++; + } + bank_second_E() = 0.0; + + // Subtract secondary particle energy from interim pulse-height results. + // In shared secondary mode, this subtraction was already done on the parent + // particle during create_secondary(), so skip it here. + if (!settings::use_shared_secondary_bank && + !model::active_pulse_height_tallies.empty() && this->type().is_photon()) { + // Since the birth cell of the particle has not been set we + // have to determine it before the energy of the secondary particle can be + // removed from the pulse-height of this cell. + if (lowest_coord().cell() == C_NONE) { + bool verbose = settings::verbosity >= 10 || trace(); + if (!exhaustive_find_cell(*this, verbose)) { + mark_as_lost("Could not find the cell containing particle " + + std::to_string(id())); + return; + } + // Set birth cell attribute + if (cell_born() == C_NONE) + cell_born() = lowest_coord().cell(); + + // Initialize last cells from current cell + for (int j = 0; j < n_coord(); ++j) { + cell_last(j) = coord(j).cell(); + } + n_coord_last() = n_coord(); + } + pht_secondary_particles(); + } + + // Enter new particle in particle track file + if (write_track()) + add_particle_track(*this); +} + +void Particle::event_check_limit_and_revive() { // If particle has too many events, display warning and kill it - ++n_event(); + n_event()++; if (n_event() == settings::max_particle_events) { warning("Particle " + std::to_string(id()) + " underwent maximum number of events."); wgt() = 0.0; } - // Check for secondary particles if this particle is dead - if (!alive()) { - // Write final position for this particle - if (write_track()) { - write_particle_track(*this); - } - - // If no secondary particles, break out of event loop - if (secondary_bank().empty()) - return; - - from_source(&secondary_bank().back()); - secondary_bank().pop_back(); - n_event() = 0; - bank_second_E() = 0.0; - - // Subtract secondary particle energy from interim pulse-height results - if (!model::active_pulse_height_tallies.empty() && - this->type().is_photon()) { - // Since the birth cell of the particle has not been set we - // have to determine it before the energy of the secondary particle can be - // removed from the pulse-height of this cell. - if (lowest_coord().cell() == C_NONE) { - bool verbose = settings::verbosity >= 10 || trace(); - if (!exhaustive_find_cell(*this, verbose)) { - mark_as_lost("Could not find the cell containing particle " + - std::to_string(id())); - return; - } - // Set birth cell attribute - if (cell_born() == C_NONE) - cell_born() = lowest_coord().cell(); - - // Initialize last cells from current cell - for (int j = 0; j < n_coord(); ++j) { - cell_last(j) = coord(j).cell(); - } - n_coord_last() = n_coord(); - } - pht_secondary_particles(); - } - - // Enter new particle in particle track file - if (write_track()) - add_particle_track(*this); + // In non-shared-secondary mode, revive from local secondary bank + if (!alive() && !settings::use_shared_secondary_bank && + !local_secondary_bank().empty()) { + SourceSite& site = local_secondary_bank().back(); + event_revive_from_secondary(site); + local_secondary_bank().pop_back(); } } @@ -515,6 +549,7 @@ void Particle::event_death() // Finish particle track output. if (write_track()) { + write_particle_track(*this); finalize_particle_track(*this); } @@ -538,11 +573,17 @@ void Particle::event_death() score_pulse_height_tally(*this, model::active_pulse_height_tallies); } + // Accumulate track count for this particle history + if (!settings::use_shared_secondary_bank) { +#pragma omp atomic + simulation::simulation_tracks_completed += n_tracks(); + } + // Record the number of progeny created by this particle. // This data will be used to efficiently sort the fission bank. - if (settings::run_mode == RunMode::EIGENVALUE) { - int64_t offset = id() - 1 - simulation::work_index[mpi::rank]; - simulation::progeny_per_particle[offset] = n_progeny(); + if (settings::run_mode == RunMode::EIGENVALUE || + settings::use_shared_secondary_bank) { + simulation::progeny_per_particle[current_work()] = n_progeny(); } } @@ -863,28 +904,27 @@ void Particle::write_restart() const write_dataset(file_id, "id", id()); write_dataset(file_id, "type", type().pdg_number()); + // Get source site data for the particle that got lost int64_t i = current_work(); + SourceSite site; if (settings::run_mode == RunMode::EIGENVALUE) { - // take source data from primary bank for eigenvalue simulation - write_dataset(file_id, "weight", simulation::source_bank[i - 1].wgt); - write_dataset(file_id, "energy", simulation::source_bank[i - 1].E); - write_dataset(file_id, "xyz", simulation::source_bank[i - 1].r); - write_dataset(file_id, "uvw", simulation::source_bank[i - 1].u); - write_dataset(file_id, "time", simulation::source_bank[i - 1].time); + site = simulation::source_bank[i]; + } else if (settings::run_mode == RunMode::FIXED_SOURCE && + settings::use_shared_secondary_bank && + i < simulation::shared_secondary_bank_read.size()) { + site = simulation::shared_secondary_bank_read[i]; } else if (settings::run_mode == RunMode::FIXED_SOURCE) { - // re-sample using rng random number seed used to generate source particle - int64_t id = (simulation::total_gen + overall_generation() - 1) * - settings::n_particles + - simulation::work_index[mpi::rank] + i; + // Re-sample using the same seed used to generate the source particle. + // current_work() is 0-indexed, compute_particle_id expects 1-indexed. + int64_t id = compute_transport_seed(compute_particle_id(i + 1)); uint64_t seed = init_seed(id, STREAM_SOURCE); - // re-sample source site - auto site = sample_external_source(&seed); - write_dataset(file_id, "weight", site.wgt); - write_dataset(file_id, "energy", site.E); - write_dataset(file_id, "xyz", site.r); - write_dataset(file_id, "uvw", site.u); - write_dataset(file_id, "time", site.time); + site = sample_external_source(&seed); } + write_dataset(file_id, "weight", site.wgt); + write_dataset(file_id, "energy", site.E); + write_dataset(file_id, "xyz", site.r); + write_dataset(file_id, "uvw", site.u); + write_dataset(file_id, "time", site.time); // Close file file_close(file_id); diff --git a/src/particle_restart.cpp b/src/particle_restart.cpp index f02fcb94b..c226d51ec 100644 --- a/src/particle_restart.cpp +++ b/src/particle_restart.cpp @@ -1,6 +1,7 @@ #include "openmc/particle_restart.h" #include "openmc/array.h" +#include "openmc/bank.h" #include "openmc/constants.h" #include "openmc/hdf5_interface.h" #include "openmc/mgxs_interface.h" @@ -106,20 +107,16 @@ void run_particle_restart() // Set all tallies to 0 for now (just tracking errors) model::tallies.clear(); - // Compute random number seed - int64_t particle_seed; - switch (previous_run_mode) { - case RunMode::EIGENVALUE: - case RunMode::FIXED_SOURCE: - particle_seed = (simulation::total_gen + overall_generation() - 1) * - settings::n_particles + - p.id(); - break; - default: - throw std::runtime_error { - "Unexpected run mode: " + - std::to_string(static_cast(previous_run_mode))}; + // Allocate progeny_per_particle if needed for shared secondary mode + // (event_death() writes to this array). Set current_work to 0 since we + // only have one particle being restarted. + if (settings::use_shared_secondary_bank) { + p.current_work() = 0; + simulation::progeny_per_particle.resize(1, 0); } + + // Compute random number seed + int64_t particle_seed = compute_transport_seed(p.id()); init_particle_seeds(particle_seed, p.seeds()); // Force calculation of cross-sections by setting last energy to zero diff --git a/src/physics.cpp b/src/physics.cpp index 106bd1aa2..892b39479 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -44,7 +44,7 @@ void collision(Particle& p) { // Add to collision counter for particle ++(p.n_collision()); - p.secondary_bank_index() = p.secondary_bank().size(); + p.secondary_bank_index() = p.local_secondary_bank().size(); // Sample reaction for the material the particle is in switch (p.type().pdg_number()) { @@ -127,7 +127,8 @@ void sample_neutron_reaction(Particle& p) // Make sure particle population doesn't grow out of control for // subcritical multiplication problems. - if (p.secondary_bank().size() >= settings::max_secondaries) { + if (p.local_secondary_bank().size() >= settings::max_secondaries && + !settings::use_shared_secondary_bank) { fatal_error( "The secondary particle bank appears to be growing without " "bound. You are likely running a subcritical multiplication problem " @@ -228,7 +229,7 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) } // Set parent and progeny IDs - site.parent_id = p.id(); + site.parent_id = p.current_work(); site.progeny_id = p.n_progeny()++; // Store fission site in bank @@ -253,7 +254,10 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) ifp(p, idx); } } else { - p.secondary_bank().push_back(site); + site.wgt_born = p.wgt_born(); + site.wgt_ww_born = p.wgt_ww_born(); + site.n_split = p.n_split(); + p.local_secondary_bank().push_back(site); p.n_secondaries()++; } @@ -1225,7 +1229,7 @@ void sample_secondary_photons(Particle& p, int i_nuclide) // Tag secondary particle with parent nuclide if (created_photon && settings::use_decay_photons) { - p.secondary_bank().back().parent_nuclide = + p.local_secondary_bank().back().parent_nuclide = rx->products_[i_product].parent_nuclide_; } } diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index 6991c7365..212ba765c 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -27,7 +27,7 @@ void collision_mg(Particle& p) { // Add to the collision counter for the particle p.n_collision()++; - p.secondary_bank_index() = p.secondary_bank().size(); + p.secondary_bank_index() = p.local_secondary_bank().size(); // Sample the reaction type sample_reaction(p); @@ -179,7 +179,7 @@ void create_fission_sites(Particle& p) } // Set parent and progeny ID - site.parent_id = p.id(); + site.parent_id = p.current_work(); site.progeny_id = p.n_progeny()++; // Store fission site in bank @@ -200,7 +200,10 @@ void create_fission_sites(Particle& p) break; } } else { - p.secondary_bank().push_back(site); + site.wgt_born = p.wgt_born(); + site.wgt_ww_born = p.wgt_ww_born(); + site.n_split = p.n_split(); + p.local_secondary_bank().push_back(site); p.n_secondaries()++; } diff --git a/src/settings.cpp b/src/settings.cpp index bb6991da9..8ae252ae1 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -83,6 +83,7 @@ bool uniform_source_sampling {false}; bool ufs_on {false}; bool urr_ptables_on {true}; bool use_decay_photons {false}; +bool use_shared_secondary_bank {false}; bool weight_windows_on {false}; bool weight_window_checkpoint_surface {false}; bool weight_window_checkpoint_collision {true}; @@ -1320,6 +1321,27 @@ void read_settings_xml(pugi::xml_node root) settings::use_decay_photons = get_node_value_bool(root, "use_decay_photons"); } + + // If weight windows are on, also enable shared secondary bank (unless + // explicitly disabled by user). + if (check_for_node(root, "shared_secondary_bank")) { + bool val = get_node_value_bool(root, "shared_secondary_bank"); + if (val && run_mode == RunMode::EIGENVALUE) { + warning( + "Shared secondary bank is not supported in eigenvalue calculations. " + "Setting will be ignored."); + } else { + settings::use_shared_secondary_bank = val; + } + } else if (settings::weight_windows_on) { + if (run_mode == RunMode::EIGENVALUE) { + warning( + "Shared secondary bank is not supported in eigenvalue calculations. " + "Particle local secondary banks will be used instead."); + } else if (run_mode == RunMode::FIXED_SOURCE) { + settings::use_shared_secondary_bank = true; + } + } } void free_memory_settings() diff --git a/src/simulation.cpp b/src/simulation.cpp index 4fad196a6..89aa9ca0e 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -86,7 +86,7 @@ int openmc_simulation_init() } // Determine how much work each process should do - calculate_work(); + calculate_work(settings::n_particles); // Allocate source, fission and surface source banks. allocate_banks(); @@ -214,6 +214,20 @@ int openmc_simulation_finalize() // Stop timers and show timing statistics simulation::time_finalize.stop(); simulation::time_total.stop(); + +#ifdef OPENMC_MPI + // Reduce track count across ranks for correct reporting. In shared secondary + // bank mode, all ranks already have the global count; in non-shared mode, + // each rank only has its own count. + if (settings::weight_windows_on && !settings::use_shared_secondary_bank) { + int64_t total_tracks; + MPI_Reduce(&simulation::simulation_tracks_completed, &total_tracks, 1, + MPI_INT64_T, MPI_SUM, 0, mpi::intracomm); + if (mpi::master) + simulation::simulation_tracks_completed = total_tracks; + } +#endif + if (mpi::master) { if (settings::solver_type != SolverType::RANDOM_RAY) { if (settings::verbosity >= 6) @@ -254,9 +268,17 @@ int openmc_next_batch(int* status) // Transport loop if (settings::event_based) { - transport_event_based(); + if (settings::use_shared_secondary_bank) { + transport_event_based_shared_secondary(); + } else { + transport_event_based(); + } } else { - transport_history_based(); + if (settings::use_shared_secondary_bank) { + transport_history_based_shared_secondary(); + } else { + transport_history_based(); + } } // Accumulate time for transport @@ -324,6 +346,8 @@ const RegularMesh* ufs_mesh {nullptr}; vector k_generation; vector work_index; +int64_t simulation_tracks_completed {0}; + } // namespace simulation //============================================================================== @@ -547,7 +571,7 @@ void finalize_generation() // If using shared memory, stable sort the fission bank (by parent IDs) // so as to allow for reproducibility regardless of which order particles // are run in. - sort_fission_bank(); + sort_bank(simulation::fission_bank, true); // Distribute fission bank across processors evenly synchronize_bank(); @@ -571,26 +595,35 @@ void finalize_generation() } } -void initialize_history(Particle& p, int64_t index_source) +void sample_source_particle(Particle& p, int64_t index_source) { - // set defaults + // Sample a particle from the source bank if (settings::run_mode == RunMode::EIGENVALUE) { - // set defaults for eigenvalue simulations from primary bank p.from_source(&simulation::source_bank[index_source - 1]); } else if (settings::run_mode == RunMode::FIXED_SOURCE) { // initialize random number seed - int64_t id = (simulation::total_gen + overall_generation() - 1) * - settings::n_particles + - simulation::work_index[mpi::rank] + index_source; + int64_t id = compute_transport_seed(compute_particle_id(index_source)); uint64_t seed = init_seed(id, STREAM_SOURCE); // sample from external source distribution or custom library then set auto site = sample_external_source(&seed); p.from_source(&site); } - p.current_work() = index_source; +} + +void initialize_particle_track( + Particle& p, int64_t index_source, bool is_secondary) +{ + // Note: index_source is 1-based (first particle = 1), but current_work() is + // stored as 0-based for direct use as an array index into + // progeny_per_particle, source_bank, ifp banks, etc. + if (!is_secondary) { + sample_source_particle(p, index_source); + } + + p.current_work() = index_source - 1; // set identifier for particle - p.id() = simulation::work_index[mpi::rank] + index_source; + p.id() = compute_particle_id(index_source); // set progeny count to zero p.n_progeny() = 0; @@ -598,6 +631,9 @@ void initialize_history(Particle& p, int64_t index_source) // Reset particle event counter p.n_event() = 0; + // Initialize track counter (1 for this primary/secondary track) + p.n_tracks() = 1; + // Reset split counter p.n_split() = 0; @@ -611,9 +647,7 @@ void initialize_history(Particle& p, int64_t index_source) std::fill(p.pht_storage().begin(), p.pht_storage().end(), 0); // set random number seed - int64_t particle_seed = - (simulation::total_gen + overall_generation() - 1) * settings::n_particles + - p.id(); + int64_t particle_seed = compute_transport_seed(p.id()); init_particle_seeds(particle_seed, p.seeds()); // set particle trace @@ -627,17 +661,21 @@ void initialize_history(Particle& p, int64_t index_source) p.write_track() = check_track_criteria(p); // Set the particle's initial weight window value. - p.wgt_ww_born() = -1.0; - apply_weight_windows(p); + if (!is_secondary) { + p.wgt_ww_born() = -1.0; + apply_weight_windows(p); + } // Display message if high verbosity or trace is on if (settings::verbosity >= 9 || p.trace()) { write_message("Simulating Particle {}", p.id()); } -// Add particle's starting weight to count for normalizing tallies later + // Add particle's starting weight to count for normalizing tallies later + if (!is_secondary) { #pragma omp atomic - simulation::total_weight += p.wgt(); + simulation::total_weight += p.wgt(); + } // Force calculation of cross-sections by setting last energy to zero if (settings::run_CE) { @@ -655,13 +693,34 @@ int overall_generation() return settings::gen_per_batch * (current_batch - 1) + current_gen; } -void calculate_work() +int64_t compute_particle_id(int64_t index_source) +{ + if (settings::use_shared_secondary_bank) { + return simulation::work_index[mpi::rank] + index_source + + simulation::simulation_tracks_completed; + } else { + return simulation::work_index[mpi::rank] + index_source; + } +} + +int64_t compute_transport_seed(int64_t particle_id) +{ + if (settings::use_shared_secondary_bank) { + return particle_id; + } else { + return (simulation::total_gen + overall_generation() - 1) * + settings::n_particles + + particle_id; + } +} + +void calculate_work(int64_t n_particles) { // Determine minimum amount of particles to simulate on each processor - int64_t min_work = settings::n_particles / mpi::n_procs; + int64_t min_work = n_particles / mpi::n_procs; // Determine number of processors that have one extra particle - int64_t remainder = settings::n_particles % mpi::n_procs; + int64_t remainder = n_particles % mpi::n_procs; int64_t i_bank = 0; simulation::work_index.resize(mpi::n_procs + 1); @@ -816,7 +875,7 @@ void transport_history_based_single_particle(Particle& p) p.event_collide(); } } - p.event_revive_from_secondary(); + p.event_check_limit_and_revive(); } p.event_death(); } @@ -826,11 +885,137 @@ void transport_history_based() #pragma omp parallel for schedule(runtime) for (int64_t i_work = 1; i_work <= simulation::work_per_rank; ++i_work) { Particle p; - initialize_history(p, i_work); + initialize_particle_track(p, i_work, false); transport_history_based_single_particle(p); } } +// The shared secondary bank transport algorithm works in two phases. In the +// first phase, all primary particles are sampled then transported, and their +// secondary particles are deposited into a shared secondary bank. The second +// phase occurs in a loop, where all secondary tracks in the shared secondary +// bank are transported. Any secondary particles generated during this phase are +// deposited back into the shared secondary bank. The shared secondary bank is +// sorted for consistent ordering and load balanced across MPI ranks. This loop +// continues until there are no more secondary tracks left to transport. +void transport_history_based_shared_secondary() +{ + // Clear shared secondary banks from any prior use + simulation::shared_secondary_bank_read.clear(); + simulation::shared_secondary_bank_write.clear(); + + if (mpi::master) { + write_message(fmt::format(" Primary source particles: {}", + settings::n_particles), + 6); + } + + simulation::progeny_per_particle.resize(simulation::work_per_rank); + std::fill(simulation::progeny_per_particle.begin(), + simulation::progeny_per_particle.end(), 0); + + // Phase 1: Transport primary particles and deposit first generation of + // secondaries in the shared secondary bank +#pragma omp parallel + { + vector thread_bank; + +#pragma omp for schedule(runtime) + for (int64_t i = 1; i <= simulation::work_per_rank; i++) { + Particle p; + initialize_particle_track(p, i, false); + transport_history_based_single_particle(p); + for (auto& site : p.local_secondary_bank()) { + thread_bank.push_back(site); + } + } + + // Drain thread-local bank into the shared secondary bank (once per thread) +#pragma omp critical(SharedSecondaryBank) + { + for (auto& site : thread_bank) { + simulation::shared_secondary_bank_write.thread_unsafe_append(site); + } + } + } + + simulation::simulation_tracks_completed += settings::n_particles; + + // Phase 2: Now that the secondary bank has been populated, enter loop over + // all secondary generations + int n_generation_depth = 1; + int64_t alive_secondary = 1; + while (alive_secondary) { + + // Sort the shared secondary bank by parent ID then progeny ID to + // ensure reproducibility. + sort_bank(simulation::shared_secondary_bank_write, false); + + // Synchronize the shared secondary bank amongst all MPI ranks, such + // that each MPI rank has an approximately equal number of secondary + // tracks. Also reports the total number of secondaries alive across + // all MPI ranks. + alive_secondary = synchronize_global_secondary_bank( + simulation::shared_secondary_bank_write); + + // Recalculate work for each MPI rank based on number of alive secondary + // tracks + calculate_work(alive_secondary); + + // Display the number of secondary tracks in this generation. This + // is useful for user monitoring so as to see if the secondary population is + // exploding and to determine how many generations of secondaries are being + // transported. + if (mpi::master) { + write_message(fmt::format(" Secondary generation {:<2} tracks: {}", + n_generation_depth, alive_secondary), + 6); + } + + simulation::shared_secondary_bank_read = + std::move(simulation::shared_secondary_bank_write); + simulation::shared_secondary_bank_write = SharedArray(); + simulation::progeny_per_particle.resize( + simulation::shared_secondary_bank_read.size()); + std::fill(simulation::progeny_per_particle.begin(), + simulation::progeny_per_particle.end(), 0); + + // Transport all secondary tracks from the shared secondary bank +#pragma omp parallel + { + vector thread_bank; + +#pragma omp for schedule(runtime) + for (int64_t i = 1; i <= simulation::shared_secondary_bank_read.size(); + i++) { + Particle p; + initialize_particle_track(p, i, true); + SourceSite& site = simulation::shared_secondary_bank_read[i - 1]; + p.event_revive_from_secondary(site); + transport_history_based_single_particle(p); + for (auto& secondary_site : p.local_secondary_bank()) { + thread_bank.push_back(secondary_site); + } + } + + // Drain thread-local bank into the shared secondary bank (once per + // thread) +#pragma omp critical(SharedSecondaryBank) + { + for (auto& secondary_site : thread_bank) { + simulation::shared_secondary_bank_write.thread_unsafe_append( + secondary_site); + } + } + } // End of transport loop over tracks in shared secondary bank + n_generation_depth++; + simulation::simulation_tracks_completed += alive_secondary; + } // End of loop over secondary generations + + // Reset work so that fission bank etc works correctly + calculate_work(settings::n_particles); +} + void transport_event_based() { int64_t remaining_work = simulation::work_per_rank; @@ -848,33 +1033,7 @@ void transport_event_based() // Initialize all particle histories for this subiteration process_init_events(n_particles, source_offset); - - // Event-based transport loop - while (true) { - // Determine which event kernel has the longest queue - int64_t max = std::max({simulation::calculate_fuel_xs_queue.size(), - simulation::calculate_nonfuel_xs_queue.size(), - simulation::advance_particle_queue.size(), - simulation::surface_crossing_queue.size(), - simulation::collision_queue.size()}); - - // Execute event with the longest queue - if (max == 0) { - break; - } else if (max == simulation::calculate_fuel_xs_queue.size()) { - process_calculate_xs_events(simulation::calculate_fuel_xs_queue); - } else if (max == simulation::calculate_nonfuel_xs_queue.size()) { - process_calculate_xs_events(simulation::calculate_nonfuel_xs_queue); - } else if (max == simulation::advance_particle_queue.size()) { - process_advance_particle_events(); - } else if (max == simulation::surface_crossing_queue.size()) { - process_surface_crossing_events(); - } else if (max == simulation::collision_queue.size()) { - process_collision_events(); - } - } - - // Execute death event for all particles + process_transport_events(); process_death_events(n_particles); // Adjust remaining work and source offset variables @@ -883,4 +1042,122 @@ void transport_event_based() } } +void transport_event_based_shared_secondary() +{ + // Clear shared secondary banks from any prior use + simulation::shared_secondary_bank_read.clear(); + simulation::shared_secondary_bank_write.clear(); + + if (mpi::master) { + write_message(fmt::format(" Primary source particles: {}", + settings::n_particles), + 6); + } + + simulation::progeny_per_particle.resize(simulation::work_per_rank); + std::fill(simulation::progeny_per_particle.begin(), + simulation::progeny_per_particle.end(), 0); + + // Phase 1: Transport primary particles using event-based processing and + // deposit first generation of secondaries in the shared secondary bank + int64_t remaining_work = simulation::work_per_rank; + int64_t source_offset = 0; + + while (remaining_work > 0) { + int64_t n_particles = + std::min(remaining_work, settings::max_particles_in_flight); + + process_init_events(n_particles, source_offset); + process_transport_events(); + process_death_events(n_particles); + + // Collect secondaries from all particle buffers into shared bank + for (int64_t i = 0; i < n_particles; i++) { + for (auto& site : simulation::particles[i].local_secondary_bank()) { + simulation::shared_secondary_bank_write.thread_unsafe_append(site); + } + simulation::particles[i].local_secondary_bank().clear(); + } + + remaining_work -= n_particles; + source_offset += n_particles; + } + + simulation::simulation_tracks_completed += settings::n_particles; + + // Phase 2: Now that the secondary bank has been populated, enter loop over + // all secondary generations + int n_generation_depth = 1; + int64_t alive_secondary = 1; + while (alive_secondary) { + + // Sort the shared secondary bank by parent ID then progeny ID to + // ensure reproducibility. + sort_bank(simulation::shared_secondary_bank_write, false); + + // Synchronize the shared secondary bank amongst all MPI ranks, such + // that each MPI rank has an approximately equal number of secondary + // tracks. + alive_secondary = synchronize_global_secondary_bank( + simulation::shared_secondary_bank_write); + + // Recalculate work for each MPI rank based on number of alive secondary + // tracks + calculate_work(alive_secondary); + + if (mpi::master) { + write_message(fmt::format(" Secondary generation {:<2} tracks: {}", + n_generation_depth, alive_secondary), + 6); + } + + simulation::shared_secondary_bank_read = + std::move(simulation::shared_secondary_bank_write); + simulation::shared_secondary_bank_write = SharedArray(); + simulation::progeny_per_particle.resize( + simulation::shared_secondary_bank_read.size()); + std::fill(simulation::progeny_per_particle.begin(), + simulation::progeny_per_particle.end(), 0); + + // Ensure particle buffer is large enough for this secondary generation + int64_t sec_buffer_length = std::min( + static_cast(simulation::shared_secondary_bank_read.size()), + settings::max_particles_in_flight); + if (sec_buffer_length > + static_cast(simulation::particles.size())) { + init_event_queues(sec_buffer_length); + } + + // Transport secondary tracks using event-based processing + int64_t sec_remaining = simulation::shared_secondary_bank_read.size(); + int64_t sec_offset = 0; + + while (sec_remaining > 0) { + int64_t n_particles = + std::min(sec_remaining, settings::max_particles_in_flight); + + process_init_secondary_events( + n_particles, sec_offset, simulation::shared_secondary_bank_read); + process_transport_events(); + process_death_events(n_particles); + + // Collect secondaries from all particle buffers into shared bank + for (int64_t i = 0; i < n_particles; i++) { + for (auto& site : simulation::particles[i].local_secondary_bank()) { + simulation::shared_secondary_bank_write.thread_unsafe_append(site); + } + simulation::particles[i].local_secondary_bank().clear(); + } + + sec_remaining -= n_particles; + sec_offset += n_particles; + } // End of subiteration loop over secondary tracks + n_generation_depth++; + simulation::simulation_tracks_completed += alive_secondary; + } // End of loop over secondary generations + + // Reset work so that fission bank etc works correctly + calculate_work(settings::n_particles); +} + } // namespace openmc diff --git a/src/tallies/filter_particle_production.cpp b/src/tallies/filter_particle_production.cpp index f8d82fdd9..899809679 100644 --- a/src/tallies/filter_particle_production.cpp +++ b/src/tallies/filter_particle_production.cpp @@ -19,7 +19,7 @@ void ParticleProductionFilter::get_all_bins( // Loop over secondary bank entries for (int bank_idx = start_idx; bank_idx < end_idx; bank_idx++) { - const auto& site = p.secondary_bank(bank_idx); + const auto& site = p.local_secondary_bank(bank_idx); // Find which particle-type slot this secondary belongs to auto it = type_to_index_.find(site.particle.pdg_number()); diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 5dcd07331..b7947afab 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -945,7 +945,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, if (p.type().is_neutron() && p.fission()) { if (is_generation_time_or_both()) { const auto& lifetimes = - simulation::ifp_source_lifetime_bank[p.current_work() - 1]; + simulation::ifp_source_lifetime_bank[p.current_work()]; if (lifetimes.size() == settings::ifp_n_generation) { score = lifetimes[0] * p.wgt_last(); } @@ -959,7 +959,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, if (p.type().is_neutron() && p.fission()) { if (is_beta_effective_or_both()) { const auto& delayed_groups = - simulation::ifp_source_delayed_group_bank[p.current_work() - 1]; + simulation::ifp_source_delayed_group_bank[p.current_work()]; if (delayed_groups.size() == settings::ifp_n_generation) { if (delayed_groups[0] > 0) { score = p.wgt_last(); @@ -985,12 +985,11 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, int ifp_data_size; if (is_beta_effective_or_both()) { ifp_data_size = static_cast( - simulation::ifp_source_delayed_group_bank[p.current_work() - 1] + simulation::ifp_source_delayed_group_bank[p.current_work()] .size()); } else { ifp_data_size = static_cast( - simulation::ifp_source_lifetime_bank[p.current_work() - 1] - .size()); + simulation::ifp_source_lifetime_bank[p.current_work()].size()); } if (ifp_data_size == settings::ifp_n_generation) { score = p.wgt_last(); diff --git a/tests/regression_tests/mg_fixed_source_ww_fission_shared_secondary/__init__.py b/tests/regression_tests/mg_fixed_source_ww_fission_shared_secondary/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/regression_tests/mg_fixed_source_ww_fission_shared_secondary/inputs_true.dat b/tests/regression_tests/mg_fixed_source_ww_fission_shared_secondary/inputs_true.dat new file mode 100644 index 000000000..318fdc7b9 --- /dev/null +++ b/tests/regression_tests/mg_fixed_source_ww_fission_shared_secondary/inputs_true.dat @@ -0,0 +1,65 @@ + + + + 2g.h5 + + + + + + + + + + + + fixed source + 100 + 2 + 0 + + + 0.0 -1000.0 -1000.0 929.45 1000.0 1000.0 + + + + false + + multi-group + + false + + true + + 1 + neutron + 0.0 0.625 20000000.0 + 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 + 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 + 3.0 + 10 + 1e-38 + + + 5 1 1 + 0.0 -1000.0 -1000.0 + 929.45 1000.0 1000.0 + + true + 100 + + + + 5 1 1 + 0.0 -1000.0 -1000.0 + 929.45 1000.0 1000.0 + + + 2 + + + 1 + flux + + + diff --git a/tests/regression_tests/mg_fixed_source_ww_fission_shared_secondary/results_true.dat b/tests/regression_tests/mg_fixed_source_ww_fission_shared_secondary/results_true.dat new file mode 100644 index 000000000..38ff94417 --- /dev/null +++ b/tests/regression_tests/mg_fixed_source_ww_fission_shared_secondary/results_true.dat @@ -0,0 +1,11 @@ +tally 1: +1.765369E+04 +3.087787E+08 +1.708316E+04 +2.842002E+08 +9.444106E+03 +7.488341E+07 +2.066528E+03 +2.142445E+06 +8.689619E+02 +5.652099E+05 diff --git a/tests/regression_tests/mg_fixed_source_ww_fission_shared_secondary/test.py b/tests/regression_tests/mg_fixed_source_ww_fission_shared_secondary/test.py new file mode 100644 index 000000000..336c8461f --- /dev/null +++ b/tests/regression_tests/mg_fixed_source_ww_fission_shared_secondary/test.py @@ -0,0 +1,91 @@ +import os + +import numpy as np +import openmc +from openmc.examples import slab_mg + +from tests.testing_harness import PyAPITestHarness + + +def create_library(): + # Instantiate the energy group data and file object + groups = openmc.mgxs.EnergyGroups([0.0, 0.625, 20.0e6]) + + mg_cross_sections_file = openmc.MGXSLibrary(groups) + + # Make the base, isotropic data + nu = [2.50, 2.50] + fiss = np.array([0.002817, 0.097]) + capture = [0.008708, 0.02518] + absorption = np.add(capture, fiss) + scatter = np.array( + [[[0.31980, 0.06694], [0.004555, -0.0003972]], + [[0.00000, 0.00000], [0.424100, 0.05439000]]]) + total = [0.33588, 0.54628] + chi = [1., 0.] + + mat_1 = openmc.XSdata('mat_1', groups) + mat_1.order = 1 + mat_1.set_nu_fission(np.multiply(nu, fiss)) + mat_1.set_absorption(absorption) + mat_1.set_scatter_matrix(scatter) + mat_1.set_total(total) + mat_1.set_chi(chi) + mg_cross_sections_file.add_xsdata(mat_1) + + # Write the file + mg_cross_sections_file.export_to_hdf5('2g.h5') + + +class MGXSTestHarness(PyAPITestHarness): + def _cleanup(self): + super()._cleanup() + f = '2g.h5' + if os.path.exists(f): + os.remove(f) + + +def test_mg_fixed_source_ww_fission_shared_secondary(): + create_library() + model = slab_mg() + + # Override settings for fixed-source mode with shared secondary bank + model.settings.run_mode = 'fixed source' + model.settings.inactive = 0 + model.settings.batches = 2 + model.settings.particles = 100 + model.settings.create_fission_neutrons = True + model.settings.shared_secondary_bank = True + model.settings.max_history_splits = 100 + + # Add weight windows on a simple 1D mesh + ww_mesh = openmc.RegularMesh() + ww_mesh.lower_left = (0.0, -1000.0, -1000.0) + ww_mesh.upper_right = (929.45, 1000.0, 1000.0) + ww_mesh.dimension = (5, 1, 1) + + # Uniform lower bounds for 2 energy groups, 5 spatial bins + lower_bounds = np.full((2, 5, 1, 1), 0.5) + ww = openmc.WeightWindows( + ww_mesh, + lower_bounds.flatten(), + None, + 5.0, + [0.0, 0.625, 20.0e6], + 'neutron' + ) + model.settings.weight_windows = [ww] + + # Add a flux tally + mesh = openmc.RegularMesh() + mesh.lower_left = (0.0, -1000.0, -1000.0) + mesh.upper_right = (929.45, 1000.0, 1000.0) + mesh.dimension = (5, 1, 1) + + tally = openmc.Tally() + tally.filters = [openmc.MeshFilter(mesh)] + tally.scores = ['flux'] + model.tallies = [tally] + + harness = MGXSTestHarness('statepoint.2.h5', model) + harness.main() diff --git a/tests/regression_tests/particle_production_fission/__init__.py b/tests/regression_tests/particle_production_fission/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/regression_tests/particle_production_fission/local/inputs_true.dat b/tests/regression_tests/particle_production_fission/local/inputs_true.dat new file mode 100644 index 000000000..e93dfc36e --- /dev/null +++ b/tests/regression_tests/particle_production_fission/local/inputs_true.dat @@ -0,0 +1,38 @@ + + + + + + + + + + + + + + fixed source + 100 + 2 + + + 1000000.0 1.0 + + + true + false + + + + neutron + + + neutron + + + 2 1 + events + analog + + + diff --git a/tests/regression_tests/particle_production_fission/local/results_true.dat b/tests/regression_tests/particle_production_fission/local/results_true.dat new file mode 100644 index 000000000..2f5288eba --- /dev/null +++ b/tests/regression_tests/particle_production_fission/local/results_true.dat @@ -0,0 +1,3 @@ +tally 1: +4.570000E+00 +1.047890E+01 diff --git a/tests/regression_tests/particle_production_fission/shared/inputs_true.dat b/tests/regression_tests/particle_production_fission/shared/inputs_true.dat new file mode 100644 index 000000000..b50a39cd0 --- /dev/null +++ b/tests/regression_tests/particle_production_fission/shared/inputs_true.dat @@ -0,0 +1,38 @@ + + + + + + + + + + + + + + fixed source + 100 + 2 + + + 1000000.0 1.0 + + + true + true + + + + neutron + + + neutron + + + 2 1 + events + analog + + + diff --git a/tests/regression_tests/particle_production_fission/shared/results_true.dat b/tests/regression_tests/particle_production_fission/shared/results_true.dat new file mode 100644 index 000000000..5a4dfc516 --- /dev/null +++ b/tests/regression_tests/particle_production_fission/shared/results_true.dat @@ -0,0 +1,3 @@ +tally 1: +5.180000E+00 +1.355140E+01 diff --git a/tests/regression_tests/particle_production_fission/test.py b/tests/regression_tests/particle_production_fission/test.py new file mode 100644 index 000000000..2ad6ca515 --- /dev/null +++ b/tests/regression_tests/particle_production_fission/test.py @@ -0,0 +1,49 @@ +import openmc +import pytest +from openmc.utility_funcs import change_directory + +from tests.testing_harness import PyAPITestHarness + + +@pytest.mark.parametrize("shared_secondary,subdir", [ + (False, "local"), + (True, "shared"), +]) +def test_particle_production_fission(shared_secondary, subdir): + """Fixed-source model with fissionable material to test that + ParticleProductionFilter correctly counts fission-born neutrons, + with both local and shared secondary bank modes.""" + with change_directory(subdir): + openmc.reset_auto_ids() + model = openmc.Model() + + mat = openmc.Material() + mat.set_density('g/cm3', 18.0) + mat.add_nuclide('U235', 1.0) + model.materials.append(mat) + + sph = openmc.Sphere(r=5.0, boundary_type='vacuum') + cell = openmc.Cell(fill=mat, region=-sph) + model.geometry = openmc.Geometry([cell]) + + source = openmc.IndependentSource() + source.energy = openmc.stats.delta_function(1.0e6) + + model.settings.particles = 100 + model.settings.run_mode = 'fixed source' + model.settings.batches = 2 + model.settings.source = source + model.settings.create_fission_neutrons = True + model.settings.shared_secondary_bank = shared_secondary + + # ParticleProductionFilter tracking fission neutron production + ppf = openmc.ParticleProductionFilter(['neutron']) + neutron_filter = openmc.ParticleFilter(['neutron']) + tally = openmc.Tally() + tally.filters = [neutron_filter, ppf] + tally.scores = ['events'] + tally.estimator = 'analog' + model.tallies = [tally] + + harness = PyAPITestHarness('statepoint.2.h5', model) + harness.main() diff --git a/tests/regression_tests/particle_restart_fixed_shared_secondary/__init__.py b/tests/regression_tests/particle_restart_fixed_shared_secondary/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/regression_tests/particle_restart_fixed_shared_secondary/geometry.xml b/tests/regression_tests/particle_restart_fixed_shared_secondary/geometry.xml new file mode 100644 index 000000000..c86e016c6 --- /dev/null +++ b/tests/regression_tests/particle_restart_fixed_shared_secondary/geometry.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/regression_tests/particle_restart_fixed_shared_secondary/materials.xml b/tests/regression_tests/particle_restart_fixed_shared_secondary/materials.xml new file mode 100644 index 000000000..f3851d7ef --- /dev/null +++ b/tests/regression_tests/particle_restart_fixed_shared_secondary/materials.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/regression_tests/particle_restart_fixed_shared_secondary/results_true.dat b/tests/regression_tests/particle_restart_fixed_shared_secondary/results_true.dat new file mode 100644 index 000000000..0c84541f9 --- /dev/null +++ b/tests/regression_tests/particle_restart_fixed_shared_secondary/results_true.dat @@ -0,0 +1,16 @@ +current batch: +4.000000E+00 +current generation: +1.000000E+00 +particle id: +3.241000E+03 +run mode: +fixed source +particle weight: +1.000000E+00 +particle energy: +3.896365E+06 +particle xyz: +8.710681E-01 3.698823E+00 -2.286229E+00 +particle uvw: +-5.882735E-01 4.665422E-01 -6.605093E-01 diff --git a/tests/regression_tests/particle_restart_fixed_shared_secondary/settings.xml b/tests/regression_tests/particle_restart_fixed_shared_secondary/settings.xml new file mode 100644 index 000000000..3a80fb17b --- /dev/null +++ b/tests/regression_tests/particle_restart_fixed_shared_secondary/settings.xml @@ -0,0 +1,15 @@ + + + + fixed source + 12 + 1000 + true + + + + -10 -10 -5 10 10 5 + + + + diff --git a/tests/regression_tests/particle_restart_fixed_shared_secondary/test.py b/tests/regression_tests/particle_restart_fixed_shared_secondary/test.py new file mode 100644 index 000000000..770010900 --- /dev/null +++ b/tests/regression_tests/particle_restart_fixed_shared_secondary/test.py @@ -0,0 +1,6 @@ +from tests.testing_harness import ParticleRestartTestHarness + + +def test_particle_restart_fixed_shared_secondary(): + harness = ParticleRestartTestHarness('particle_4_3241.h5') + harness.main() diff --git a/tests/regression_tests/pulse_height/inputs_true.dat b/tests/regression_tests/pulse_height/local/inputs_true.dat similarity index 95% rename from tests/regression_tests/pulse_height/inputs_true.dat rename to tests/regression_tests/pulse_height/local/inputs_true.dat index 590928e43..50e026bfb 100644 --- a/tests/regression_tests/pulse_height/inputs_true.dat +++ b/tests/regression_tests/pulse_height/local/inputs_true.dat @@ -18,14 +18,12 @@ 100 5 - - 0.0 0.0 0.0 - 1000000.0 1.0 true + false diff --git a/tests/regression_tests/pulse_height/results_true.dat b/tests/regression_tests/pulse_height/local/results_true.dat similarity index 100% rename from tests/regression_tests/pulse_height/results_true.dat rename to tests/regression_tests/pulse_height/local/results_true.dat diff --git a/tests/regression_tests/pulse_height/shared/inputs_true.dat b/tests/regression_tests/pulse_height/shared/inputs_true.dat new file mode 100644 index 000000000..c2aef829c --- /dev/null +++ b/tests/regression_tests/pulse_height/shared/inputs_true.dat @@ -0,0 +1,40 @@ + + + + + + + + + + + + + + + + + fixed source + 100 + 5 + + + 1000000.0 1.0 + + + true + true + + + + 1 + + + 0.0 10000.0 20000.0 30000.0 40000.0 50000.0 60000.0 70000.0 80000.0 90000.0 100000.0 110000.0 120000.0 130000.0 140000.0 150000.0 160000.0 170000.0 180000.0 190000.0 200000.0 210000.0 220000.0 230000.0 240000.0 250000.0 260000.0 270000.0 280000.0 290000.0 300000.0 310000.0 320000.0 330000.0 340000.0 350000.0 360000.0 370000.0 380000.0 390000.0 400000.0 410000.0 420000.0 430000.0 440000.0 450000.0 460000.0 470000.0 480000.0 490000.0 500000.0 510000.0 520000.0 530000.0 540000.0 550000.0 560000.0 570000.0 580000.0 590000.0 600000.0 610000.0 620000.0 630000.0 640000.0 650000.0 660000.0 670000.0 680000.0 690000.0 700000.0 710000.0 720000.0 730000.0 740000.0 750000.0 760000.0 770000.0 780000.0 790000.0 800000.0 810000.0 820000.0 830000.0 840000.0 850000.0 860000.0 870000.0 880000.0 890000.0 900000.0 910000.0 920000.0 930000.0 940000.0 950000.0 960000.0 970000.0 980000.0 990000.0 1000000.0 + + + 1 2 + pulse-height + + + diff --git a/tests/regression_tests/pulse_height/shared/results_true.dat b/tests/regression_tests/pulse_height/shared/results_true.dat new file mode 100644 index 000000000..c57e8ff1c --- /dev/null +++ b/tests/regression_tests/pulse_height/shared/results_true.dat @@ -0,0 +1,201 @@ +tally 1: +4.140000E+00 +3.443000E+00 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +3.000000E-02 +5.000000E-04 +2.000000E-02 +4.000000E-04 +2.000000E-02 +4.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +3.000000E-02 +3.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +4.000000E-02 +6.000000E-04 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +2.000000E-02 +2.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +2.000000E-02 +4.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +2.000000E-02 +2.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +3.000000E-02 +5.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +2.000000E-02 +2.000000E-04 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +2.000000E-02 +2.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +2.000000E-02 +2.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +3.000000E-02 +5.000000E-04 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +2.000000E-02 +2.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +2.000000E-01 +8.600000E-03 diff --git a/tests/regression_tests/pulse_height/test.py b/tests/regression_tests/pulse_height/test.py index 90d960f66..2cbb9152e 100644 --- a/tests/regression_tests/pulse_height/test.py +++ b/tests/regression_tests/pulse_height/test.py @@ -1,53 +1,53 @@ import numpy as np import openmc import pytest +from openmc.utility_funcs import change_directory from tests.testing_harness import PyAPITestHarness -@pytest.fixture -def sphere_model(): - - model = openmc.model.Model() +@pytest.mark.parametrize("shared_secondary,subdir", [ + (False, "local"), + (True, "shared"), +]) +def test_pulse_height(shared_secondary, subdir): + with change_directory(subdir): + openmc.reset_auto_ids() + model = openmc.Model() - # Define materials - NaI = openmc.Material() - NaI.set_density('g/cc', 3.7) - NaI.add_element('Na', 1.0) - NaI.add_element('I', 1.0) + # Define materials + NaI = openmc.Material() + NaI.set_density('g/cc', 3.7) + NaI.add_element('Na', 1.0) + NaI.add_element('I', 1.0) - model.materials = openmc.Materials([NaI]) + model.materials = openmc.Materials([NaI]) - # Define geometry: two spheres in each other - s1 = openmc.Sphere(r=1) - s2 = openmc.Sphere(r=2, boundary_type='vacuum') - inner_sphere = openmc.Cell(name='inner sphere', fill=NaI, region=-s1) - outer_sphere = openmc.Cell(name='outer sphere', region=+s1 & -s2) - model.geometry = openmc.Geometry([inner_sphere, outer_sphere]) + # Define geometry: two spheres in each other + s1 = openmc.Sphere(r=1) + s2 = openmc.Sphere(r=2, boundary_type='vacuum') + inner_sphere = openmc.Cell(name='inner sphere', fill=NaI, region=-s1) + outer_sphere = openmc.Cell(name='outer sphere', region=+s1 & -s2) + model.geometry = openmc.Geometry([inner_sphere, outer_sphere]) - # Define settings - model.settings.run_mode = 'fixed source' - model.settings.batches = 5 - model.settings.particles = 100 - model.settings.photon_transport = True - model.settings.source = openmc.IndependentSource( - space=openmc.stats.Point(), - energy=openmc.stats.Discrete([1e6], [1]), - particle='photon' - ) + # Define settings + model.settings.run_mode = 'fixed source' + model.settings.batches = 5 + model.settings.particles = 100 + model.settings.photon_transport = True + model.settings.shared_secondary_bank = shared_secondary + model.settings.source = openmc.IndependentSource( + energy=openmc.stats.delta_function(1e6), + particle='photon' + ) - # Define tallies - tally = openmc.Tally(name="pht tally") - tally.scores = ['pulse-height'] - cell_filter = openmc.CellFilter(inner_sphere) - energy_filter = openmc.EnergyFilter(np.linspace(0, 1_000_000, 101)) - tally.filters = [cell_filter, energy_filter] - model.tallies = [tally] + # Define tallies + tally = openmc.Tally(name="pht tally") + tally.scores = ['pulse-height'] + cell_filter = openmc.CellFilter(inner_sphere) + energy_filter = openmc.EnergyFilter(np.linspace(0, 1_000_000, 101)) + tally.filters = [cell_filter, energy_filter] + model.tallies = [tally] - return model - - - -def test_pulse_height(sphere_model): - harness = PyAPITestHarness('statepoint.5.h5', sphere_model) - harness.main() + harness = PyAPITestHarness('statepoint.5.h5', model) + harness.main() diff --git a/tests/regression_tests/weightwindows/local/inputs_true.dat b/tests/regression_tests/weightwindows/local/inputs_true.dat new file mode 100644 index 000000000..57dc42ca4 --- /dev/null +++ b/tests/regression_tests/weightwindows/local/inputs_true.dat @@ -0,0 +1,94 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + fixed source + 500 + 2 + + + 0.001 0.001 0.001 + + + 14000000.0 1.0 + + + true + + 2 + neutron + 0.0 0.5 20000000.0 + -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.0008758251780046591 -1.0 -1.0 -1.0 -1.0 0.00044494017319042853 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 1.6620271125067025e-05 0.0006278777700923334 3.3816651814344154e-05 -1.0 -1.0 0.0024363004681119066 0.04404124277352227 0.002389900091734746 -1.0 -1.0 0.001096572213298872 0.04862206884590391 0.0011530054432113332 -1.0 -1.0 -1.0 0.00029204950777272335 2.2332845991385424e-05 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.0001298973206651331 0.0028826281529980655 0.00018941326535850932 -1.0 5.4657521927731274e-05 0.014670839729146354 0.49999999999999994 0.011271060592765664 -1.0 -1.0 0.014846491082523524 0.4897211700090108 0.013284874451810723 -1.0 -1.0 5.14657989105535e-05 0.0010115278438204965 0.0008429411802845685 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.0009344129479768359 -1.0 -1.0 -1.0 0.0024828273390431962 0.04299740304329489 0.0020539079252555113 -1.0 -1.0 0.003034165905819096 0.04870927636605937 0.00313101668086297 -1.0 -1.0 -1.0 6.339910999436737e-05 7.442176086066386e-05 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.00011276372995589871 0.0002236100157024887 9.56304298913265e-08 -1.0 -1.0 0.0001596955110781239 9.960576598335084e-06 5.3833030623153676e-05 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 4.0453018186074215e-05 4.6013290110121894e-05 2.709738801641897e-05 -1.0 -1.0 -1.0 9.538410811938703e-05 1.992978141244964e-05 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 1.0723231851298364e-05 -1.0 -1.0 -1.0 -1.0 0.0008030100296698905 0.017386220476187222 0.0005027758935317528 -1.0 -1.0 0.00103568411326969 0.015609071124009234 0.0007473731339616588 -1.0 -1.0 -1.0 6.979537549963374e-05 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.0014257756927429828 5.9098070539559466e-05 -1.0 5.0168071459366625e-06 0.005588863190166347 0.5 0.003930588100414563 -1.0 -1.0 0.005163345217476071 0.48250750993887326 0.003510432129522241 -1.0 -1.0 3.371977841720992e-05 0.0006640103276501794 9.554899988419713e-05 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.000187872910697387 -1.0 -1.0 -1.0 0.0009134587866163143 0.017081858078684467 0.0002972747357542354 -1.0 -1.0 0.0007973472495507653 0.019300903254809747 0.000902203032280694 -1.0 -1.0 6.170015233787292e-05 1.898984300674969e-05 9.85411583595722e-07 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 1.0173913765490095e-05 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 + -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 0.008758251780046591 -10.0 -10.0 -10.0 -10.0 0.004449401731904285 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 0.00016620271125067024 0.006278777700923334 0.00033816651814344155 -10.0 -10.0 0.024363004681119065 0.4404124277352227 0.02389900091734746 -10.0 -10.0 0.01096572213298872 0.4862206884590391 0.011530054432113333 -10.0 -10.0 -10.0 0.0029204950777272335 0.00022332845991385425 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 0.0012989732066513312 0.028826281529980655 0.0018941326535850931 -10.0 0.0005465752192773128 0.14670839729146354 4.999999999999999 0.11271060592765664 -10.0 -10.0 0.14846491082523525 4.897211700090108 0.13284874451810724 -10.0 -10.0 0.000514657989105535 0.010115278438204964 0.008429411802845685 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 0.00934412947976836 -10.0 -10.0 -10.0 0.024828273390431962 0.4299740304329489 0.020539079252555114 -10.0 -10.0 0.03034165905819096 0.48709276366059373 0.0313101668086297 -10.0 -10.0 -10.0 0.0006339910999436737 0.0007442176086066385 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 0.0011276372995589871 0.002236100157024887 9.56304298913265e-07 -10.0 -10.0 0.001596955110781239 9.960576598335084e-05 0.0005383303062315367 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 0.00040453018186074213 0.00046013290110121896 0.0002709738801641897 -10.0 -10.0 -10.0 0.0009538410811938703 0.00019929781412449641 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 0.00010723231851298364 -10.0 -10.0 -10.0 -10.0 0.008030100296698905 0.17386220476187222 0.0050277589353175285 -10.0 -10.0 0.010356841132696899 0.15609071124009233 0.007473731339616587 -10.0 -10.0 -10.0 0.0006979537549963374 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 0.014257756927429827 0.0005909807053955946 -10.0 5.016807145936663e-05 0.055888631901663474 5.0 0.039305881004145636 -10.0 -10.0 0.05163345217476071 4.825075099388733 0.03510432129522241 -10.0 -10.0 0.00033719778417209923 0.006640103276501793 0.0009554899988419713 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 0.00187872910697387 -10.0 -10.0 -10.0 0.009134587866163143 0.17081858078684467 0.002972747357542354 -10.0 -10.0 0.007973472495507653 0.19300903254809748 0.009022030322806941 -10.0 -10.0 0.0006170015233787292 0.0001898984300674969 9.85411583595722e-06 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 0.00010173913765490096 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 + 3.0 + 1.5 + 10 + 1e-38 + + + 5 6 7 + -240 -240 -240 + 240 240 240 + + + 2 + photon + 0.0 0.5 20000000.0 + -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 2.0445691096074744e-05 -1.0 -1.0 -1.0 0.00013281328509288383 0.0006267128082339883 -1.0 -1.0 -1.0 -1.0 0.0004209808495056742 5.340221214300681e-05 -1.0 -1.0 -1.0 1.553693492042344e-05 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 3.7096847855980484e-05 -1.0 -1.0 5.41969614042385e-05 0.0007918090215756442 9.763572147337743e-05 -1.0 -1.0 0.0026398319229115233 0.04039871468258457 0.002806654735003014 -1.0 -1.0 0.0027608366482897244 0.040190978087723005 0.0025084416103979975 -1.0 -1.0 9.774868433298412e-05 0.0006126404916879651 0.00014841730774317252 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 9.75342408717107e-06 1.7610849405640072e-05 -1.0 -1.0 0.0001160370359598608 0.002822592083222468 0.0003897177301239376 -1.0 -1.0 0.014965179733172532 0.5 0.011723491704290401 -1.0 2.260493623875525e-05 0.015157364795884922 0.49376993953402387 0.013111419473204967 -1.0 -1.0 0.00014919915366377564 0.002197964829133137 0.0003209711829219673 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 1.4806422867046436e-05 -1.0 -1.0 -1.0 9.207268175745281e-05 0.0007066033561638983 -1.0 -1.0 -1.0 0.003246064903858183 0.03905493028065908 0.0025380574501073605 -1.0 -1.0 0.0032499260211917933 0.04335567738779479 0.0031577214557017056 4.521352809838806e-05 -1.0 0.00018268090756001903 0.0004353654810741932 0.00011273259573092372 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 9.650194311662424e-06 0.00029403373970835075 0.00013506203419326438 -1.0 -1.0 4.5065302725431816e-06 0.00027725323638744237 3.547290864255937e-05 -1.0 -1.0 -1.0 8.786172732576295e-05 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 + -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 0.00020445691096074745 -10.0 -10.0 -10.0 0.0013281328509288383 0.006267128082339883 -10.0 -10.0 -10.0 -10.0 0.004209808495056742 0.0005340221214300681 -10.0 -10.0 -10.0 0.00015536934920423439 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 0.00037096847855980485 -10.0 -10.0 0.000541969614042385 0.007918090215756443 0.0009763572147337743 -10.0 -10.0 0.026398319229115234 0.4039871468258457 0.02806654735003014 -10.0 -10.0 0.027608366482897245 0.40190978087723006 0.025084416103979976 -10.0 -10.0 0.0009774868433298411 0.0061264049168796506 0.0014841730774317252 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 9.753424087171071e-05 0.00017610849405640073 -10.0 -10.0 0.001160370359598608 0.02822592083222468 0.003897177301239376 -10.0 -10.0 0.14965179733172532 5.0 0.11723491704290401 -10.0 0.0002260493623875525 0.15157364795884923 4.937699395340239 0.13111419473204966 -10.0 -10.0 0.0014919915366377564 0.02197964829133137 0.003209711829219673 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 0.00014806422867046437 -10.0 -10.0 -10.0 0.0009207268175745281 0.007066033561638983 -10.0 -10.0 -10.0 0.03246064903858183 0.39054930280659084 0.025380574501073606 -10.0 -10.0 0.03249926021191793 0.4335567738779479 0.03157721455701706 0.0004521352809838806 -10.0 0.0018268090756001904 0.004353654810741932 0.001127325957309237 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 9.650194311662424e-05 0.0029403373970835075 0.0013506203419326437 -10.0 -10.0 4.5065302725431816e-05 0.002772532363874424 0.0003547290864255937 -10.0 -10.0 -10.0 0.0008786172732576295 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 -10.0 + 3.0 + 1.5 + 10 + 1e-38 + + false + + true + true + + 200 + + + + 5 10 15 + -240 -240 -240 + 240 240 240 + + + 1 + + + 0.0 0.5 20000000.0 + + + neutron photon + + + 1 2 3 + flux + + + diff --git a/tests/regression_tests/weightwindows/local/results_true.dat b/tests/regression_tests/weightwindows/local/results_true.dat new file mode 100644 index 000000000..d3dfa119d --- /dev/null +++ b/tests/regression_tests/weightwindows/local/results_true.dat @@ -0,0 +1 @@ +a78972fe3c0dadfc256cfb139c5ca8c7b634738e1895ad2d6286ed5e2567c6dc518e499363908e9058432684148a706a179a39110f703d5322491184a1d0c3e4 \ No newline at end of file diff --git a/tests/regression_tests/weightwindows/results_true.dat b/tests/regression_tests/weightwindows/results_true.dat deleted file mode 100644 index 897089e05..000000000 --- a/tests/regression_tests/weightwindows/results_true.dat +++ /dev/null @@ -1 +0,0 @@ -a4a3ccca43666e2ca1e71800201b152cca20c387b93d67522c5339807348dcee5cada9acbed3238f37e2e86e76b374b06988742f07d4ea1b413e4e75d0c180b1 \ No newline at end of file diff --git a/tests/regression_tests/weightwindows/inputs_true.dat b/tests/regression_tests/weightwindows/shared/inputs_true.dat similarity index 99% rename from tests/regression_tests/weightwindows/inputs_true.dat rename to tests/regression_tests/weightwindows/shared/inputs_true.dat index aa7a0379d..6ec7d12e3 100644 --- a/tests/regression_tests/weightwindows/inputs_true.dat +++ b/tests/regression_tests/weightwindows/shared/inputs_true.dat @@ -64,6 +64,7 @@ 10 1e-38 + true true true diff --git a/tests/regression_tests/weightwindows/shared/results_true.dat b/tests/regression_tests/weightwindows/shared/results_true.dat new file mode 100644 index 000000000..01381db32 --- /dev/null +++ b/tests/regression_tests/weightwindows/shared/results_true.dat @@ -0,0 +1 @@ +27c2d218ecf46161088003fbb39ccd63495f80d6e0b27dbc3b3e845b49b282dabf0eaba4fb619be15c6116fc963ec5ce1fe4197efbe4084972ec761e02e70eb5 \ No newline at end of file diff --git a/tests/regression_tests/weightwindows/survival_biasing/local/inputs_true.dat b/tests/regression_tests/weightwindows/survival_biasing/local/inputs_true.dat new file mode 100644 index 000000000..fe144dc62 --- /dev/null +++ b/tests/regression_tests/weightwindows/survival_biasing/local/inputs_true.dat @@ -0,0 +1,69 @@ + + + + + + + + + + + + + + + + + + + + + + + fixed source + 50 + 5 + + + + 0.01 1.0 + + + + + + + 14100000.0 1.0 + + + true + + 1 + neutron + 0.46135961568957096 0.2874485390202603 0.13256013950494605 0.052765209609807295 0.020958440832685638 0.006317250035749876 0.002627037175682506 0.00030032343592437284 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.4611178493390162 0.28624862560842024 0.13999870622621136 0.05508479408959756 0.018281241958254993 0.0055577764872361555 0.0015265432226293397 0.00024298772352636966 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.46294957913063317 0.2857734367546051 0.13365623190336945 0.05034742166227103 0.017525851812913266 0.005096725451851077 0.0026297640951531403 0.00033732424392026144 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.45772598750472554 0.281494921471495 0.13604397962762246 0.054792881472295364 0.019802376898755525 0.0073351745579128495 0.002067392946066426 9.529474085926143e-05 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.4739762983490318 0.28419061537046764 0.12757509901507888 0.0516920293019733 0.01919167874787235 0.007593035964707848 0.0017488176314107277 9.678267909681988e-05 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.47560199098558276 0.27728389146956994 0.12760802572535623 0.05251965811713552 0.022498960644951632 0.01063372459325151 0.004242071054914633 0.00045603112202665183 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.47659199471934693 0.28886666944252215 0.1258637573398161 0.05818862375955657 0.020476313720744762 0.007143193244018263 0.0022364083022515273 0.0003953123781408474 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.4693487369606823 0.27180862139490125 0.12833455570423483 0.055146743559938344 0.020667403529470333 0.007891508723137146 0.001643551705407358 0.0006501416781353278 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.4538724931023337 0.2824491421084296 0.1275341706351566 0.055027501141893705 0.02305114640508087 0.008599303158607198 0.0025082611194161804 0.000518530311231696 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.45258444753220123 0.26839102466484255 0.12991633918797083 0.05442263709947767 0.019695895223471486 0.005353204961887518 0.0015074698293840157 0.0001598680898180058 8.872047880811405e-05 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.4482467985982804 0.2813276470183359 0.13449860790768647 0.056103278190533436 0.026274320334196962 0.006884071908429303 0.0033099390738735458 0.0005337454397674922 0.0001386380399465575 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.4875748950139487 0.2908586705316248 0.14212060658392278 0.06068108009637272 0.019467970468789477 0.00637102427946399 0.0028510425266191687 0.0012184421778115488 0.00029404541567146187 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.5 0.2852101551297824 0.13622307172321643 0.058302519884339946 0.023124734915152625 0.007401527839638796 0.002527635652743992 0.0008509612315162482 0.00018213334574554768 7.815977984795461e-06 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.47706397688984176 0.28313436859495966 0.13170624744221976 0.0525453512766549 0.021111753158319105 0.0067451713955609645 0.003204270539718037 0.0010308107242263192 0.0002759210144794204 0.00013727805365387318 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.48641988659284513 0.2816173231471242 0.131541850061199 0.054793200248327956 0.016517701691156576 0.006757591781856257 0.0025425350750908644 0.0006383278085839443 2.279421641064226e-05 0.0003597432472224371 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.4307755159245372 0.2652713784080849 0.12753857957535653 0.05470396852881577 0.02049817309420563 0.00566559741247352 0.0007213846765465147 6.848024591311009e-05 -1.0 8.609308814924014e-05 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.43631324024956025 0.2639077825530567 0.13368737962742414 0.05213531603720763 0.020628860469743833 0.008268384844678654 0.0028599221013934375 0.00013059757129982714 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.45926287238219 0.2635824148883888 0.13104280996799478 0.055379812186041905 0.019255042561941074 0.007252095690246872 0.0018026593488140996 0.0003215836458542421 9.14840176000331e-06 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.4574378998400604 0.28714971179927723 0.1363507911051337 0.05010315954654347 0.017963278989571074 0.0058998176297971605 0.0021393135666168 0.00015677575033083482 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.4821442782978812 0.27704132233083506 0.14149074035662307 0.05546042815834048 0.016800649382269998 0.004096804603054233 0.0019161108398578026 0.0005289600253155307 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 + 2.306798078447855 1.4372426951013015 0.6628006975247303 0.2638260480490365 0.10479220416342819 0.03158625017874938 0.013135185878412529 0.0015016171796218643 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.3055892466950807 1.4312431280421012 0.6999935311310568 0.2754239704479878 0.09140620979127496 0.027788882436180776 0.007632716113146698 0.0012149386176318483 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.314747895653166 1.4288671837730256 0.6682811595168472 0.25173710831135515 0.08762925906456634 0.025483627259255386 0.013148820475765701 0.0016866212196013073 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.288629937523628 1.407474607357475 0.6802198981381123 0.2739644073614768 0.09901188449377762 0.036675872789564246 0.01033696473033213 0.00047647370429630714 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.369881491745159 1.4209530768523382 0.6378754950753944 0.2584601465098665 0.09595839373936176 0.03796517982353924 0.008744088157053638 0.00048391339548409945 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.3780099549279137 1.3864194573478497 0.6380401286267812 0.26259829058567763 0.11249480322475816 0.053168622966257545 0.021210355274573163 0.0022801556101332593 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.382959973596735 1.4443333472126108 0.6293187866990806 0.29094311879778284 0.10238156860372381 0.035715966220091315 0.011182041511257637 0.001976561890704237 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.3467436848034113 1.3590431069745064 0.6416727785211741 0.2757337177996917 0.10333701764735166 0.03945754361568573 0.00821775852703679 0.0032507083906766388 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.2693624655116684 1.412245710542148 0.637670853175783 0.27513750570946854 0.11525573202540434 0.04299651579303599 0.012541305597080901 0.00259265155615848 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.262922237661006 1.3419551233242126 0.6495816959398542 0.2721131854973884 0.09847947611735743 0.02676602480943759 0.007537349146920078 0.000799340449090029 0.0004436023940405703 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.2412339929914022 1.4066382350916795 0.6724930395384323 0.28051639095266717 0.1313716016709848 0.03442035954214652 0.016549695369367727 0.0026687271988374613 0.0006931901997327875 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.4378744750697434 1.4542933526581239 0.7106030329196139 0.3034054004818636 0.09733985234394739 0.03185512139731995 0.014255212633095843 0.006092210889057744 0.0014702270783573093 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.5 1.4260507756489118 0.6811153586160822 0.2915125994216997 0.11562367457576313 0.03700763919819398 0.012638178263719959 0.004254806157581241 0.0009106667287277384 3.9079889923977307e-05 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.385319884449209 1.4156718429747983 0.6585312372110987 0.2627267563832745 0.10555876579159552 0.03372585697780482 0.016021352698590185 0.005154053621131596 0.001379605072397102 0.000686390268269366 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.432099432964226 1.4080866157356209 0.657709250305995 0.2739660012416398 0.08258850845578287 0.03378795890928128 0.012712675375454322 0.0031916390429197216 0.0001139710820532113 0.0017987162361121855 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.1538775796226863 1.3263568920404245 0.6376928978767826 0.27351984264407886 0.10249086547102815 0.028327987062367603 0.0036069233827325737 0.0003424012295655504 -5.0 0.0004304654407462007 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.181566201247801 1.3195389127652835 0.6684368981371207 0.2606765801860381 0.10314430234871916 0.041341924223393264 0.014299610506967188 0.0006529878564991358 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.29631436191095 1.317912074441944 0.6552140498399739 0.27689906093020955 0.09627521280970537 0.03626047845123436 0.009013296744070498 0.0016079182292712106 4.574200880001655e-05 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.287189499200302 1.435748558996386 0.6817539555256685 0.25051579773271737 0.08981639494785537 0.029499088148985803 0.010696567833083998 0.0007838787516541741 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.410721391489406 1.3852066116541752 0.7074537017831153 0.2773021407917024 0.08400324691134999 0.020484023015271167 0.009580554199289014 0.0026448001265776534 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 + 3.0 + 10 + 1e-38 + + + 20 20 1 + 0.0 0.0 0.0 + 160.0 160.0 160.0 + + false + + true + true + + + + + 1 + + + 1 + flux + + + diff --git a/tests/regression_tests/weightwindows/survival_biasing/results_true.dat b/tests/regression_tests/weightwindows/survival_biasing/local/results_true.dat similarity index 100% rename from tests/regression_tests/weightwindows/survival_biasing/results_true.dat rename to tests/regression_tests/weightwindows/survival_biasing/local/results_true.dat diff --git a/tests/regression_tests/weightwindows/survival_biasing/inputs_true.dat b/tests/regression_tests/weightwindows/survival_biasing/shared/inputs_true.dat similarity index 99% rename from tests/regression_tests/weightwindows/survival_biasing/inputs_true.dat rename to tests/regression_tests/weightwindows/survival_biasing/shared/inputs_true.dat index d5aa135d4..bafff89c1 100644 --- a/tests/regression_tests/weightwindows/survival_biasing/inputs_true.dat +++ b/tests/regression_tests/weightwindows/survival_biasing/shared/inputs_true.dat @@ -51,6 +51,7 @@ 0.0 0.0 0.0 160.0 160.0 160.0 + true true true diff --git a/tests/regression_tests/weightwindows/survival_biasing/shared/results_true.dat b/tests/regression_tests/weightwindows/survival_biasing/shared/results_true.dat new file mode 100644 index 000000000..11e5ffa77 --- /dev/null +++ b/tests/regression_tests/weightwindows/survival_biasing/shared/results_true.dat @@ -0,0 +1 @@ +d17a437262d3316985fba4b48e21a7fcd5f32ac2d96ef0e66849f567deaeadc1e0f18d5d0bf5e9cab4246dbf823e7f43f249a1f67e928b27ae8c70b89f3cefbf \ No newline at end of file diff --git a/tests/regression_tests/weightwindows/survival_biasing/test.py b/tests/regression_tests/weightwindows/survival_biasing/test.py index 92604e338..c3a3e91a2 100644 --- a/tests/regression_tests/weightwindows/survival_biasing/test.py +++ b/tests/regression_tests/weightwindows/survival_biasing/test.py @@ -1,14 +1,16 @@ +from pathlib import Path + import pytest import numpy as np - import openmc -from openmc.stats import Discrete, Point +from openmc.utility_funcs import change_directory from tests.testing_harness import HashedPyAPITestHarness -@pytest.fixture -def model(): +def build_model(shared_secondary): + openmc.reset_auto_ids() + # Material w = openmc.Material(name='Tungsten') w.add_element('W', 1.0) @@ -38,13 +40,14 @@ def model(): angle = openmc.stats.Monodirectional((1.0, 0.0, 0.0)) energy = openmc.stats.Discrete([14.1e6], [1.0]) - source = openmc.Source(space=space, angle=angle, energy=energy) + source = openmc.IndependentSource(space=space, angle=angle, energy=energy) settings = openmc.Settings() settings.run_mode = 'fixed source' settings.batches = 5 settings.particles = 50 settings.source = source + settings.shared_secondary_bank = shared_secondary model = openmc.Model(geometry=geometry, materials=materials, settings=settings) @@ -61,7 +64,8 @@ def model(): tallies = openmc.Tallies([flux_tally]) model.tallies = tallies - lower_ww_bounds = np.loadtxt('ww_n.txt') + parent_dir = Path(__file__).parent + lower_ww_bounds = np.loadtxt(parent_dir / 'ww_n.txt') weight_windows = openmc.WeightWindows(mesh, lower_ww_bounds, @@ -76,6 +80,12 @@ def model(): return model -def test_weight_windows_with_survival_biasing(model): - harness = HashedPyAPITestHarness('statepoint.5.h5', model) - harness.main() +@pytest.mark.parametrize("shared_secondary,subdir", [ + (False, "local"), + (True, "shared"), +]) +def test_weight_windows_with_survival_biasing(shared_secondary, subdir): + with change_directory(subdir): + model = build_model(shared_secondary) + harness = HashedPyAPITestHarness('statepoint.5.h5', model) + harness.main() diff --git a/tests/regression_tests/weightwindows/test.py b/tests/regression_tests/weightwindows/test.py index c22ca1b75..1d3b063bd 100644 --- a/tests/regression_tests/weightwindows/test.py +++ b/tests/regression_tests/weightwindows/test.py @@ -1,14 +1,17 @@ +from pathlib import Path + import pytest import numpy as np import openmc from openmc.stats import Discrete, Point +from openmc.utility_funcs import change_directory from tests.testing_harness import HashedPyAPITestHarness -@pytest.fixture -def model(): +def build_model(shared_secondary): + openmc.reset_auto_ids() model = openmc.Model() # materials (M4 steel alloy) @@ -43,6 +46,7 @@ def model(): settings.batches = 2 settings.max_history_splits = 200 settings.photon_transport = True + settings.shared_secondary_bank = shared_secondary settings.weight_window_checkpoints = {'surface': True, 'collision': True} space = Point((0.001, 0.001, 0.001)) @@ -71,10 +75,10 @@ def model(): # weight windows - # load pre-generated weight windows - # (created using the same tally as above) - ww_n_lower_bnds = np.loadtxt('ww_n.txt') - ww_p_lower_bnds = np.loadtxt('ww_p.txt') + # load pre-generated weight windows from parent directory + parent_dir = Path(__file__).parent + ww_n_lower_bnds = np.loadtxt(parent_dir / 'ww_n.txt') + ww_p_lower_bnds = np.loadtxt(parent_dir / 'ww_p.txt') # create a mesh matching the one used # to generate the weight windows @@ -104,9 +108,15 @@ def model(): return model -def test_weightwindows(model): - test = HashedPyAPITestHarness('statepoint.2.h5', model) - test.main() +@pytest.mark.parametrize("shared_secondary,subdir", [ + (False, "local"), + (True, "shared"), +]) +def test_weightwindows(shared_secondary, subdir): + with change_directory(subdir): + model = build_model(shared_secondary) + test = HashedPyAPITestHarness('statepoint.2.h5', model) + test.main() def test_wwinp_cylindrical(): diff --git a/tests/regression_tests/weightwindows_pulse_height/__init__.py b/tests/regression_tests/weightwindows_pulse_height/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/regression_tests/weightwindows_pulse_height/local/inputs_true.dat b/tests/regression_tests/weightwindows_pulse_height/local/inputs_true.dat new file mode 100644 index 000000000..4cd24d529 --- /dev/null +++ b/tests/regression_tests/weightwindows_pulse_height/local/inputs_true.dat @@ -0,0 +1,60 @@ + + + + + + + + + + + + + + + + + fixed source + 100 + 5 + + + 1000000.0 1.0 + + + true + + 1 + photon + 0.0 2000000.0 + 0.01 + 0.05 + 3.0 + 10 + 1e-38 + + + 1 1 1 + -2 -2 -2 + 2 2 2 + + false + + true + true + + 50 + + + + 1 + + + 0.0 10000.0 20000.0 30000.0 40000.0 50000.0 60000.0 70000.0 80000.0 90000.0 100000.0 110000.0 120000.0 130000.0 140000.0 150000.0 160000.0 170000.0 180000.0 190000.0 200000.0 210000.0 220000.0 230000.0 240000.0 250000.0 260000.0 270000.0 280000.0 290000.0 300000.0 310000.0 320000.0 330000.0 340000.0 350000.0 360000.0 370000.0 380000.0 390000.0 400000.0 410000.0 420000.0 430000.0 440000.0 450000.0 460000.0 470000.0 480000.0 490000.0 500000.0 510000.0 520000.0 530000.0 540000.0 550000.0 560000.0 570000.0 580000.0 590000.0 600000.0 610000.0 620000.0 630000.0 640000.0 650000.0 660000.0 670000.0 680000.0 690000.0 700000.0 710000.0 720000.0 730000.0 740000.0 750000.0 760000.0 770000.0 780000.0 790000.0 800000.0 810000.0 820000.0 830000.0 840000.0 850000.0 860000.0 870000.0 880000.0 890000.0 900000.0 910000.0 920000.0 930000.0 940000.0 950000.0 960000.0 970000.0 980000.0 990000.0 1000000.0 + + + 1 2 + pulse-height + + + diff --git a/tests/regression_tests/weightwindows_pulse_height/local/results_true.dat b/tests/regression_tests/weightwindows_pulse_height/local/results_true.dat new file mode 100644 index 000000000..c57e8ff1c --- /dev/null +++ b/tests/regression_tests/weightwindows_pulse_height/local/results_true.dat @@ -0,0 +1,201 @@ +tally 1: +4.140000E+00 +3.443000E+00 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +3.000000E-02 +5.000000E-04 +2.000000E-02 +4.000000E-04 +2.000000E-02 +4.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +3.000000E-02 +3.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +4.000000E-02 +6.000000E-04 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +2.000000E-02 +2.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +2.000000E-02 +4.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +2.000000E-02 +2.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +3.000000E-02 +5.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +2.000000E-02 +2.000000E-04 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +2.000000E-02 +2.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +2.000000E-02 +2.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +3.000000E-02 +5.000000E-04 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +2.000000E-02 +2.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +2.000000E-01 +8.600000E-03 diff --git a/tests/regression_tests/weightwindows_pulse_height/shared/inputs_true.dat b/tests/regression_tests/weightwindows_pulse_height/shared/inputs_true.dat new file mode 100644 index 000000000..8bb499a73 --- /dev/null +++ b/tests/regression_tests/weightwindows_pulse_height/shared/inputs_true.dat @@ -0,0 +1,60 @@ + + + + + + + + + + + + + + + + + fixed source + 100 + 5 + + + 1000000.0 1.0 + + + true + + 1 + photon + 0.0 2000000.0 + 0.01 + 0.05 + 3.0 + 10 + 1e-38 + + + 1 1 1 + -2 -2 -2 + 2 2 2 + + true + + true + true + + 50 + + + + 1 + + + 0.0 10000.0 20000.0 30000.0 40000.0 50000.0 60000.0 70000.0 80000.0 90000.0 100000.0 110000.0 120000.0 130000.0 140000.0 150000.0 160000.0 170000.0 180000.0 190000.0 200000.0 210000.0 220000.0 230000.0 240000.0 250000.0 260000.0 270000.0 280000.0 290000.0 300000.0 310000.0 320000.0 330000.0 340000.0 350000.0 360000.0 370000.0 380000.0 390000.0 400000.0 410000.0 420000.0 430000.0 440000.0 450000.0 460000.0 470000.0 480000.0 490000.0 500000.0 510000.0 520000.0 530000.0 540000.0 550000.0 560000.0 570000.0 580000.0 590000.0 600000.0 610000.0 620000.0 630000.0 640000.0 650000.0 660000.0 670000.0 680000.0 690000.0 700000.0 710000.0 720000.0 730000.0 740000.0 750000.0 760000.0 770000.0 780000.0 790000.0 800000.0 810000.0 820000.0 830000.0 840000.0 850000.0 860000.0 870000.0 880000.0 890000.0 900000.0 910000.0 920000.0 930000.0 940000.0 950000.0 960000.0 970000.0 980000.0 990000.0 1000000.0 + + + 1 2 + pulse-height + + + diff --git a/tests/regression_tests/weightwindows_pulse_height/shared/results_true.dat b/tests/regression_tests/weightwindows_pulse_height/shared/results_true.dat new file mode 100644 index 000000000..c57e8ff1c --- /dev/null +++ b/tests/regression_tests/weightwindows_pulse_height/shared/results_true.dat @@ -0,0 +1,201 @@ +tally 1: +4.140000E+00 +3.443000E+00 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +3.000000E-02 +5.000000E-04 +2.000000E-02 +4.000000E-04 +2.000000E-02 +4.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +3.000000E-02 +3.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +4.000000E-02 +6.000000E-04 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +2.000000E-02 +2.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +2.000000E-02 +4.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +2.000000E-02 +2.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +3.000000E-02 +5.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +2.000000E-02 +2.000000E-04 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +2.000000E-02 +2.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +2.000000E-02 +2.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +3.000000E-02 +5.000000E-04 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +2.000000E-02 +2.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +2.000000E-01 +8.600000E-03 diff --git a/tests/regression_tests/weightwindows_pulse_height/test.py b/tests/regression_tests/weightwindows_pulse_height/test.py new file mode 100644 index 000000000..b6b49be26 --- /dev/null +++ b/tests/regression_tests/weightwindows_pulse_height/test.py @@ -0,0 +1,73 @@ +import numpy as np +import openmc +import pytest +from openmc.utility_funcs import change_directory + +from tests.testing_harness import PyAPITestHarness + + +@pytest.mark.parametrize("shared_secondary,subdir", [ + (False, "local"), + (True, "shared"), +]) +def test_weightwindows_pulse_height(shared_secondary, subdir): + with change_directory(subdir): + openmc.reset_auto_ids() + model = openmc.Model() + + # Define materials (NaI scintillator) + NaI = openmc.Material() + NaI.set_density('g/cc', 3.7) + NaI.add_element('Na', 1.0) + NaI.add_element('I', 1.0) + + model.materials = openmc.Materials([NaI]) + + # Define geometry: NaI sphere inside vacuum sphere + s1 = openmc.Sphere(r=1) + s2 = openmc.Sphere(r=2, boundary_type='vacuum') + inner_sphere = openmc.Cell(name='inner sphere', fill=NaI, region=-s1) + outer_sphere = openmc.Cell(name='outer sphere', region=+s1 & -s2) + model.geometry = openmc.Geometry([inner_sphere, outer_sphere]) + + # Define settings + model.settings.run_mode = 'fixed source' + model.settings.batches = 5 + model.settings.particles = 100 + model.settings.photon_transport = True + model.settings.shared_secondary_bank = shared_secondary + model.settings.max_history_splits = 50 + model.settings.weight_window_checkpoints = { + 'surface': True, + 'collision': True, + } + model.settings.source = openmc.IndependentSource( + energy=openmc.stats.delta_function(1e6), + particle='photon' + ) + + # Define pulse-height tally + tally = openmc.Tally(name="pht tally") + tally.scores = ['pulse-height'] + cell_filter = openmc.CellFilter(inner_sphere) + energy_filter = openmc.EnergyFilter(np.linspace(0, 1_000_000, 101)) + tally.filters = [cell_filter, energy_filter] + model.tallies = [tally] + + # Define weight windows on a simple mesh covering the geometry + ww_mesh = openmc.RegularMesh() + ww_mesh.lower_left = (-2, -2, -2) + ww_mesh.upper_right = (2, 2, 2) + ww_mesh.dimension = (1, 1, 1) + + # Single energy bin for photons + e_bnds = [0.0, 2e6] + + # Uniform weight window bounds (low enough to trigger some splitting) + lower_bounds = np.array([0.01]) + + ww = openmc.WeightWindows(ww_mesh, lower_bounds, None, 5.0, e_bnds, 'photon') + model.settings.weight_windows = [ww] + + harness = PyAPITestHarness('statepoint.5.h5', model) + harness.main() diff --git a/tests/unit_tests/weightwindows/test.py b/tests/unit_tests/weightwindows/test.py index d6e509522..efa203b51 100644 --- a/tests/unit_tests/weightwindows/test.py +++ b/tests/unit_tests/weightwindows/test.py @@ -122,7 +122,8 @@ def model(): return model -def test_weightwindows(model, wws): +@pytest.mark.parametrize("shared_secondary", [False, True]) +def test_weightwindows(model, wws, shared_secondary): ww_files = ('ww_n.txt', 'ww_p.txt') cwd = Path(__file__).parent.absolute() @@ -131,6 +132,7 @@ def test_weightwindows(model, wws): with cdtemp(filepaths): # run once with variance reduction off model.settings.weight_windows_on = False + model.settings.shared_secondary_bank = shared_secondary analog_sp = model.run() os.rename(analog_sp, 'statepoint.analog.h5') @@ -223,7 +225,8 @@ def test_lower_ww_bounds_shape(): assert ww.lower_ww_bounds.shape == (2, 3, 4, 1) -def test_photon_heating(run_in_tmpdir): +@pytest.mark.parametrize("shared_secondary", [False, True]) +def test_photon_heating(run_in_tmpdir, shared_secondary): water = openmc.Material() water.add_nuclide('H1', 1.0) water.add_nuclide('O16', 2.0) @@ -246,7 +249,8 @@ def test_photon_heating(run_in_tmpdir): model.settings.run_mode = 'fixed source' model.settings.batches = 5 - model.settings.particles = 100 + model.settings.particles = 101 + model.settings.shared_secondary_bank = shared_secondary tally = openmc.Tally() tally.scores = ['heating'] @@ -260,7 +264,11 @@ def test_photon_heating(run_in_tmpdir): with openmc.StatePoint(sp_file) as sp: tally_mean = sp.tallies[tally.id].mean - # these values should be nearly identical + # Note: Our current physics model actually does allow this tally to + # occasionally go slightly negative. However, larger bugs can + # make this more common. We have selected a particle count for + # this test that happens to produce no negative tallies for both + # the shared and non-shared secondary PRNG streams. assert np.all(tally_mean >= 0) From 3c7a030d43d79a9784982bc37a1c9c87cba5a375 Mon Sep 17 00:00:00 2001 From: Hridoy Kabiraj Date: Thu, 21 May 2026 00:49:00 +0600 Subject: [PATCH 07/10] Guard average molar mass and validate materials depletion inputs (#3941) --- openmc/material.py | 13 +++++++- tests/unit_tests/test_materials.py | 49 ++++++++++++++++++++++++++++++ 2 files changed, 61 insertions(+), 1 deletion(-) diff --git a/openmc/material.py b/openmc/material.py index 86cc3d793..eded9c3fa 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -296,6 +296,8 @@ class Material(IDManagerMixin): mass += nuc.percent # Compute and return the molar mass + if moles == 0.0: + raise ValueError("Material has no nuclides; cannot compute molar mass") return mass / moles @property @@ -2287,7 +2289,7 @@ class Materials(cv.CheckedList): multigroup_fluxes: Sequence[Sequence[float]] Energy-dependent multigroup flux values, where each sublist corresponds to a specific material. Will be normalized so that it sums to 1. - energy_group_structures': Sequence[Sequence[float] | str] + energy_group_structures: Sequence[Sequence[float] | str] Energy group boundaries in [eV] or the name of the group structure. timesteps : iterable of float or iterable of tuple Array of timesteps. Note that values are not cumulative. The units are @@ -2322,6 +2324,11 @@ class Materials(cv.CheckedList): for mat in self: mat.depletable = True + if len(multigroup_fluxes) != len(self): + raise ValueError("multigroup_fluxes length must match number of materials") + if len(energy_group_structures) != len(self): + raise ValueError("energy_group_structures length must match number of materials") + chain = _get_chain(chain_file) # Create MicroXS objects for all materials @@ -2332,6 +2339,10 @@ class Materials(cv.CheckedList): for material, flux, energy in zip( self, multigroup_fluxes, energy_group_structures ): + if material.volume is None: + raise ValueError( + f"Material {material.id} has no volume; cannot deplete" + ) temperature = material.temperature or 293.6 micro_xs = openmc.deplete.MicroXS.from_multigroup_flux( energies=energy, diff --git a/tests/unit_tests/test_materials.py b/tests/unit_tests/test_materials.py index 5a382b777..6620402fd 100644 --- a/tests/unit_tests/test_materials.py +++ b/tests/unit_tests/test_materials.py @@ -1,5 +1,7 @@ from pathlib import Path +import pytest + import openmc from openmc.deplete import Chain @@ -78,3 +80,50 @@ def test_export_duplicate_materials_to_xml(run_in_tmpdir): materials_in = openmc.Materials.from_xml("materials.xml") assert len(materials_in) == 2 + + +def test_materials_deplete_length_mismatch(): + mats = openmc.Materials([openmc.Material()]) + + with pytest.raises(ValueError, match="multigroup_fluxes length"): + mats.deplete( + multigroup_fluxes=[], + energy_group_structures=["VITAMIN-J-42"], + timesteps=[1.0], + source_rates=1.0, + ) + + with pytest.raises(ValueError, match="energy_group_structures length"): + mats.deplete( + multigroup_fluxes=[[1.0]], + energy_group_structures=[], + timesteps=[1.0], + source_rates=1.0, + ) + + +def test_materials_deplete_missing_volume(monkeypatch): + mat = openmc.Material() + mat.add_nuclide("Ni58", 1.0) + mat.set_density("g/cm3", 7.87) + + mats = openmc.Materials([mat]) + + class DummySession: + def __enter__(self): + return self + + def __exit__(self, exc_type, exc, tb): + return False + + monkeypatch.setattr(openmc.lib, "TemporarySession", DummySession) + + chain = Path(__file__).parents[1] / "chain_ni.xml" + with pytest.raises(ValueError, match="has no volume"): + mats.deplete( + multigroup_fluxes=[[1.0]], + energy_group_structures=["VITAMIN-J-42"], + timesteps=[1.0], + source_rates=1.0, + chain_file=chain, + ) From 66497e76b11653554676c2ac305a14f764d82c12 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Wed, 20 May 2026 14:01:21 -0500 Subject: [PATCH 08/10] Support writing of hex elements to VTKHDF format. (#3623) Co-authored-by: Paul Romano --- docs/source/usersguide/processing.rst | 10 ++- openmc/mesh.py | 107 ++++++++++++-------------- tests/unit_tests/test_mesh.py | 45 ++++++++++- tests/unit_tests/test_mesh_hexes.exo | 1 + 4 files changed, 102 insertions(+), 61 deletions(-) create mode 120000 tests/unit_tests/test_mesh_hexes.exo diff --git a/docs/source/usersguide/processing.rst b/docs/source/usersguide/processing.rst index fe6ab0826..ed3e8564a 100644 --- a/docs/source/usersguide/processing.rst +++ b/docs/source/usersguide/processing.rst @@ -97,7 +97,15 @@ VTK Mesh File Generation ------------------------ VTK files of OpenMC meshes can be created using the -:meth:`openmc.Mesh.write_data_to_vtk` method. Data can be applied to the +:meth:`openmc.Mesh.write_data_to_vtk` method. This method supports several VTK +formats depending on the mesh type. Structured meshes +(:class:`~openmc.RegularMesh`, :class:`~openmc.RectilinearMesh`, +:class:`~openmc.CylindricalMesh`, and :class:`~openmc.SphericalMesh`) can be +exported to legacy VTK format (``.vtk``). The :class:`~openmc.UnstructuredMesh` +class supports VTK unstructured grid formats (``.vtu``) as well as an HDF5-based +format (``.vtkhdf``) that does not require the ``vtk`` module to write. + +Data can be applied to the elements of the resulting mesh from mesh filter objects. This data can be provided either as a flat array or, in the case of structured meshes (:class:`~openmc.RegularMesh`, :class:`~openmc.RectilinearMesh`, diff --git a/openmc/mesh.py b/openmc/mesh.py index 4129913d8..61effd850 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -2738,7 +2738,8 @@ class UnstructuredMesh(MeshBase): _UNSUPPORTED_ELEM = -1 _LINEAR_TET = 0 _LINEAR_HEX = 1 - _VTK_TETRA = 10 + _VTK_TET = 10 + _VTK_HEX = 12 def __init__(self, filename: PathLike, library: str, mesh_id: int | None = None, name: str = '', length_multiplier: float = 1.0, @@ -3049,7 +3050,9 @@ class UnstructuredMesh(MeshBase): n_skipped += 1 continue else: - raise RuntimeError(f"Invalid element type {elem_type} found") + raise RuntimeError( + f"Invalid element type {elem_type} found in mesh {self.id}" + ) for i, c in enumerate(conn): if c == -1: @@ -3106,35 +3109,42 @@ class UnstructuredMesh(MeshBase): datasets: dict | None = None, volume_normalization: bool = True, ): - def append_dataset(dset, array): - """Convenience function to append data to an HDF5 dataset""" - origLen = dset.shape[0] - dset.resize(origLen + array.shape[0], axis=0) - dset[origLen:] = array + # This writer supports linear tetrahedra and linear hexahedra elements + conn_list = [] # flattened connectivity ids + cell_sizes = [] # number of points per cell + vtk_types = [] # VTK cell types per cell (uint8) + n_skipped = 0 - if self.library != "moab": - raise NotImplementedError("VTKHDF output is only supported for MOAB meshes") - - # the self.connectivity contains arrays of length 8 to support hex - # elements as well, in the case of tetrahedra mesh elements, the - # last 4 values are -1 and are removed - trimmed_connectivity = [] - for cell in self.connectivity: - # Find the index of the first -1 value, if any - first_negative_index = np.where(cell == -1)[0] - if first_negative_index.size > 0: - # Slice the array up to the first -1 value - trimmed_connectivity.append(cell[: first_negative_index[0]]) + for conn, etype in zip(self.connectivity, self.element_types): + if etype == self._LINEAR_TET: + ids = conn[:4] + vtk_types.append(self._VTK_TET) + elif etype == self._LINEAR_HEX: + ids = conn[:8] + vtk_types.append(self._VTK_HEX) + elif etype == self._UNSUPPORTED_ELEM: + n_skipped += 1 + continue else: - # No -1 values, append the whole cell - trimmed_connectivity.append(cell) - trimmed_connectivity = np.array(trimmed_connectivity, dtype="int32").flatten() + raise RuntimeError( + f"Invalid element type {etype} found in mesh {self.id}" + ) + conn_list.extend(ids.tolist()) + cell_sizes.append(len(ids)) - # MOAB meshes supports tet elements only so we know it has 4 points per cell - points_per_cell = 4 + if n_skipped > 0: + warnings.warn( + f"{n_skipped} elements were not written because " + "they are not of type linear tet/hex" + ) - # offsets are the indices of the first point of each cell in the array of points - offsets = np.arange(0, self.n_elements * points_per_cell + 1, points_per_cell) + connectivity = np.asarray(conn_list, dtype=np.int64) + + # Offsets must be length (numCells + 1) with a leading 0 and + # cumulative end-indices thereafter, per VTK's layout + cell_sizes_arr = np.asarray(cell_sizes, dtype=np.int64) + offsets = np.zeros(cell_sizes_arr.size + 1, dtype=np.int64) + np.cumsum(cell_sizes_arr, out=offsets[1:]) for name, data in datasets.items(): if data.shape != self.dimension: @@ -3156,42 +3166,27 @@ class UnstructuredMesh(MeshBase): dtype=h5py.string_dtype("ascii", len(ascii_type)), ) - # create hdf5 file structure - root.create_dataset("NumberOfPoints", (0,), maxshape=(None,), dtype="i8") - root.create_dataset("Types", (0,), maxshape=(None,), dtype="uint8") - root.create_dataset("Points", (0, 3), maxshape=(None, 3), dtype="f") - root.create_dataset( - "NumberOfConnectivityIds", (0,), maxshape=(None,), dtype="i8" - ) - root.create_dataset("NumberOfCells", (0,), maxshape=(None,), dtype="i8") - root.create_dataset("Offsets", (0,), maxshape=(None,), dtype="i8") - root.create_dataset("Connectivity", (0,), maxshape=(None,), dtype="i8") + # Create HDF5 file structure compliant with VTKHDF UnstructuredGrid + n_points = int(len(self.vertices)) + n_cells = int(len(cell_sizes)) + n_conn_ids = int(len(connectivity)) - append_dataset(root["NumberOfPoints"], np.array([len(self.vertices)])) - append_dataset(root["Points"], self.vertices) - append_dataset( - root["NumberOfConnectivityIds"], - np.array([len(trimmed_connectivity)]), - ) - append_dataset(root["Connectivity"], trimmed_connectivity) - append_dataset(root["NumberOfCells"], np.array([self.n_elements])) - append_dataset(root["Offsets"], offsets) - - append_dataset( - root["Types"], np.full(self.n_elements, self._VTK_TETRA, dtype="uint8") - ) + root.create_dataset("NumberOfPoints", data=(n_points,), dtype="i8") + root.create_dataset("NumberOfCells", data=(n_cells,), dtype="i8") + root.create_dataset("NumberOfConnectivityIds", data=(n_conn_ids,), dtype="i8") + root.create_dataset("Points", data=self.vertices.astype(np.float64, copy=False), dtype="f8") + root.create_dataset("Types", data=np.asarray(vtk_types, dtype=np.uint8), dtype="uint8") + root.create_dataset("Offsets", data=offsets.astype("i8"), dtype="i8") + root.create_dataset("Connectivity", data=connectivity.astype("i8"), dtype="i8") cell_data_group = root.create_group("CellData") for name, data in datasets.items(): - - cell_data_group.create_dataset( - name, (0,), maxshape=(None,), dtype="float64", chunks=True - ) - if volume_normalization: data /= self.volumes - append_dataset(cell_data_group[name], data) + cell_data_group.create_dataset( + name, data=data, dtype="float64", chunks=True + ) @classmethod def from_hdf5(cls, group: h5py.Group, mesh_id: int, name: str): diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py index 0b28bdfbe..85cb68a83 100644 --- a/tests/unit_tests/test_mesh.py +++ b/tests/unit_tests/test_mesh.py @@ -490,13 +490,28 @@ def test_umesh(run_in_tmpdir, simple_umesh, export_type): simple_umesh.write_data_to_vtk(datasets={'mean': ref_data[:-2]}, filename=filename) -@pytest.mark.skipif(not openmc.lib._dagmc_enabled(), reason="DAGMC not enabled.") -def test_write_vtkhdf(request, run_in_tmpdir): +vtkhdf_tests = [ + ( + Path("test_mesh_dagmc_tets.vtk"), + "moab" + ), + ( + Path("test_mesh_hexes.exo"), + "libmesh" + ) +] +@pytest.mark.parametrize('mesh_file, mesh_library', vtkhdf_tests) +def test_write_vtkhdf(mesh_file, mesh_library, request, run_in_tmpdir): """Performs a minimal UnstructuredMesh simulation, reads in the resulting statepoint file and writes the mesh data to vtk and vtkhdf files. It is necessary to read in the unstructured mesh from a statepoint file to ensure it has all the required attributes """ + if mesh_library == 'moab' and not openmc.lib._dagmc_enabled(): + pytest.skip("DAGMC not enabled.") + if mesh_library == 'libmesh' and not openmc.lib._libmesh_enabled(): + pytest.skip("LibMesh not enabled.") + model = openmc.Model() surf1 = openmc.Sphere(r=1000.0, boundary_type="vacuum") @@ -504,8 +519,8 @@ def test_write_vtkhdf(request, run_in_tmpdir): model.geometry = openmc.Geometry([cell1]) umesh = openmc.UnstructuredMesh( - request.path.parent / "test_mesh_dagmc_tets.vtk", - "moab", + request.path.parent / mesh_file, + mesh_library, mesh_id = 1 ) mesh_filter = openmc.MeshFilter(umesh) @@ -514,6 +529,8 @@ def test_write_vtkhdf(request, run_in_tmpdir): mesh_tally = openmc.Tally(name="test_tally") mesh_tally.filters = [mesh_filter] mesh_tally.scores = ["flux"] + if mesh_library == "libmesh": + mesh_tally.estimator = "collision" model.tallies = [mesh_tally] @@ -556,6 +573,26 @@ def test_write_vtkhdf(request, run_in_tmpdir): with h5py.File("test_mesh.vtkhdf", "r"): ... + import vtk + reader = vtk.vtkHDFReader() + reader.SetFileName("test_mesh.vtkhdf") + reader.Update() + + # Get mean from file and make sure it matches original data + num_elements = reader.GetOutput().GetNumberOfCells() + assert num_elements == umesh_from_sp.n_elements + + num_vertices = reader.GetOutput().GetNumberOfPoints() + assert num_vertices == umesh_from_sp.vertices.shape[0] + + arr = reader.GetOutput().GetCellData().GetArray("mean") + mean = np.array([arr.GetTuple1(i) for i in range(my_tally.mean.size)]) + np.testing.assert_almost_equal(mean, my_tally.mean.flatten()/umesh_from_sp.volumes) + + arr = reader.GetOutput().GetCellData().GetArray("std_dev") + std_dev = np.array([arr.GetTuple1(i) for i in range(my_tally.std_dev.size)]) + np.testing.assert_almost_equal(std_dev, my_tally.std_dev.flatten()/umesh_from_sp.volumes) + def test_mesh_get_homogenized_materials(): """Test the get_homogenized_materials method""" diff --git a/tests/unit_tests/test_mesh_hexes.exo b/tests/unit_tests/test_mesh_hexes.exo new file mode 120000 index 000000000..272da92a4 --- /dev/null +++ b/tests/unit_tests/test_mesh_hexes.exo @@ -0,0 +1 @@ +../regression_tests/unstructured_mesh/test_mesh_hexes.exo \ No newline at end of file From 7d09a12606850b37d57cb61372ad2edc4ed327c6 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Wed, 20 May 2026 16:43:21 -0500 Subject: [PATCH 09/10] DAGMC Cell Override Updates (#3888) Co-authored-by: Claude Sonnet 4.6 Co-authored-by: Paul Romano --- docs/source/io_formats/geometry.rst | 61 +++++-- include/openmc/cell.h | 24 +++ include/openmc/dagmc.h | 19 +- openmc/dagmc.py | 230 ++++++++++++++---------- openmc/model/model.py | 2 + src/cell.cpp | 94 ++++++---- src/dagmc.cpp | 227 +++++++++++++++++++----- tests/unit_tests/dagmc/test_model.py | 253 +++++++++++++++++++-------- tests/unit_tests/test_cell.py | 43 +++++ 9 files changed, 689 insertions(+), 264 deletions(-) diff --git a/docs/source/io_formats/geometry.rst b/docs/source/io_formats/geometry.rst index cef3bb79c..60ff0afae 100644 --- a/docs/source/io_formats/geometry.rst +++ b/docs/source/io_formats/geometry.rst @@ -406,24 +406,55 @@ Each ```` element can have the following attributes or sub-eleme *Default*: None - :material_overrides: - This element contains information on material overrides to be applied to the - DAGMC universe. It has the following attributes and sub-elements: + :cell: + Zero or more ```` sub-elements may appear to override properties of + individual DAGMC volumes. Each ```` element supports the following + attributes and sub-elements: - :cell: - Material override information for a single cell. It contains the following - attributes and sub-elements: + :id: + The integer cell ID in the DAGMC geometry to override. Required. - :id: - The cell ID in the DAGMC geometry for which the material override will - apply. + :name: + An optional string label for the cell. - :materials: - A list of material IDs that will apply to instances of the cell. If the - list contains only one ID, it will replace the original material - assignment of all instances of the DAGMC cell. If the list contains more - than one material, each material ID of the list will be assigned to the - various instances of the DAGMC cell. + *Default*: None + + :material: + The material ID to assign to this cell. Use ``void`` for vacuum. Multiple + space-separated IDs may be given to specify a distribmat (distributed + material) assignment. Required. + + :temperature: + Temperature(s) in [K] to assign to the cell. Must be ≥ 0. Multiple + space-separated values may be given. + + *Default*: None + + :density: + Density in [g/cm³] to assign to the cell. Must be > 0. Requires a non-void + material fill. Multiple space-separated values may be given. + + *Default*: None + + :volume: + Volume of the cell in [cm³]. + + .. note:: DAGMC can compute cell volumes exactly from the triangulated + mesh surfaces. Specifying a manual volume risks inconsistency + with that capability. + + *Default*: None + + The following standard ```` attributes are **not** supported inside + ```` and will raise an error if present: ``region``, + ``fill``, ``universe``, ``translation``, ``rotation``. + + .. deprecated:: + The ```` sub-element (containing ```` + children with ````) is deprecated. A deprecation warning is + emitted and the overrides are converted to the ```` format at parse + time. It is an error to specify both ```` and + ```` sub-elements on the same ````. *Default*: None diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 26c34ebd4..1a27ad516 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -145,6 +145,30 @@ private: bool simple_; //!< Does the region contain only intersections? }; +//============================================================================== +// XML parsing helpers for nodes +//============================================================================== + +//! Parse material IDs from a XML node. +//! \param node XML node containing a "material" attribute or child element +//! \param cell_id Cell ID used in error messages +//! \return Vector of material IDs (MATERIAL_VOID for "void") +vector parse_cell_material_xml(pugi::xml_node node, int32_t cell_id); + +//! Parse temperatures in [K] from a XML node. +//! Validates that all values are non-negative and the list is non-empty. +//! \param node XML node containing a "temperature" attribute or child element +//! \param cell_id Cell ID used in error messages +//! \return Vector of temperatures in [K] +vector parse_cell_temperature_xml(pugi::xml_node node, int32_t cell_id); + +//! Parse densities in [g/cm³] from a XML node. +//! Validates that all values are positive and the list is non-empty. +//! \param node XML node containing a "density" attribute or child element +//! \param cell_id Cell ID used in error messages +//! \return Vector of densities in [g/cm³] +vector parse_cell_density_xml(pugi::xml_node node, int32_t cell_id); + //============================================================================== class Cell { diff --git a/include/openmc/dagmc.h b/include/openmc/dagmc.h index 0e27402a1..858636c0b 100644 --- a/include/openmc/dagmc.h +++ b/include/openmc/dagmc.h @@ -94,6 +94,10 @@ private: class DAGUniverse : public Universe { public: + using MaterialOverrides = std::unordered_map>; + using TemperatureOverrides = std::unordered_map>; + using DensityOverrides = std::unordered_map>; + explicit DAGUniverse(pugi::xml_node node); //! Create a new DAGMC universe @@ -112,6 +116,9 @@ public: //! Initialize the DAGMC accel. data structures, indices, material //! assignments, etc. void initialize(); + void initialize(const MaterialOverrides& material_overrides, + const TemperatureOverrides& temperature_overrides, + const DensityOverrides& density_overrides = {}); //! Reads UWUW materials and returns an ID map void read_uwuw_materials(); @@ -146,7 +153,8 @@ public: //! Assign a material overriding normal assignement to a cell //! \param[in] c The OpenMC cell to which the material is assigned - void override_assign_material(std::unique_ptr& c) const; + void override_assign_material(std::unique_ptr& c, + const MaterialOverrides& material_overrides) const; //! Return the index into the model cells vector for a given DAGMC volume //! handle in the universe @@ -187,7 +195,9 @@ private: void set_id(); //!< Deduce the universe id from model::universes void init_dagmc(); //!< Create and initialise DAGMC pointer void init_metadata(); //!< Create and initialise dagmcMetaData pointer - void init_geometry(); //!< Create cells and surfaces from DAGMC entities + void init_geometry(const MaterialOverrides& material_overrides, + const TemperatureOverrides& temperature_overrides, + const DensityOverrides& density_overrides); std::string filename_; //!< Name of the DAGMC file used to create this universe @@ -201,11 +211,6 @@ private: //!< generate new material IDs for the universe bool has_graveyard_; //!< Indicates if the DAGMC geometry has a "graveyard" //!< volume - std::unordered_map> - material_overrides_; //!< Map of material overrides - //!< keys correspond to the DAGMCCell id - //!< values are a list of material ids used - //!< for the override }; //============================================================================== diff --git a/openmc/dagmc.py b/openmc/dagmc.py index fd4258225..75ead5e43 100644 --- a/openmc/dagmc.py +++ b/openmc/dagmc.py @@ -36,12 +36,6 @@ class DAGMCUniverse(openmc.UniverseBase): auto_mat_ids : bool Set IDs automatically on initialization (True) or report overlaps in ID space between OpenMC and UWUW materials (False) - material_overrides : dict, optional - A dictionary of material overrides. The keys are material name strings - and the values are Iterables of openmc.Material objects. If a material - name is found in the DAGMC file, the material will be replaced with the - openmc.Material object in the value. - Attributes ---------- id : int @@ -78,15 +72,6 @@ class DAGMCUniverse(openmc.UniverseBase): The number of surfaces in the model. .. versionadded:: 0.13.2 - material_overrides : dict - A dictionary of material overrides. Keys are cell IDs; values are - iterables of :class:`openmc.Material` objects. The material assignment - of each DAGMC cell ID key will be replaced with the - :class:`~openmc.Material` object in the value. If the value contains - multiple :class:`~openmc.Material` objects, each Material in the list - will be assigned to the corresponding instance of the cell. - - .. versionadded:: 0.15.1 """ def __init__(self, @@ -94,16 +79,12 @@ class DAGMCUniverse(openmc.UniverseBase): universe_id=None, name='', auto_geom_ids=False, - auto_mat_ids=False, - material_overrides=None): + auto_mat_ids=False): super().__init__(universe_id, name) # Initialize class attributes self.filename = filename self.auto_geom_ids = auto_geom_ids self.auto_mat_ids = auto_mat_ids - self._material_overrides = {} - if material_overrides is not None: - self.material_overrides = material_overrides def __repr__(self): string = super().__repr__() @@ -130,47 +111,17 @@ class DAGMCUniverse(openmc.UniverseBase): @property def material_overrides(self): - return self._material_overrides + raise AttributeError( + "DAGMCUniverse.material_overrides has been removed. Use " + "DAGMCCell objects added via add_cell() to manage per-cell " + "material assignments.") @material_overrides.setter def material_overrides(self, val): - cv.check_type('material overrides', val, Mapping) - for key, value in val.items(): - self.add_material_override(key, value) - - def replace_material_assignment(self, material_name: str, material: openmc.Material): - """Replace a material assignment within the DAGMC universe. - - Replace the material assignment of all cells filled with a material in - the DAGMC universe. The universe must be synchronized in an initialized - Model (see :meth:`~openmc.DAGMCUniverse.sync_dagmc_cells`) before - calling this method. - - .. versionadded:: 0.15.1 - - Parameters - ---------- - material_name : str - Material name to replace - material : openmc.Material - Material to replace the material_name with - - """ - if material_name not in self.material_names: - raise ValueError( - f"No material with name '{material_name}' found in the DAGMC universe") - - if not self.cells: - raise RuntimeError("This DAGMC universe has not been synchronized " - "on an initialized Model.") - - for cell in self.cells.values(): - if cell.fill is None: - continue - if isinstance(cell.fill, openmc.Iterable): - cell.fill = list(map(lambda x: material if x.name == material_name else x, cell.fill)) - else: - cell.fill = material if cell.fill.name == material_name else cell.fill + raise AttributeError( + "DAGMCUniverse.material_overrides has been removed. Use " + "DAGMCCell objects added via add_cell() to manage per-cell " + "material assignments.") def add_material_override(self, key, overrides=None): """Add a material override to the universe. @@ -201,7 +152,10 @@ class DAGMCUniverse(openmc.UniverseBase): if key not in self.cells: raise ValueError(f"Cell ID '{key}' not found in DAGMC universe") - self._material_overrides[key] = overrides + if len(overrides) == 1: + self.cells[key].fill = overrides[0] + else: + self.cells[key].fill = list(overrides) @property def auto_geom_ids(self): @@ -290,12 +244,6 @@ class DAGMCUniverse(openmc.UniverseBase): memo.add(self) - # Ensure that the material overrides are up-to-date - for cell in self.cells.values(): - if cell.fill is None: - continue - self.add_material_override(cell, cell.fill) - # Set xml element values dagmc_element = ET.Element('dagmc_universe') dagmc_element.set('id', str(self.id)) @@ -307,17 +255,10 @@ class DAGMCUniverse(openmc.UniverseBase): if self.auto_mat_ids: dagmc_element.set('auto_mat_ids', 'true') dagmc_element.set('filename', str(self.filename)) - if self._material_overrides: - mat_element = ET.Element('material_overrides') - for key in self._material_overrides: - cell_overrides = ET.Element('cell_override') - cell_overrides.set("id", str(key)) - material_element = ET.Element('material_ids') - material_element.text = ' '.join( - str(t.id) for t in self._material_overrides[key]) - cell_overrides.append(material_element) - mat_element.append(cell_overrides) - dagmc_element.append(mat_element) + if self.cells: + for cell in self.cells.values(): + cell_element = cell.create_xml_subelement(xml_element, memo) + dagmc_element.append(cell_element) xml_element.append(dagmc_element) def bounding_region( @@ -442,7 +383,7 @@ class DAGMCUniverse(openmc.UniverseBase): return out @classmethod - def from_xml_element(cls, elem, mats = None): + def from_xml_element(cls, elem, mats=None): """Generate DAGMC universe from XML element Parameters @@ -471,21 +412,56 @@ class DAGMCUniverse(openmc.UniverseBase): out.auto_geom_ids = bool(get_text(elem, "auto_geom_ids")) out.auto_mat_ids = bool(get_text(elem, "auto_mat_ids")) - el_mat_override = elem.find('material_overrides') - if el_mat_override is not None: - if mats is None: - raise ValueError("Material overrides found in DAGMC universe " - "but no materials were provided to populate " - "the mapping.") - out._material_overrides = {} - for elem in el_mat_override.findall('cell_override'): - cell_id = int(get_text(elem, 'id')) - mat_ids = get_elem_list(elem, "material_ids", str) or [] - mat_objs = [mats[mat_id] for mat_id in mat_ids] - out._material_overrides[cell_id] = mat_objs + has_overrides = elem.find('material_overrides') is not None + has_cells = elem.find('cell') is not None + + if has_overrides and has_cells: + raise ValueError( + "DAGMCUniverse cannot specify both and " + " sub-elements. Use elements only.") + + if has_overrides: + warnings.warn( + "DAGMCUniverse is deprecated and will be " + "removed in a future version. Use nested elements " + "instead.", DeprecationWarning, stacklevel=2) + out._parse_legacy_material_overrides(elem, mats) + elif has_cells: + out._parse_cell_overrides(elem, mats) return out + def _parse_legacy_material_overrides(self, elem, mats): + """Parse the deprecated XML format and populate + the universe with equivalent DAGMCCell objects.""" + if mats is None: + raise ValueError( + "DAGMC material overrides found but no materials were " + "provided to populate the mapping.") + mo_elem = elem.find('material_overrides') + for co_elem in mo_elem.findall('cell_override'): + cell_id = int(get_text(co_elem, 'id')) + mat_ids = co_elem.find('material_ids').text.split() + fill_objs = [mats[mid] for mid in mat_ids] + fill = fill_objs[0] if len(fill_objs) == 1 else fill_objs + if cell_id in self.cells: + raise ValueError( + f"Duplicate DAGMC cell override specified for cell {cell_id}.") + self.add_cell(DAGMCCell(cell_id=cell_id, fill=fill)) + + def _parse_cell_overrides(self, elem, mats): + if mats is None: + raise ValueError("DAGMC cell overrides found in DAGMC universe but " + "no materials were provided to populate the " + "mapping.") + + for cell_elem in elem.findall('cell'): + cell_id = int(get_text(cell_elem, 'id')) + if cell_id in self.cells: + raise ValueError( + f"Duplicate DAGMC cell override specified for cell {cell_id}.") + DAGMCCell.from_xml_element(cell_elem, mats, self) + def _partial_deepcopy(self): """Clone all of the openmc.DAGMCUniverse object's attributes except for its cells, as they are copied within the clone function. This should @@ -565,7 +541,13 @@ class DAGMCUniverse(openmc.UniverseBase): fill = [mats_per_id[mat.id] for mat in dag_cell.fill if mat] else: fill = mats_per_id[dag_cell.fill.id] if dag_cell.fill else None - self.add_cell(openmc.DAGMCCell(cell_id=dag_cell_id, fill=fill)) + name = dag_cell.name + if dag_cell_id in self._cells: + self._cells[dag_cell_id].name = name + self._cells[dag_cell_id].fill = fill + else: + self.add_cell( + openmc.DAGMCCell(cell_id=dag_cell_id, name=name, fill=fill)) @add_plot_params def plot(self, *args, **kwargs): @@ -594,6 +576,14 @@ class DAGMCCell(openmc.Cell): DAG_parent_universe : int The parent universe of the cell. + Notes + ----- + DAGMC geometries are composed of triangulated surfaces, which means cell + volumes can in principle be computed exactly (e.g. via mesh-based + integration). Manually specifying :attr:`volume` overrides any such + calculation and may introduce inconsistencies if the value does not + accurately reflect the true geometric volume. + """ def __init__(self, cell_id=None, name='', fill=None): super().__init__(cell_id, name, fill, None) @@ -625,8 +615,62 @@ class DAGMCCell(openmc.Cell): raise TypeError("plot is not available for DAGMC cells.") def create_xml_subelement(self, xml_element, memo=None): - raise TypeError("create_xml_subelement is not available for DAGMC cells.") + if self.fill_type not in ('void', 'material', 'distribmat'): + raise TypeError("DAGMC cell overrides currently only support " + "material fills.") + if self.temperature is not None and self.fill_type not in ( + 'material', 'distribmat' + ): + raise TypeError("DAGMC cell temperature overrides require a " + "material fill.") + if self.density is not None and self.fill_type not in ('material', 'distribmat'): + raise TypeError("DAGMC cell density overrides require a " + "material fill.") + if any(getattr(self, attr) is not None for attr in ('translation', 'rotation')): + raise TypeError("DAGMC cell overrides do not support translation " + "or rotation.") + return super().create_xml_subelement(xml_element, memo) @classmethod - def from_xml_element(cls, elem, surfaces, materials, get_universe): - raise TypeError("from_xml_element is not available for DAGMC cells.") + def from_xml_element(cls, elem, mats, universe): + """Generate a DAGMCCell from an XML override element. + + Parameters + ---------- + elem : lxml.etree._Element + `` element containing a DAGMC cell property override + mats : dict + Dictionary mapping material ID strings to :class:`openmc.Material` + instances + universe : DAGMCUniverse + Universe to add the parsed cell to. + + Returns + ------- + DAGMCCell + DAGMCCell instance + """ + if not isinstance(universe, DAGMCUniverse): + raise TypeError( + f"universe must be a DAGMCUniverse instance, " + f"got {type(universe).__name__}.") + + cell_id = int(get_text(elem, 'id')) + + # Validate attributes that are unsupported for DAGMC cell overrides + for tag in ('region', 'fill', 'universe'): + if get_text(elem, tag) is not None: + raise ValueError( + f"DAGMC cell {cell_id} override cannot specify '{tag}'.") + for tag in ('translation', 'rotation'): + if get_text(elem, tag) is not None: + raise ValueError( + f"DAGMC cell {cell_id} override does not support " + f"'{tag}'.") + if get_elem_list(elem, 'material', str) is None: + raise ValueError( + f"DAGMC cell {cell_id} must specify a material override.") + + return super().from_xml_element( + elem, surfaces={}, materials=mats, + get_universe=lambda _: universe) diff --git a/openmc/model/model.py b/openmc/model/model.py index c9a24b8b3..11c26d4e7 100644 --- a/openmc/model/model.py +++ b/openmc/model/model.py @@ -467,6 +467,8 @@ class Model: This method iterates over all DAGMC universes in the geometry and synchronizes their cells with the current material assignments. Requires that the model has been initialized via :meth:`Model.init_lib`. + Synchronized DAGMC cells can then be edited and exported as nested + `` overrides inside each `` element. .. versionadded:: 0.15.1 diff --git a/src/cell.cpp b/src/cell.cpp index ebe28c3d2..06bc29a66 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -340,6 +340,63 @@ void Cell::to_hdf5(hid_t cell_group) const close_group(group); } +//============================================================================== +// XML parsing helpers for nodes +//============================================================================== + +vector parse_cell_material_xml(pugi::xml_node node, int32_t cell_id) +{ + vector mats { + get_node_array(node, "material", true)}; + if (mats.empty()) { + fatal_error(fmt::format( + "An empty material element was specified for cell {}", cell_id)); + } + vector material; + material.reserve(mats.size()); + for (const auto& mat : mats) { + if (mat == "void") { + material.push_back(MATERIAL_VOID); + } else { + material.push_back(std::stoi(mat)); + } + } + return material; +} + +vector parse_cell_temperature_xml(pugi::xml_node node, int32_t cell_id) +{ + auto temperatures = get_node_array(node, "temperature"); + if (temperatures.empty()) { + fatal_error(fmt::format( + "An empty temperature element was specified for cell {}", cell_id)); + } + for (auto T : temperatures) { + if (T < 0) { + fatal_error(fmt::format( + "Cell {} was specified with a negative temperature", cell_id)); + } + } + return temperatures; +} + +vector parse_cell_density_xml(pugi::xml_node node, int32_t cell_id) +{ + auto densities = get_node_array(node, "density"); + if (densities.empty()) { + fatal_error(fmt::format( + "An empty density element was specified for cell {}", cell_id)); + } + for (auto rho : densities) { + if (rho <= 0) { + fatal_error(fmt::format( + "Cell {} was specified with a density less than or equal to zero", + cell_id)); + } + } + return densities; +} + //============================================================================== // CSGCell implementation //============================================================================== @@ -390,26 +447,12 @@ CSGCell::CSGCell(pugi::xml_node cell_node) // universe), more than one material (distribmats), and some materials may // be "void". if (material_present) { - vector mats { - get_node_array(cell_node, "material", true)}; - if (mats.size() > 0) { - material_.reserve(mats.size()); - for (std::string mat : mats) { - if (mat.compare("void") == 0) { - material_.push_back(MATERIAL_VOID); - } else { - material_.push_back(std::stoi(mat)); - } - } - } else { - fatal_error(fmt::format( - "An empty material element was specified for cell {}", id_)); - } + material_ = parse_cell_material_xml(cell_node, id_); } // Read the temperature element which may be distributed like materials. if (check_for_node(cell_node, "temperature")) { - sqrtkT_ = get_node_array(cell_node, "temperature"); + sqrtkT_ = parse_cell_temperature_xml(cell_node, id_); sqrtkT_.shrink_to_fit(); // Make sure this is a material-filled cell. @@ -420,14 +463,6 @@ CSGCell::CSGCell(pugi::xml_node cell_node) id_)); } - // Make sure all temperatures are non-negative. - for (auto T : sqrtkT_) { - if (T < 0) { - fatal_error(fmt::format( - "Cell {} was specified with a negative temperature", id_)); - } - } - // Convert to sqrt(k*T). for (auto& T : sqrtkT_) { T = std::sqrt(K_BOLTZMANN * T); @@ -440,7 +475,7 @@ CSGCell::CSGCell(pugi::xml_node cell_node) // Note: calculating the actual density multiplier is deferred until materials // are finalized. density_mult_ contains the true density in the meantime. if (check_for_node(cell_node, "density")) { - density_mult_ = get_node_array(cell_node, "density"); + density_mult_ = parse_cell_density_xml(cell_node, id_); density_mult_.shrink_to_fit(); // Make sure this is a material-filled cell. @@ -461,15 +496,6 @@ CSGCell::CSGCell(pugi::xml_node cell_node) id_)); } } - - // Make sure all densities are non-negative and greater than zero. - for (auto rho : density_mult_) { - if (rho <= 0) { - fatal_error(fmt::format( - "Cell {} was specified with a density less than or equal to zero", - id_)); - } - } } // Read the region specification. diff --git a/src/dagmc.cpp b/src/dagmc.cpp index 571182fa6..a5cc71b52 100644 --- a/src/dagmc.cpp +++ b/src/dagmc.cpp @@ -50,6 +50,10 @@ namespace openmc { DAGUniverse::DAGUniverse(pugi::xml_node node) { + MaterialOverrides material_overrides; + TemperatureOverrides temperature_overrides; + DensityOverrides density_overrides; + if (check_for_node(node, "id")) { id_ = std::stoi(get_node_value(node, "id")); } else { @@ -76,24 +80,77 @@ DAGUniverse::DAGUniverse(pugi::xml_node node) adjust_material_ids_ = get_node_value_bool(node, "auto_mat_ids"); } - // get material assignment overloading - if (check_for_node(node, "material_overrides")) { - auto mat_node = node.child("material_overrides"); - // loop over all subelements (each subelement corresponds to a material) - for (pugi::xml_node cell_node : mat_node.children("cell_override")) { - // Store assignment reference name - int32_t ref_assignment = std::stoi(get_node_value(cell_node, "id")); + // Get material assignment overrides from nested DAGMC cell elements. + if (node.child("cell")) { + for (pugi::xml_node cell_node : node.children("cell")) { + if (!check_for_node(cell_node, "id")) { + fatal_error( + "Must specify id for each DAGMC cell override in ."); + } - // Get mat name for each assignement instances - vector instance_mats = - get_node_array(cell_node, "material_ids"); + int32_t cell_id = std::stoi(get_node_value(cell_node, "id")); - // Store mat name for each instances - material_overrides_.emplace(ref_assignment, instance_mats); + if (check_for_node(cell_node, "region")) { + fatal_error(fmt::format( + "DAGMC cell {} override cannot specify a region.", cell_id)); + } + if (check_for_node(cell_node, "fill")) { + fatal_error(fmt::format( + "DAGMC cell {} override currently only supports material fills.", + cell_id)); + } + if (check_for_node(cell_node, "universe")) { + fatal_error(fmt::format( + "DAGMC cell {} override cannot specify a universe.", cell_id)); + } + if (check_for_node(cell_node, "translation") || + check_for_node(cell_node, "rotation")) { + fatal_error(fmt::format( + "DAGMC cell {} override does not support translation or rotation.", + cell_id)); + } + if (!check_for_node(cell_node, "material")) { + fatal_error(fmt::format( + "DAGMC cell {} override must specify material.", cell_id)); + } + + auto inserted = material_overrides.emplace( + cell_id, parse_cell_material_xml(cell_node, cell_id)); + if (!inserted.second) { + fatal_error(fmt::format( + "Duplicate DAGMC cell override specified for cell {}", cell_id)); + } + + if (check_for_node(cell_node, "temperature")) { + temperature_overrides.emplace( + cell_id, parse_cell_temperature_xml(cell_node, cell_id)); + } + + if (check_for_node(cell_node, "density")) { + density_overrides.emplace( + cell_id, parse_cell_density_xml(cell_node, cell_id)); + } + } + } else if (check_for_node(node, "material_overrides")) { + if (node.child("cell")) { + fatal_error("DAGMCUniverse cannot specify both and " + " sub-elements. Use elements only."); + } + warning("DAGMCUniverse is deprecated. Use nested " + " elements under instead."); + for (pugi::xml_node co : + node.child("material_overrides").children("cell_override")) { + int32_t cell_id = std::stoi(get_node_value(co, "id")); + std::istringstream iss(co.child("material_ids").text().get()); + vector mats; + for (std::string s; iss >> s;) { + mats.push_back(s == "void" ? MATERIAL_VOID : std::stoi(s)); + } + material_overrides.emplace(cell_id, mats); } } - initialize(); + initialize(material_overrides, temperature_overrides, density_overrides); } DAGUniverse::DAGUniverse( @@ -110,9 +167,12 @@ DAGUniverse::DAGUniverse(std::shared_ptr dagmc_ptr, : dagmc_instance_(dagmc_ptr), filename_(filename), adjust_geometry_ids_(auto_geom_ids), adjust_material_ids_(auto_mat_ids) { + MaterialOverrides material_overrides; + TemperatureOverrides temperature_overrides; + DensityOverrides density_overrides; set_id(); init_metadata(); - init_geometry(); + init_geometry(material_overrides, temperature_overrides, density_overrides); } void DAGUniverse::set_id() @@ -130,6 +190,15 @@ void DAGUniverse::set_id() } void DAGUniverse::initialize() +{ + MaterialOverrides material_overrides; + TemperatureOverrides temperature_overrides; + initialize(material_overrides, temperature_overrides); +} + +void DAGUniverse::initialize(const MaterialOverrides& material_overrides, + const TemperatureOverrides& temperature_overrides, + const DensityOverrides& density_overrides) { #ifdef OPENMC_UWUW_ENABLED // read uwuw materials from the .h5m file if present @@ -140,7 +209,7 @@ void DAGUniverse::initialize() init_metadata(); - init_geometry(); + init_geometry(material_overrides, temperature_overrides, density_overrides); } void DAGUniverse::init_dagmc() @@ -176,7 +245,9 @@ void DAGUniverse::init_metadata() MB_CHK_ERR_CONT(rval); } -void DAGUniverse::init_geometry() +void DAGUniverse::init_geometry(const MaterialOverrides& material_overrides, + const TemperatureOverrides& temperature_overrides, + const DensityOverrides& density_overrides) { moab::ErrorCode rval; @@ -202,6 +273,9 @@ void DAGUniverse::init_geometry() : dagmc_instance_->id_by_index(3, c->dag_index()); c->universe_ = this->id_; c->fill_ = C_NONE; // no fill, single universe + if (dagmc_instance_->is_implicit_complement(vol_handle)) { + c->name_ = "implicit complement"; + } auto in_map = model::cell_map.find(c->id_); if (in_map == model::cell_map.end()) { @@ -230,16 +304,68 @@ void DAGUniverse::init_geometry() if (mat_str == "graveyard") { graveyard = vol_handle; } - // material void checks - if (mat_str == "void" || mat_str == "vacuum" || mat_str == "graveyard") { + if (material_overrides.count(c->id_)) { + override_assign_material(c, material_overrides); + } else if (mat_str == "void" || mat_str == "vacuum" || + mat_str == "graveyard") { c->material_.push_back(MATERIAL_VOID); + } else if (uses_uwuw()) { + uwuw_assign_material(vol_handle, c); } else { - if (material_overrides_.count(c->id_)) { - override_assign_material(c); - } else if (uses_uwuw()) { - uwuw_assign_material(vol_handle, c); - } else { - legacy_assign_material(mat_str, c); + legacy_assign_material(mat_str, c); + } + + if (temperature_overrides.count(c->id_)) { + if (c->material_.empty() || c->material_[0] == MATERIAL_VOID) { + fatal_error(fmt::format("DAGMC cell {} was specified with a " + "temperature but no non-void material.", + c->id_)); + } + + c->sqrtkT_.clear(); + const auto& temp_overrides = temperature_overrides.at(c->id_); + c->sqrtkT_.reserve(temp_overrides.size()); + for (auto T : temp_overrides) { + c->sqrtkT_.push_back(std::sqrt(K_BOLTZMANN * T)); + } + + if (settings::verbosity >= 10) { + std::stringstream override_values; + for (size_t i = 0; i < temp_overrides.size(); ++i) { + if (i > 0) { + override_values << " "; + } + override_values << temp_overrides[i]; + } + auto msg = fmt::format("Overriding DAGMC cell {} property " + "'temperature [K]' with value(s): {}", + c->id_, override_values.str()); + write_message(msg, 10); + } + } + + if (density_overrides.count(c->id_)) { + if (c->material_.empty() || c->material_[0] == MATERIAL_VOID) { + fatal_error(fmt::format("DAGMC cell {} was specified with a density " + "but no non-void material.", + c->id_)); + } + // density_mult_ holds the true density until materials are finalized, + // at which point it is converted to a proper multiplier (same as CSG). + c->density_mult_ = density_overrides.at(c->id_); + + if (settings::verbosity >= 10) { + const auto& dens = density_overrides.at(c->id_); + std::stringstream override_values; + for (size_t i = 0; i < dens.size(); ++i) { + if (i > 0) + override_values << " "; + override_values << dens[i]; + } + write_message(fmt::format("Overriding DAGMC cell {} property " + "'density [g/cm³]' with value(s): {}", + c->id_, override_values.str()), + 10); } } @@ -252,18 +378,21 @@ void DAGUniverse::init_geometry() continue; } - // assign cell temperature - const auto& mat = model::materials[model::material_map.at(c->material_[0])]; - if (dagmc_instance_->has_prop(vol_handle, "temp")) { - rval = dagmc_instance_->prop_value(vol_handle, "temp", temp_value); - MB_CHK_ERR_CONT(rval); - double temp = std::stod(temp_value); - c->sqrtkT_.push_back(std::sqrt(K_BOLTZMANN * temp)); - } else if (mat->temperature() > 0.0) { - c->sqrtkT_.push_back(std::sqrt(K_BOLTZMANN * mat->temperature())); - } else { - c->sqrtkT_.push_back( - std::sqrt(K_BOLTZMANN * settings::temperature_default)); + // assign cell temperature if not explicitly overridden + if (c->sqrtkT_.empty()) { + const auto& mat = + model::materials[model::material_map.at(c->material_[0])]; + if (dagmc_instance_->has_prop(vol_handle, "temp")) { + rval = dagmc_instance_->prop_value(vol_handle, "temp", temp_value); + MB_CHK_ERR_CONT(rval); + double temp = std::stod(temp_value); + c->sqrtkT_.push_back(std::sqrt(K_BOLTZMANN * temp)); + } else if (mat->temperature() > 0.0) { + c->sqrtkT_.push_back(std::sqrt(K_BOLTZMANN * mat->temperature())); + } else { + c->sqrtkT_.push_back( + std::sqrt(K_BOLTZMANN * settings::temperature_default)); + } } model::cells.emplace_back(std::move(c)); @@ -630,7 +759,8 @@ void DAGUniverse::uwuw_assign_material( #endif // OPENMC_UWUW_ENABLED } -void DAGUniverse::override_assign_material(std::unique_ptr& c) const +void DAGUniverse::override_assign_material(std::unique_ptr& c, + const MaterialOverrides& material_overrides) const { // if Cell ID matches an override key, use it to override the material // assignment else if UWUW is used, get the material assignment from the DAGMC @@ -638,17 +768,30 @@ void DAGUniverse::override_assign_material(std::unique_ptr& c) const // Notify User that an override is being applied on a DAGMCCell write_message(fmt::format("Applying override for DAGMCCell {}", c->id_), 8); + const auto& mat_overrides = material_overrides.at(c->id_); if (settings::verbosity >= 10) { - auto msg = fmt::format("Assigning DAGMC cell {} material(s) based on " - "override information (see input XML).", - c->id_); + std::stringstream override_values; + for (size_t i = 0; i < mat_overrides.size(); ++i) { + if (i > 0) { + override_values << " "; + } + if (mat_overrides[i] == MATERIAL_VOID) { + override_values << "void"; + } else { + override_values << mat_overrides[i]; + } + } + auto msg = fmt::format("Overriding DAGMC cell {} property 'material' " + "with value(s): {}", + c->id_, override_values.str()); write_message(msg, 10); } // Override the material assignment for each cell instance using the legacy // assignement - for (auto mat_id : material_overrides_.at(c->id_)) { - if (model::material_map.find(mat_id) == model::material_map.end()) { + for (auto mat_id : mat_overrides) { + if (mat_id != MATERIAL_VOID && + model::material_map.find(mat_id) == model::material_map.end()) { fatal_error(fmt::format( "Material with ID '{}' not found for DAGMC cell {}", mat_id, c->id_)); } diff --git a/tests/unit_tests/dagmc/test_model.py b/tests/unit_tests/dagmc/test_model.py index 0917f5b23..5889fa542 100644 --- a/tests/unit_tests/dagmc/test_model.py +++ b/tests/unit_tests/dagmc/test_model.py @@ -70,28 +70,20 @@ def model(request): openmc.reset_auto_ids() -def test_dagmc_replace_material_assignment(model): - mats = {} - - mats["foo"] = openmc.Material(name="foo") - mats["foo"].add_nuclide("H1", 2.0) - mats["foo"].add_element("O", 1.0) - mats["foo"].set_density("g/cm3", 1.0) - mats["foo"].add_s_alpha_beta("c_H_in_H2O") - +def test_dagmc_sync_cell_names(model): + dag_univ = None for univ in model.geometry.get_all_universes().values(): - if not isinstance(univ, openmc.DAGMCUniverse): + if isinstance(univ, openmc.DAGMCUniverse): + dag_univ = univ break - cells_with_41 = [] - for cell in univ.cells.values(): - if cell.fill is None: - continue - if cell.fill.name == "41": - cells_with_41.append(cell.id) - univ.replace_material_assignment("41", mats["foo"]) - for cell_id in cells_with_41: - assert univ.cells[cell_id] == mats["foo"] + assert dag_univ is not None + + for cell_id, cell in dag_univ.cells.items(): + assert cell.name == openmc.lib.cells[cell_id].name + + assert any(cell.name == "implicit complement" + for cell in dag_univ.cells.values()) def test_dagmc_add_material_override_with_id(model): @@ -114,7 +106,7 @@ def test_dagmc_add_material_override_with_id(model): cells_with_41.append(cell.id) univ.add_material_override(cell.id, mats["foo"]) for cell_id in cells_with_41: - assert univ.cells[cell_id] == mats["foo"] + assert univ.cells[cell_id].fill == mats["foo"] def test_dagmc_add_material_override_with_cell(model): @@ -137,7 +129,7 @@ def test_dagmc_add_material_override_with_cell(model): cells_with_41.append(cell.id) univ.add_material_override(cell, mats["foo"]) for cell_id in cells_with_41: - assert univ.cells[cell_id] == mats["foo"] + assert univ.cells[cell_id].fill == mats["foo"] def test_model_differentiate_depletable_with_dagmc(model, run_in_tmpdir): @@ -174,57 +166,20 @@ def test_model_differentiate_with_dagmc(model): assert len(model.materials) == 4*2 + 4 -def test_bad_override_cell_id(model): - for univ in model.geometry.get_all_universes().values(): - if isinstance(univ, openmc.DAGMCUniverse): - break - with pytest.raises(ValueError, match="Cell ID '1' not found in DAGMC universe"): - univ.material_overrides = {1: model.materials[0]} - - -def test_bad_override_type(model): - not_a_dag_cell = openmc.Cell() - for univ in model.geometry.get_all_universes().values(): - if isinstance(univ, openmc.DAGMCUniverse): - break - with pytest.raises(ValueError, match="Unrecognized key type. Must be an integer or openmc.DAGMCCell object"): - univ.material_overrides = {not_a_dag_cell: model.materials[0]} - - -def test_bad_replacement_mat_name(model): - for univ in model.geometry.get_all_universes().values(): - if isinstance(univ, openmc.DAGMCUniverse): - break - with pytest.raises(ValueError, match="No material with name 'not_a_mat' found in the DAGMC universe"): - univ.replace_material_assignment("not_a_mat", model.materials[0]) - - def test_dagmc_xml(model): - # Set the environment - mats = {} - mats["no-void fuel"] = openmc.Material(1, name="no-void fuel") - mats["no-void fuel"].add_nuclide("U235", 0.03) - mats["no-void fuel"].add_nuclide("U238", 0.97) - mats["no-void fuel"].add_nuclide("O16", 2.0) - mats["no-void fuel"].set_density("g/cm3", 10.0) - - mats[5] = openmc.Material(name="41") - mats[5].add_nuclide("H1", 2.0) - mats[5].add_element("O", 1.0) - mats[5].set_density("g/cm3", 1.0) - mats[5].add_s_alpha_beta("c_H_in_H2O") + override_mat = openmc.Material(name="41") + override_mat.add_nuclide("H1", 2.0) + override_mat.add_element("O", 1.0) + override_mat.set_density("g/cm3", 1.0) + override_mat.add_s_alpha_beta("c_H_in_H2O") + model.materials.append(override_mat) for univ in model.geometry.get_all_universes().values(): if isinstance(univ, openmc.DAGMCUniverse): dag_univ = univ break - for k, v in mats.items(): - if isinstance(k, int): - dag_univ.add_material_override(k, v) - model.materials.append(v) - elif isinstance(k, str): - dag_univ.replace_material_assignment(k, v) + dag_univ.add_material_override(5, override_mat) # Tesing the XML subelement generation root = ET.Element('dagmc_universe') @@ -236,12 +191,24 @@ def test_dagmc_xml(model): assert dagmc_ele.get('filename') == str(dag_univ.filename) assert dagmc_ele.get('auto_geom_ids') == str(dag_univ.auto_geom_ids).lower() - override_eles = dagmc_ele.find('material_overrides').findall('cell_override') - assert len(override_eles) == 4 + assert dagmc_ele.find('material_overrides') is None - for i, override_ele in enumerate(override_eles): - cell_id = override_ele.get('id') - assert dag_univ.material_overrides[int(cell_id)][0].id == int(override_ele.find('material_ids').text) + override_elements = dagmc_ele.findall('cell') + assert len(override_elements) == len(dag_univ.cells) + xml_cells = {int(elem.get('id')): elem for elem in override_elements} + for cell_id, cell in dag_univ.cells.items(): + assert cell_id in xml_cells + xml_cell = xml_cells[cell_id] + if cell.fill_type == 'void': + assert xml_cell.get('material') == 'void' + elif cell.fill_type == 'material': + assert xml_cell.get('material') == str(cell.fill.id) + elif cell.fill_type == 'distribmat': + mat_list = xml_cell.find('material').text.split() + expected = ["void" if m is None else str(m.id) for m in cell.fill] + assert mat_list == expected + else: + pytest.fail(f"Unexpected DAGMC cell fill type: {cell.fill_type}") model.export_to_model_xml() @@ -252,7 +219,147 @@ def test_dagmc_xml(model): xml_dagmc_univ = univ break - assert xml_dagmc_univ._material_overrides.keys() == dag_univ._material_overrides.keys() + assert xml_dagmc_univ.cells.keys() == dag_univ.cells.keys() - for xml_mats, model_mats in zip(xml_dagmc_univ._material_overrides.values(), dag_univ._material_overrides.values()): - assert all([xml_mat.id == orig_mat.id for xml_mat, orig_mat in zip(xml_mats, model_mats)]) + for cell_id, cell in dag_univ.cells.items(): + xml_cell = xml_dagmc_univ.cells[cell_id] + assert xml_cell.fill_type == cell.fill_type + if cell.fill_type == 'void': + assert xml_cell.fill is None + elif cell.fill_type == 'material': + assert xml_cell.fill.id == cell.fill.id + elif cell.fill_type == 'distribmat': + xml_ids = [m.id if m is not None else None for m in xml_cell.fill] + model_ids = [m.id if m is not None else None for m in cell.fill] + assert xml_ids == model_ids + else: + pytest.fail(f"Unexpected DAGMC cell fill type: {cell.fill_type}") + + +def test_dagmc_xml_reject_fill_override(): + mats = {'1': openmc.Material(1), 'void': None} + elem = ET.fromstring( + '' + '' + '' + ) + with pytest.raises(ValueError, match="cannot specify 'fill'"): + openmc.DAGMCUniverse.from_xml_element(elem, mats) + + +def test_dagmc_xml_reject_region_override(): + mats = {'1': openmc.Material(1), 'void': None} + elem = ET.fromstring( + '' + '' + '' + ) + with pytest.raises(ValueError, match="cannot specify 'region'"): + openmc.DAGMCUniverse.from_xml_element(elem, mats) + + +def _legacy_xml(cell_overrides): + """Helper to build a with old-format .""" + inner = ''.join( + f'{mids}' + for cid, mids in cell_overrides.items() + ) + return ET.fromstring( + f'' + f'{inner}' + f'' + ) + + +def test_dagmc_xml_legacy_single_material_compat(): + mat = openmc.Material(1) + mats = {'1': mat, 'void': None} + elem = _legacy_xml({3: '1'}) + with pytest.warns(DeprecationWarning, match="deprecated"): + univ = openmc.DAGMCUniverse.from_xml_element(elem, mats) + assert 3 in univ.cells + assert univ.cells[3].fill is mat + + +def test_dagmc_xml_legacy_distribmat_compat(): + mat1, mat2 = openmc.Material(2), openmc.Material(3) + mats = {'2': mat1, '3': mat2, 'void': None} + elem = _legacy_xml({5: '2 3'}) + with pytest.warns(DeprecationWarning): + univ = openmc.DAGMCUniverse.from_xml_element(elem, mats) + assert univ.cells[5].fill_type == 'distribmat' + assert list(univ.cells[5].fill) == [mat1, mat2] + + +def test_dagmc_xml_legacy_void_compat(): + mats = {'void': None} + elem = _legacy_xml({7: 'void'}) + with pytest.warns(DeprecationWarning): + univ = openmc.DAGMCUniverse.from_xml_element(elem, mats) + assert univ.cells[7].fill_type == 'void' + + +def test_dagmc_xml_legacy_both_raises(): + mat = openmc.Material(1) + mats = {'1': mat, 'void': None} + elem = ET.fromstring( + '' + '' + '1' + '' + '' + '' + ) + with pytest.raises(ValueError, match="both"): + openmc.DAGMCUniverse.from_xml_element(elem, mats) + + +def test_dagmc_xml_legacy_deprecation_warning(): + mats = {'1': openmc.Material(1), 'void': None} + elem = _legacy_xml({3: '1'}) + with pytest.warns(DeprecationWarning): + openmc.DAGMCUniverse.from_xml_element(elem, mats) + + +def test_dagmc_xml_legacy_roundtrip(): + """Old-format XML loads correctly and re-exports using the new format.""" + mat = openmc.Material(1) + mats = {'1': mat, 'void': None} + elem = _legacy_xml({3: '1'}) + with pytest.warns(DeprecationWarning): + univ = openmc.DAGMCUniverse.from_xml_element(elem, mats) + + root = ET.Element('geometry') + univ.create_xml_subelement(root) + dagmc_elem = root.find('dagmc_universe') + + assert dagmc_elem.find('material_overrides') is None + cell_elems = dagmc_elem.findall('cell') + assert len(cell_elems) == 1 + assert int(cell_elems[0].get('id')) == 3 + assert cell_elems[0].get('material') == '1' + + +def test_dagmc_xml_temperature_roundtrip(): + mat = openmc.Material(1) + mats = {'1': mat, 'void': None} + + elem = ET.fromstring( + '' + '' + '' + ) + + dag_univ = openmc.DAGMCUniverse.from_xml_element(elem, mats) + assert dag_univ.cells[7].fill.id == 1 + assert dag_univ.cells[7].temperature == pytest.approx(825.0) + + root = ET.Element('geometry') + dag_univ.create_xml_subelement(root) + dagmc_elem = root.find('dagmc_universe') + xml_cell = dagmc_elem.find('cell') + assert xml_cell.get('temperature') == '825.0' + + dag_univ_roundtrip = openmc.DAGMCUniverse.from_xml_element(dagmc_elem, mats) + assert dag_univ_roundtrip.cells[7].fill.id == 1 + assert dag_univ_roundtrip.cells[7].temperature == pytest.approx(825.0) diff --git a/tests/unit_tests/test_cell.py b/tests/unit_tests/test_cell.py index 2baa59f1b..8d39c071b 100644 --- a/tests/unit_tests/test_cell.py +++ b/tests/unit_tests/test_cell.py @@ -374,6 +374,49 @@ def test_rotation_from_xml(rotation): np.testing.assert_allclose(new_cell.rotation, cell.rotation) +def test_dagmccell_from_xml_element(): + """DAGMCCell.from_xml_element parses material, temperature, density, + and volume; rejects unsupported attributes.""" + mat = openmc.Material(1) + mat.add_nuclide('U235', 1.0) + mats = {'1': mat} + # In practice, from_xml_element is always called from DAGMCUniverse + # during XML parsing. A placeholder universe is used here so the test + # can exercise the method directly without a real DAGMC model file. + placeholder_univ = openmc.DAGMCUniverse('model.h5m') + + # material + temperature + density round-trip + xml = '' + cell = openmc.DAGMCCell.from_xml_element(ET.fromstring(xml), mats, placeholder_univ) + assert cell.id == 5 + assert cell.name == 'fuel' + assert cell.fill is mat + assert cell.density == 10.5 + assert cell.temperature == 900.0 + + # volume round-trip + xml = '' + cell = openmc.DAGMCCell.from_xml_element(ET.fromstring(xml), mats, placeholder_univ) + assert cell.volume == 42.0 + + # forbidden: region, fill, universe + for tag, val in [('region', '-1'), ('fill', '2'), ('universe', '0')]: + xml = f'' + with pytest.raises(ValueError, match=tag): + openmc.DAGMCCell.from_xml_element(ET.fromstring(xml), mats, placeholder_univ) + + # forbidden: translation, rotation + for tag, val in [('translation', '1 0 0'), ('rotation', '0 0 90')]: + xml = f'' + with pytest.raises(ValueError, match=tag): + openmc.DAGMCCell.from_xml_element(ET.fromstring(xml), mats, placeholder_univ) + + # missing material raises + xml = '' + with pytest.raises(ValueError, match='material'): + openmc.DAGMCCell.from_xml_element(ET.fromstring(xml), mats, placeholder_univ) + + def test_plot(run_in_tmpdir): zcyl = openmc.ZCylinder() c = openmc.Cell(region=-zcyl) From dfb6c5699c6fb966ca4dfd396d3ccb22f4e28f60 Mon Sep 17 00:00:00 2001 From: Ahnaf Tahmid Chowdhury Date: Fri, 22 May 2026 02:33:56 +0600 Subject: [PATCH 10/10] Modernize CMake packaging (#3653) Co-authored-by: Paul Romano --- CMakeLists.txt | 24 +++++++++++++++------ cmake/OpenMCConfig.cmake.in | 34 +++++++++++++++++++----------- cmake/OpenMCConfigVersion.cmake.in | 11 ---------- 3 files changed, 40 insertions(+), 29 deletions(-) delete mode 100644 cmake/OpenMCConfigVersion.cmake.in diff --git a/CMakeLists.txt b/CMakeLists.txt index 9fe133a22..d5af01610 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -611,9 +611,7 @@ add_custom_command(TARGET libopenmc POST_BUILD #=============================================================================== # Install executable, scripts, manpage, license #=============================================================================== - -configure_file(cmake/OpenMCConfig.cmake.in "${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/OpenMCConfig.cmake" @ONLY) -configure_file(cmake/OpenMCConfigVersion.cmake.in "${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/OpenMCConfigVersion.cmake" @ONLY) +include(CMakePackageConfigHelpers) set(INSTALL_CONFIGDIR ${CMAKE_INSTALL_LIBDIR}/cmake/OpenMC) install(TARGETS openmc libopenmc @@ -627,10 +625,24 @@ install(EXPORT openmc-targets NAMESPACE OpenMC:: DESTINATION ${INSTALL_CONFIGDIR}) +configure_package_config_file( + "cmake/OpenMCConfig.cmake.in" + "${CMAKE_BINARY_DIR}/${CMAKE_FILES_DIRECTORY}/OpenMCConfig.cmake" + INSTALL_DESTINATION ${INSTALL_CONFIGDIR} +) + +write_basic_package_version_file( + "${CMAKE_BINARY_DIR}/${CMAKE_FILES_DIRECTORY}/OpenMCConfigVersion.cmake" + VERSION ${OPENMC_VERSION} + COMPATIBILITY AnyNewerVersion +) + install(FILES - "${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/OpenMCConfig.cmake" - "${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/OpenMCConfigVersion.cmake" - DESTINATION ${INSTALL_CONFIGDIR}) + "${CMAKE_BINARY_DIR}/${CMAKE_FILES_DIRECTORY}/OpenMCConfig.cmake" + "${CMAKE_BINARY_DIR}/${CMAKE_FILES_DIRECTORY}/OpenMCConfigVersion.cmake" + DESTINATION "${INSTALL_CONFIGDIR}" +) + install(FILES man/man1/openmc.1 DESTINATION ${CMAKE_INSTALL_MANDIR}/man1) install(FILES LICENSE DESTINATION "${CMAKE_INSTALL_DOCDIR}" RENAME copyright) install(DIRECTORY include/ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}) diff --git a/cmake/OpenMCConfig.cmake.in b/cmake/OpenMCConfig.cmake.in index cc837bba7..0e15060a3 100644 --- a/cmake/OpenMCConfig.cmake.in +++ b/cmake/OpenMCConfig.cmake.in @@ -1,12 +1,18 @@ -get_filename_component(OpenMC_CMAKE_DIR "${CMAKE_CURRENT_LIST_FILE}" DIRECTORY) +@PACKAGE_INIT@ -# Compute the install prefix from this file's location -get_filename_component(_OPENMC_PREFIX "${OpenMC_CMAKE_DIR}/../../.." ABSOLUTE) +include("${CMAKE_CURRENT_LIST_DIR}/OpenMCConfigVersion.cmake") +include(CMakeFindDependencyMacro) + +# Explicitly calculate prefix if it was not generated above +if(NOT DEFINED PACKAGE_PREFIX_DIR) + get_filename_component(PACKAGE_PREFIX_DIR "${CMAKE_CURRENT_LIST_DIR}/../../.." ABSOLUTE) +endif() + +find_dependency(fmt CONFIG REQUIRED HINTS ${PACKAGE_PREFIX_DIR}) +find_dependency(pugixml CONFIG REQUIRED HINTS ${PACKAGE_PREFIX_DIR}) -find_package(fmt CONFIG REQUIRED HINTS ${_OPENMC_PREFIX}) -find_package(pugixml CONFIG REQUIRED HINTS ${_OPENMC_PREFIX}) if(@OPENMC_USE_DAGMC@) - find_package(DAGMC REQUIRED HINTS @DAGMC_DIR@) + find_dependency(DAGMC REQUIRED HINTS @DAGMC_DIR@) endif() if(@OPENMC_USE_LIBMESH@) @@ -16,20 +22,24 @@ if(@OPENMC_USE_LIBMESH@) pkg_check_modules(LIBMESH REQUIRED @LIBMESH_PC_FILE@>=1.7.0 IMPORTED_TARGET) endif() -find_package(PNG) - -if(NOT TARGET OpenMC::libopenmc) - include("${OpenMC_CMAKE_DIR}/OpenMCTargets.cmake") +if("@PNG_FOUND@") + find_dependency(PNG) endif() if(@OPENMC_USE_MPI@) - find_package(MPI REQUIRED) + find_dependency(MPI REQUIRED) endif() if(@OPENMC_USE_OPENMP@) - find_package(OpenMP REQUIRED) + find_dependency(OpenMP REQUIRED) endif() if(@OPENMC_USE_UWUW@ AND NOT ${DAGMC_BUILD_UWUW}) message(FATAL_ERROR "UWUW is enabled in OpenMC but the DAGMC installation discovered was not configured with UWUW.") endif() + +include("${CMAKE_CURRENT_LIST_DIR}/OpenMCTargets.cmake") + +if(NOT OpenMC_FIND_QUIETLY) + message(STATUS "Found OpenMC: ${PACKAGE_VERSION} (found in ${PACKAGE_PREFIX_DIR})") +endif() diff --git a/cmake/OpenMCConfigVersion.cmake.in b/cmake/OpenMCConfigVersion.cmake.in deleted file mode 100644 index 90d345de4..000000000 --- a/cmake/OpenMCConfigVersion.cmake.in +++ /dev/null @@ -1,11 +0,0 @@ -set(PACKAGE_VERSION "@OPENMC_VERSION@") - -# Check whether the requested PACKAGE_FIND_VERSION is compatible -if("${PACKAGE_VERSION}" VERSION_LESS "${PACKAGE_FIND_VERSION}") - set(PACKAGE_VERSION_COMPATIBLE FALSE) -else() - set(PACKAGE_VERSION_COMPATIBLE TRUE) - if ("${PACKAGE_VERSION}" VERSION_EQUAL "${PACKAGE_FIND_VERSION}") - set(PACKAGE_VERSION_EXACT TRUE) - endif() -endif()