Add user setting for source rejection fraction (#3433)

docs/source/io_formats/settings.rst

@@ -487,7 +487,7 @@ found in the :ref:`random ray user guide <random_ray>`.
 
         :type:
           The type of the domain. Can be ``material``, ``cell``, or ``universe``.
-  
+
   :diagonal_stabilization_rho:
     The rho factor for use with diagonal stabilization. This technique is
     applied when negative diagonal (in-group) elements are detected in
@@ -917,6 +917,16 @@ variable and whose sub-elements/attributes are as follows:
   :dist:
     This sub-element of a ``pair`` element provides information on the corresponding univariate distribution.
 
+---------------------------------------
+``<source_rejection_fraction>`` Element
+---------------------------------------
+
+The ``<source_rejection_fraction>`` element specifies the minimum fraction of
+external source sites that must be accepted when applying rejection sampling
+based on constraints.
+
+   *Default*: 0.05
+
 -------------------------
 ``<state_point>`` Element
 -------------------------

include/openmc/settings.h

@@ -145,6 +145,8 @@ extern std::unordered_set<int>
   statepoint_batch; //!< Batches when state should be written
 extern std::unordered_set<int>
   source_write_surf_id; //!< Surface ids where sources will be written
+extern double source_rejection_fraction; //!< Minimum fraction of source sites
+                                         //!< that must be accepted
 
 extern int
   max_history_splits; //!< maximum number of particle splits for weight windows

include/openmc/source.h

@@ -21,10 +21,9 @@ namespace openmc {
 // Constants
 //==============================================================================
 
-// Maximum number of external source spatial resamples to encounter before an
-// error is thrown.
+// Minimum number of external source sites rejected before checking againts the
+// source_rejection_fraction
 constexpr int EXTSRC_REJECT_THRESHOLD {10000};
-constexpr double EXTSRC_REJECT_FRACTION {0.05};
 
 //==============================================================================
 // Global variables

openmc/settings.py

@@ -222,6 +222,10 @@ class Settings:
         Number of random numbers allocated for each source particle history
     source : Iterable of openmc.SourceBase
         Distribution of source sites in space, angle, and energy
+    source_rejection_fraction : float
+        Minimum fraction of source sites that must be accepted when applying
+        rejection sampling based on constraints. If not specified, the default
+        value is 0.05.
     sourcepoint : dict
         Options for writing source points. Acceptable keys are:
 
@@ -357,6 +361,7 @@ class Settings:
 
         # Source subelement
         self._source = cv.CheckedList(SourceBase, 'source distributions')
+        self._source_rejection_fraction = None
 
         self._confidence_intervals = None
         self._electron_treatment = None
@@ -1224,6 +1229,17 @@ class Settings:
         cv.check_type('use decay photons', value, bool)
         self._use_decay_photons = value
 
+    @property
+    def source_rejection_fraction(self) -> float:
+        return self._source_rejection_fraction
+
+    @source_rejection_fraction.setter
+    def source_rejection_fraction(self, source_rejection_fraction: float):
+        cv.check_type('source_rejection_fraction', source_rejection_fraction, Real)
+        cv.check_greater_than('source_rejection_fraction', source_rejection_fraction, 0)
+        cv.check_less_than('source_rejection_fraction', source_rejection_fraction, 1)
+        self._source_rejection_fraction = source_rejection_fraction
+
     def _create_run_mode_subelement(self, root):
         elem = ET.SubElement(root, "run_mode")
         elem.text = self._run_mode.value
@@ -1685,6 +1701,11 @@ class Settings:
                     subelement = ET.SubElement(element, key)
                     subelement.text = str(value)
 
+    def _create_source_rejection_fraction_subelement(self, root):
+        if self._source_rejection_fraction is not None:
+            element = ET.SubElement(root, "source_rejection_fraction")
+            element.text = str(self._source_rejection_fraction)
+
     def _eigenvalue_from_xml_element(self, root):
         elem = root.find('eigenvalue')
         if elem is not None:
@@ -2104,6 +2125,11 @@ class Settings:
         if text is not None:
             self.use_decay_photons = text in ('true', '1')
 
+    def _source_rejection_fraction_from_xml_element(self, root):
+        text = get_text(root, 'source_rejection_fraction')
+        if text is not None:
+            self.source_rejection_fraction = float(text)
+
     def to_xml_element(self, mesh_memo=None):
         """Create a 'settings' element to be written to an XML file.
 
@@ -2171,6 +2197,7 @@ class Settings:
         self._create_max_tracks_subelement(element)
         self._create_random_ray_subelement(element, mesh_memo)
         self._create_use_decay_photons_subelement(element)
+        self._create_source_rejection_fraction_subelement(element)
 
         # Clean the indentation in the file to be user-readable
         clean_indentation(element)
@@ -2279,6 +2306,7 @@ class Settings:
         settings._max_tracks_from_xml_element(elem)
         settings._random_ray_from_xml_element(elem)
         settings._use_decay_photons_from_xml_element(elem)
+        settings._source_rejection_fraction_from_xml_element(elem)
 
         return settings
 

src/error.cpp

@@ -110,23 +110,26 @@ void write_message(const std::string& message, int level)
 
 void fatal_error(const std::string& message, int err)
 {
+#pragma omp critical(FatalError)
+  {
 #ifdef _POSIX_VERSION
-  // Make output red if user is in a terminal
-  if (isatty(STDERR_FILENO)) {
-    std::cerr << "\033[0;31m";
-  }
+    // Make output red if user is in a terminal
+    if (isatty(STDERR_FILENO)) {
+      std::cerr << "\033[0;31m";
+    }
 #endif
 
-  // Write error message
-  std::cerr << " ERROR: ";
-  output(message, std::cerr, 8);
+    // Write error message
+    std::cerr << " ERROR: ";
+    output(message, std::cerr, 8);
 
 #ifdef _POSIX_VERSION
-  // Reset color for terminal
-  if (isatty(STDERR_FILENO)) {
-    std::cerr << "\033[0m";
-  }
+    // Reset color for terminal
+    if (isatty(STDERR_FILENO)) {
+      std::cerr << "\033[0m";
+    }
 #endif
+  }
 
 #ifdef OPENMC_MPI
   MPI_Abort(mpi::intracomm, err);

src/finalize.cpp

@@ -119,6 +119,7 @@ int openmc_finalize()
   settings::run_CE = true;
   settings::run_mode = RunMode::UNSET;
   settings::source_latest = false;
+  settings::source_rejection_fraction = 0.05;
   settings::source_separate = false;
   settings::source_write = true;
   settings::ssw_cell_id = C_NONE;

src/settings.cpp

@@ -124,6 +124,7 @@ RunMode run_mode {RunMode::UNSET};
 SolverType solver_type {SolverType::MONTE_CARLO};
 std::unordered_set<int> sourcepoint_batch;
 std::unordered_set<int> statepoint_batch;
+double source_rejection_fraction {0.05};
 std::unordered_set<int> source_write_surf_id;
 int64_t ssw_max_particles;
 int64_t ssw_max_files;
@@ -644,6 +645,12 @@ void read_settings_xml(pugi::xml_node root)
     write_initial_source = get_node_value_bool(root, "write_initial_source");
   }
 
+  // Get relative number of lost particles
+  if (check_for_node(root, "source_rejection_fraction")) {
+    source_rejection_fraction =
+      std::stod(get_node_value(root, "source_rejection_fraction"));
+  }
+
   // Survival biasing
   if (check_for_node(root, "survival_biasing")) {
     survival_biasing = get_node_value_bool(root, "survival_biasing");

src/source.cpp

@@ -157,11 +157,28 @@ void Source::read_constraints(pugi::xml_node node)
   }
 }
 
+void check_rejection_fraction(int64_t n_reject, int64_t n_accept)
+{
+  // Don't check unless we've hit a minimum number of total sites rejected
+  if (n_reject < EXTSRC_REJECT_THRESHOLD)
+    return;
+
+  // Compute fraction of accepted sites and compare against minimum
+  double fraction = static_cast<double>(n_accept) / n_reject;
+  if (fraction <= settings::source_rejection_fraction) {
+    fatal_error(fmt::format(
+      "Too few source sites satisfied the constraints (minimum source "
+      "rejection fraction = {}). Please check your source definition or "
+      "set a lower value of Settings.source_rejection_fraction.",
+      settings::source_rejection_fraction));
+  }
+}
+
 SourceSite Source::sample_with_constraints(uint64_t* seed) const
 {
   bool accepted = false;
-  static int n_reject = 0;
-  static int n_accept = 0;
+  static int64_t n_reject = 0;
+  static int64_t n_accept = 0;
   SourceSite site;
 
   while (!accepted) {
@@ -176,13 +193,9 @@ SourceSite Source::sample_with_constraints(uint64_t* seed) const
                  satisfies_energy_constraints(site.E) &&
                  satisfies_time_constraints(site.time);
       if (!accepted) {
+        // Increment number of rejections and check against minimum fraction
         ++n_reject;
-        if (n_reject >= EXTSRC_REJECT_THRESHOLD &&
-            static_cast<double>(n_accept) / n_reject <=
-              EXTSRC_REJECT_FRACTION) {
-          fatal_error("More than 95% of external source sites sampled were "
-                      "rejected. Please check your source definition.");
-        }
+        check_rejection_fraction(n_reject, n_accept);
 
         // For the "kill" strategy, accept particle but set weight to 0 so that
         // it is terminated immediately
@@ -337,8 +350,8 @@ SourceSite IndependentSource::sample(uint64_t* seed) const
 
   // Repeat sampling source location until a good site has been accepted
   bool accepted = false;
-  static int n_reject = 0;
-  static int n_accept = 0;
+  static int64_t n_reject = 0;
+  static int64_t n_accept = 0;
 
   while (!accepted) {
 
@@ -351,12 +364,7 @@ SourceSite IndependentSource::sample(uint64_t* seed) const
     // Check for rejection
     if (!accepted) {
       ++n_reject;
-      if (n_reject >= EXTSRC_REJECT_THRESHOLD &&
-          static_cast<double>(n_accept) / n_reject <= EXTSRC_REJECT_FRACTION) {
-        fatal_error("More than 95% of external source sites sampled were "
-                    "rejected. Please check your external source's spatial "
-                    "definition.");
-      }
+      check_rejection_fraction(n_reject, n_accept);
     }
   }
 
@@ -381,18 +389,12 @@ SourceSite IndependentSource::sample(uint64_t* seed) const
       site.E = energy_->sample(seed);
 
       // Resample if energy falls above maximum particle energy
-      if (site.E < data::energy_max[p] and
+      if (site.E < data::energy_max[p] &&
           (satisfies_energy_constraints(site.E)))
         break;
 
       n_reject++;
-      if (n_reject >= EXTSRC_REJECT_THRESHOLD &&
-          static_cast<double>(n_accept) / n_reject <= EXTSRC_REJECT_FRACTION) {
-        fatal_error(
-          "More than 95% of external source sites sampled were "
-          "rejected. Please check your external source energy spectrum "
-          "definition.");
-      }
+      check_rejection_fraction(n_reject, n_accept);
     }
 
     // Sample particle creation time

tests/unit_tests/test_settings.py

@@ -66,8 +66,8 @@ def test_export_to_xml(run_in_tmpdir):
             space=openmc.stats.Box((-1., -1., -1.), (1., 1., 1.))
         )
     }
-
     s.max_particle_events = 100
+    s.source_rejection_fraction = 0.01
 
     # Make sure exporting XML works
     s.export_to_xml()
@@ -130,7 +130,7 @@ def test_export_to_xml(run_in_tmpdir):
     assert s.log_grid_bins == 2000
     assert not s.photon_transport
     assert s.electron_treatment == 'led'
-    assert s.write_initial_source == True
+    assert s.write_initial_source
     assert len(s.volume_calculations) == 1
     vol = s.volume_calculations[0]
     assert vol.domain_type == 'cell'
@@ -144,3 +144,4 @@ def test_export_to_xml(run_in_tmpdir):
     assert s.random_ray['distance_active'] == 100.0
     assert s.random_ray['ray_source'].space.lower_left == [-1., -1., -1.]
     assert s.random_ray['ray_source'].space.upper_right == [1., 1., 1.]
+    assert s.source_rejection_fraction == 0.01

tests/unit_tests/test_source.py

@@ -264,21 +264,36 @@ def test_constraints_file(sphere_box_model, run_in_tmpdir):
 
 
 @pytest.mark.skipif(config['mpi'], reason='Not compatible with MPI')
-def test_rejection_limit(sphere_box_model, run_in_tmpdir):
-    model, cell1 = sphere_box_model[:2]
+def test_rejection_fraction(run_in_tmpdir):
+    mat = openmc.Material()
+    mat.add_nuclide('H1', 1.0)
+    w = 0.25
+    rpp1 = openmc.model.RectangularParallelepiped(
+        -w/2, w/2, -w/2, w/2, -w/2, w/2)
+    rpp2 = openmc.model.RectangularParallelepiped(
+        -0.5, 0.5, -0.5, 0.5, -0.5, 0.5, boundary_type='vacuum')
+    cell1 = openmc.Cell(fill=mat, region=-rpp1)
+    cell2 = openmc.Cell(region=+rpp1 & -rpp2)
+    model = openmc.Model()
+    model.geometry = openmc.Geometry([cell1, cell2])
 
-    # Define a point source that will get rejected 100% of the time
+    # Create a box source over a 1 cm³ volume that is constrained to the source
+    # cell of volume (0.25 cm)³ = 0.0125 cm³, which means the default rejection
+    # fraction of 0.05 won't work
+    model.settings.particles = 1000
+    model.settings.batches = 1
+    model.settings.run_mode = 'fixed source'
     model.settings.source = openmc.IndependentSource(
-        space=openmc.stats.Point((-3., 0., 0.)),
+        space=openmc.stats.Box(*(-rpp2).bounding_box),
         constraints={'domains': [cell1]}
     )
-
-    # Confirm that OpenMC doesn't run in an infinite loop. Note that this may
-    # work when running with MPI since it won't necessarily capture the error
-    # message correctly
-    with pytest.raises(RuntimeError, match="rejected"):
+    with pytest.raises(RuntimeError, match='Too few source sites'):
         model.run(openmc_exec=config['exe'])
 
+    # With a source rejection fraction below 0.0125, the simulation should run
+    model.settings.source_rejection_fraction = 0.005
+    model.run(openmc_exec=config['exe'])
+
 
 def test_exceptions():