Add fission bank birth-mesh tagging

docs/source/io_formats/settings.rst

@@ -33,6 +33,43 @@ source neutrons.
 
   *Default*: None
 
+------------------------
+``<birth_mesh>`` Element
+------------------------
+
+The ``<birth_mesh>`` 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.<batch>.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.
+
 -----------------------------
 ``<collision_track>`` Element
 -----------------------------

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 <io_statepoint>`.
+
 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 ``<birth_mesh>`` 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 <io_statepoint>`. Only present in source
+             point files of an eigenvalue calculation for which the
+             ``<birth_mesh>`` 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``.

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 ``<birth_mesh>`` 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
+             ``<birth_mesh>`` settings element was specified for an eigenvalue
+             calculation.
 
 **/tallies/**
 

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.<batch>.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
 --------------------
 

include/openmc/constants.h

@@ -26,7 +26,7 @@ using double_4dvec = vector<vector<vector<vector<double>>>>;
 constexpr int HDF5_VERSION[] {3, 1};
 
 // Version numbers for binary files
-constexpr array<int, 2> VERSION_STATEPOINT {18, 2};
+constexpr array<int, 2> VERSION_STATEPOINT {18, 3};
 constexpr array<int, 2> VERSION_PARTICLE_RESTART {2, 1};
 constexpr array<int, 2> VERSION_TRACK {3, 1};
 constexpr array<int, 2> VERSION_SUMMARY {6, 1};

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 <birth_mesh> 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 {

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?

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 <birth_mesh> (see settings::birth_mesh).
 void write_h5_source_point(const char* filename, span<SourceSite> source_bank,
-  const vector<int64_t>& bank_index);
+  const vector<int64_t>& bank_index, bool birth_mesh_tags = false);
 
 void write_source_point(std::string, span<SourceSite> source_bank,
-  const vector<int64_t>& bank_index, bool use_mcpl);
+  const vector<int64_t>& 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<SourceSite> source_bank,
-  const vector<int64_t>& bank_index);
+  const vector<int64_t>& bank_index, bool birth_mesh_tags = false);
 
 void read_source_bank(
   hid_t group_id, vector<SourceSite>& sites, bool distribute);

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

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 <mesh> 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)

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)
 
 

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