diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index dcb63d0fc10..18ddf79edcf 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -33,6 +33,43 @@ 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 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 +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..72ad199d836 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 and with the +random ray solver, 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..853aea41a73 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -46,8 +46,15 @@ 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 and with the +// random ray solver +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..07b53f37c48 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 or with the random ray solver. + + .. 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..590f3bfdaa5 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,36 @@ 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, + // 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(); + birth_mesh = m.get(); + } + } + // Check if the user has specified to write state points if (check_for_node(root, "state_point")) { @@ -942,6 +973,18 @@ 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. 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."); + } + // Check is the user specified to convert strength to statistical weight if (check_for_node(root, "uniform_source_sampling")) { uniform_source_sampling = @@ -1379,6 +1422,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/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/__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..56b5a1cd3f0 --- /dev/null +++ b/tests/regression_tests/birth_mesh/test.py @@ -0,0 +1,676 @@ +"""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. 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() +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 + +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 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']}) + 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) + + # 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): + # 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_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]} + 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, 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'])) + + # 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 + + +# 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. 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}, + '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] + 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 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') + 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 +# 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..be42f75905d --- /dev/null +++ b/tests/regression_tests/birth_mesh_ce/test.py @@ -0,0 +1,213 @@ +"""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. + +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 +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 + + +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) 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', '