From 9efa07a96cf8476457d47c495f1afbabffcb366f Mon Sep 17 00:00:00 2001 From: Ondrej Chvala Date: Fri, 2 Oct 2026 22:26:21 -0400 Subject: [PATCH 1/7] Add fission bank birth-mesh tagging Add a settings element (Settings.birth_mesh) that tags every fission site banked in an eigenvalue calculation with the bin of a mesh that contains the birth position, Particle::r_born(), of the neutron that produced the site, or with -1 if that position is outside the mesh. Any type of mesh can be used. The tag is the 0-based bin that a MeshFilter on the mesh gives for the birth position, which is the binning of MeshBornFilter. No translation or rotation is applied: these are attributes of a filter, not of a mesh, so a MeshFilter with a translation corresponds to a birth mesh that is itself moved. The tags of a saved source bank can thus be combined with the site positions, energies, times and delayed groups to form fission-matrix estimates in post-processing. Continuous-energy and multigroup fission sites are both tagged; the bin is computed once per fission event. Tagging consumes no random numbers and changes no scoring, so transport and tally results are unchanged. In fixed-source mode, where no fission sites are banked, the element is ignored with a warning. The tag is a new SourceSite member, birth_mesh_bin, which fills the existing tail padding of the struct (sizeof(SourceSite) is unchanged) and is added to the MPI datatype and to the openmc.lib ctypes mirror. It is not part of the source_bank compound type. Instead, the tags are written with the existing write_bank_dataset() template as a separate int32 birth_mesh_bin dataset, aligned by index with source_bank, together with a birth_mesh_id attribute, in state points that contain the source bank and in the source point files source..h5 and source.h5. A state point of a tagged run records birth_mesh_id even when the source bank is written separately. Surface source files, the initial source file and MCPL files carry no tags. read_source_bank() reads the tags back for a restart or a file source. The dataset must be one-dimensional with the extent of source_bank, otherwise the run stops with a fatal error. The tags are read through a memory dataspace sized to the buffer, selecting this process's share of the sites when the bank is distributed. The state point and source file format version is bumped to 18.3. Files written by earlier versions still load in Python and as file sources: StatePoint checks only the major version and returns None for birth_mesh_bin and birth_mesh_id when a file has no tags, and read_source_bank() handles the particle codes of any file as before. As for earlier minor version bumps, a restart still requires a state point of the current version. ParticleList.export_to_hdf5() now writes version 18.3 too. Tests: the data-independent multigroup regression test birth_mesh uses two fissile slabs separated by an absorber that no neutron crosses except with probability <= e^-100 (and 2^-64 per sampled flight distance, for a random number of exactly 0), so that every tag must equal the bin of its own site on a two-bin mesh. Finer regular and nonuniform rectilinear meshes show that the tags come from the birth position rather than from the site. The test also checks which files carry tags, a restart (through openmc.lib), the rejection of tag datasets that are one entry too short or too long on restart and for a file source, and fixed-source mode. Unit tests cover the Settings XML round trip, including a cylindrical birth mesh, and StatePoint.birth_mesh_bin and StatePoint.birth_mesh_id. --- docs/source/io_formats/settings.rst | 36 ++ docs/source/io_formats/source.rst | 15 + docs/source/io_formats/statepoint.rst | 16 +- docs/source/usersguide/settings.rst | 14 + include/openmc/constants.h | 2 +- include/openmc/particle_data.h | 4 + include/openmc/settings.h | 6 + include/openmc/state_point.h | 15 +- openmc/lib/core.py | 3 +- openmc/settings.py | 55 +++ openmc/source.py | 2 +- openmc/statepoint.py | 28 ++ src/initialize.cpp | 12 +- src/physics.cpp | 7 + src/physics_mg.cpp | 7 + src/settings.cpp | 28 ++ src/simulation.cpp | 5 +- src/state_point.cpp | 102 +++++- tests/regression_tests/birth_mesh/__init__.py | 0 tests/regression_tests/birth_mesh/test.py | 318 ++++++++++++++++++ tests/unit_tests/test_settings.py | 57 ++++ tests/unit_tests/test_statepoint.py | 44 +++ 22 files changed, 755 insertions(+), 21 deletions(-) create mode 100644 tests/regression_tests/birth_mesh/__init__.py create mode 100644 tests/regression_tests/birth_mesh/test.py diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index dcb63d0fc10..6fbbc1fd51c 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -33,6 +33,42 @@ source neutrons. *Default*: None +------------------------ +```` Element +------------------------ + +The ```` element indicates the ID of a mesh used to tag the fission +sites banked in an eigenvalue calculation, for the post-processing of +fission-matrix data. When this element is present, every site banked in the +fission bank, in continuous-energy and in multigroup mode, is tagged with the +bin of the mesh that contains the *birth position of the neutron that produced +the site*, or with ``-1`` if that position is outside the mesh. Only an +eigenvalue calculation banks fission sites, so in fixed-source mode the element +has no effect and is ignored with a warning. Tagging consumes no random numbers +and changes no transport or tally result. + +The birth position is the position at which the neutron was started from a +site: its fission site for a neutron started from the source bank, the +collision point for the secondary neutrons of :math:`(n,xn)` reactions, and the +point of the split for particles created by weight-window splitting. This is +the same notion of birth as used by :class:`openmc.MeshBornFilter`. + +The mesh can be of any type and must be defined in the same settings file using +a :ref:`mesh_element`. The tags are 0-based bin indices numbered in the same +order as the bins of a :class:`openmc.MeshFilter` on the mesh; for example, for +a regular mesh with :math:`n_x \times n_y \times n_z` elements and 0-based +element indices :math:`(i_x, i_y, i_z)`, the bin is +:math:`i_x + n_x (i_y + n_y i_z)`. No translation or rotation is applied: these +are attributes of a mesh filter, not of the mesh. + +The tags are written as an int32 ``birth_mesh_bin`` dataset, aligned by index +with the ``source_bank`` dataset, in state point files that contain the source +bank and in the source point files ``source..h5`` and ``source.h5``, +and the mesh ID as a ``birth_mesh_id`` attribute (see :ref:`io_statepoint` and +:ref:`io_source`). Surface source files, the initial source file and MCPL files +carry no tags. The tags are read back with the sites for a restart or a file +source. + ----------------------------- ```` Element ----------------------------- diff --git a/docs/source/io_formats/source.rst b/docs/source/io_formats/source.rst index 2e81e0a3f2f..417ed573a28 100644 --- a/docs/source/io_formats/source.rst +++ b/docs/source/io_formats/source.rst @@ -4,6 +4,9 @@ Source File Format ================== +The current version of the source file format is 18.3. Source files share their +version number with the :ref:`state point file format `. + Normally, source data is stored in a state point file. However, it is possible to request that the source be written separately, in which case the format used is that documented here. @@ -17,6 +20,9 @@ following the same format. :Attributes: - **filetype** (*char[]*) -- String indicating the type of file. - **version** (*int[2]*) -- Major and minor version of the source file format. + - **birth_mesh_id** (*int*) -- ID of the mesh used to tag fission + sites with the ```` settings element. Only present + together with the ``birth_mesh_bin`` dataset. :Datasets: @@ -26,3 +32,12 @@ following the same format. which represent the position, direction, energy, time, weight, delayed group, surface ID, and particle type (PDG number), respectively. + - **birth_mesh_bin** (*int[]*) -- Fission bank birth tags, aligned + by index with ``source_bank``, as described for the + :ref:`state point file `. Only present in source + point files of an eigenvalue calculation for which the + ```` settings element was specified; surface source + files and the initial source file never contain this dataset. + MCPL source files do not carry the tags. When a file containing + this dataset is read, the dataset must have the same length as + ``source_bank``. diff --git a/docs/source/io_formats/statepoint.rst b/docs/source/io_formats/statepoint.rst index 6e0561e2481..fd058e4dbe4 100644 --- a/docs/source/io_formats/statepoint.rst +++ b/docs/source/io_formats/statepoint.rst @@ -4,7 +4,7 @@ State Point File Format ======================= -The current version of the statepoint file format is 18.2. +The current version of the statepoint file format is 18.3. **/** @@ -21,6 +21,12 @@ The current version of the statepoint file format is 18.2. are present (1) or not (0). - **source_present** (*int*) -- Flag indicating whether the source bank is present (1) or not (0). + - **birth_mesh_id** (*int*) -- ID of the mesh used to tag fission + sites with the ```` settings element. Only present + when `run_mode` is 'eigenvalue' and fission bank tagging is + enabled; it is written even if the source bank is not stored in + the file (e.g., because it is written to a separate source + file). :Datasets: - **seed** (*int8_t*) -- Pseudo-random number generator seed. - **stride** (*uint64_t*) -- Pseudo-random number generator stride. @@ -58,6 +64,14 @@ The current version of the statepoint file format is 18.2. time, weight, delayed group, surface ID, and particle type (PDG number), respectively. Only present when `run_mode` is 'eigenvalue'. + - **birth_mesh_bin** (*int[]*) -- Fission bank birth tags, aligned + by index with ``source_bank``: for each site, the 0-based bin, + numbered as for a mesh filter, of the birth position of the + neutron that produced the site in the mesh referenced by the + ``birth_mesh_id`` attribute, or ``-1`` if that position is outside + the mesh. Only present when ``source_bank`` is present and the + ```` settings element was specified for an eigenvalue + calculation. **/tallies/** diff --git a/docs/source/usersguide/settings.rst b/docs/source/usersguide/settings.rst index 5fad64f3c10..8363e34e7a3 100644 --- a/docs/source/usersguide/settings.rst +++ b/docs/source/usersguide/settings.rst @@ -718,6 +718,20 @@ As an example, to write a statepoint file every five batches:: settings.batches = n settings.statepoint = {'batches': range(5, n + 5, 5)} +In an eigenvalue calculation, the fission sites of the source bank can be +tagged for fission-matrix post-processing by setting +:attr:`Settings.birth_mesh` to a mesh, which can be of any type. Each site is +tagged with the 0-based bin that a :class:`openmc.MeshFilter` on that mesh +gives for the birth position of the neutron that produced the site (the binning +of :class:`openmc.MeshBornFilter`), or with -1 if that position is outside the +mesh. No translation or rotation is applied, since these are attributes of a +mesh filter, not of the mesh. The tags are written as a ``birth_mesh_bin`` +dataset to the statepoint files that contain the source bank and to the source +files ``source..h5`` and ``source.h5``. In fixed-source mode, the +attribute has no effect. For example:: + + settings.birth_mesh = mesh + Particle Track Files -------------------- diff --git a/include/openmc/constants.h b/include/openmc/constants.h index a210083de32..6ff0d0dbf51 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -26,7 +26,7 @@ using double_4dvec = vector>>>; constexpr int HDF5_VERSION[] {3, 1}; // Version numbers for binary files -constexpr array VERSION_STATEPOINT {18, 2}; +constexpr array VERSION_STATEPOINT {18, 3}; constexpr array VERSION_PARTICLE_RESTART {2, 1}; constexpr array VERSION_TRACK {3, 1}; constexpr array VERSION_SUMMARY {6, 1}; diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index 24b7eb53c0a..af36a24b7cf 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -56,6 +56,10 @@ struct SourceSite { double wgt_ww_born {-1.0}; int64_t n_split {0}; int n_collision {0}; + // Bin of the mesh that contains the birth position of the + // neutron that produced this fission site; -1 if it is outside the mesh or + // the site is not tagged. Not part of the HDF5 "source_bank" compound type. + int birth_mesh_bin {-1}; }; struct CollisionTrackSite { diff --git a/include/openmc/settings.h b/include/openmc/settings.h index 3ca7e9c3a82..a4340969f2c 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -46,8 +46,14 @@ struct CollisionTrackConfig { // Global variable declarations //============================================================================== +class Mesh; + namespace settings { +// Mesh used to tag banked fission sites (see SourceSite::birth_mesh_bin); +// nullptr when no tagging is done, including in fixed-source mode +extern const Mesh* birth_mesh; + // Boolean flags extern bool assume_separate; //!< assume tallies are spatially separate? extern bool check_overlaps; //!< check overlaps in geometry? diff --git a/include/openmc/state_point.h b/include/openmc/state_point.h index fb1aaf7b985..4caf755b82c 100644 --- a/include/openmc/state_point.h +++ b/include/openmc/state_point.h @@ -36,16 +36,23 @@ void load_state_point(); // // The source_bank variable is used as work space if MPI is used, // so it cannot be given as a const span. +// +// If birth_mesh_tags is true, the fission bank birth tags of the sites are +// written as a "birth_mesh_bin" dataset together with a "birth_mesh_id" +// attribute. This should only be requested for a bank of fission sites banked +// in an eigenvalue calculation with a (see settings::birth_mesh). void write_h5_source_point(const char* filename, span source_bank, - const vector& bank_index); + const vector& bank_index, bool birth_mesh_tags = false); void write_source_point(std::string, span source_bank, - const vector& bank_index, bool use_mcpl); + const vector& bank_index, bool use_mcpl, + bool birth_mesh_tags = false); // This appends a source bank specification to an HDF5 file -// that's already open. It is used internally by write_source_point. +// that's already open. It is used internally by write_source_point. If +// birth_mesh_tags is true, the "birth_mesh_bin" dataset is also written. void write_source_bank(hid_t group_id, span source_bank, - const vector& bank_index); + const vector& bank_index, bool birth_mesh_tags = false); void read_source_bank( hid_t group_id, vector& sites, bool distribute); diff --git a/openmc/lib/core.py b/openmc/lib/core.py index 72901008712..11131797d14 100644 --- a/openmc/lib/core.py +++ b/openmc/lib/core.py @@ -35,7 +35,8 @@ class _SourceSite(Structure): ('wgt_born', c_double), ('wgt_ww_born', c_double), ('n_split', c_int64), - ('n_collision', c_int)] + ('n_collision', c_int), + ('birth_mesh_bin', c_int)] # 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 5d71df091c5..6ae0a20cec9 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -47,6 +47,21 @@ class Settings: interactions. batches : int Number of batches to simulate + birth_mesh : openmc.MeshBase + Mesh used to tag the fission sites banked in an eigenvalue calculation + with the 0-based bin, numbered as for :class:`openmc.MeshFilter`, of + the birth position of the neutron that produced each site, or -1 + outside the mesh. Any type of mesh can be used. No translation or + rotation is applied, since these are attributes of a mesh filter, not + of the mesh. The birth position is where that neutron was started + from a site: its fission site, or the collision or split point for + (n,xn) secondaries and split particles, as for + :class:`openmc.MeshBornFilter`. The tags are written to state point and + source point files as a ``birth_mesh_bin`` dataset aligned with the + source bank (see :ref:`io_statepoint`). Has no effect in fixed-source + mode. + + .. versionadded:: 0.16.1 confidence_intervals : bool If True, uncertainties on tally results will be reported as the half-width of the 95% two-sided confidence interval. If False, @@ -485,6 +500,9 @@ def __init__(self, **kwargs): # Uniform fission source subelement self._ufs_mesh = None + # Fission bank birth tagging mesh + self._birth_mesh = None + self._resonance_scattering = {} self._volume_calculations = cv.CheckedList( VolumeCalculation, 'volume calculations') @@ -1189,6 +1207,15 @@ def ufs_mesh(self, ufs_mesh: RegularMesh): cv.check_length('UFS mesh upper-right corner', ufs_mesh.upper_right, 3) self._ufs_mesh = ufs_mesh + @property + def birth_mesh(self) -> MeshBase: + return self._birth_mesh + + @birth_mesh.setter + def birth_mesh(self, mesh: MeshBase): + cv.check_type('birth mesh', mesh, MeshBase) + self._birth_mesh = mesh + @property def resonance_scattering(self) -> dict: return self._resonance_scattering @@ -1873,6 +1900,23 @@ def _create_ufs_mesh_subelement(self, root, mesh_memo=None): if mesh_memo is not None: mesh_memo.add(self.ufs_mesh.id) + def _create_birth_mesh_subelement(self, root, mesh_memo=None): + if self.birth_mesh is None: + return + + subelement = ET.SubElement(root, "birth_mesh") + subelement.text = str(self.birth_mesh.id) + + if mesh_memo and self.birth_mesh.id in mesh_memo: + return + + # See if a element already exists -- if not, add it + path = f"./mesh[@id='{self.birth_mesh.id}']" + if root.find(path) is None: + root.append(self.birth_mesh.to_xml_element()) + if mesh_memo is not None: + mesh_memo.add(self.birth_mesh.id) + def _create_use_decay_photons_subelement(self, root): if self._use_decay_photons is not None: element = ET.SubElement(root, "use_decay_photons") @@ -2400,6 +2444,15 @@ def _ufs_mesh_from_xml_element(self, root, meshes): raise ValueError(f'Could not locate mesh with ID "{mesh_id}"') self.ufs_mesh = meshes[mesh_id] + def _birth_mesh_from_xml_element(self, root, meshes): + text = get_text(root, 'birth_mesh') + if text is None: + return + mesh_id = int(text) + if mesh_id not in meshes: + raise ValueError(f'Could not locate mesh with ID "{mesh_id}"') + self.birth_mesh = meshes[mesh_id] + def _resonance_scattering_from_xml_element(self, root): elem = root.find('resonance_scattering') if elem is not None: @@ -2639,6 +2692,7 @@ def to_xml_element(self, mesh_memo=None): self._create_trace_subelement(element) self._create_track_subelement(element) self._create_ufs_mesh_subelement(element, mesh_memo) + self._create_birth_mesh_subelement(element, mesh_memo) self._create_resonance_scattering_subelement(element) self._create_volume_calcs_subelement(element) self._create_create_fission_neutrons_subelement(element) @@ -2758,6 +2812,7 @@ def from_xml_element(cls, elem, meshes=None): settings._trace_from_xml_element(elem) settings._track_from_xml_element(elem) settings._ufs_mesh_from_xml_element(elem, meshes) + settings._birth_mesh_from_xml_element(elem, meshes) settings._resonance_scattering_from_xml_element(elem) settings._create_fission_neutrons_from_xml_element(elem) settings._create_delayed_neutrons_from_xml_element(elem) diff --git a/openmc/source.py b/openmc/source.py index e68bdd5cf25..3a04035b3c9 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -1582,7 +1582,7 @@ def export_to_hdf5(self, filename: PathLike, **kwargs): kwargs.setdefault('mode', 'w') with h5py.File(filename, **kwargs) as fh: fh.attrs['filetype'] = np.bytes_("source") - fh.attrs['version'] = np.array([_VERSION_STATEPOINT, 2]) + fh.attrs['version'] = np.array([_VERSION_STATEPOINT, 3]) fh.create_dataset('source_bank', data=arr, dtype=source_dtype) diff --git a/openmc/statepoint.py b/openmc/statepoint.py index e22c55d40f0..96efe07c801 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -36,6 +36,22 @@ class StatePoint: Attributes ---------- + birth_mesh_bin : numpy.ndarray or None + Fission bank birth tag of each source site, aligned by index with + :attr:`source`: the 0-based bin, numbered as the bins of a + :class:`openmc.MeshFilter`, of the birth position of the neutron that + produced the site in the mesh given by + :attr:`openmc.Settings.birth_mesh`, or -1 outside the mesh. None if + the file contains no tags, either because tagging was not enabled or + because the source bank is not stored in this file. + + .. versionadded:: 0.16.1 + birth_mesh_id : int or None + ID of the mesh used to tag fission sites (see :attr:`birth_mesh_bin`), + or None if tagging was not enabled. It is also given when the source + bank, and with it the tags, is written to a separate source file. + + .. versionadded:: 0.16.1 cmfd_on : bool Indicate whether CMFD is active cmfd_balance : numpy.ndarray @@ -371,6 +387,18 @@ def stride(self): def source(self): return self._f['source_bank'][()] if self.source_present else None + @property + def birth_mesh_bin(self): + if self.source_present and 'birth_mesh_bin' in self._f: + return self._f['birth_mesh_bin'][()] + return None + + @property + def birth_mesh_id(self): + if 'birth_mesh_id' in self._f.attrs: + return int(self._f.attrs['birth_mesh_id']) + return None + @property def source_present(self): return self._f.attrs['source_present'] > 0 diff --git a/src/initialize.cpp b/src/initialize.cpp index fc259246f7a..66e2b470e9d 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -173,7 +173,7 @@ void initialize_mpi(MPI_Comm intracomm) // Create bank datatype SourceSite b; - MPI_Aint disp[15]; + MPI_Aint disp[16]; MPI_Get_address(&b.r, &disp[0]); MPI_Get_address(&b.u, &disp[1]); MPI_Get_address(&b.E, &disp[2]); @@ -189,12 +189,13 @@ void initialize_mpi(MPI_Comm intracomm) MPI_Get_address(&b.wgt_ww_born, &disp[12]); MPI_Get_address(&b.n_split, &disp[13]); MPI_Get_address(&b.n_collision, &disp[14]); - for (int i = 14; i >= 0; --i) { + MPI_Get_address(&b.birth_mesh_bin, &disp[15]); + for (int i = 15; i >= 0; --i) { disp[i] -= disp[0]; } // Block counts for each field - int blocks[] = {3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}; + int blocks[] = {3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}; // Types for each field MPI_Datatype types[] = { @@ -212,10 +213,11 @@ void initialize_mpi(MPI_Comm intracomm) MPI_DOUBLE, // wgt_born MPI_DOUBLE, // wgt_ww_born MPI_INT64_T, // n_split - MPI_INT // n_collision + MPI_INT, // n_collision + MPI_INT // birth_mesh_bin }; - MPI_Type_create_struct(15, blocks, disp, types, &mpi::source_site); + MPI_Type_create_struct(16, blocks, disp, types, &mpi::source_site); MPI_Type_commit(&mpi::source_site); CollisionTrackSite bc; diff --git a/src/physics.cpp b/src/physics.cpp index 342d4144685..8acdfba63ca 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -11,6 +11,7 @@ #include "openmc/ifp.h" #include "openmc/material.h" #include "openmc/math_functions.h" +#include "openmc/mesh.h" #include "openmc/message_passing.h" #include "openmc/ncrystal_interface.h" #include "openmc/nuclide.h" @@ -208,6 +209,11 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) // fission bank or the secondary particle bank int n_sites_stored; + // Tag for the sites: the bin of this neutron's birth position in the birth + // mesh, if one was specified (no random numbers are consumed) + int birth_mesh_bin = + settings::birth_mesh ? settings::birth_mesh->get_bin(p.r_born()) : -1; + for (n_sites_stored = 0; n_sites_stored < nu; n_sites_stored++) { // Initialize fission site object with particle data SourceSite site; @@ -216,6 +222,7 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) site.time = p.time(); site.wgt = 1. / weight; site.surf_id = 0; + site.birth_mesh_bin = birth_mesh_bin; // Sample delayed group and angle/energy for fission reaction sample_fission_neutron(i_nuclide, rx, &site, p); diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index 212ba765ca0..6c7da4cb7d4 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -11,6 +11,7 @@ #include "openmc/error.h" #include "openmc/material.h" #include "openmc/math_functions.h" +#include "openmc/mesh.h" #include "openmc/message_passing.h" #include "openmc/mgxs_interface.h" #include "openmc/particle.h" @@ -133,6 +134,11 @@ void create_fission_sites(Particle& p) // fission bank or the secondary particle bank int n_sites_stored; + // Tag for the sites: the bin of this neutron's birth position in the birth + // mesh, if one was specified (no random numbers are consumed) + int birth_mesh_bin = + settings::birth_mesh ? settings::birth_mesh->get_bin(p.r_born()) : -1; + for (n_sites_stored = 0; n_sites_stored < nu; n_sites_stored++) { // Initialize fission site object with particle data SourceSite site; @@ -140,6 +146,7 @@ void create_fission_sites(Particle& p) site.particle = ParticleType::neutron(); site.time = p.time(); site.wgt = 1. / weight; + site.birth_mesh_bin = birth_mesh_bin; // Sample the cosine of the angle, assuming fission neutrons are emitted // isotropically diff --git a/src/settings.cpp b/src/settings.cpp index 5e2da3a52e4..1f0809198b0 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -113,6 +113,7 @@ int max_particle_events {1000000}; ElectronTreatment electron_treatment {ElectronTreatment::TTB}; array energy_cutoff {0.0, 1000.0, 0.0, 0.0}; array time_cutoff {INFTY, INFTY, INFTY, INFTY}; +const Mesh* birth_mesh {nullptr}; int ifp_n_generation {-1}; int legendre_to_tabular_points {C_NONE}; int max_order {0}; @@ -881,6 +882,32 @@ void read_settings_xml(pugi::xml_node root) "it by specifying its ID in a element."); } + // Fission bank birth-tag mesh: every banked fission site is tagged with the + // 0-based bin that a MeshFilter on this mesh gives for the birth position of + // the creating neutron (the binning of MeshBornFilter), or with -1 if that + // position is outside the mesh. Any type of mesh can be used. No translation + // or rotation is applied: these are attributes of a filter, not of a mesh. + // The mesh must be defined in settings.xml (read above). + if (check_for_node(root, "birth_mesh")) { + auto temp = std::stoi(get_node_value(root, "birth_mesh")); + if (model::mesh_map.find(temp) == model::mesh_map.end()) { + fatal_error(fmt::format("Mesh {} specified for fission bank birth " + "tagging does not exist.", + temp)); + } + + // Only an eigenvalue calculation banks fission sites in the fission bank, + // so birth tagging has no effect in fixed-source mode + if (run_mode == RunMode::FIXED_SOURCE) { + warning("The element has no effect in fixed-source mode " + "and is ignored."); + } else { + auto& m = model::meshes[model::mesh_map.at(temp)]; + m->prepare_for_point_location(); + birth_mesh = m.get(); + } + } + // Check if the user has specified to write state points if (check_for_node(root, "state_point")) { @@ -1379,6 +1406,7 @@ void free_memory_settings() settings::track_identifiers.clear(); settings::ifp_delayed_group_on = false; settings::ifp_lifetime_on = false; + settings::birth_mesh = nullptr; } //============================================================================== diff --git a/src/simulation.cpp b/src/simulation.cpp index 4e74193f39a..505bdeeece7 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -573,7 +573,8 @@ void finalize_batch() settings::path_output, simulation::current_batch, w); span bankspan(simulation::source_bank); write_source_point(source_point_filename, bankspan, - simulation::work_index, settings::source_mcpl_write); + simulation::work_index, settings::source_mcpl_write, + settings::birth_mesh != nullptr); } // Write a continously-overwritten source point if requested. @@ -581,7 +582,7 @@ void finalize_batch() auto filename = settings::path_output + "source"; span bankspan(simulation::source_bank); write_source_point(filename, bankspan, simulation::work_index, - settings::source_mcpl_write); + settings::source_mcpl_write, settings::birth_mesh != nullptr); } } diff --git a/src/state_point.cpp b/src/state_point.cpp index 06d376fe25b..978e0c92353 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -68,6 +68,11 @@ extern "C" int openmc_statepoint_write(const char* filename, bool* write_source) // Determine whether or not to write the source bank bool write_source_ = write_source ? *write_source : true; + // Fission bank birth tags are only defined for the source bank of an + // eigenvalue calculation, which consists of banked fission sites + bool birth_mesh_tags = + settings::birth_mesh && settings::run_mode == RunMode::EIGENVALUE; + // Write message write_message("Creating state point " + filename_ + "...", 5); @@ -123,6 +128,10 @@ extern "C" int openmc_statepoint_write(const char* filename, bool* write_source) // Indicate whether source bank is stored in statepoint write_attribute(file_id, "source_present", write_source_); + // Record the mesh used for fission bank birth tags, if any + if (birth_mesh_tags) + write_attribute(file_id, "birth_mesh_id", settings::birth_mesh->id()); + // Write out information for eigenvalue run if (settings::run_mode == RunMode::EIGENVALUE) write_eigenvalue_hdf5(file_id); @@ -350,7 +359,8 @@ extern "C" int openmc_statepoint_write(const char* filename, bool* write_source) if (write_source_) { if (mpi::master || parallel) file_id = file_open(filename_, 'a', true); - write_source_bank(file_id, simulation::source_bank, simulation::work_index); + write_source_bank(file_id, simulation::source_bank, simulation::work_index, + birth_mesh_tags); if (mpi::master || parallel) file_close(file_id); } @@ -596,7 +606,7 @@ hid_t h5banktype(bool memory) } void write_source_point(std::string filename, span source_bank, - const vector& bank_index, bool use_mcpl) + const vector& bank_index, bool use_mcpl, bool birth_mesh_tags) { std::string ext = use_mcpl ? "mcpl" : "h5"; @@ -616,12 +626,13 @@ void write_source_point(std::string filename, span source_bank, write_mcpl_source_point(filename.c_str(), source_bank, bank_index); } else { filename.append(".h5"); - write_h5_source_point(filename.c_str(), source_bank, bank_index); + write_h5_source_point( + filename.c_str(), source_bank, bank_index, birth_mesh_tags); } } void write_h5_source_point(const char* filename, span source_bank, - const vector& bank_index) + const vector& bank_index, bool birth_mesh_tags) { // When using parallel HDF5, the file is written to collectively by all // processes. With MPI-only, the file is opened and written by the master @@ -648,17 +659,19 @@ void write_h5_source_point(const char* filename, span source_bank, file_id = file_open(filename_.c_str(), 'w', true); write_attribute(file_id, "filetype", "source"); write_attribute(file_id, "version", VERSION_STATEPOINT); + if (birth_mesh_tags) + write_attribute(file_id, "birth_mesh_id", settings::birth_mesh->id()); } // Get pointer to source bank and write to file - write_source_bank(file_id, source_bank, bank_index); + write_source_bank(file_id, source_bank, bank_index, birth_mesh_tags); if (mpi::master || parallel) file_close(file_id); } void write_source_bank(hid_t group_id, span source_bank, - const vector& bank_index) + const vector& bank_index, bool birth_mesh_tags) { hid_t membanktype = h5banktype(true); hid_t filebanktype = h5banktype(false); @@ -673,6 +686,23 @@ void write_source_bank(hid_t group_id, span source_bank, H5Tclose(membanktype); H5Tclose(filebanktype); + + // Write the fission bank birth tags as an int32 dataset aligned by index with + // "source_bank". The buffer has the size of the local source bank, which + // write_bank_dataset also uses as a receive buffer on the master process. + if (birth_mesh_tags) { + vector bins(source_bank.size()); + for (size_t i = 0; i < source_bank.size(); ++i) + bins[i] = source_bank[i].birth_mesh_bin; + span bins_span(bins); +#ifdef OPENMC_MPI + write_bank_dataset("birth_mesh_bin", group_id, bins_span, bank_index, + H5T_NATIVE_INT32, H5T_NATIVE_INT32, MPI_INT32_T); +#else + write_bank_dataset("birth_mesh_bin", group_id, bins_span, bank_index, + H5T_NATIVE_INT32, H5T_NATIVE_INT32); +#endif + } } // Determine member names of a compound HDF5 datatype @@ -765,6 +795,66 @@ void read_source_bank( H5Dclose(dset); H5Tclose(banktype); + // Read the birth mesh bin tag dataset, when present, to keep the tags + // aligned with the source sites (e.g. for restart runs). The dataset must + // have exactly one entry per site of the "source_bank" dataset. + if (object_exists(group_id, "birth_mesh_bin")) { + hid_t dset_bin = H5Dopen(group_id, "birth_mesh_bin", H5P_DEFAULT); + hid_t dspace_bin = H5Dget_space(dset_bin); + int rank_bin = H5Sget_simple_extent_ndims(dspace_bin); + hsize_t n_bins = 0; + if (rank_bin == 1) + H5Sget_simple_extent_dims(dspace_bin, &n_bins, nullptr); + if (rank_bin != 1 || n_bins != n_sites) { + H5Sclose(dspace_bin); + H5Dclose(dset_bin); + fatal_error(fmt::format( + "The birth_mesh_bin dataset in the source file must be " + "one-dimensional with the same length as the source_bank dataset " + "({} sites).", + n_sites)); + } + + // Number of tags to read: this process's share of the sites when + // distributing them, otherwise all of them + hsize_t n_read = + distribute ? static_cast(simulation::work_per_rank) : n_sites; + vector birth_bins(n_read); + + // Memory dataspace sized to the buffer, and the matching selection of + // the file dataspace (the whole dataset unless distributing) + hid_t memspace_bin = H5Screate_simple(1, &n_read, nullptr); + if (distribute) { + hsize_t offset = simulation::work_index[mpi::rank]; + H5Sselect_hyperslab( + dspace_bin, H5S_SELECT_SET, &offset, nullptr, &n_read, nullptr); + } + +#ifdef PHDF5 + // Read data in parallel + hid_t plist = H5Pcreate(H5P_DATASET_XFER); + H5Pset_dxpl_mpio(plist, H5FD_MPIO_COLLECTIVE); + herr_t status = H5Dread(dset_bin, H5T_NATIVE_INT32, memspace_bin, + dspace_bin, plist, birth_bins.data()); + H5Pclose(plist); +#else + herr_t status = H5Dread(dset_bin, H5T_NATIVE_INT32, memspace_bin, + dspace_bin, H5P_DEFAULT, birth_bins.data()); +#endif + + H5Sclose(memspace_bin); + H5Sclose(dspace_bin); + H5Dclose(dset_bin); + if (status < 0) { + fatal_error("Failed to read the birth_mesh_bin dataset from the source " + "file."); + } + + for (hsize_t i = 0; i < n_read; ++i) { + sites[i].birth_mesh_bin = birth_bins[i]; + } + } + if (legacy_particle_codes) { for (auto& site : sites) { site.particle = legacy_particle_index_to_type(site.particle.pdg_number()); diff --git a/tests/regression_tests/birth_mesh/__init__.py b/tests/regression_tests/birth_mesh/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/birth_mesh/test.py b/tests/regression_tests/birth_mesh/test.py new file mode 100644 index 00000000000..46675810def --- /dev/null +++ b/tests/regression_tests/birth_mesh/test.py @@ -0,0 +1,318 @@ +"""Test of the fission bank birth-mesh tags in a decoupled two-slab problem. + +Two one-group fissile slabs are separated by a 10 cm pure absorber with +Sigma_t = 10/cm, in which a neutron is absorbed at its first collision. A +neutron needs at least 10 cm of uncollided path in the absorber to cross it, +so no neutron crosses it except with probability <= e^-100 (and 2^-64 per +sampled flight distance, for a random number of exactly 0, which gives an +infinite flight). The birth position of the neutron that produced a banked +fission site is therefore in the same slab as the site itself. With a birth +mesh that has one bin per slab, every tag must equal the bin of the position +of its own site, so the test does not need reference results. A finer mesh, +with two bins per slab, checks that the tags are taken from the birth position +of the neutron that produced the site rather than from the position of the +site, for a regular and for a nonuniform rectilinear mesh. + +The tests also check which output files carry the tags, the restart path that +reads them back, and the rejection of a tag dataset whose length differs from +that of the source bank. + +""" +import shutil + +import h5py +import numpy as np +import openmc +import openmc.lib +import pytest + +from tests.regression_tests import config + +# One-group data: subcritical fuel, so that the fixed-source variant with +# fission neutrons is well defined, and a pure absorber +FUEL = {'total': 1.0, 'absorption': 0.5, 'fission': 0.15, 'nu': 2.5} +ABSORBER_SIGMA = 10.0 # 1/cm, total = absorption +HALF_WIDTH = 5.0 # cm, half-thickness of the absorber +FUEL_WIDTH = 5.0 # cm, thickness of each fuel slab + + +def make_model(library): + groups = openmc.mgxs.EnergyGroups([0.0, 20.0e6]) + fuel_xs = openmc.XSdata('fuel', groups) + fuel_xs.order = 0 + fuel_xs.set_total([FUEL['total']]) + fuel_xs.set_absorption([FUEL['absorption']]) + fuel_xs.set_fission([FUEL['fission']]) + fuel_xs.set_nu_fission([FUEL['nu'] * FUEL['fission']]) + fuel_xs.set_scatter_matrix([[[FUEL['total'] - FUEL['absorption']]]]) + fuel_xs.set_chi([1.0]) + absorber_xs = openmc.XSdata('absorber', groups) + absorber_xs.order = 0 + absorber_xs.set_total([ABSORBER_SIGMA]) + absorber_xs.set_absorption([ABSORBER_SIGMA]) + absorber_xs.set_scatter_matrix([[[0.0]]]) + mgxs_library = openmc.MGXSLibrary(groups) + mgxs_library.add_xsdata(fuel_xs) + mgxs_library.add_xsdata(absorber_xs) + mgxs_library.export_to_hdf5(library) + + fuel = openmc.Material(name='fuel') + fuel.set_density('macro', 1.0) + fuel.add_macroscopic('fuel') + absorber = openmc.Material(name='absorber') + absorber.set_density('macro', 1.0) + absorber.add_macroscopic('absorber') + + x_max = HALF_WIDTH + FUEL_WIDTH + x0 = openmc.XPlane(-x_max, boundary_type='reflective') + x1 = openmc.XPlane(-HALF_WIDTH) + x2 = openmc.XPlane(HALF_WIDTH) + x3 = openmc.XPlane(x_max, boundary_type='reflective') + y0 = openmc.YPlane(-5.0, boundary_type='reflective') + y1 = openmc.YPlane(5.0, boundary_type='reflective') + z0 = openmc.ZPlane(-5.0, boundary_type='reflective') + z1 = openmc.ZPlane(5.0, boundary_type='reflective') + yz = +y0 & -y1 & +z0 & -z1 + cells = [openmc.Cell(fill=fuel, region=+x0 & -x1 & yz), + openmc.Cell(fill=absorber, region=+x1 & -x2 & yz), + openmc.Cell(fill=fuel, region=+x2 & -x3 & yz)] + + model = openmc.Model() + model.materials = openmc.Materials([fuel, absorber]) + model.materials.cross_sections = str(library) + model.geometry = openmc.Geometry(cells) + model.settings.energy_mode = 'multi-group' + model.settings.particles = 2000 + model.settings.source = openmc.IndependentSource( + space=openmc.stats.Box((-x_max, -5.0, -5.0), (x_max, 5.0, 5.0)), + constraints={'fissionable': True}) + return model, x1 + + +def slab_mesh(nx): + """Regular mesh over the problem with nx bins of equal width along x.""" + x_max = HALF_WIDTH + FUEL_WIDTH + mesh = openmc.RegularMesh() + mesh.lower_left = (-x_max, -5.0, -5.0) + mesh.upper_right = (x_max, 5.0, 5.0) + mesh.dimension = (nx, 1, 1) + return mesh + + +def rectilinear_slab_mesh(): + """Nonuniform rectilinear mesh with two bins along x in each fuel slab, one + for the absorber, and two unequal bins along y.""" + x_max = HALF_WIDTH + FUEL_WIDTH + mesh = openmc.RectilinearMesh() + mesh.x_grid = [-x_max, -7.0, -HALF_WIDTH, HALF_WIDTH, 7.5, x_max] + mesh.y_grid = [-5.0, 1.0, 5.0] + mesh.z_grid = [-5.0, 5.0] + return mesh + + +def bin_of(mesh, r): + """0-based bin (x fastest) of each position of r in a Cartesian mesh.""" + xyz = np.column_stack([r['x'], r['y'], r['z']]) + if isinstance(mesh, openmc.RegularMesh): + ll = np.asarray(mesh.lower_left) + width = (np.asarray(mesh.upper_right) - ll) / mesh.dimension + ijk = np.floor((xyz - ll) / width).astype(int) + else: + grids = (mesh.x_grid, mesh.y_grid, mesh.z_grid) + ijk = np.column_stack([np.searchsorted(g, xyz[:, i], side='right') - 1 + for i, g in enumerate(grids)]) + nx, ny, _ = mesh.dimension + return ijk[:, 0] + nx * (ijk[:, 1] + ny * ijk[:, 2]) + + +def run_model(model, path, **kwargs): + kwargs.update({'cwd': path, 'openmc_exec': config['exe'], + 'event_based': config['event']}) + if config['mpi']: + kwargs['mpi_args'] = [config['mpiexec'], '-n', config['mpi_np']] + return model.run(**kwargs) + + +def eigenvalue_model(tmp_path, mesh): + model, surface = make_model(tmp_path / 'two_slab.h5') + model.settings.run_mode = 'eigenvalue' + model.settings.batches = 6 + model.settings.inactive = 2 + model.settings.birth_mesh = mesh + return model, surface + + +def read_source(path): + with h5py.File(path, 'r') as f: + attrs = dict(f.attrs) + source = f['source_bank'][()] + tags = f['birth_mesh_bin'][()] if 'birth_mesh_bin' in f else None + return attrs, source, tags + + +def test_birth_mesh_eigenvalue(tmp_path): + # One bin per slab: bin 0 holds the left slab and bin 1 the right one + mesh = slab_mesh(2) + model, surface = eigenvalue_model(tmp_path, mesh) + # Also write a surface source file and the initial source file + model.settings.surf_source_write = {'surface_ids': [surface.id], + 'max_particles': 1000} + model.settings.write_initial_source = True + sp_path = run_model(model, tmp_path) + + attrs, source, tags = read_source(sp_path) + with openmc.StatePoint(sp_path) as sp: + assert sp.source_present + assert sp.birth_mesh_id == mesh.id + np.testing.assert_array_equal(sp.birth_mesh_bin, tags) + + # The tags are int32, aligned with the source bank, and refer to the mesh + assert attrs['birth_mesh_id'] == mesh.id + assert tags is not None + assert tags.dtype == np.int32 + assert tags.shape == source.shape == (model.settings.particles,) + + # Every site is tagged with the bin of its own slab, and both slabs hold + # sites + np.testing.assert_array_equal(tags, bin_of(mesh, source['r'])) + assert set(np.unique(tags)) == {0, 1} + + # Surface source and initial source sites are not fission bank sites, so + # these files carry no tags + for name in ('surface_source.h5', 'initial_source.h5'): + file_attrs, _, file_tags = read_source(tmp_path / name) + assert 'birth_mesh_id' not in file_attrs + assert file_tags is None + + +def test_birth_mesh_source_files(tmp_path): + mesh = slab_mesh(2) + model, _ = eigenvalue_model(tmp_path, mesh) + # Write the source bank of the last batch to source.6.h5 and the latest + # source bank to source.h5 instead of to the state point + model.settings.sourcepoint = {'separate': True, 'overwrite': True} + sp_path = run_model(model, tmp_path) + + # The state point records the mesh but holds no source bank or tags + with h5py.File(sp_path, 'r') as f: + assert f.attrs['birth_mesh_id'] == mesh.id + assert 'source_bank' not in f + assert 'birth_mesh_bin' not in f + with openmc.StatePoint(sp_path) as sp: + assert not sp.source_present + assert sp.birth_mesh_id == mesh.id + assert sp.birth_mesh_bin is None + + # Both source point files carry the tags of the final source bank + _, source, tags = read_source(tmp_path / 'source.6.h5') + np.testing.assert_array_equal(tags, bin_of(mesh, source['r'])) + for name in ('source.6.h5', 'source.h5'): + file_attrs, file_source, file_tags = read_source(tmp_path / name) + assert file_attrs['birth_mesh_id'] == mesh.id + assert file_tags.dtype == np.int32 + np.testing.assert_array_equal(file_source, source) + np.testing.assert_array_equal(file_tags, tags) + + +@pytest.mark.parametrize('rectilinear', [False, True]) +def test_birth_mesh_birth_position(tmp_path, rectilinear): + # With two bins along x in each fuel slab (the first two and the last two), + # a neutron is often born in a different bin from the one in which it + # causes a fission, but always in the same slab + mesh = rectilinear_slab_mesh() if rectilinear else slab_mesh(8) + model, _ = eigenvalue_model(tmp_path, mesh) + _, source, tags = read_source(run_model(model, tmp_path)) + + nx = mesh.dimension[0] + own_bin = bin_of(mesh, source['r']) + assert tags.min() >= 0 + left, right = {0, 1}, {nx - 2, nx - 1} + assert set(np.unique(own_bin % nx)) == left | right + for slab in (left, right): + in_slab = np.isin(own_bin % nx, list(slab)) + assert set(np.unique(tags[in_slab] % nx)) <= slab + assert set(np.unique(tags // nx)) == set(np.unique(own_bin // nx)) + assert np.count_nonzero(tags != own_bin) > 0 + + +def test_birth_mesh_restart(tmp_path, monkeypatch): + mesh = slab_mesh(2) + model, _ = eigenvalue_model(tmp_path, mesh) + model.settings.statepoint = {'batches': [3, 6]} + run_model(model, tmp_path) + restart_file = tmp_path / 'statepoint.3.h5' + _, source, tags = read_source(restart_file) + + # A restart run reads the tagged source bank back and tags the new sites + sp_path = run_model(model, tmp_path, restart_file=restart_file) + _, new_source, new_tags = read_source(sp_path) + np.testing.assert_array_equal(new_tags, bin_of(mesh, new_source['r'])) + + # The tags read from the restart file are those of the source sites + if not config['mpi']: + monkeypatch.chdir(tmp_path) + model.init_lib(restart_file=restart_file, output=False) + try: + openmc.lib.simulation_init() + bank = openmc.lib.source_bank() + np.testing.assert_array_equal(bank['birth_mesh_bin'], tags) + np.testing.assert_array_equal(bank['r'][:, 0], source['r']['x']) + openmc.lib.simulation_finalize() + finally: + model.finalize_lib() + + # A tag dataset whose length differs from that of the source bank is + # rejected, both on restart and for a file source. Without the length + # check, a dataset that is one entry too long would be read silently on + # restart, where each process reads only its share of the sites. + length_error = ('birth_mesh_bin dataset in the source file must be ' + 'one-dimensional with the same length as the source_bank') + initial_source = model.settings.source + for n_tags in (tags.size - 1, tags.size + 1): + bad_restart = tmp_path / f'statepoint.bad{n_tags}.h5' + shutil.copy(restart_file, bad_restart) + with h5py.File(bad_restart, 'r+') as f: + del f['birth_mesh_bin'] + f.create_dataset('birth_mesh_bin', + data=np.zeros(n_tags, dtype=np.int32)) + with pytest.raises(RuntimeError, match=length_error): + run_model(model, tmp_path, restart_file=bad_restart) + + bad_source = tmp_path / f'source.bad{n_tags}.h5' + with h5py.File(bad_source, 'w') as f: + f.attrs['filetype'] = np.bytes_('source') + with h5py.File(restart_file, 'r') as sp: + f.attrs['version'] = sp.attrs['version'] + f.create_dataset('source_bank', data=source) + f.create_dataset('birth_mesh_bin', + data=np.zeros(n_tags, dtype=np.int32)) + model.settings.source = openmc.FileSource(bad_source) + with pytest.raises(RuntimeError, match=length_error): + run_model(model, tmp_path) + + # Restore the initial source so that the next restart check reads + # only the restart file + model.settings.source = initial_source + + +def test_birth_mesh_fixed_source(tmp_path, capsys): + mesh = slab_mesh(2) + model, surface = make_model(tmp_path / 'two_slab.h5') + model.settings.run_mode = 'fixed source' + model.settings.batches = 2 + model.settings.birth_mesh = mesh + model.settings.surf_source_write = {'surface_ids': [surface.id], + 'max_particles': 1000} + sp_path = run_model(model, tmp_path) + + # The element is ignored with a warning, and no file carries tags + assert 'no effect in fixed-source mode' in capsys.readouterr().out + with h5py.File(sp_path, 'r') as f: + assert 'birth_mesh_id' not in f.attrs + assert 'birth_mesh_bin' not in f + with openmc.StatePoint(sp_path) as sp: + assert sp.birth_mesh_id is None + assert sp.birth_mesh_bin is None + file_attrs, _, file_tags = read_source(tmp_path / 'surface_source.h5') + assert 'birth_mesh_id' not in file_attrs + assert file_tags is None diff --git a/tests/unit_tests/test_settings.py b/tests/unit_tests/test_settings.py index 12195b876a1..3a0b856f761 100644 --- a/tests/unit_tests/test_settings.py +++ b/tests/unit_tests/test_settings.py @@ -58,6 +58,11 @@ def test_export_to_xml(run_in_tmpdir): s.trace = (10, 1, 20) s.track = [(1, 1, 1), (2, 1, 1)] s.ufs_mesh = mesh + birth_mesh = openmc.RegularMesh() + birth_mesh.lower_left = (-20., -20., -5.) + birth_mesh.upper_right = (20., 20., 5.) + birth_mesh.dimension = (4, 2, 1) + s.birth_mesh = birth_mesh s.resonance_scattering = {'enable': True, 'method': 'rvs', 'energy_min': 1.0, 'energy_max': 1000.0, 'nuclides': ['U235', 'U238', 'Pu239']} @@ -99,6 +104,10 @@ def test_export_to_xml(run_in_tmpdir): # Make sure exporting XML works s.export_to_xml() + # The birth mesh is referenced by ID and written to settings.xml + root = openmc.Settings.from_xml().to_xml_element() + assert root.find('birth_mesh').text == str(birth_mesh.id) + assert root.find(f"mesh[@id='{birth_mesh.id}']") is not None # Generate settings from XML s = openmc.Settings.from_xml() @@ -153,6 +162,11 @@ def test_export_to_xml(run_in_tmpdir): assert s.ufs_mesh.lower_left == [-10., -10., -10.] assert s.ufs_mesh.upper_right == [10., 10., 10.] assert s.ufs_mesh.dimension == (5, 5, 5) + assert isinstance(s.birth_mesh, openmc.RegularMesh) + assert s.birth_mesh.id == birth_mesh.id + assert s.birth_mesh.lower_left == [-20., -20., -5.] + assert s.birth_mesh.upper_right == [20., 20., 5.] + assert s.birth_mesh.dimension == (4, 2, 1) assert s.resonance_scattering == {'enable': True, 'method': 'rvs', 'energy_min': 1.0, 'energy_max': 1000.0, 'nuclides': ['U235', 'U238', 'Pu239']} @@ -194,6 +208,49 @@ def test_export_to_xml(run_in_tmpdir): assert s.free_gas_threshold == 800.0 +def test_birth_mesh(run_in_tmpdir): + mesh = openmc.RegularMesh() + mesh.lower_left = (-1., -1., -1.) + mesh.upper_right = (1., 1., 1.) + mesh.dimension = (2, 1, 1) + + s = openmc.Settings() + assert s.birth_mesh is None + with pytest.raises(TypeError): + s.birth_mesh = mesh.id + + # No element is written without a birth mesh + assert s.to_xml_element().find('birth_mesh') is None + + # A mesh shared with another setting is written to settings.xml once + s.birth_mesh = mesh + s.ufs_mesh = mesh + s.export_to_xml() + root = openmc.Settings.from_xml().to_xml_element() + assert len(root.findall('mesh')) == 1 + s = openmc.Settings.from_xml() + assert s.birth_mesh.id == mesh.id + assert s.birth_mesh is s.ufs_mesh + + # Any type of mesh can be used + cyl_mesh = openmc.CylindricalMesh([0., 1., 2.], [-1., 1.]) + s = openmc.Settings() + s.birth_mesh = cyl_mesh + s.export_to_xml() + s = openmc.Settings.from_xml() + assert isinstance(s.birth_mesh, openmc.CylindricalMesh) + assert s.birth_mesh.id == cyl_mesh.id + assert list(s.birth_mesh.r_grid) == [0., 1., 2.] + + # A reference to a mesh that is not defined is an error + elem = openmc.Settings().to_xml_element() + birth = elem.makeelement('birth_mesh', {}) + birth.text = '1234' + elem.append(birth) + with pytest.raises(ValueError, match='1234'): + openmc.Settings.from_xml_element(elem) + + def test_properties_file_load(tmp_path, mpi_intracomm): model = openmc.examples.pwr_assembly() diff --git a/tests/unit_tests/test_statepoint.py b/tests/unit_tests/test_statepoint.py index 51f5378c055..3c94e5aaebd 100644 --- a/tests/unit_tests/test_statepoint.py +++ b/tests/unit_tests/test_statepoint.py @@ -1,3 +1,5 @@ +import h5py +import numpy as np import openmc @@ -65,3 +67,45 @@ def test_get_tally_filter_type(run_in_tmpdir): tally_found = sp.get_tally(id=2) assert tally_found.id == 2 + + +def _write_tagged_file(path, source_present, tags=None): + """Write a minimal state point file with an optional birth tag dataset.""" + with h5py.File(path, 'w') as f: + f.attrs['filetype'] = np.bytes_('statepoint') + f.attrs['version'] = np.array( + [openmc.statepoint._VERSION_STATEPOINT, 0], dtype=np.int32) + f.attrs['source_present'] = np.int32(source_present) + if tags is not None: + f.attrs['birth_mesh_id'] = np.int32(3) + if source_present: + source = np.zeros(4, dtype=[('E', ' Date: Tue, 13 Oct 2026 23:05:18 -0400 Subject: [PATCH 2/7] Test birth-mesh tags with openmc.lib, MPI, any mesh and in CE Add tests of the fission bank birth-mesh tags: - test_source_site (Catch2), also run on two MPI ranks as test_source_site_mpi: the MPI datatype of SourceSite holds every byte of the struct, and sites sent from one rank to another keep every member. - test_birth_mesh_lib_source_bank: every field of openmc.lib.source_bank() equals the source bank and the tags written to the state point, and each tag is the bin of the source site from which the neutron that produced the site was started. - test_birth_mesh_any_mesh: in a problem of eight decoupled fuel islands, the tag of every site is the bin that a MeshFilter gives in its island, for a regular, a translated regular, a rectilinear, a cylindrical and a spherical mesh, and for a regular mesh with a translated filter. - birth_mesh_ce: the two-slab checks of birth_mesh in continuous-energy transport, with U-235 and H-1 slabs separated by B-10, and the comparison of the tags with the bins of the source sites from which the histories were started. test_birth_mesh_restart also checks that a restart run reads the tags without errors from HDF5, and test_birth_mesh_source_files that ParticleList.export_to_hdf5() writes the format version of the source files written by OpenMC. --- tests/cpp_unit_tests/CMakeLists.txt | 11 + tests/cpp_unit_tests/test_source_site.cpp | 177 +++++++++++++ tests/regression_tests/birth_mesh/test.py | 248 +++++++++++++++++- .../birth_mesh_ce/__init__.py | 0 tests/regression_tests/birth_mesh_ce/test.py | 177 +++++++++++++ 5 files changed, 611 insertions(+), 2 deletions(-) create mode 100644 tests/cpp_unit_tests/test_source_site.cpp create mode 100644 tests/regression_tests/birth_mesh_ce/__init__.py create mode 100644 tests/regression_tests/birth_mesh_ce/test.py diff --git a/tests/cpp_unit_tests/CMakeLists.txt b/tests/cpp_unit_tests/CMakeLists.txt index 82bcf319ca0..717f0df66b3 100644 --- a/tests/cpp_unit_tests/CMakeLists.txt +++ b/tests/cpp_unit_tests/CMakeLists.txt @@ -13,6 +13,7 @@ set(TEST_NAMES test_surface test_tensor test_geometry + test_source_site # Add additional unit test files here ) @@ -27,6 +28,16 @@ foreach(test ${TEST_NAMES}) set_tests_properties(${test} PROPERTIES SKIP_RETURN_CODE 4) endforeach() +# Also run the SourceSite test on two MPI ranks, so that sites are sent from +# one rank to another +if(OPENMC_USE_MPI) + # The ranks exchange sites, so they must run the test cases in the same order + add_test(NAME test_source_site_mpi + COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} 2 ${MPIEXEC_PREFLAGS} + $ ${MPIEXEC_POSTFLAGS} --order decl + WORKING_DIRECTORY ${UNIT_TEST_BIN_OUTPUT_DIR}) +endif() + if(OPENMC_USE_XDG) target_compile_definitions(test_mesh PRIVATE OPENMC_XDG_ENABLED) endif() diff --git a/tests/cpp_unit_tests/test_source_site.cpp b/tests/cpp_unit_tests/test_source_site.cpp new file mode 100644 index 00000000000..653a11700fd --- /dev/null +++ b/tests/cpp_unit_tests/test_source_site.cpp @@ -0,0 +1,177 @@ +#include +#include +#include + +#include + +#include "openmc/particle_data.h" +#include "openmc/vector.h" + +#ifdef OPENMC_MPI +#include +#include + +#include "openmc/initialize.h" +#include "openmc/message_passing.h" +#include "openmc/particle_type.h" +#endif + +using namespace openmc; + +#ifdef OPENMC_MPI + +namespace { + +// Set up MPI and the SourceSite datatype as openmc_init() does, once for the +// whole test run +class MPIListener : public Catch::EventListenerBase { +public: + using Catch::EventListenerBase::EventListenerBase; + + void testRunStarting(const Catch::TestRunInfo&) override + { + initialize_mpi(MPI_COMM_WORLD); + } + + void testRunEnded(const Catch::TestRunStats&) override + { + MPI_Type_free(&mpi::source_site); + MPI_Type_free(&mpi::collision_track_site); + MPI_Finalize(); + } +}; + +// Number of sites that each rank sends, and number of extra sites after them +// in each buffer, which must not be written to +constexpr int N_SITES {4}; +constexpr int N_EXTRA {2}; + +// Site i of the given rank. Every byte of the site is first given a value +// that depends on the rank, the site and the position of the byte, so that +// members not set below also differ from site to site. Every member is then +// set to a value that differs from its default value and from the value of +// the same member in every other site. +SourceSite make_site(int rank, int i) +{ + unsigned char bytes[sizeof(SourceSite)]; + for (std::size_t k = 0; k < sizeof(SourceSite); ++k) { + bytes[k] = static_cast(1 + k + 37 * i + 101 * rank); + } + SourceSite site; + std::memcpy(static_cast(&site), bytes, sizeof(SourceSite)); + + int n = 100 * rank + 10 * i; + double x = n; + site.r = {x + 0.125, x + 0.25, x + 0.375}; + site.u = {x + 0.5, x + 0.625, x + 0.75}; + site.E = 1.0e6 + x; + site.time = 1.0e-6 * (x + 1.0); + site.wgt = x + 1.5; + site.delayed_group = n + 1; + site.surf_id = n + 2; + site.particle = ParticleType {static_cast(n + 3)}; + site.parent_nuclide = n + 4; + site.parent_id = (int64_t {1} << 40) + n + 5; + site.progeny_id = (int64_t {1} << 41) + n + 6; + site.wgt_born = x + 2.5; + site.wgt_ww_born = x + 3.5; + site.n_split = (int64_t {1} << 42) + n + 7; + site.n_collision = n + 8; + site.birth_mesh_bin = n + 9; + return site; +} + +// Fill every byte of a buffer of sites with the same value +void fill_bytes(vector& sites, unsigned char value) +{ + std::memset( + static_cast(sites.data()), value, sites.size() * sizeof(SourceSite)); +} + +} // namespace + +CATCH_REGISTER_LISTENER(MPIListener) + +TEST_CASE("The MPI datatype of SourceSite covers every byte of the struct") +{ + // Banks are sent as arrays of SourceSite, so the datatype must have the + // extent of the struct + MPI_Aint lb, extent; + MPI_Type_get_extent(mpi::source_site, &lb, &extent); + REQUIRE(lb == 0); + REQUIRE(extent == static_cast(sizeof(SourceSite))); + + // SourceSite has no padding: every byte belongs to a member. A datatype that + // leaves out a member or part of one, or that gives a member a type of + // another size, therefore does not hold as many bytes as the struct. + MPI_Aint true_lb, true_extent; + MPI_Type_get_true_extent(mpi::source_site, &true_lb, &true_extent); + REQUIRE(true_lb == 0); + REQUIRE(true_extent == static_cast(sizeof(SourceSite))); + int size; + MPI_Type_size(mpi::source_site, &size); + REQUIRE(size == static_cast(sizeof(SourceSite))); +} + +TEST_CASE("Source sites sent between ranks keep every member") +{ + // Each rank sends its sites to the next rank and receives the sites of the + // previous one. With one rank, a rank sends its sites to itself. + int next = (mpi::rank + 1) % mpi::n_procs; + int previous = (mpi::rank + mpi::n_procs - 1) % mpi::n_procs; + + vector sent(N_SITES + N_EXTRA); + vector received(N_SITES + N_EXTRA); + fill_bytes(sent, 0xA5); + fill_bytes(received, 0x5A); + for (int i = 0; i < N_SITES; ++i) { + sent[i] = make_site(mpi::rank, i); + } + + MPI_Sendrecv(sent.data(), N_SITES, mpi::source_site, next, 0, received.data(), + N_SITES, mpi::source_site, previous, 0, mpi::intracomm, MPI_STATUS_IGNORE); + + for (int i = 0; i < N_SITES; ++i) { + const SourceSite expected = make_site(previous, i); + const SourceSite& site = received[i]; + REQUIRE(site.r.x == expected.r.x); + REQUIRE(site.r.y == expected.r.y); + REQUIRE(site.r.z == expected.r.z); + REQUIRE(site.u.x == expected.u.x); + REQUIRE(site.u.y == expected.u.y); + REQUIRE(site.u.z == expected.u.z); + REQUIRE(site.E == expected.E); + REQUIRE(site.time == expected.time); + REQUIRE(site.wgt == expected.wgt); + REQUIRE(site.delayed_group == expected.delayed_group); + REQUIRE(site.surf_id == expected.surf_id); + REQUIRE(site.particle.pdg_number() == expected.particle.pdg_number()); + REQUIRE(site.parent_nuclide == expected.parent_nuclide); + REQUIRE(site.parent_id == expected.parent_id); + REQUIRE(site.progeny_id == expected.progeny_id); + REQUIRE(site.wgt_born == expected.wgt_born); + REQUIRE(site.wgt_ww_born == expected.wgt_ww_born); + REQUIRE(site.n_split == expected.n_split); + REQUIRE(site.n_collision == expected.n_collision); + REQUIRE(site.birth_mesh_bin == expected.birth_mesh_bin); + + // Since SourceSite has no padding, equal bytes mean equal members, also + // for any member that is not compared above + REQUIRE(std::memcmp(&site, &expected, sizeof(SourceSite)) == 0); + } + + // Nothing was written past the received sites + vector untouched(N_EXTRA); + fill_bytes(untouched, 0x5A); + REQUIRE(std::memcmp(&received[N_SITES], untouched.data(), + N_EXTRA * sizeof(SourceSite)) == 0); +} + +#else + +TEST_CASE("Source sites sent between ranks keep every member") +{ + SKIP("OpenMC was built without MPI"); +} + +#endif diff --git a/tests/regression_tests/birth_mesh/test.py b/tests/regression_tests/birth_mesh/test.py index 46675810def..8ea35da6760 100644 --- a/tests/regression_tests/birth_mesh/test.py +++ b/tests/regression_tests/birth_mesh/test.py @@ -17,6 +17,22 @@ reads them back, and the rejection of a tag dataset whose length differs from that of the source bank. +Two more tests compare the tags with references. The first runs the two-slab +problem through openmc.lib and compares every field of openmc.lib.source_bank() +with the source bank and the tags written to the state point, which checks the +ctypes mirror of SourceSite. It also compares each tag with the bin of the +source site from which the neutron that produced the site was started. These +data create no secondary neutrons, so that source site is the birth position +of the neutron. The second test places eight fuel islands, at least 18 cm (180 +mean free paths of the absorber) apart, in the absorber, so that the birth +position of the neutron that produced a site lies in the island of the site. +Each island lies within one element of each birth mesh, or outside the mesh, +so the tag of a site must be the bin that a MeshFilter on the mesh gives in +that island. A tally with a CellFilter on the islands and that MeshFilter +provides this bin, for a regular mesh, a translated regular mesh, a regular +mesh with a translated filter, a rectilinear mesh, a cylindrical mesh and a +spherical mesh. + """ import shutil @@ -125,6 +141,18 @@ def bin_of(mesh, r): return ijk[:, 0] + nx * (ijk[:, 1] + ny * ijk[:, 2]) +def regular_bin(mesh, xyz): + """0-based bin (x fastest) of each position of xyz, an array of shape + (n, 3), in a regular mesh, or -1 for a position outside the mesh.""" + ll = np.asarray(mesh.lower_left) + dimension = np.asarray(mesh.dimension) + width = (np.asarray(mesh.upper_right) - ll) / dimension + ijk = np.floor((xyz - ll) / width).astype(int) + inside = np.all((ijk >= 0) & (ijk < dimension), axis=1) + bins = ijk[:, 0] + dimension[0] * (ijk[:, 1] + dimension[1] * ijk[:, 2]) + return np.where(inside, bins, -1) + + def run_model(model, path, **kwargs): kwargs.update({'cwd': path, 'openmc_exec': config['exe'], 'event_based': config['event']}) @@ -213,6 +241,15 @@ def test_birth_mesh_source_files(tmp_path): np.testing.assert_array_equal(file_source, source) np.testing.assert_array_equal(file_tags, tags) + # A source file written by ParticleList has the format version of the + # source files written by OpenMC + particles = openmc.ParticleList([openmc.SourceParticle()]) + particles.export_to_hdf5(tmp_path / 'particles.h5') + particle_attrs, _, _ = read_source(tmp_path / 'particles.h5') + file_attrs, _, _ = read_source(tmp_path / 'source.h5') + np.testing.assert_array_equal(particle_attrs['version'], + file_attrs['version']) + @pytest.mark.parametrize('rectilinear', [False, True]) def test_birth_mesh_birth_position(tmp_path, rectilinear): @@ -235,7 +272,62 @@ def test_birth_mesh_birth_position(tmp_path, rectilinear): assert np.count_nonzero(tags != own_bin) > 0 -def test_birth_mesh_restart(tmp_path, monkeypatch): +def test_birth_mesh_lib_source_bank(tmp_path, monkeypatch): + # openmc.lib reads the source bank through a ctypes mirror of SourceSite, + # so a member at a wrong offset or of a wrong type shows up as a mismatch + # with the source bank and the tags written to the state point + if config['mpi']: + pytest.skip('openmc.lib runs in this process, without MPI') + + # Bins of 2.5 cm along x and 5 cm along y. The outer half of the right slab + # is outside the mesh. + x_max = HALF_WIDTH + FUEL_WIDTH + mesh = openmc.RegularMesh() + mesh.lower_left = (-x_max, -5.0, -5.0) + mesh.upper_right = (HALF_WIDTH + FUEL_WIDTH / 2, 5.0, 5.0) + mesh.dimension = (7, 2, 1) + model, _ = eigenvalue_model(tmp_path, mesh) + batches = model.settings.batches + + monkeypatch.chdir(tmp_path) + model.init_lib(output=False) + try: + openmc.lib.simulation_init() + for _ in range(batches - 1): + openmc.lib.next_batch() + # Source sites of the last batch, then the sites that they produced + parents = openmc.lib.source_bank().copy() + openmc.lib.next_batch() + bank = openmc.lib.source_bank().copy() + openmc.lib.simulation_finalize() + finally: + model.finalize_lib() + _, source, tags = read_source(tmp_path / f'statepoint.{batches}.h5') + + # The source bank after the last batch is the one in the state point + assert bank.shape == source.shape == tags.shape + for name in ('r', 'u'): + xyz = np.column_stack([source[name][c] for c in 'xyz']) + np.testing.assert_array_equal(bank[name], xyz) + for name in ('E', 'time', 'wgt', 'delayed_group', 'surf_id', 'particle'): + np.testing.assert_array_equal(bank[name], source[name]) + np.testing.assert_array_equal(bank['birth_mesh_bin'], tags) + # Sites outside the mesh and in several bins + assert -1 in tags + assert np.count_nonzero(np.unique(tags) >= 0) >= 2 + + # parent_id is the index of the source site from which the neutron that + # produced a site was started. That neutron was born there, since no + # secondary neutrons are created with these data, so the tag is the bin + # of that source site, which often differs from the bin of the site. + parent = bank['parent_id'] + assert parent.min() >= 0 and parent.max() < parents.size + birth_bin = regular_bin(mesh, parents['r'][parent]) + np.testing.assert_array_equal(tags, birth_bin) + assert np.count_nonzero(tags != regular_bin(mesh, bank['r'])) > 0 + + +def test_birth_mesh_restart(tmp_path, monkeypatch, capfd): mesh = slab_mesh(2) model, _ = eigenvalue_model(tmp_path, mesh) model.settings.statepoint = {'batches': [3, 6]} @@ -243,8 +335,11 @@ def test_birth_mesh_restart(tmp_path, monkeypatch): restart_file = tmp_path / 'statepoint.3.h5' _, source, tags = read_source(restart_file) - # A restart run reads the tagged source bank back and tags the new sites + # A restart run reads the tagged source bank back, without errors from + # HDF5, and tags the new sites + capfd.readouterr() sp_path = run_model(model, tmp_path, restart_file=restart_file) + assert 'HDF5-DIAG' not in capfd.readouterr().out _, new_source, new_tags = read_source(sp_path) np.testing.assert_array_equal(new_tags, bin_of(mesh, new_source['r'])) @@ -316,3 +411,152 @@ def test_birth_mesh_fixed_source(tmp_path, capsys): file_attrs, _, file_tags = read_source(tmp_path / 'surface_source.h5') assert 'birth_mesh_id' not in file_attrs assert file_tags is None + + +# Eight fuel islands, x in s_x [a, a + 10] and y in s_y [a, a + 10] cm for a in +# ISLAND_OFFSETS and s_x, s_y = +1 or -1, in the absorber of the two-slab +# problem, which fills the rest of [-60, 60]^2 x [-5, 5] cm. Any two islands +# are at least 18 cm apart. +ISLAND_OFFSETS = (12.0, 40.0) +ISLAND_WIDTH = 10.0 +ISLAND_SPAN = 60.0 +TRANSLATION = (30.0, -30.0, 0.0) +ANY_MESH_CASES = ('regular', 'translated', 'translated_filter', 'rectilinear', + 'cylindrical', 'spherical') + + +def island_boxes(): + """x and y ranges of the islands.""" + boxes = [] + for a in ISLAND_OFFSETS: + for sx in (1, -1): + for sy in (1, -1): + x = sorted((sx * a, sx * (a + ISLAND_WIDTH))) + y = sorted((sy * a, sy * (a + ISLAND_WIDTH))) + boxes.append((x, y)) + return boxes + + +def any_mesh_case(case): + """Birth mesh, mesh of the reference MeshFilter, the translation of that + filter, and the islands outside the birth mesh.""" + span = ISLAND_SPAN + t = np.asarray(TRANSLATION) + ll = np.array((-span, -span, -5.0)) + ur = np.array((span, span, 5.0)) + phi = np.linspace(0.0, 2.0 * np.pi, 5) + + def regular(shift): + mesh = openmc.RegularMesh() + mesh.lower_left = ll + shift + mesh.upper_right = ur + shift + mesh.dimension = (4, 4, 1) + return mesh + + if case == 'regular': + mesh = regular(0.0) + elif case in ('translated', 'translated_filter'): + # The regular mesh moved by TRANSLATION, which leaves the outer islands + # with x < 0 or y > 0 outside + mesh = regular(t) + outside = {k for k, (x, y) in enumerate(island_boxes()) + if min(abs(x[0]), abs(x[1])) == ISLAND_OFFSETS[-1] + and (x[1] < 0.0 or y[0] > 0.0)} + if case == 'translated': + return mesh, mesh, None, outside + # A MeshFilter with a translation on the unmoved mesh bins positions + # as the moved mesh does + return mesh, regular(0.0), TRANSLATION, outside + elif case == 'rectilinear': + mesh = openmc.RectilinearMesh() + mesh.x_grid = [-60.0, -31.0, -5.0, 7.0, 33.0, 60.0] + mesh.y_grid = [-60.0, -25.0, 3.0, 60.0] + mesh.z_grid = [-5.0, 5.0] + elif case == 'cylindrical': + mesh = openmc.CylindricalMesh( + r_grid=[0.0, 36.0, 80.0], z_grid=[-5.0, 5.0], phi_grid=phi, + origin=(3.0, -2.0, 0.0)) + elif case == 'spherical': + mesh = openmc.SphericalMesh( + r_grid=[0.0, 36.0, 80.0], phi_grid=phi, theta_grid=[0.0, np.pi], + origin=(3.0, -2.0, 0.0)) + return mesh, mesh, None, set() + + +def island_model(tmp_path, birth_mesh, filter_mesh, translation): + """Eigenvalue model of the islands with a reference tally that gives the + MeshFilter bins of the collisions in each island.""" + model, _ = make_model(tmp_path / 'islands.h5') + fuel, absorber = model.materials + span = ISLAND_SPAN + z0 = openmc.ZPlane(-5.0, boundary_type='reflective') + z1 = openmc.ZPlane(5.0, boundary_type='reflective') + outer = (+openmc.XPlane(-span, boundary_type='vacuum') + & -openmc.XPlane(span, boundary_type='vacuum') + & +openmc.YPlane(-span, boundary_type='vacuum') + & -openmc.YPlane(span, boundary_type='vacuum') & +z0 & -z1) + islands = [] + absorber_region = outer + for x, y in island_boxes(): + region = (+openmc.XPlane(x[0]) & -openmc.XPlane(x[1]) + & +openmc.YPlane(y[0]) & -openmc.YPlane(y[1]) & +z0 & -z1) + islands.append(openmc.Cell(fill=fuel, region=region)) + absorber_region &= ~region + model.geometry = openmc.Geometry( + islands + [openmc.Cell(fill=absorber, region=absorber_region)]) + + model.settings.run_mode = 'eigenvalue' + model.settings.batches = 6 + model.settings.inactive = 2 + model.settings.source = openmc.IndependentSource( + space=openmc.stats.Box((-span, -span, -5.0), (span, span, 5.0)), + constraints={'fissionable': True}) + # The islands fill only 1/18 of the source box + model.settings.source_rejection_fraction = 0.01 + model.settings.birth_mesh = birth_mesh + + mesh_filter = openmc.MeshFilter(filter_mesh) + if translation is not None: + mesh_filter.translation = translation + tally = openmc.Tally() + tally.filters = [openmc.CellFilter(islands), mesh_filter] + tally.scores = ['flux'] + tally.estimator = 'collision' + model.tallies = [tally] + return model, tally + + +@pytest.mark.parametrize('case', ANY_MESH_CASES) +def test_birth_mesh_any_mesh(tmp_path, case): + birth_mesh, filter_mesh, translation, outside = any_mesh_case(case) + model, tally = island_model(tmp_path, birth_mesh, filter_mesh, translation) + sp_path = run_model(model, tmp_path) + boxes = island_boxes() + with openmc.StatePoint(sp_path) as sp: + flux = sp.tallies[tally.id].mean.reshape(len(boxes), -1) + source = sp.source + tags = sp.birth_mesh_bin + + # The collisions in an island are all in one bin of the filter, which is + # the bin of the island, or in none for an island outside the mesh. The + # islands inside the mesh are in distinct bins. + island_bin = np.full(len(boxes), -1) + for k in range(len(boxes)): + nonzero = np.flatnonzero(flux[k]) + assert nonzero.size == (0 if k in outside else 1) + if nonzero.size == 1: + island_bin[k] = nonzero[0] + inside = island_bin[island_bin >= 0] + assert np.unique(inside).size == inside.size == len(boxes) - len(outside) + + # Every site lies in one island, every island holds sites, and the tag of + # each site is the bin of its island + assert tags.dtype == np.int32 + assert tags.shape == source.shape + x, y = source['r']['x'], source['r']['y'] + in_box = np.array([(x >= bx[0]) & (x <= bx[1]) + & (y >= by[0]) & (y <= by[1]) for bx, by in boxes]) + assert np.all(in_box.sum(axis=0) == 1) + island = in_box.argmax(axis=0) + assert np.unique(island).size == len(boxes) + np.testing.assert_array_equal(tags, island_bin[island]) diff --git a/tests/regression_tests/birth_mesh_ce/__init__.py b/tests/regression_tests/birth_mesh_ce/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/birth_mesh_ce/test.py b/tests/regression_tests/birth_mesh_ce/test.py new file mode 100644 index 00000000000..35e1e68c030 --- /dev/null +++ b/tests/regression_tests/birth_mesh_ce/test.py @@ -0,0 +1,177 @@ +"""Test of the fission bank birth-mesh tags in continuous-energy transport. + +Two slabs of a mixture of U-235 and H-1, 10 cm thick, are separated by 10 cm +of B-10 at 50 g/cm^3. The total cross section of B-10 is at least 1.4 b below +20 MeV, so the layer is more than 40 mean free paths thick at every energy. A +neutron that scatters in the layer slows down and is absorbed: at 1 keV, its +mean free path for absorption is about 0.02 cm. A neutron born in one slab +therefore causes fissions only in that slab, and the tags of the banked +fission sites can be checked without reference results, as in the multigroup +birth_mesh test: + +- with one mesh bin per slab, every tag is the bin of its own site; +- with two bins per slab along x (and two along y), some tags differ from the + bin of the site, but none leaves the slab of the site; +- with a mesh over the left half of the problem, the sites in the right slab + are tagged -1. + +A last test runs the problem through openmc.lib and compares each tag with the +bin of the source site from which the neutron that produced the site was +started, its birth position. + +""" +import numpy as np +import openmc +import openmc.lib +import pytest + +from tests.regression_tests import config + +FUEL_WIDTH = 10.0 # cm, thickness of each fuel slab +HALF_GAP = 5.0 # cm, half-thickness of the B-10 layer +X_MAX = HALF_GAP + FUEL_WIDTH + + +def make_model(mesh): + fuel = openmc.Material(name='fuel') + fuel.add_nuclide('U235', 2.2e-4) + fuel.add_nuclide('H1', 6.6e-2) + fuel.set_density('sum') + absorber = openmc.Material(name='absorber') + absorber.add_nuclide('B10', 1.0) + absorber.set_density('g/cm3', 50.0) + + x0 = openmc.XPlane(-X_MAX, boundary_type='reflective') + x1 = openmc.XPlane(-HALF_GAP) + x2 = openmc.XPlane(HALF_GAP) + x3 = openmc.XPlane(X_MAX, boundary_type='reflective') + y0 = openmc.YPlane(-5.0, boundary_type='reflective') + y1 = openmc.YPlane(5.0, boundary_type='reflective') + z0 = openmc.ZPlane(-5.0, boundary_type='reflective') + z1 = openmc.ZPlane(5.0, boundary_type='reflective') + yz = +y0 & -y1 & +z0 & -z1 + cells = [openmc.Cell(fill=fuel, region=+x0 & -x1 & yz), + openmc.Cell(fill=absorber, region=+x1 & -x2 & yz), + openmc.Cell(fill=fuel, region=+x2 & -x3 & yz)] + + model = openmc.Model() + model.materials = openmc.Materials([fuel, absorber]) + model.geometry = openmc.Geometry(cells) + model.settings.run_mode = 'eigenvalue' + model.settings.particles = 1000 + model.settings.batches = 6 + model.settings.inactive = 2 + model.settings.source = openmc.IndependentSource( + space=openmc.stats.Box((-X_MAX, -5.0, -5.0), (X_MAX, 5.0, 5.0)), + constraints={'fissionable': True}) + model.settings.birth_mesh = mesh + return model + + +def regular_mesh(x_max, dimension): + mesh = openmc.RegularMesh() + mesh.lower_left = (-X_MAX, -5.0, -5.0) + mesh.upper_right = (x_max, 5.0, 5.0) + mesh.dimension = dimension + return mesh + + +def regular_bin(mesh, xyz): + """0-based bin (x fastest) of each position of xyz, an array of shape + (n, 3), in a regular mesh, or -1 for a position outside the mesh.""" + ll = np.asarray(mesh.lower_left) + dimension = np.asarray(mesh.dimension) + width = (np.asarray(mesh.upper_right) - ll) / dimension + ijk = np.floor((xyz - ll) / width).astype(int) + inside = np.all((ijk >= 0) & (ijk < dimension), axis=1) + bins = ijk[:, 0] + dimension[0] * (ijk[:, 1] + dimension[1] * ijk[:, 2]) + return np.where(inside, bins, -1) + + +def run_tagged(tmp_path, mesh): + """Run the eigenvalue problem and return the final source bank, the tags, + the bin of each site, and whether each site is in the left slab.""" + model = make_model(mesh) + kwargs = {'cwd': tmp_path, 'openmc_exec': config['exe'], + 'event_based': config['event']} + if config['mpi']: + kwargs['mpi_args'] = [config['mpiexec'], '-n', config['mpi_np']] + sp_path = model.run(**kwargs) + with openmc.StatePoint(sp_path) as sp: + source = sp.source + tags = sp.birth_mesh_bin + assert tags.dtype == np.int32 + assert tags.shape == source.shape == (model.settings.particles,) + + xyz = np.column_stack([source['r'][c] for c in 'xyz']) + left = xyz[:, 0] < -HALF_GAP + right = xyz[:, 0] > HALF_GAP + # All sites are in the fuel, and both slabs hold sites + assert np.all(left | right) + assert np.any(left) and np.any(right) + return tags, regular_bin(mesh, xyz), left + + +def test_birth_mesh_ce_one_bin_per_slab(tmp_path): + # Bin 0 holds the left slab and bin 1 the right one + mesh = regular_mesh(X_MAX, (2, 1, 1)) + tags, own_bin, left = run_tagged(tmp_path, mesh) + np.testing.assert_array_equal(own_bin, np.where(left, 0, 1)) + np.testing.assert_array_equal(tags, own_bin) + + +def test_birth_mesh_ce_two_bins_per_slab(tmp_path): + # Bins of 5 cm along x and y: x indices 0 and 1 in the left slab, 4 and 5 + # in the right slab + nx = 6 + mesh = regular_mesh(X_MAX, (nx, 2, 1)) + tags, own_bin, left = run_tagged(tmp_path, mesh) + assert tags.min() >= 0 + for in_slab, slab in ((left, {0, 1}), (~left, {nx - 2, nx - 1})): + assert set(np.unique(own_bin[in_slab] % nx)) == slab + assert set(np.unique(tags[in_slab] % nx)) <= slab + assert set(np.unique(tags // nx)) == {0, 1} + assert np.count_nonzero(tags != own_bin) > 0 + + +def test_birth_mesh_ce_left_half(tmp_path): + # One bin over the left half of the problem; the right slab is outside + mesh = regular_mesh(0.0, (1, 1, 1)) + tags, _, left = run_tagged(tmp_path, mesh) + np.testing.assert_array_equal(tags, np.where(left, 0, -1)) + + +def test_birth_mesh_ce_birth_position(tmp_path, monkeypatch): + if config['mpi']: + pytest.skip('openmc.lib runs in this process, without MPI') + mesh = regular_mesh(X_MAX, (6, 2, 1)) + model = make_model(mesh) + monkeypatch.chdir(tmp_path) + model.init_lib(output=False) + try: + openmc.lib.simulation_init() + for _ in range(model.settings.batches - 1): + openmc.lib.next_batch() + # Source sites of the last batch, then the sites that they produced + parents = openmc.lib.source_bank().copy() + openmc.lib.next_batch() + bank = openmc.lib.source_bank().copy() + openmc.lib.simulation_finalize() + finally: + model.finalize_lib() + + # parent_id is the index of the source site from which the history of the + # neutron that produced a site was started. That neutron was born at the + # source site, unless it was created by an (n,2n) or similar reaction in + # that history, where it was born at the point of the reaction. Such + # neutrons cause few of the fissions here, of the order of one in a + # thousand, so at least 99% of the tags are the bin of the source site, + # while the bin of the site itself differs in many of them. + parent = bank['parent_id'] + assert parent.min() >= 0 and parent.max() < parents.size + tags = bank['birth_mesh_bin'] + n_sites = tags.size + birth_bin = regular_bin(mesh, parents['r'][parent]) + own_bin = regular_bin(mesh, bank['r']) + assert np.count_nonzero(tags != birth_bin) <= 0.01 * n_sites + assert np.count_nonzero(tags != own_bin) > 0.1 * n_sites From 35a001201e58fd7c1b58a23f01acffe8154bd20e Mon Sep 17 00:00:00 2001 From: Ondrej Chvala Date: Wed, 14 Oct 2026 02:26:15 -0400 Subject: [PATCH 3/7] Test the birth-mesh tags of sites read from a source file Add test_birth_mesh_ce_file_source: a tagged continuous-energy run writes source.h5, and an eigenvalue run with a FileSource on that file is initialized through openmc.lib. Each site sampled for the first batch must have the tag of the file site that it copies, which checks that read_source_bank() copies the tags into the sites when the whole bank is read. Like the other openmc.lib tests, it is skipped under MPI. --- tests/regression_tests/birth_mesh_ce/test.py | 42 ++++++++++++++++++-- 1 file changed, 39 insertions(+), 3 deletions(-) diff --git a/tests/regression_tests/birth_mesh_ce/test.py b/tests/regression_tests/birth_mesh_ce/test.py index 35e1e68c030..be42f75905d 100644 --- a/tests/regression_tests/birth_mesh_ce/test.py +++ b/tests/regression_tests/birth_mesh_ce/test.py @@ -15,11 +15,13 @@ - with a mesh over the left half of the problem, the sites in the right slab are tagged -1. -A last test runs the problem through openmc.lib and compares each tag with the -bin of the source site from which the neutron that produced the site was -started, its birth position. +Two last tests run the problem through openmc.lib. The first compares each tag +with the bin of the source site from which the neutron that produced the site +was started, its birth position. The second checks that the sites sampled from +a tagged source file for the first batch carry the tags of the file. """ +import h5py import numpy as np import openmc import openmc.lib @@ -175,3 +177,37 @@ def test_birth_mesh_ce_birth_position(tmp_path, monkeypatch): own_bin = regular_bin(mesh, bank['r']) assert np.count_nonzero(tags != birth_bin) <= 0.01 * n_sites assert np.count_nonzero(tags != own_bin) > 0.1 * n_sites + + +def test_birth_mesh_ce_file_source(tmp_path, monkeypatch): + if config['mpi']: + pytest.skip('openmc.lib runs in this process, without MPI') + # A tagged run writes its final source bank and tags to source.h5 + mesh = regular_mesh(X_MAX, (2, 1, 1)) + model = make_model(mesh) + model.settings.sourcepoint = {'separate': True, 'overwrite': True} + model.run(cwd=tmp_path, openmc_exec=config['exe'], + event_based=config['event']) + source_file = tmp_path / 'source.h5' + with h5py.File(source_file, 'r') as f: + file_sites = f['source_bank'][()] + file_tags = f['birth_mesh_bin'][()] + assert set(np.unique(file_tags)) == {0, 1} + file_tag = {tuple(r): tag for r, tag in zip(file_sites['r'], file_tags)} + + # A FileSource on that file reads the tags with the sites, so each site + # sampled for the first batch has the tag of the file site that it copies + model = make_model(mesh) + model.settings.source = openmc.FileSource(source_file) + run_dir = tmp_path / 'file_source' + run_dir.mkdir() + monkeypatch.chdir(run_dir) + model.init_lib(output=False) + try: + openmc.lib.simulation_init() + bank = openmc.lib.source_bank().copy() + openmc.lib.simulation_finalize() + finally: + model.finalize_lib() + expected = [file_tag[tuple(r)] for r in bank['r']] + np.testing.assert_array_equal(bank['birth_mesh_bin'], expected) From e9ece1953b589eb41cbe307a368ecd11e2bf078d Mon Sep 17 00:00:00 2001 From: Ondrej Chvala Date: Wed, 14 Oct 2026 08:10:38 -0400 Subject: [PATCH 4/7] Ignore the birth mesh with random ray and warn for MCPL source points The random ray solver banks no fission sites, so a element tagged nothing in a random ray eigenvalue run, while its state points still recorded a birth_mesh_id. The element is now ignored with a warning, as in fixed-source mode, and no birth_mesh_id is written. MCPL source point files have no place for the tags, which were dropped without notice. When an eigenvalue run with a birth mesh writes its source points to separate MCPL files, a warning now says so. It is printed once, by the master process, when the settings are read, and not at every source point written. --- docs/source/io_formats/settings.rst | 7 ++++--- docs/source/usersguide/settings.rst | 4 ++-- include/openmc/settings.h | 3 ++- openmc/settings.py | 2 +- src/settings.cpp | 14 +++++++++++++- 5 files changed, 22 insertions(+), 8 deletions(-) diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 6fbbc1fd51c..18ddf79edcf 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -43,9 +43,10 @@ fission-matrix data. When this element is present, every site banked in the fission bank, in continuous-energy and in multigroup mode, is tagged with the bin of the mesh that contains the *birth position of the neutron that produced the site*, or with ``-1`` if that position is outside the mesh. Only an -eigenvalue calculation banks fission sites, so in fixed-source mode the element -has no effect and is ignored with a warning. Tagging consumes no random numbers -and changes no transport or tally result. +eigenvalue calculation with the Monte Carlo solver banks fission sites, so in +fixed-source mode and with the random ray solver the element has no effect and +is ignored with a warning. Tagging consumes no random numbers and changes no +transport or tally result. The birth position is the position at which the neutron was started from a site: its fission site for a neutron started from the source bank, the diff --git a/docs/source/usersguide/settings.rst b/docs/source/usersguide/settings.rst index 8363e34e7a3..72ad199d836 100644 --- a/docs/source/usersguide/settings.rst +++ b/docs/source/usersguide/settings.rst @@ -727,8 +727,8 @@ of :class:`openmc.MeshBornFilter`), or with -1 if that position is outside the mesh. No translation or rotation is applied, since these are attributes of a mesh filter, not of the mesh. The tags are written as a ``birth_mesh_bin`` dataset to the statepoint files that contain the source bank and to the source -files ``source..h5`` and ``source.h5``. In fixed-source mode, the -attribute has no effect. For example:: +files ``source..h5`` and ``source.h5``. In fixed-source mode and with the +random ray solver, the attribute has no effect. For example:: settings.birth_mesh = mesh diff --git a/include/openmc/settings.h b/include/openmc/settings.h index a4340969f2c..853aea41a73 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -51,7 +51,8 @@ class Mesh; namespace settings { // Mesh used to tag banked fission sites (see SourceSite::birth_mesh_bin); -// nullptr when no tagging is done, including in fixed-source mode +// nullptr when no tagging is done, including in fixed-source mode and with the +// random ray solver extern const Mesh* birth_mesh; // Boolean flags diff --git a/openmc/settings.py b/openmc/settings.py index 6ae0a20cec9..07b53f37c48 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -59,7 +59,7 @@ class Settings: :class:`openmc.MeshBornFilter`. The tags are written to state point and source point files as a ``birth_mesh_bin`` dataset aligned with the source bank (see :ref:`io_statepoint`). Has no effect in fixed-source - mode. + mode or with the random ray solver. .. versionadded:: 0.16.1 confidence_intervals : bool diff --git a/src/settings.cpp b/src/settings.cpp index 1f0809198b0..9e296640dfc 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -897,10 +897,14 @@ void read_settings_xml(pugi::xml_node root) } // Only an eigenvalue calculation banks fission sites in the fission bank, - // so birth tagging has no effect in fixed-source mode + // and the random ray solver banks none, so birth tagging has no effect in + // fixed-source mode or with the random ray solver if (run_mode == RunMode::FIXED_SOURCE) { warning("The element has no effect in fixed-source mode " "and is ignored."); + } else if (solver_type == SolverType::RANDOM_RAY) { + warning("The element has no effect with the random ray " + "solver and is ignored."); } else { auto& m = model::meshes[model::mesh_map.at(temp)]; m->prepare_for_point_location(); @@ -969,6 +973,14 @@ void read_settings_xml(pugi::xml_node root) sourcepoint_batch = statepoint_batch; } + // MCPL source point files have no place for the fission bank birth tags. + // Say so once, rather than at every source point written. + if (birth_mesh && source_write && source_separate && source_mcpl_write && + mpi::master) { + warning("The tags are not written to MCPL source point " + "files."); + } + // Check is the user specified to convert strength to statistical weight if (check_for_node(root, "uniform_source_sampling")) { uniform_source_sampling = From b912e7a7aee0701eea0c69b9f505e82dac25c4f7 Mon Sep 17 00:00:00 2001 From: Ondrej Chvala Date: Wed, 14 Oct 2026 08:28:34 -0400 Subject: [PATCH 5/7] Test the birth mesh with random ray and the MCPL source point warning Add test_birth_mesh_random_ray: with the random ray solver, a element is ignored with a warning, printed once by each process, and the state point records no birth_mesh_id and no tags. Add test_birth_mesh_mcpl_warning: the warning that MCPL source point files carry no tags is given exactly once, by the master process, for an eigenvalue run that writes its source points to separate MCPL files, and not for HDF5 source points, for source points kept in the state point, when source points are not written, in fixed-source mode or with the random ray solver. Source points are requested at a batch that the run does not reach, so that these cases need no MCPL library. A case that writes MCPL source points at two batches is skipped without mcpl-config, as the other MCPL tests are. --- tests/regression_tests/birth_mesh/test.py | 98 ++++++++++++++++++++++- 1 file changed, 97 insertions(+), 1 deletion(-) diff --git a/tests/regression_tests/birth_mesh/test.py b/tests/regression_tests/birth_mesh/test.py index 8ea35da6760..d1c6b02f080 100644 --- a/tests/regression_tests/birth_mesh/test.py +++ b/tests/regression_tests/birth_mesh/test.py @@ -15,7 +15,9 @@ The tests also check which output files carry the tags, the restart path that reads them back, and the rejection of a tag dataset whose length differs from -that of the source bank. +that of the source bank. With the random ray solver, which banks no fission +sites, the element is ignored with a warning, and a warning says that source +point files written in MCPL format carry no tags. Two more tests compare the tags with references. The first runs the two-slab problem through openmc.lib and compares every field of openmc.lib.source_bank() @@ -413,6 +415,100 @@ def test_birth_mesh_fixed_source(tmp_path, capsys): assert file_tags is None +# Fragments of the warnings about the element. Warnings are +# wrapped at 80 columns, and each fragment lies within the first line. +FIXED_SOURCE_WARNING = 'no effect in fixed-source mode' +RANDOM_RAY_WARNING = 'no effect with the random ray solver' +MCPL_WARNING = 'tags are not written to MCPL source point files' + + +def random_ray_model(tmp_path, mesh): + """The two-slab eigenvalue problem solved with the random ray solver.""" + model, _ = eigenvalue_model(tmp_path, mesh) + x_max = HALF_WIDTH + FUEL_WIDTH + model.settings.particles = 200 + model.settings.random_ray = { + 'distance_inactive': 10.0, + 'distance_active': 50.0, + 'ray_source': openmc.IndependentSource( + space=openmc.stats.Box((-x_max, -5.0, -5.0), (x_max, 5.0, 5.0))) + } + return model + + +def n_processes(): + return int(config['mpi_np']) if config['mpi'] else 1 + + +def test_birth_mesh_random_ray(tmp_path, capfd): + # The random ray solver banks no fission sites, so the element is ignored + # with a warning, which each process prints, and no file records the mesh + model = random_ray_model(tmp_path, slab_mesh(2)) + sp_path = run_model(model, tmp_path) + out = capfd.readouterr().out + assert out.count(RANDOM_RAY_WARNING) == n_processes() + assert FIXED_SOURCE_WARNING not in out + assert MCPL_WARNING not in out + with h5py.File(sp_path, 'r') as f: + assert 'birth_mesh_id' not in f.attrs + assert 'birth_mesh_bin' not in f + with openmc.StatePoint(sp_path) as sp: + assert sp.birth_mesh_id is None + assert sp.birth_mesh_bin is None + + +# Source point settings of test_birth_mesh_mcpl_warning. In the cases 'mcpl', +# 'fixed_source' and 'random_ray', source points are requested in MCPL format +# at a batch that the run does not reach, so no MCPL file is written and the +# MCPL library is not needed. +MCPL_CASES = { + 'mcpl': {'batches': [7], 'separate': True, 'mcpl': True}, + 'mcpl_written': {'batches': [3, 6], 'separate': True, 'mcpl': True}, + 'hdf5': {'batches': [3, 6], 'separate': True}, + 'not_separate': {'batches': [6], 'separate': False, 'mcpl': True}, + 'not_written': {'write': False, 'separate': True, 'mcpl': True}, + 'fixed_source': {'batches': [7], 'separate': True, 'mcpl': True}, + 'random_ray': {'batches': [7], 'separate': True, 'mcpl': True}, +} + + +@pytest.mark.parametrize('case', MCPL_CASES) +def test_birth_mesh_mcpl_warning(tmp_path, capfd, case): + if case == 'mcpl_written' and shutil.which('mcpl-config') is None: + pytest.skip('MCPL is not available.') + mesh = slab_mesh(2) + if case == 'random_ray': + model = random_ray_model(tmp_path, mesh) + else: + model, _ = eigenvalue_model(tmp_path, mesh) + if case == 'fixed_source': + model.settings.run_mode = 'fixed source' + model.settings.batches = 2 + model.settings.sourcepoint = MCPL_CASES[case] + sp_path = run_model(model, tmp_path) + out = capfd.readouterr().out + + # The warning is given once, by the master process, when the settings are + # read, and only for an eigenvalue run with the Monte Carlo solver that + # writes its source points to separate MCPL files + n_warnings = 1 if case in ('mcpl', 'mcpl_written') else 0 + assert out.count(MCPL_WARNING) == n_warnings + if case == 'mcpl_written': + for batch in (3, 6): + assert (tmp_path / f'source.{batch}.mcpl').exists() + elif case == 'hdf5': + for batch in (3, 6): + _, _, tags = read_source(tmp_path / f'source.{batch}.h5') + assert tags is not None + elif case == 'not_separate': + _, _, tags = read_source(sp_path) + assert tags is not None + elif case == 'fixed_source': + assert FIXED_SOURCE_WARNING in out + elif case == 'random_ray': + assert RANDOM_RAY_WARNING in out + + # Eight fuel islands, x in s_x [a, a + 10] and y in s_y [a, a + 10] cm for a in # ISLAND_OFFSETS and s_x, s_y = +1 or -1, in the absorber of the two-slab # problem, which fills the rest of [-60, 60]^2 x [-5, 5] cm. Any two islands From 7c8d393feb32d8a85a5b96c6c7bfecd5703c1e8c Mon Sep 17 00:00:00 2001 From: Ondrej Chvala Date: Wed, 14 Oct 2026 08:40:39 -0400 Subject: [PATCH 6/7] Warn for MCPL source points written with overwrite_latest The warning that MCPL source point files carry no birth-mesh tags required the source to be written to separate files. With overwrite_latest, finalize_batch() writes the latest source at every batch even when the source is otherwise not written, so a run with write set to false and MCPL output dropped the tags without the warning. The warning now follows both writes: the source points of the given batches when the source is written to separate files, and the latest source with overwrite_latest. --- src/settings.cpp | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/settings.cpp b/src/settings.cpp index 9e296640dfc..590f3bfdaa5 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -974,9 +974,13 @@ void read_settings_xml(pugi::xml_node root) } // MCPL source point files have no place for the fission bank birth tags. - // Say so once, rather than at every source point written. - if (birth_mesh && source_write && source_separate && source_mcpl_write && - mpi::master) { + // Say so once, rather than at every source point written. finalize_batch() + // writes the source points of the given batches to separate files, and with + // overwrite_latest the latest source at every batch, even when the source is + // otherwise not written. + bool source_point_files = source_latest || (source_write && source_separate && + !sourcepoint_batch.empty()); + if (birth_mesh && source_mcpl_write && source_point_files && mpi::master) { warning("The tags are not written to MCPL source point " "files."); } From 71590c0a5933073adf8ea28ffb804800586e67ee Mon Sep 17 00:00:00 2001 From: Ondrej Chvala Date: Wed, 14 Oct 2026 08:49:26 -0400 Subject: [PATCH 7/7] Test the MCPL source point warning with overwrite_latest Extend test_birth_mesh_mcpl_warning. With overwrite_latest the latest source is written to an MCPL file at every batch, also when write is False, so that run gives the warning once. An empty list of source point batches and HDF5 output with overwrite_latest give none. Without the MCPL library, the overwrite_latest run stops where it writes the first MCPL file, after the warning; the test accepts that stop, and with the library it checks that source.mcpl is written. --- tests/regression_tests/birth_mesh/test.py | 26 +++++++++++++++++++---- 1 file changed, 22 insertions(+), 4 deletions(-) diff --git a/tests/regression_tests/birth_mesh/test.py b/tests/regression_tests/birth_mesh/test.py index d1c6b02f080..56b5a1cd3f0 100644 --- a/tests/regression_tests/birth_mesh/test.py +++ b/tests/regression_tests/birth_mesh/test.py @@ -460,11 +460,15 @@ def test_birth_mesh_random_ray(tmp_path, capfd): # Source point settings of test_birth_mesh_mcpl_warning. In the cases 'mcpl', # 'fixed_source' and 'random_ray', source points are requested in MCPL format # at a batch that the run does not reach, so no MCPL file is written and the -# MCPL library is not needed. +# MCPL library is not needed. With 'overwrite', the latest source is written +# at every batch, also when 'write' is False. MCPL_CASES = { 'mcpl': {'batches': [7], 'separate': True, 'mcpl': True}, 'mcpl_written': {'batches': [3, 6], 'separate': True, 'mcpl': True}, + 'overwrite_not_written': {'write': False, 'overwrite': True, 'mcpl': True}, 'hdf5': {'batches': [3, 6], 'separate': True}, + 'overwrite_hdf5': {'overwrite': True}, + 'empty_batches': {'batches': [], 'separate': True, 'mcpl': True}, 'not_separate': {'batches': [6], 'separate': False, 'mcpl': True}, 'not_written': {'write': False, 'separate': True, 'mcpl': True}, 'fixed_source': {'batches': [7], 'separate': True, 'mcpl': True}, @@ -485,17 +489,31 @@ def test_birth_mesh_mcpl_warning(tmp_path, capfd, case): model.settings.run_mode = 'fixed source' model.settings.batches = 2 model.settings.sourcepoint = MCPL_CASES[case] - sp_path = run_model(model, tmp_path) + try: + sp_path = run_model(model, tmp_path) + except RuntimeError as err: + # Without the MCPL library, the run of 'overwrite_not_written' stops + # where it writes the first MCPL file, after the settings are read + if (case != 'overwrite_not_written' or shutil.which('mcpl-config') + or 'the MCPL library is not available' not in str(err)): + raise + sp_path = None out = capfd.readouterr().out # The warning is given once, by the master process, when the settings are # read, and only for an eigenvalue run with the Monte Carlo solver that - # writes its source points to separate MCPL files - n_warnings = 1 if case in ('mcpl', 'mcpl_written') else 0 + # writes source points to MCPL files + mcpl_files = ('mcpl', 'mcpl_written', 'overwrite_not_written') + n_warnings = 1 if case in mcpl_files else 0 assert out.count(MCPL_WARNING) == n_warnings if case == 'mcpl_written': for batch in (3, 6): assert (tmp_path / f'source.{batch}.mcpl').exists() + elif case == 'overwrite_not_written' and sp_path is not None: + assert (tmp_path / 'source.mcpl').exists() + elif case == 'overwrite_hdf5': + _, _, tags = read_source(tmp_path / 'source.h5') + assert tags is not None elif case == 'hdf5': for batch in (3, 6): _, _, tags = read_source(tmp_path / f'source.{batch}.h5')