Add a max_events setting. (#2945)

Co-authored-by: Paul Romano <paul.k.romano@gmail.com>

docs/source/io_formats/settings.rst

@@ -252,19 +252,27 @@ to false.
 
   *Default*: true
 
-----------------------------------------
+-------------------------------------
 ``<max_particles_in_flight>`` Element
-----------------------------------------
+-------------------------------------
 
 This element indicates the number of neutrons to run in flight concurrently
 when using event-based parallelism. A higher value uses more memory, but
 may be more efficient computationally.
 
   *Default*: 100000
 
----------------------------
+---------------------------------
+``<max_particle_events>`` Element
+---------------------------------
+
+This element indicates the maximum number of events a particle can undergo.
+
+  *Default*: 1000000
+
+-----------------------
 ``<max_order>`` Element
----------------------------
+-----------------------
 
 The ``<max_order>`` element allows the user to set a maximum scattering order
 to apply to every nuclide/material in the problem.  That is, if the data
@@ -276,11 +284,11 @@ then, OpenMC will only use up to the :math:`P_1` data.
   .. note:: This element is not used in the continuous-energy
     :ref:`energy_mode`.
 
----------------------------
+------------------------
 ``<max_splits>`` Element
----------------------------
+------------------------
 
-The ``<max_splits>`` element indicates the number of times a particle can split during a history. 
+The ``<max_splits>`` element indicates the number of times a particle can split during a history.
 
   *Default*: 1000
 

include/openmc/particle_data.h

@@ -27,9 +27,6 @@ constexpr int MAX_DELAYED_GROUPS {8};
 
 constexpr double CACHE_INVALID {-1.0};
 
-// Maximum number of collisions/crossings
-constexpr int MAX_EVENTS {1000000};
-
 //==========================================================================
 // Aliases and type definitions
 

include/openmc/settings.h

@@ -93,8 +93,8 @@ extern "C" int32_t gen_per_batch; //!< number of generations per batch
 extern "C" int64_t n_particles;   //!< number of particles per generation
 
 extern int64_t
-  max_particles_in_flight; //!< Max num. event-based particles in flight
-
+  max_particles_in_flight;      //!< Max num. event-based particles in flight
+extern int max_particle_events; //!< Maximum number of particle events
 extern ElectronTreatment
   electron_treatment; //!< how to treat secondary electrons
 extern array<double, 4>

openmc/settings.py

@@ -109,6 +109,10 @@ class Settings:
         parallelism.
 
         .. versionadded:: 0.12
+    max_particle_events : int
+        Maximum number of allowed particle events per source particle.
+
+        .. versionadded:: 0.14.1
     max_order : None or int
         Maximum scattering order to apply globally when in multi-group mode.
     max_splits : int
@@ -334,6 +338,7 @@ def __init__(self, **kwargs):
 
         self._event_based = None
         self._max_particles_in_flight = None
+        self._max_particle_events = None
         self._write_initial_source = None
         self._weight_windows = cv.CheckedList(WeightWindows, 'weight windows')
         self._weight_window_generators = cv.CheckedList(WeightWindowGenerator, 'weight window generators')
@@ -938,6 +943,16 @@ def max_particles_in_flight(self, value: int):
         cv.check_greater_than('max particles in flight', value, 0)
         self._max_particles_in_flight = value
 
+    @property
+    def max_particle_events(self) -> int:
+        return self._max_particle_events
+
+    @max_particle_events.setter
+    def max_particle_events(self, value: int):
+        cv.check_type('max particle events', value, Integral)
+        cv.check_greater_than('max particle events', value, 0)
+        self._max_particle_events = value
+
     @property
     def write_initial_source(self) -> bool:
         return self._write_initial_source
@@ -1341,6 +1356,11 @@ def _create_max_particles_in_flight_subelement(self, root):
             elem = ET.SubElement(root, "max_particles_in_flight")
             elem.text = str(self._max_particles_in_flight).lower()
 
+    def _create_max_events_subelement(self, root):
+        if self._max_particle_events is not None:
+            elem = ET.SubElement(root, "max_particle_events")
+            elem.text = str(self._max_particle_events).lower()
+
     def _create_material_cell_offsets_subelement(self, root):
         if self._material_cell_offsets is not None:
             elem = ET.SubElement(root, "material_cell_offsets")
@@ -1719,6 +1739,11 @@ def _max_particles_in_flight_from_xml_element(self, root):
         if text is not None:
             self.max_particles_in_flight = int(text)
 
+    def _max_particle_events_from_xml_element(self, root):
+        text = get_text(root, 'max_particle_events')
+        if text is not None:
+            self.max_particle_events = int(text)
+
     def _material_cell_offsets_from_xml_element(self, root):
         text = get_text(root, 'material_cell_offsets')
         if text is not None:
@@ -1820,6 +1845,7 @@ def to_xml_element(self, mesh_memo=None):
         self._create_delayed_photon_scaling_subelement(element)
         self._create_event_based_subelement(element)
         self._create_max_particles_in_flight_subelement(element)
+        self._create_max_events_subelement(element)
         self._create_material_cell_offsets_subelement(element)
         self._create_log_grid_bins_subelement(element)
         self._create_write_initial_source_subelement(element)
@@ -1923,6 +1949,7 @@ def from_xml_element(cls, elem, meshes=None):
         settings._delayed_photon_scaling_from_xml_element(elem)
         settings._event_based_from_xml_element(elem)
         settings._max_particles_in_flight_from_xml_element(elem)
+        settings._max_particle_events_from_xml_element(elem)
         settings._material_cell_offsets_from_xml_element(elem)
         settings._log_grid_bins_from_xml_element(elem)
         settings._write_initial_source_from_xml_element(elem)

src/finalize.cpp

@@ -91,6 +91,7 @@ int openmc_finalize()
   settings::max_lost_particles = 10;
   settings::max_order = 0;
   settings::max_particles_in_flight = 100000;
+  settings::max_particle_events = 1000000;
   settings::max_splits = 1000;
   settings::max_tracks = 1000;
   settings::max_write_lost_particles = -1;

src/particle.cpp

@@ -381,7 +381,7 @@ void Particle::event_revive_from_secondary()
 {
   // If particle has too many events, display warning and kill it
   ++n_event();
-  if (n_event() == MAX_EVENTS) {
+  if (n_event() == settings::max_particle_events) {
     warning("Particle " + std::to_string(id()) +
             " underwent maximum number of events.");
     wgt() = 0.0;

src/settings.cpp

@@ -96,6 +96,7 @@ int32_t gen_per_batch {1};
 int64_t n_particles {-1};
 
 int64_t max_particles_in_flight {100000};
+int max_particle_events {1000000};
 
 ElectronTreatment electron_treatment {ElectronTreatment::TTB};
 array<double, 4> energy_cutoff {0.0, 1000.0, 0.0, 0.0};
@@ -156,6 +157,12 @@ void get_run_parameters(pugi::xml_node node_base)
       std::stoll(get_node_value(node_base, "max_particles_in_flight"));
   }
 
+  // Get maximum number of events allowed per particle
+  if (check_for_node(node_base, "max_particle_events")) {
+    max_particle_events =
+      std::stoll(get_node_value(node_base, "max_particle_events"));
+  }
+
   // Get number of basic batches
   if (check_for_node(node_base, "batches")) {
     n_batches = std::stoi(get_node_value(node_base, "batches"));

tests/unit_tests/test_settings.py

@@ -27,8 +27,8 @@ def test_export_to_xml(run_in_tmpdir):
     s.survival_biasing = True
     s.cutoff = {'weight': 0.25, 'weight_avg': 0.5, 'energy_neutron': 1.0e-5,
                 'energy_photon': 1000.0, 'energy_electron': 1.0e-5,
-                'energy_positron': 1.0e-5, 'time_neutron': 1.0e-5, 
-                'time_photon': 1.0e-5, 'time_electron': 1.0e-5, 
+                'energy_positron': 1.0e-5, 'time_neutron': 1.0e-5,
+                'time_photon': 1.0e-5, 'time_electron': 1.0e-5,
                 'time_positron': 1.0e-5}
     mesh = openmc.RegularMesh()
     mesh.lower_left = (-10., -10., -10.)
@@ -59,6 +59,8 @@ def test_export_to_xml(run_in_tmpdir):
     s.write_initial_source = True
     s.weight_window_checkpoints = {'surface': True, 'collision': False}
 
+    s.max_particle_events = 100
+
     # Make sure exporting XML works
     s.export_to_xml()
 
@@ -92,7 +94,7 @@ def test_export_to_xml(run_in_tmpdir):
     assert s.cutoff == {'weight': 0.25, 'weight_avg': 0.5,
                         'energy_neutron': 1.0e-5, 'energy_photon': 1000.0,
                         'energy_electron': 1.0e-5, 'energy_positron': 1.0e-5,
-                        'time_neutron': 1.0e-5, 'time_photon': 1.0e-5, 
+                        'time_neutron': 1.0e-5, 'time_photon': 1.0e-5,
                         'time_electron': 1.0e-5, 'time_positron': 1.0e-5}
     assert isinstance(s.entropy_mesh, openmc.RegularMesh)
     assert s.entropy_mesh.lower_left == [-10., -10., -10.]
@@ -128,3 +130,4 @@ def test_export_to_xml(run_in_tmpdir):
     assert vol.lower_left == (-10., -10., -10.)
     assert vol.upper_right == (10., 10., 10.)
     assert s.weight_window_checkpoints == {'surface': True, 'collision': False}
+    assert s.max_particle_events == 100