Implement user-configurable random number stride

docs/source/io_formats/settings.rst

@@ -514,6 +514,15 @@ pseudo-random number generator.
 
   *Default*: 1
 
+------------------
+``<stride>`` Element
+------------------
+
+The ``stride`` element is used to specify how many random numbers are allocated
+for each particle history.
+
+  *Default*: 152,917
+
 .. _source_element:
 
 --------------------

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();

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

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
 

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)

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

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;
 }

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

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

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);

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);

tests/regression_tests/stride/geometry.xml

@@ -0,0 +1,8 @@
+<?xml version="1.0"?>
+<geometry>
+
+  <!-- Sphere with radius 10 -->
+  <surface id="1" type="sphere" coeffs="0 0 0 10" boundary="vacuum"/>
+  <cell id="1" material="1" region="-1" />
+
+</geometry>

tests/regression_tests/stride/materials.xml

@@ -0,0 +1,9 @@
+<?xml version="1.0"?>
+<materials>
+
+  <material id="1">
+    <density value="4.5" units="g/cc" />
+    <nuclide name="U235" ao="1.0" />
+  </material>
+
+</materials>

tests/regression_tests/stride/results_true.dat

@@ -0,0 +1,2 @@
+k-combined:
+2.978080E-01 6.106774E-03

tests/regression_tests/stride/settings.xml

@@ -0,0 +1,17 @@
+<?xml version="1.0"?>
+<settings>
+
+  <stride>1529170</stride>
+
+  <run_mode>eigenvalue</run_mode>
+  <batches>10</batches>
+  <inactive>5</inactive>
+  <particles>1000</particles>
+
+  <source>
+    <space type="box">
+      <parameters>-4 -4 -4  4  4  4</parameters>
+    </space>
+  </source>
+
+</settings>

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()