crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 05:35:49 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

added particle based seed initialization. Now for sampling functions to use the new seeds.

Browse source
This commit is contained in:
John Tramm 2019-11-20 22:55:40 +00:00
parent 5c7e7ff13f
commit 0ebfae4ef5
4 changed files with 24 additions and 21 deletions
Download patch file Download diff file Expand all files Collapse all files

2
include/openmc/particle.h
View file

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

6
include/openmc/random_lcg.h
View file

@ -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.

31
src/random_lcg.cpp
View file

@ -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<uint64_t>(n), prn_seed[stream]) * prn_norm;
return future_seed(static_cast<uint64_t>(n), *prn_seed) * prn_norm;
return future_seed(static_cast<uint64_t>(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<uint64_t>(id) * prn_stride, seed + i);
prn_seeds[i] = future_seed(static_cast<uint64_t>(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<uint64_t>(n), prn_seed[stream]);
*prn_seed = future_seed(static_cast<uint64_t>(n), *prn_seed);
prn_seeds[stream] = future_seed(static_cast<uint64_t>(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

6
src/simulation.cpp
View file

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