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

removed some comments. Also attempted to fix the math unit tests that now require direct passing of seeds and stream variables

Browse source
This commit is contained in:
John Tramm 2019-11-22 16:17:46 +00:00
parent 3b985f1976
commit ecb5f828d2
5 changed files with 29 additions and 68 deletions
Download patch file Download diff file Expand all files Collapse all files

23
include/openmc/random_lcg.h
View file

@ -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
//==============================================================================

8
openmc/lib/math.py
View file

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

29
src/random_lcg.cpp
View file

@ -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<uint64_t>(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

4
src/simulation.cpp
View file

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

33
tests/unit_tests/test_math.py
View file

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