crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 22:26:08 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Use unsigned integers in the RNG

Browse source
This commit is contained in:
Sterling Harper 2017-12-01 17:32:27 -05:00
parent cc52bfc435
commit 6811a19732
2 changed files with 28 additions and 31 deletions
Download patch file Download diff file Expand all files Collapse all files

49
src/random_lcg.C
View file

@ -6,12 +6,12 @@
int64_t seed = 1;
// LCG parameters
const int64_t prn_mult = 2806196910506780709LL; // multiplication factor, g
const int64_t prn_add = 1; // additive factor, c
const int64_t prn_mod = -0x8000000000000000; // -2^63
const int64_t prn_mask = 0x7fffffffffffffff; // 2^63 - 1
const int64_t prn_stride = 152917LL; // stride between particles
const double prn_norm = pow(2, -63); // 2^-63
const uint64_t prn_mult = 2806196910506780709LL; // multiplication factor, g
const uint64_t prn_add = 1; // additive factor, c
const uint64_t prn_mod = 0x8000000000000000; // 2^63
const uint64_t prn_mask = 0x7fffffffffffffff; // 2^63 - 1
const uint64_t prn_stride = 152917LL; // stride between particles
const double prn_norm = pow(2, -63); // 2^-63
// Module constants
const int N_STREAMS = 5;
@ -22,8 +22,8 @@ const int STREAM_URR_PTABLE = 3;
const int STREAM_VOLUME = 4;
// Current PRNG state
int64_t prn_seed[N_STREAMS]; // current seed
int stream; // current RNG stream
uint64_t prn_seed[N_STREAMS]; // current seed
int stream; // current RNG stream
#pragma omp threadprivate(prn_seed, stream)
@ -40,8 +40,7 @@ prn()
// Once the integer is calculated, we just need to divide by 2**m,
// represented here as multiplying by a pre-calculated factor
double pseudo_rn = prn_seed[stream] * prn_norm;
return pseudo_rn;
return prn_seed[stream] * prn_norm;
}
//==============================================================================
@ -49,10 +48,9 @@ prn()
//==============================================================================
extern "C" double
future_prn(int64_t n)
future_prn(uint64_t n)
{
double pseudo_rn = future_seed(n, prn_seed[stream]) * prn_norm;
return pseudo_rn;
return future_seed(n, prn_seed[stream]) * prn_norm;
}
//==============================================================================
@ -60,7 +58,7 @@ future_prn(int64_t n)
//==============================================================================
extern "C" void
set_particle_seed(int64_t id)
set_particle_seed(uint64_t id)
{
for (int i = 0; i < N_STREAMS; i++) {
prn_seed[i] = future_seed(id * prn_stride, seed + i);
@ -72,7 +70,7 @@ set_particle_seed(int64_t id)
//==============================================================================
extern "C" void
advance_prn_seed(int64_t n)
advance_prn_seed(uint64_t n)
{
prn_seed[stream] = future_seed(n, prn_seed[stream]);
}
@ -81,16 +79,16 @@ advance_prn_seed(int64_t n)
// FUTURE_SEED
//==============================================================================
extern "C" int64_t
future_seed(int64_t n, int64_t seed)
extern "C" uint64_t
future_seed(uint64_t n, uint64_t seed)
{
// In cases where we want to skip backwards, we add the period of the random
// number generator until the number of PRNs to skip is positive since
// skipping ahead that much is the same as skipping backwards by the original
// amount.
int64_t nskip = n;
while (nskip < 0) nskip += prn_mod;
uint64_t nskip = n;
while (nskip > prn_mod) nskip += prn_mod;
// Make sure nskip is less than 2^M.
nskip &= prn_mask;
@ -102,10 +100,10 @@ future_seed(int64_t n, int64_t seed)
// and C which can then be used to find x_N = G*x_0 + C mod 2^M.
// Initialize constants
int64_t g = prn_mult;
int64_t c = prn_add;
int64_t g_new = 1;
int64_t c_new = 0;
uint64_t g = prn_mult;
uint64_t c = prn_add;
uint64_t g_new = 1;
uint64_t c_new = 0;
while (nskip > 0) {
// Check if the least significant bit is 1.
@ -121,8 +119,7 @@ future_seed(int64_t n, int64_t seed)
}
// With G and C, we can now find the new seed.
int64_t new_seed = (g_new * seed + c_new) & prn_mask;
return new_seed;
return (g_new * seed + c_new) & prn_mask;
}
//==============================================================================
@ -140,7 +137,7 @@ prn_set_stream(int i)
//==============================================================================
extern "C" int
openmc_set_seed(int64_t new_seed)
openmc_set_seed(uint64_t new_seed)
{
seed = new_seed;
#pragma omp parallel

10
src/random_lcg.h
View file

@ -14,21 +14,21 @@ extern "C" double prn();
// current seed.
//==============================================================================
extern "C" double future_prn(int64_t n);
extern "C" double future_prn(uint64_t n);
//==============================================================================
// SET_PARTICLE_SEED sets the seed to a unique value based on the ID of the
// particle.
//==============================================================================
extern "C" void set_particle_seed(int64_t id);
extern "C" void set_particle_seed(uint64_t id);
//==============================================================================
// ADVANCE_PRN_SEED advances the random number seed 'n' times from the current
// seed.
//==============================================================================
extern "C" void advance_prn_seed(int64_t n);
extern "C" void advance_prn_seed(uint64_t n);
//==============================================================================
// FUTURE_SEED advances the random number seed 'skip' times. This is usually
@ -37,7 +37,7 @@ extern "C" void advance_prn_seed(int64_t n);
// are used.
//==============================================================================
extern "C" int64_t future_seed(int64_t n, int64_t seed);
extern "C" uint64_t future_seed(uint64_t n, uint64_t seed);
//==============================================================================
// PRN_SET_STREAM changes the random number stream. If random numbers are needed
@ -51,6 +51,6 @@ extern "C" void prn_set_stream(int n);
// API FUNCTIONS
//==============================================================================
extern "C" int openmc_set_seed(int64_t new_seed);
extern "C" int openmc_set_seed(uint64_t new_seed);
#endif // RANDOM_LCG_H