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

changed basic prototypes

Browse source
This commit is contained in:
John Tramm 2019-11-20 19:20:03 +00:00
parent 50a271bfda
commit 5c7e7ff13f
3 changed files with 32 additions and 17 deletions
Download patch file Download diff file Expand all files Collapse all files

5
include/openmc/particle.h
View file

@ -13,6 +13,7 @@
#include "openmc/constants.h"
#include "openmc/position.h"
#include "openmc/random_lcg.h"
namespace openmc {
@ -285,6 +286,10 @@ public:
// Track output
bool write_track_ {false};
// Current PRNG state
uint64_t prn_seed[N_STREAMS]; // current seed
int stream; // current RNG stream
};
} // namespace openmc

12
include/openmc/random_lcg.h
View file

@ -24,7 +24,8 @@ constexpr int64_t DEFAULT_SEED = 1;
//! @return A random number between 0 and 1
//==============================================================================
extern "C" double prn();
//extern "C" double prn();
extern "C" double prn(uint64_t * seed);
//==============================================================================
//! Generate a random number which is 'n' times ahead from the current seed.
@ -35,21 +36,22 @@ extern "C" double prn();
//! @return A random number between 0 and 1
//==============================================================================
extern "C" double future_prn(int64_t n);
//extern "C" double future_prn(int64_t n);
extern "C" double future_prn(int64_t n, uint64_t seed);
//==============================================================================
//! Set the RNG seed to a unique value based on the ID of the particle.
//! @param id The particle ID
//==============================================================================
extern "C" void set_particle_seed(int64_t id);
extern "C" void set_particle_seed(int64_t id, uint64_t * prn_seed );
//==============================================================================
//! 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);
extern "C" void advance_prn_seed(int64_t n, uint64_t * seed);
//==============================================================================
//! Advance a random number seed 'n' times.
@ -73,7 +75,7 @@ uint64_t future_seed(uint64_t n, uint64_t seed);
//! `STREAM_TRACKING` and `STREAM_TALLIES` for this argument.
//==============================================================================
extern "C" void prn_set_stream(int n);
//extern "C" void prn_set_stream(int n);
//==============================================================================
// API FUNCTIONS

32
src/random_lcg.cpp
View file

@ -29,9 +29,9 @@ constexpr uint64_t prn_stride {152917LL}; // stride between
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)
//uint64_t prn_seed[N_STREAMS]; // current seed
//int stream; // current RNG stream
//#pragma omp threadprivate(prn_seed, stream)
//==============================================================================
@ -39,15 +39,17 @@ int stream; // current RNG stream
//==============================================================================
extern "C" double
prn()
prn(uint64_t * prn_seed)
{
// 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[stream] = (prn_mult*prn_seed[stream] + prn_add) & prn_mask;
*prn_seed = (prn_mult * (*prn_seed) + 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[stream] * prn_norm;
return (*prn_seed) * prn_norm;
}
//==============================================================================
@ -55,9 +57,10 @@ prn()
//==============================================================================
extern "C" double
future_prn(int64_t n)
future_prn(int64_t n, uint64_t prn_seed)
{
return future_seed(static_cast<uint64_t>(n), prn_seed[stream]) * prn_norm;
//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;
}
//==============================================================================
@ -65,7 +68,8 @@ future_prn(int64_t n)
//==============================================================================
extern "C" void
set_particle_seed(int64_t id)
//set_particle_seed(int64_t id)
set_particle_seed(int64_t id, uint64_t * prn_seed)
{
for (int i = 0; i < N_STREAMS; i++) {
prn_seed[i] = future_seed(static_cast<uint64_t>(id) * prn_stride, seed + i);
@ -77,9 +81,11 @@ set_particle_seed(int64_t id)
//==============================================================================
extern "C" void
advance_prn_seed(int64_t n)
//advance_prn_seed(int64_t n)
advance_prn_seed(int64_t n, uint64_t * prn_seed)
{
prn_seed[stream] = future_seed(static_cast<uint64_t>(n), prn_seed[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);
}
//==============================================================================
@ -125,11 +131,13 @@ future_seed(uint64_t n, uint64_t seed)
// 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
@ -141,7 +149,7 @@ extern "C" void
openmc_set_seed(int64_t new_seed)
{
seed = new_seed;
#pragma omp parallel
//#pragma omp parallel
{
for (int i = 0; i < N_STREAMS; i++) {
prn_seed[i] = seed + i;