diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index e395ada6826..55cd3b4f74a 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -514,6 +514,15 @@ pseudo-random number generator. *Default*: 1 +------------------ +```` Element +------------------ + +The ``stride`` element is used to specify how many random numbers are allocated +for each particle history. + + *Default*: 152,917 + .. _source_element: -------------------- diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 9401156a64f..3d44a9c6bae 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -68,6 +68,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); @@ -134,6 +135,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 4157b7cfe7a..fd2c4cc497b 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 {152917}; //============================================================================== //! 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 062670ef843..4fba8d48b6e 100644 --- a/openmc/lib/settings.py +++ b/openmc/lib/settings.py @@ -12,6 +12,8 @@ _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 @@ def seed(self): 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 fae6f638364..746bb25bd04 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -174,6 +174,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 particle history source : Iterable of openmc.SourceBase Distribution of source sites in space, angle, and energy sourcepoint : dict @@ -312,6 +314,7 @@ def __init__(self, **kwargs): self._plot_seed = None self._ptables = None self._seed = None + self._stride = None self._survival_biasing = None # Shannon entropy mesh @@ -578,6 +581,16 @@ def seed(self, seed: int): 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 @@ -1266,6 +1279,11 @@ def _create_seed_subelement(self, root): 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") @@ -1681,6 +1699,11 @@ def _seed_from_xml_element(self, root): 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: @@ -1912,6 +1935,7 @@ def to_xml_element(self, mesh_memo=None): 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) @@ -2018,6 +2042,7 @@ def from_xml_element(cls, elem, meshes=None): 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 13aac0caf7e..fae291d655c 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -98,6 +98,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 @@ -355,6 +357,10 @@ def runtime(self): 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 2bde0e3387c..01135f89701 100644 --- a/src/finalize.cpp +++ b/src/finalize.cpp @@ -153,6 +153,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(); @@ -215,5 +216,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 cc1eac9cf35..33a02cdfd77 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -98,9 +98,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 581d696176f..50d93a7bd80 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 {152917LL}; //============================================================================== // PRN @@ -132,4 +132,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 6c3226b04a6..2e74435062e 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -465,6 +465,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 c7b7d6ad85c..ee59a606092 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -90,6 +90,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"); @@ -400,6 +403,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 00000000000..e69de29bb2d diff --git a/tests/regression_tests/stride/geometry.xml b/tests/regression_tests/stride/geometry.xml new file mode 100644 index 00000000000..bc56030e18b --- /dev/null +++ b/tests/regression_tests/stride/geometry.xml @@ -0,0 +1,8 @@ + + + + + + + + diff --git a/tests/regression_tests/stride/materials.xml b/tests/regression_tests/stride/materials.xml new file mode 100644 index 00000000000..2472a747174 --- /dev/null +++ b/tests/regression_tests/stride/materials.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/regression_tests/stride/results_true.dat b/tests/regression_tests/stride/results_true.dat new file mode 100644 index 00000000000..a654111500a --- /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/settings.xml b/tests/regression_tests/stride/settings.xml new file mode 100644 index 00000000000..5c79943330a --- /dev/null +++ b/tests/regression_tests/stride/settings.xml @@ -0,0 +1,17 @@ + + + + 1529170 + + eigenvalue + 10 + 5 + 1000 + + + + -4 -4 -4 4 4 4 + + + + diff --git a/tests/regression_tests/stride/test.py b/tests/regression_tests/stride/test.py new file mode 100644 index 00000000000..8b2f031b3e0 --- /dev/null +++ b/tests/regression_tests/stride/test.py @@ -0,0 +1,6 @@ +from tests.testing_harness import TestHarness + + +def test_seed(): + harness = TestHarness('statepoint.10.h5') + harness.main()