diff --git a/include/openmc/random_lcg.h b/include/openmc/random_lcg.h index 4a6e13a83..b99f914a4 100644 --- a/include/openmc/random_lcg.h +++ b/include/openmc/random_lcg.h @@ -21,10 +21,11 @@ constexpr int64_t DEFAULT_SEED = 1; //============================================================================== //! Generate a pseudo-random number using a linear congruential generator. +//! @param prn_seeds Pseudorandom number seed array +//! @param stream Pseudorandom number stream index //! @return A random number between 0 and 1 //============================================================================== -//extern "C" double prn(); extern "C" double prn(uint64_t * seeds, int stream); //============================================================================== @@ -33,14 +34,16 @@ extern "C" double prn(uint64_t * seeds, int stream); //! The result of this function will be the same as the result from calling //! `prn()` 'n' times. //! @param n The number of RNG seeds to skip ahead by +//! @param prn_seeds Pseudorandom number seed array +//! @param stream Pseudorandom number stream index //! @return A random number between 0 and 1 //============================================================================== -//extern "C" double future_prn(int64_t n); extern "C" double future_prn(int64_t n, uint64_t * prn_seeds, int stream); //============================================================================== -//! Set the RNG seed to a unique value based on the ID of the particle. +//! Set the RNG seeds to unique values based on the ID of the particle. +//! @param prn_seeds Pseudorandom number seed array //! @param id The particle ID //============================================================================== @@ -48,6 +51,8 @@ extern "C" void set_particle_seed(int64_t id, uint64_t * prn_seeds ); //============================================================================== //! Advance the random number seed 'n' times from the current seed. +//! @param prn_seeds Pseudorandom number seed array +//! @param stream Pseudorandom number stream index //! @param n The number of RNG seeds to skip ahead by //============================================================================== @@ -65,18 +70,6 @@ extern "C" void advance_prn_seed(int64_t n, uint64_t * prn_seeds, int stream); uint64_t future_seed(uint64_t n, uint64_t seed); -//============================================================================== -//! Switch the RNG to a different stream of random numbers. -//! -//! If random numbers are needed in routines not used directly for tracking -//! (e.g. physics), this allows the numbers to be generated without affecting -//! reproducibility of the physics. -//! @param n The RNG stream to switch to. Use the constants such as -//! `STREAM_TRACKING` and `STREAM_TALLIES` for this argument. -//============================================================================== - -//extern "C" void prn_set_stream(int n); - //============================================================================== // API FUNCTIONS //============================================================================== diff --git a/openmc/lib/math.py b/openmc/lib/math.py index 1baa9cf64..781402c46 100644 --- a/openmc/lib/math.py +++ b/openmc/lib/math.py @@ -1,4 +1,4 @@ -from ctypes import (c_int, c_double, POINTER) +from ctypes import (c_int, c_double, POINTER, c_uint64) import numpy as np from numpy.ctypeslib import ndpointer @@ -26,9 +26,9 @@ _dll.calc_zn_rad.argtypes = [c_int, c_double, ndpointer(c_double)] _dll.rotate_angle_c.restype = None _dll.rotate_angle_c.argtypes = [ndpointer(c_double), c_double, - POINTER(c_double)] + POINTER(c_double), POINTER(c_uint64), c_int] _dll.maxwell_spectrum.restype = c_double -_dll.maxwell_spectrum.argtypes = [c_double] +_dll.maxwell_spectrum.argtypes = [c_double, POINTER(c_uint64), c_int] _dll.watt_spectrum.restype = c_double _dll.watt_spectrum.argtypes = [c_double, c_double] @@ -38,7 +38,7 @@ _dll.broaden_wmp_polynomials.argtypes = [c_double, c_double, c_int, ndpointer(c_double)] _dll.normal_variate.restype = c_double -_dll.normal_variate.argtypes = [c_double, c_double] +_dll.normal_variate.argtypes = [c_double, c_double, POINTER(c_uint64), c_int] def t_percentile(p, df): """ Calculate the percentile of the Student's t distribution with a diff --git a/src/random_lcg.cpp b/src/random_lcg.cpp index c3b0c6a62..8c6445ccb 100644 --- a/src/random_lcg.cpp +++ b/src/random_lcg.cpp @@ -27,12 +27,6 @@ constexpr uint64_t prn_stride {152917LL}; // stride between // particles constexpr double prn_norm {1.0 / prn_mod}; // 2^-63 -// Current PRNG state -//uint64_t prn_seed[N_STREAMS]; // current seed -//int stream; // current RNG stream -//#pragma omp threadprivate(prn_seed, stream) - - //============================================================================== // PRN //============================================================================== @@ -64,7 +58,6 @@ future_prn(int64_t n, uint64_t * prn_seeds, int stream) //============================================================================== extern "C" void -//set_particle_seed(int64_t id) set_particle_seed(int64_t id, uint64_t * prn_seeds) { for (int i = 0; i < N_STREAMS; i++) { @@ -77,7 +70,6 @@ set_particle_seed(int64_t id, uint64_t * prn_seeds) //============================================================================== extern "C" void -//advance_prn_seed(int64_t n) advance_prn_seed(int64_t n, uint64_t * prn_seeds, int stream) { prn_seeds[stream] = future_seed(static_cast(n), prn_seeds[stream]); @@ -122,18 +114,6 @@ future_seed(uint64_t n, uint64_t seed) return (g_new * seed + c_new) & prn_mask; } -//============================================================================== -// PRN_SET_STREAM -//============================================================================== - -/* -extern "C" void -prn_set_stream(int i) -{ - stream = i; // Shift by one to move from Fortran to C indexing. -} -*/ - //============================================================================== // API FUNCTIONS //============================================================================== @@ -144,15 +124,6 @@ extern "C" void openmc_set_seed(int64_t new_seed) { master_seed = new_seed; - /* - //#pragma omp parallel - { - for (int i = 0; i < N_STREAMS; i++) { - prn_seed[i] = seed + i; - } - prn_set_stream(STREAM_TRACKING); - } - */ } } // namespace openmc diff --git a/src/simulation.cpp b/src/simulation.cpp index 90b67ab8b..592878426 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -473,13 +473,9 @@ void initialize_history(Particle* p, int64_t index_source) // set identifier for particle p->id_ = simulation::work_index[mpi::rank] + index_source; - // TODO: Why doesn't this include the particle skip offset? - // ANSWER: Stride is included in the set seed function, which actually forwards based - // on the master_seed. // set random number seed int64_t particle_seed = (simulation::total_gen + overall_generation() - 1) * settings::n_particles + p->id_; - //set_particle_seed(particle_seed); set_particle_seed(particle_seed, p->prn_seeds); // set particle trace diff --git a/tests/unit_tests/test_math.py b/tests/unit_tests/test_math.py index 7c0218fbd..359181d23 100644 --- a/tests/unit_tests/test_math.py +++ b/tests/unit_tests/test_math.py @@ -156,11 +156,13 @@ def test_rotate_angle(): uvw0 = np.array([1., 0., 0.]) phi = 0. mu = 0. + prn_seeds = np.array([1, 2, 3, 4, 5, 6], dtype=np.uint64) + stream = 0 # reference: mu of 0 pulls the vector the bottom, so: ref_uvw = np.array([0., 0., -1.]) - test_uvw = openmc.lib.math.rotate_angle(uvw0, mu, phi) + test_uvw = openmc.lib.math.rotate_angle(uvw0, mu, phi, prn_seeds, stream) assert np.array_equal(ref_uvw, test_uvw) @@ -168,59 +170,58 @@ def test_rotate_angle(): mu = 1. ref_uvw = np.array([1., 0., 0.]) - test_uvw = openmc.lib.math.rotate_angle(uvw0, mu, phi) + test_uvw = openmc.lib.math.rotate_angle(uvw0, mu, phi, prn_seeds, stream) assert np.array_equal(ref_uvw, test_uvw) # Now to test phi is None mu = 0.9 - settings = openmc.lib.settings - settings.seed = 1 # When seed = 1, phi will be sampled as 1.9116495709698769 # The resultant reference is from hand-calculations given the above ref_uvw = [0.9, 0.410813051297112, 0.1457142302040] - test_uvw = openmc.lib.math.rotate_angle(uvw0, mu) + test_uvw = openmc.lib.math.rotate_angle(uvw0, mu, phi, prn_seeds, stream) assert np.allclose(ref_uvw, test_uvw) def test_maxwell_spectrum(): - settings = openmc.lib.settings - settings.seed = 1 + prn_seeds = np.array([1, 2, 3, 4, 5, 6], dtype=np.uint64) + stream = 0 T = 0.5 ref_val = 0.6129982175261098 - test_val = openmc.lib.math.maxwell_spectrum(T) + test_val = openmc.lib.math.maxwell_spectrum(T, prn_seeds, stream) assert ref_val == test_val def test_watt_spectrum(): - settings = openmc.lib.settings - settings.seed = 1 + prn_seeds = np.array([1, 2, 3, 4, 5, 6], dtype=np.uint64) + stream = 0 a = 0.5 b = 0.75 ref_val = 0.6247242713640233 - test_val = openmc.lib.math.watt_spectrum(a, b) + test_val = openmc.lib.math.watt_spectrum(a, b, prn_seeds, stream) assert ref_val == test_val def test_normal_dist(): - settings = openmc.lib.settings - settings.seed = 1 + prn_seeds = np.array([1, 2, 3, 4, 5, 6], dtype=np.uint64) + stream = 0 a = 14.08 b = 0.0 ref_val = 14.08 - test_val = openmc.lib.math.normal_variate(a, b) + test_val = openmc.lib.math.normal_variate(a, b, prn_seeds, stream) assert ref_val == pytest.approx(test_val) - settings.seed = 1 + prn_seeds = np.array([1, 2, 3, 4, 5, 6], dtype=np.uint64) + stream = 0 a = 14.08 b = 1.0 ref_val = 16.436645416691427 - test_val = openmc.lib.math.normal_variate(a, b) + test_val = openmc.lib.math.normal_variate(a, b, prn_seeds, stream) assert ref_val == pytest.approx(test_val)