From 239f7fed5e39eadfdad4f47c08f3f7b633862ab9 Mon Sep 17 00:00:00 2001 From: ahman24 <79746189+ahman24@users.noreply.github.com> Date: Wed, 5 Mar 2025 08:26:38 +0900 Subject: [PATCH] Implement user-configurable random number stride (#3067) Co-authored-by: Paul Romano --- docs/source/io_formats/settings.rst | 9 ++++++ docs/source/io_formats/statepoint.rst | 1 + include/openmc/capi.h | 2 ++ include/openmc/random_lcg.h | 14 +++++++++ openmc/lib/settings.py | 10 +++++++ openmc/settings.py | 25 ++++++++++++++++ openmc/statepoint.py | 6 ++++ src/finalize.cpp | 2 ++ src/initialize.cpp | 5 ++-- src/random_lcg.cpp | 12 +++++++- src/settings.cpp | 6 ++++ src/state_point.cpp | 8 +++++ tests/regression_tests/stride/__init__.py | 0 tests/regression_tests/stride/inputs_true.dat | 25 ++++++++++++++++ .../regression_tests/stride/results_true.dat | 2 ++ tests/regression_tests/stride/test.py | 29 +++++++++++++++++++ 16 files changed, 153 insertions(+), 3 deletions(-) create mode 100644 tests/regression_tests/stride/__init__.py create mode 100644 tests/regression_tests/stride/inputs_true.dat create mode 100644 tests/regression_tests/stride/results_true.dat create mode 100644 tests/regression_tests/stride/test.py diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 0c225b1ed..5174cfd8b 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -538,6 +538,15 @@ pseudo-random number generator. *Default*: 1 +-------------------- +```` Element +-------------------- + +The ``stride`` element is used to specify how many random numbers are allocated +for each source particle history. + + *Default*: 152,917 + .. _source_element: -------------------- diff --git a/docs/source/io_formats/statepoint.rst b/docs/source/io_formats/statepoint.rst index fd00a3e77..3b1031769 100644 --- a/docs/source/io_formats/statepoint.rst +++ b/docs/source/io_formats/statepoint.rst @@ -23,6 +23,7 @@ The current version of the statepoint file format is 18.1. bank is present (1) or not (0). :Datasets: - **seed** (*int8_t*) -- Pseudo-random number generator seed. + - **stride** (*uint64_t*) -- Pseudo-random number generator stride. - **energy_mode** (*char[]*) -- Energy mode of the run, either 'continuous-energy' or 'multi-group'. - **run_mode** (*char[]*) -- Run mode used, either 'eigenvalue' or diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 3802692c3..54257d093 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -71,6 +71,7 @@ int openmc_get_nuclide_index(const char name[], int* index); int openmc_add_unstructured_mesh( const char filename[], const char library[], int* id); int64_t openmc_get_seed(); +uint64_t openmc_get_stride(); int openmc_get_tally_index(int32_t id, int32_t* index); void openmc_get_tally_next_id(int32_t* id); int openmc_global_tallies(double** ptr); @@ -137,6 +138,7 @@ int openmc_reset_timers(); int openmc_run(); int openmc_sample_external_source(size_t n, uint64_t* seed, void* sites); void openmc_set_seed(int64_t new_seed); +void openmc_set_stride(uint64_t new_stride); int openmc_set_n_batches( int32_t n_batches, bool set_max_batches, bool add_statepoint_batch); int openmc_simulation_finalize(); diff --git a/include/openmc/random_lcg.h b/include/openmc/random_lcg.h index 4157b7cfe..5aecafed3 100644 --- a/include/openmc/random_lcg.h +++ b/include/openmc/random_lcg.h @@ -15,6 +15,7 @@ constexpr int STREAM_SOURCE {1}; constexpr int STREAM_URR_PTABLE {2}; constexpr int STREAM_VOLUME {3}; constexpr int64_t DEFAULT_SEED {1}; +constexpr uint64_t DEFAULT_STRIDE {152917ULL}; //============================================================================== //! Generate a pseudo-random number using a linear congruential generator. @@ -98,5 +99,18 @@ extern "C" int64_t openmc_get_seed(); extern "C" void openmc_set_seed(int64_t new_seed); +//============================================================================== +//! Get OpenMC's stride. +//============================================================================== + +extern "C" uint64_t openmc_get_stride(); + +//============================================================================== +//! Set OpenMC's stride. +//! @param new_stride Stride. +//============================================================================== + +extern "C" void openmc_set_stride(uint64_t new_stride); + } // namespace openmc #endif // OPENMC_RANDOM_LCG_H diff --git a/openmc/lib/settings.py b/openmc/lib/settings.py index 062670ef8..4fba8d48b 100644 --- a/openmc/lib/settings.py +++ b/openmc/lib/settings.py @@ -12,6 +12,8 @@ _RUN_MODES = {1: 'fixed source', _dll.openmc_set_seed.argtypes = [c_int64] _dll.openmc_get_seed.restype = c_int64 +_dll.openmc_set_stride.argtypes = [c_int64] +_dll.openmc_get_stride.restype = c_int64 _dll.openmc_get_n_batches.argtypes = [POINTER(c_int), c_bool] _dll.openmc_get_n_batches.restype = c_int _dll.openmc_get_n_batches.errcheck = _error_handler @@ -68,6 +70,14 @@ class _Settings: def seed(self, seed): _dll.openmc_set_seed(seed) + @property + def stride(self): + return _dll.openmc_get_stride() + + @stride.setter + def stride(self, stride): + _dll.openmc_set_stride(stride) + def set_batches(self, n_batches, set_max_batches=True, add_sp_batch=True): """Set number of batches or maximum number of batches diff --git a/openmc/settings.py b/openmc/settings.py index 0e2a18399..882a17b68 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -193,6 +193,8 @@ class Settings: The type of calculation to perform (default is 'eigenvalue') seed : int Seed for the linear congruential pseudorandom number generator + stride : int + Number of random numbers allocated for each source particle history source : Iterable of openmc.SourceBase Distribution of source sites in space, angle, and energy sourcepoint : dict @@ -338,6 +340,7 @@ class Settings: self._ptables = None self._uniform_source_sampling = None self._seed = None + self._stride = None self._survival_biasing = None # Shannon entropy mesh @@ -614,6 +617,16 @@ class Settings: cv.check_greater_than('random number generator seed', seed, 0) self._seed = seed + @property + def stride(self) -> int: + return self._stride + + @stride.setter + def stride(self, stride: int): + cv.check_type('random number generator stride', stride, Integral) + cv.check_greater_than('random number generator stride', stride, 0) + self._stride = stride + @property def survival_biasing(self) -> bool: return self._survival_biasing @@ -1336,6 +1349,11 @@ class Settings: element = ET.SubElement(root, "seed") element.text = str(self._seed) + def _create_stride_subelement(self, root): + if self._stride is not None: + element = ET.SubElement(root, "stride") + element.text = str(self._stride) + def _create_survival_biasing_subelement(self, root): if self._survival_biasing is not None: element = ET.SubElement(root, "survival_biasing") @@ -1763,6 +1781,11 @@ class Settings: if text is not None: self.seed = int(text) + def _stride_from_xml_element(self, root): + text = get_text(root, 'stride') + if text is not None: + self.stride = int(text) + def _survival_biasing_from_xml_element(self, root): text = get_text(root, 'survival_biasing') if text is not None: @@ -2014,6 +2037,7 @@ class Settings: self._create_plot_seed_subelement(element) self._create_ptables_subelement(element) self._create_seed_subelement(element) + self._create_stride_subelement(element) self._create_survival_biasing_subelement(element) self._create_cutoff_subelement(element) self._create_entropy_mesh_subelement(element, mesh_memo) @@ -2122,6 +2146,7 @@ class Settings: settings._plot_seed_from_xml_element(elem) settings._ptables_from_xml_element(elem) settings._seed_from_xml_element(elem) + settings._stride_from_xml_element(elem) settings._survival_biasing_from_xml_element(elem) settings._cutoff_from_xml_element(elem) settings._entropy_mesh_from_xml_element(elem, meshes) diff --git a/openmc/statepoint.py b/openmc/statepoint.py index f2fd48066..715becf48 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -99,6 +99,8 @@ class StatePoint: and whose values are time values in seconds. seed : int Pseudorandom number generator seed + stride : int + Number of random numbers allocated for each particle history source : numpy.ndarray of compound datatype Array of source sites. The compound datatype has fields 'r', 'u', 'E', 'wgt', 'delayed_group', 'surf_id', and 'particle', corresponding to @@ -356,6 +358,10 @@ class StatePoint: def seed(self): return self._f['seed'][()] + @property + def stride(self): + return self._f['stride'][()] + @property def source(self): return self._f['source_bank'][()] if self.source_present else None diff --git a/src/finalize.cpp b/src/finalize.cpp index 981ec5cba..ad6f0cf62 100644 --- a/src/finalize.cpp +++ b/src/finalize.cpp @@ -159,6 +159,7 @@ int openmc_finalize() model::root_universe = -1; model::plotter_seed = 1; openmc::openmc_set_seed(DEFAULT_SEED); + openmc::openmc_set_stride(DEFAULT_STRIDE); // Deallocate arrays free_memory(); @@ -221,5 +222,6 @@ int openmc_hard_reset() // Reset the random number generator state openmc::openmc_set_seed(DEFAULT_SEED); + openmc::openmc_set_stride(DEFAULT_STRIDE); return 0; } diff --git a/src/initialize.cpp b/src/initialize.cpp index 5f717b137..4b821bee1 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -99,9 +99,10 @@ int openmc_init(int argc, char* argv[], const void* intracomm) } #endif - // Initialize random number generator -- if the user specifies a seed, it - // will be re-initialized later + // Initialize random number generator -- if the user specifies a seed and/or + // stride, it will be re-initialized later openmc::openmc_set_seed(DEFAULT_SEED); + openmc::openmc_set_stride(DEFAULT_STRIDE); // Copy previous locale and set locale to C. This is a workaround for an issue // whereby when openmc_init is called from the plotter, the Qt application diff --git a/src/random_lcg.cpp b/src/random_lcg.cpp index ca4467719..29457569b 100644 --- a/src/random_lcg.cpp +++ b/src/random_lcg.cpp @@ -10,7 +10,7 @@ int64_t master_seed {1}; // LCG parameters constexpr uint64_t prn_mult {6364136223846793005ULL}; // multiplication constexpr uint64_t prn_add {1442695040888963407ULL}; // additive factor, c -constexpr uint64_t prn_stride {152917LL}; // stride between particles +uint64_t prn_stride {DEFAULT_STRIDE}; // stride between particles //============================================================================== // PRN @@ -133,4 +133,14 @@ extern "C" void openmc_set_seed(int64_t new_seed) master_seed = new_seed; } +extern "C" uint64_t openmc_get_stride() +{ + return prn_stride; +} + +extern "C" void openmc_set_stride(uint64_t new_stride) +{ + prn_stride = new_stride; +} + } // namespace openmc diff --git a/src/settings.cpp b/src/settings.cpp index e2f5a033b..cd925be70 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -514,6 +514,12 @@ void read_settings_xml(pugi::xml_node root) openmc_set_seed(seed); } + // Copy random number stride if specified + if (check_for_node(root, "stride")) { + auto stride = std::stoull(get_node_value(root, "stride")); + openmc_set_stride(stride); + } + // Check for electron treatment if (check_for_node(root, "electron_treatment")) { auto temp_str = get_node_value(root, "electron_treatment", true, true); diff --git a/src/state_point.cpp b/src/state_point.cpp index ed8c6ed41..3b822715c 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -89,6 +89,9 @@ extern "C" int openmc_statepoint_write(const char* filename, bool* write_source) // Write out random number seed write_dataset(file_id, "seed", openmc_get_seed()); + // Write out random number stride + write_dataset(file_id, "stride", openmc_get_stride()); + // Write run information write_dataset(file_id, "energy_mode", settings::run_CE ? "continuous-energy" : "multi-group"); @@ -399,6 +402,11 @@ extern "C" int openmc_statepoint_load(const char* filename) read_dataset(file_id, "seed", seed); openmc_set_seed(seed); + // Read and overwrite random number stride + uint64_t stride; + read_dataset(file_id, "stride", stride); + openmc_set_stride(stride); + // It is not impossible for a state point to be generated from a CE run but // to be loaded in to an MG run (or vice versa), check to prevent that. read_dataset(file_id, "energy_mode", word); diff --git a/tests/regression_tests/stride/__init__.py b/tests/regression_tests/stride/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/regression_tests/stride/inputs_true.dat b/tests/regression_tests/stride/inputs_true.dat new file mode 100644 index 000000000..f93ec33d1 --- /dev/null +++ b/tests/regression_tests/stride/inputs_true.dat @@ -0,0 +1,25 @@ + + + + + + + + + + + + + + eigenvalue + 1000 + 10 + 5 + + + -4 -4 -4 4 4 4 + + + 1529170 + + diff --git a/tests/regression_tests/stride/results_true.dat b/tests/regression_tests/stride/results_true.dat new file mode 100644 index 000000000..a65411150 --- /dev/null +++ b/tests/regression_tests/stride/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +2.978080E-01 6.106774E-03 diff --git a/tests/regression_tests/stride/test.py b/tests/regression_tests/stride/test.py new file mode 100644 index 000000000..f911af1f5 --- /dev/null +++ b/tests/regression_tests/stride/test.py @@ -0,0 +1,29 @@ +import pytest +import openmc + +from tests.testing_harness import PyAPITestHarness + + +@pytest.fixture +def model(): + u = openmc.Material() + u.add_nuclide('U235', 1.0) + u.set_density('g/cm3', 4.5) + sph = openmc.Sphere(r=10.0, boundary_type='vacuum') + cell = openmc.Cell(fill=u, region=-sph) + model = openmc.Model() + model.geometry = openmc.Geometry([cell]) + + model.settings.batches = 10 + model.settings.inactive = 5 + model.settings.particles = 1000 + model.settings.stride = 1_529_170 + model.settings.source = openmc.IndependentSource( + space=openmc.stats.Box([-4, -4, -4], [4, 4, 4]) + ) + return model + + +def test_seed(model): + harness = PyAPITestHarness('statepoint.10.h5', model) + harness.main()