From 0ebfae4ef58d74ff49531e544ad7b3ed54df28f9 Mon Sep 17 00:00:00 2001 From: John Tramm Date: Wed, 20 Nov 2019 22:55:40 +0000 Subject: [PATCH] added particle based seed initialization. Now for sampling functions to use the new seeds. --- include/openmc/particle.h | 2 +- include/openmc/random_lcg.h | 6 +++--- src/random_lcg.cpp | 31 +++++++++++++++---------------- src/simulation.cpp | 6 +++++- 4 files changed, 24 insertions(+), 21 deletions(-) diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 38c694946..1e06d9474 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -288,7 +288,7 @@ public: bool write_track_ {false}; // Current PRNG state - uint64_t prn_seed[N_STREAMS]; // current seed + uint64_t prn_seeds[N_STREAMS]; // current seeds int stream; // current RNG stream }; diff --git a/include/openmc/random_lcg.h b/include/openmc/random_lcg.h index d1d3b5259..5f0fea494 100644 --- a/include/openmc/random_lcg.h +++ b/include/openmc/random_lcg.h @@ -25,7 +25,7 @@ constexpr int64_t DEFAULT_SEED = 1; //============================================================================== //extern "C" double prn(); -extern "C" double prn(uint64_t * seed); +extern "C" double prn(uint64_t * seeds, int stream); //============================================================================== //! Generate a random number which is 'n' times ahead from the current seed. @@ -44,14 +44,14 @@ extern "C" double future_prn(int64_t n, uint64_t seed); //! @param id The particle ID //============================================================================== -extern "C" void set_particle_seed(int64_t id, uint64_t * prn_seed ); +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 n The number of RNG seeds to skip ahead by //============================================================================== -extern "C" void advance_prn_seed(int64_t n, uint64_t * seed); +extern "C" void advance_prn_seed(int64_t n, uint64_t * prn_seeds, int stream); //============================================================================== //! Advance a random number seed 'n' times. diff --git a/src/random_lcg.cpp b/src/random_lcg.cpp index 497272ba7..fdd1547c6 100644 --- a/src/random_lcg.cpp +++ b/src/random_lcg.cpp @@ -16,7 +16,7 @@ extern "C" const int STREAM_VOLUME {4}; extern "C" const int STREAM_PHOTON {5}; // Starting seed -int64_t seed {1}; +int64_t master_seed {1}; // LCG parameters constexpr uint64_t prn_mult {2806196910506780709LL}; // multiplication @@ -39,17 +39,15 @@ constexpr double prn_norm {1.0 / prn_mod}; // 2^-63 //============================================================================== extern "C" double -prn(uint64_t * prn_seed) +prn(uint64_t * prn_seeds, int stream) { // This algorithm uses bit-masking to find the next integer(8) value to be // used to calculate the random number. - //prn_seed[stream] = (prn_mult*prn_seed[stream] + prn_add) & prn_mask; - *prn_seed = (prn_mult * (*prn_seed) + prn_add) & prn_mask; + prn_seeds[stream] = (prn_mult*prn_seeds[stream] + prn_add) & prn_mask; // Once the integer is calculated, we just need to divide by 2**m, // represented here as multiplying by a pre-calculated factor - //return prn_seed[stream] * prn_norm; - return (*prn_seed) * prn_norm; + return prn_seeds[stream] * prn_norm; } //============================================================================== @@ -57,10 +55,9 @@ prn(uint64_t * prn_seed) //============================================================================== extern "C" double -future_prn(int64_t n, uint64_t prn_seed) +future_prn(int64_t n, uint64_t * prn_seeds, int stream) { - //return future_seed(static_cast(n), prn_seed[stream]) * prn_norm; - return future_seed(static_cast(n), *prn_seed) * prn_norm; + return future_seed(static_cast(n), prn_seeds[stream]) * prn_norm; } //============================================================================== @@ -69,10 +66,10 @@ future_prn(int64_t n, uint64_t prn_seed) extern "C" void //set_particle_seed(int64_t id) -set_particle_seed(int64_t id, uint64_t * prn_seed) +set_particle_seed(int64_t id, uint64_t * prn_seeds) { for (int i = 0; i < N_STREAMS; i++) { - prn_seed[i] = future_seed(static_cast(id) * prn_stride, seed + i); + prn_seeds[i] = future_seed(static_cast(id) * prn_stride, master_seed + i); } } @@ -82,10 +79,9 @@ set_particle_seed(int64_t id, uint64_t * prn_seed) extern "C" void //advance_prn_seed(int64_t n) -advance_prn_seed(int64_t n, uint64_t * prn_seed) +advance_prn_seed(int64_t n, uint64_t * prn_seeds, int stream) { - //prn_seed[stream] = future_seed(static_cast(n), prn_seed[stream]); - *prn_seed = future_seed(static_cast(n), *prn_seed); + prn_seeds[stream] = future_seed(static_cast(n), prn_seeds[stream]); } //============================================================================== @@ -143,12 +139,14 @@ prn_set_stream(int i) // API FUNCTIONS //============================================================================== -extern "C" int64_t openmc_get_seed() {return seed;} +extern "C" int64_t openmc_get_seed() {return master_seed;} +// TODO: The idea of this function no longer works -- we'll need to update it somehow extern "C" void openmc_set_seed(int64_t new_seed) { - seed = new_seed; + master_seed = new_seed; + /* //#pragma omp parallel { for (int i = 0; i < N_STREAMS; i++) { @@ -156,6 +154,7 @@ openmc_set_seed(int64_t new_seed) } prn_set_stream(STREAM_TRACKING); } + */ } } // namespace openmc diff --git a/src/simulation.cpp b/src/simulation.cpp index 8de1e7fbc..90b67ab8b 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -473,10 +473,14 @@ 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); + set_particle_seed(particle_seed, p->prn_seeds); // set particle trace simulation::trace = false;