From f571be87c5e46587ce898dd5fa944512283b6ea3 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 10 Jun 2025 22:51:51 -0500 Subject: [PATCH] Add user setting for source rejection fraction (#3433) --- docs/source/io_formats/settings.rst | 12 ++++++- include/openmc/settings.h | 2 ++ include/openmc/source.h | 5 ++- openmc/settings.py | 28 ++++++++++++++++ src/error.cpp | 25 ++++++++------- src/finalize.cpp | 1 + src/settings.cpp | 7 ++++ src/source.cpp | 50 +++++++++++++++-------------- tests/unit_tests/test_settings.py | 5 +-- tests/unit_tests/test_source.py | 33 +++++++++++++------ 10 files changed, 118 insertions(+), 50 deletions(-) diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index cef9bf911f..26673faac2 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -487,7 +487,7 @@ found in the :ref:`random ray user guide `. :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. +--------------------------------------- +```` Element +--------------------------------------- + +The ```` element specifies the minimum fraction of +external source sites that must be accepted when applying rejection sampling +based on constraints. + + *Default*: 0.05 + ------------------------- ```` Element ------------------------- diff --git a/include/openmc/settings.h b/include/openmc/settings.h index 9017b2d080..069d94c3b2 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -145,6 +145,8 @@ extern std::unordered_set statepoint_batch; //!< Batches when state should be written extern std::unordered_set 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 diff --git a/include/openmc/source.h b/include/openmc/source.h index 6733eaeffc..e82c42269f 100644 --- a/include/openmc/source.h +++ b/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 diff --git a/openmc/settings.py b/openmc/settings.py index 4cacefbedb..bd3a89e11c 100644 --- a/openmc/settings.py +++ b/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 diff --git a/src/error.cpp b/src/error.cpp index 566950a973..f99f5935f0 100644 --- a/src/error.cpp +++ b/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); diff --git a/src/finalize.cpp b/src/finalize.cpp index 54aa1d1d9e..9cbebc878f 100644 --- a/src/finalize.cpp +++ b/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; diff --git a/src/settings.cpp b/src/settings.cpp index c8230a10f1..9afe743614 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -124,6 +124,7 @@ RunMode run_mode {RunMode::UNSET}; SolverType solver_type {SolverType::MONTE_CARLO}; std::unordered_set sourcepoint_batch; std::unordered_set statepoint_batch; +double source_rejection_fraction {0.05}; std::unordered_set 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"); diff --git a/src/source.cpp b/src/source.cpp index 8e6eb4f11f..6ee0b56035 100644 --- a/src/source.cpp +++ b/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(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(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(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(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 diff --git a/tests/unit_tests/test_settings.py b/tests/unit_tests/test_settings.py index 397d704342..7f202bcdc7 100644 --- a/tests/unit_tests/test_settings.py +++ b/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 diff --git a/tests/unit_tests/test_source.py b/tests/unit_tests/test_source.py index c88fbcbe62..bb8a1b7852 100644 --- a/tests/unit_tests/test_source.py +++ b/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():