Storing surface source points using a cell ID (#2888)

Co-authored-by: Paul Romano <[email protected]>

Author: JoffreyDorville

docs/source/io_formats/settings.rst

@@ -902,7 +902,12 @@ attributes/sub-elements:
 
 The ``<surf_source_write>`` element triggers OpenMC to bank particles crossing
 certain surfaces and write out the source bank in a separate file called
-``surface_source.h5``. This element has the following attributes/sub-elements:
+``surface_source.h5``. One or multiple surface IDs and one cell ID can be used
+to select the surfaces of interest. If no surface IDs are declared, every surface
+of the model is eligible to bank particles. In that case, a cell ID (using
+either the ``cell``, ``cellfrom`` or ``cellto`` attributes) can be used to select
+every surface of a specific cell. This element has the following
+attributes/sub-elements:
 
   :surface_ids:
     A list of integers separated by spaces indicating the unique IDs of surfaces
@@ -928,6 +933,34 @@ certain surfaces and write out the source bank in a separate file called
 
     .. _MCPL: https://mctools.github.io/mcpl/mcpl.pdf
 
+  :cell:
+    An integer representing the cell ID used to determine if particles crossing
+    identified surfaces are to be banked. Particles coming from or going to this
+    declared cell will be banked if they cross the identified surfaces.
+
+    *Default*: None
+
+  :cellfrom:
+    An integer representing the cell ID used to determine if particles crossing
+    identified surfaces are to be banked. Particles coming from this declared
+    cell will be banked if they cross the identified surfaces.
+
+    *Default*: None
+
+  :cellto:
+    An integer representing the cell ID used to determine if particles crossing
+    identified surfaces are to be banked. Particles going to this declared cell
+    will be banked if they cross the identified surfaces.
+
+    *Default*: None
+
+.. note:: The ``cell``, ``cellfrom`` and ``cellto`` attributes cannot be
+          used simultaneously.
+
+.. note:: Surfaces with boundary conditions that are not "transmission" or "vacuum"
+          are not eligible to store any particles when using ``cell``, ``cellfrom``
+          or ``cellto`` attributes. It is recommended to use surface IDs instead.
+
 ------------------------------
 ``<survival_biasing>`` Element
 ------------------------------

docs/source/usersguide/settings.rst

@@ -277,6 +277,9 @@ source file can be manually generated with the :func:`openmc.write_source_file`
 function. This is particularly useful for coupling OpenMC with another program
 that generates a source to be used in OpenMC.
 
+Surface Sources
++++++++++++++++
+
 A source file based on particles that cross one or more surfaces can be
 generated during a simulation using the :attr:`Settings.surf_source_write`
 attribute::
@@ -287,7 +290,51 @@ attribute::
   }
 
 In this example, at most 10,000 source particles are stored when particles cross
-surfaces with IDs of 1, 2, or 3.
+surfaces with IDs of 1, 2, or 3. If no surface IDs are declared, particles
+crossing any surface of the model will be banked::
+
+  settings.surf_source_write = {'max_particles': 10000}
+
+A cell ID can also be used to bank particles that are crossing any surface of
+a cell that particles are either coming from or going to::
+
+  settings.surf_source_write = {'cell': 1, 'max_particles': 10000}
+
+In this example, particles that are crossing any surface that bounds cell 1 will
+be banked excluding any surface that does not use a 'transmission' or 'vacuum'
+boundary condition.
+
+.. note:: Surfaces with boundary conditions that are not "transmission" or "vacuum"
+          are not eligible to store any particles when using ``cell``, ``cellfrom``
+          or ``cellto`` attributes. It is recommended to use surface IDs instead.
+
+Surface IDs can be used in combination with a cell ID::
+
+  settings.surf_source_write = {
+      'cell': 1,
+      'surfaces_ids': [1, 2, 3],
+      'max_particles': 10000
+  }
+
+In that case, only particles that are crossing the declared surfaces coming from
+cell 1 or going to cell 1 will be banked. To account specifically for particles
+leaving or entering a given cell, ``cellfrom`` and ``cellto`` are also available
+to respectively account for particles coming from a cell::
+
+  settings.surf_source_write = {
+      'cellfrom': 1,
+      'max_particles': 10000
+  }
+
+or particles going to a cell::
+
+  settings.surf_source_write = {
+      'cellto': 1,
+      'max_particles': 10000
+  }
+
+.. note:: The ``cell``, ``cellfrom`` and ``cellto`` attributes cannot be
+          used simultaneously.
 
 .. _compiled_source:
 

include/openmc/particle.h

@@ -74,7 +74,7 @@ class Particle : public ParticleData {
   void pht_secondary_particles();
 
   //! Cross a surface and handle boundary conditions
-  void cross_surface();
+  void cross_surface(const Surface& surf);
 
   //! Cross a vacuum boundary condition.
   //
@@ -127,6 +127,8 @@ std::string particle_type_to_str(ParticleType type);
 
 ParticleType str_to_particle_type(std::string str);
 
+void add_surf_source_to_bank(Particle& p, const Surface& surf);
+
 } // namespace openmc
 
 #endif // OPENMC_PARTICLE_H

include/openmc/particle_data.h

@@ -211,9 +211,15 @@ class GeometryState {
   // resets all coordinate levels for the particle
   void clear()
   {
-    for (auto& level : coord_)
+    for (auto& level : coord_) {
       level.reset();
+    }
     n_coord_ = 1;
+
+    for (auto& cell : cell_last_) {
+      cell = C_NONE;
+    }
+    n_coord_last_ = 1;
   }
 
   // Initialize all internal state from position and direction

include/openmc/settings.h

@@ -16,6 +16,14 @@
 
 namespace openmc {
 
+// Type of surface source write
+enum class SSWCellType {
+  None,
+  Both,
+  From,
+  To,
+};
+
 //==============================================================================
 // Global variable declarations
 //==============================================================================
@@ -125,6 +133,10 @@ extern int
   max_history_splits; //!< maximum number of particle splits for weight windows
 extern int64_t max_surface_particles; //!< maximum number of particles to be
                                       //!< banked on surfaces per process
+extern int64_t ssw_cell_id;           //!< Cell id for the surface source
+                                      //!< write setting
+extern SSWCellType ssw_cell_type;     //!< Type of option for the cell
+                                      //!< argument of surface source write
 extern TemperatureMethod
   temperature_method; //!< method for choosing temperatures
 extern double

openmc/settings.py

@@ -201,6 +201,15 @@ class Settings:
         :max_particles: Maximum number of particles to be banked on surfaces per
                    process (int)
         :mcpl: Output in the form of an MCPL-file (bool)
+        :cell: Cell ID used to determine if particles crossing identified
+               surfaces are to be banked. Particles coming from or going to this
+               declared cell will be banked (int)
+        :cellfrom: Cell ID used to determine if particles crossing identified
+                   surfaces are to be banked. Particles coming from this
+                   declared cell will be banked (int)
+        :cellto: Cell ID used to determine if particles crossing identified
+                 surfaces are to be banked. Particles going to this declared
+                 cell will be banked (int)
     survival_biasing : bool
         Indicate whether survival biasing is to be used
     tabular_legendre : dict
@@ -693,23 +702,30 @@ def surf_source_write(self) -> dict:
 
     @surf_source_write.setter
     def surf_source_write(self, surf_source_write: dict):
-        cv.check_type('surface source writing options', surf_source_write, Mapping)
+        cv.check_type("surface source writing options", surf_source_write, Mapping)
         for key, value in surf_source_write.items():
-            cv.check_value('surface source writing key', key,
-                           ('surface_ids', 'max_particles', 'mcpl'))
-            if key == 'surface_ids':
-                cv.check_type('surface ids for source banking', value,
-                              Iterable, Integral)
+            cv.check_value(
+                "surface source writing key",
+                key,
+                ("surface_ids", "max_particles", "mcpl", "cell", "cellfrom", "cellto"),
+            )
+            if key == "surface_ids":
+                cv.check_type(
+                    "surface ids for source banking", value, Iterable, Integral
+                )
                 for surf_id in value:
-                    cv.check_greater_than('surface id for source banking',
-                                          surf_id, 0)
-            elif key == 'max_particles':
-                cv.check_type('maximum particle banks on surfaces per process',
-                              value, Integral)
-                cv.check_greater_than('maximum particle banks on surfaces per process',
-                                      value, 0)
-            elif key == 'mcpl':
-                cv.check_type('write to an MCPL-format file', value, bool)
+                    cv.check_greater_than("surface id for source banking", surf_id, 0)
+            elif key == "mcpl":
+                cv.check_type("write to an MCPL-format file", value, bool)
+            elif key in ("max_particles", "cell", "cellfrom", "cellto"):
+                name = {
+                    "max_particles": "maximum particle banks on surfaces per process",
+                    "cell": "Cell ID for source banking (from or to)",
+                    "cellfrom": "Cell ID for source banking (from only)",
+                    "cellto": "Cell ID for source banking (to only)",
+                }[key]
+                cv.check_type(name, value, Integral)
+                cv.check_greater_than(name, value, 0)
 
         self._surf_source_write = surf_source_write
 
@@ -1207,16 +1223,18 @@ def _create_surf_source_read_subelement(self, root):
     def _create_surf_source_write_subelement(self, root):
         if self._surf_source_write:
             element = ET.SubElement(root, "surf_source_write")
-            if 'surface_ids' in self._surf_source_write:
+            if "surface_ids" in self._surf_source_write:
                 subelement = ET.SubElement(element, "surface_ids")
-                subelement.text = ' '.join(
-                    str(x) for x in self._surf_source_write['surface_ids'])
-            if 'max_particles' in self._surf_source_write:
-                subelement = ET.SubElement(element, "max_particles")
-                subelement.text = str(self._surf_source_write['max_particles'])
-            if 'mcpl' in self._surf_source_write:
+                subelement.text = " ".join(
+                    str(x) for x in self._surf_source_write["surface_ids"]
+                )
+            if "mcpl" in self._surf_source_write:
                 subelement = ET.SubElement(element, "mcpl")
-                subelement.text = str(self._surf_source_write['mcpl']).lower()
+                subelement.text = str(self._surf_source_write["mcpl"]).lower()
+            for key in ("max_particles", "cell", "cellfrom", "cellto"):
+                if key in self._surf_source_write:
+                    subelement = ET.SubElement(element, key)
+                    subelement.text = str(self._surf_source_write[key])
 
     def _create_confidence_intervals(self, root):
         if self._confidence_intervals is not None:
@@ -1381,9 +1399,9 @@ def _create_create_fission_neutrons_subelement(self, root):
             elem.text = str(self._create_fission_neutrons).lower()
 
     def _create_create_delayed_neutrons_subelement(self, root):
-       if self._create_delayed_neutrons is not None:
-           elem = ET.SubElement(root, "create_delayed_neutrons")
-           elem.text = str(self._create_delayed_neutrons).lower()
+        if self._create_delayed_neutrons is not None:
+            elem = ET.SubElement(root, "create_delayed_neutrons")
+            elem.text = str(self._create_delayed_neutrons).lower()
 
     def _create_delayed_photon_scaling_subelement(self, root):
         if self._delayed_photon_scaling is not None:
@@ -1610,17 +1628,18 @@ def _surf_source_read_from_xml_element(self, root):
 
     def _surf_source_write_from_xml_element(self, root):
         elem = root.find('surf_source_write')
-        if elem is not None:
-            for key in ('surface_ids', 'max_particles','mcpl'):
-                value = get_text(elem, key)
-                if value is not None:
-                    if key == 'surface_ids':
-                        value = [int(x) for x in value.split()]
-                    elif key in ('max_particles'):
-                        value = int(value)
-                    elif key == 'mcpl':
-                        value = value in ('true', '1')
-                    self.surf_source_write[key] = value
+        if elem is None:
+            return
+        for key in ('surface_ids', 'max_particles', 'mcpl', 'cell', 'cellto', 'cellfrom'):
+            value = get_text(elem, key)
+            if value is not None:
+                if key == 'surface_ids':
+                    value = [int(x) for x in value.split()]
+                elif key == 'mcpl':
+                    value = value in ('true', '1')
+                elif key in ('max_particles', 'cell', 'cellfrom', 'cellto'):
+                    value = int(value)
+                self.surf_source_write[key] = value
 
     def _confidence_intervals_from_xml_element(self, root):
         text = get_text(root, 'confidence_intervals')

src/dagmc.cpp

@@ -271,8 +271,10 @@ void DAGUniverse::init_geometry()
                                   : dagmc_instance_->id_by_index(2, i + 1);
 
     // set surface source attribute if needed
-    if (contains(settings::source_write_surf_id, s->id_))
+    if (contains(settings::source_write_surf_id, s->id_) ||
+        settings::source_write_surf_id.empty()) {
       s->surf_source_ = true;
+    }
 
     // set BCs
     std::string bc_value =