From f1173e82c1c5786e0ad4debdb5b0ae7d09c45d8a Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 13 May 2021 07:18:46 -0500 Subject: [PATCH] Address @jtramm comments --- src/random_dist.cpp | 32 +++++++++++++++++--------------- src/random_lcg.cpp | 2 +- tests/unit_tests/test_lib.py | 2 +- 3 files changed, 19 insertions(+), 17 deletions(-) diff --git a/src/random_dist.cpp b/src/random_dist.cpp index 3966cf5d82..c8d2380c36 100644 --- a/src/random_dist.cpp +++ b/src/random_dist.cpp @@ -9,11 +9,11 @@ namespace openmc { double uniform_distribution(double a, double b, uint64_t* seed) { - return a + (b - a)*prn(seed); + return a + (b - a) * prn(seed); } - -double maxwell_spectrum(double T, uint64_t* seed) { +double maxwell_spectrum(double T, uint64_t* seed) +{ // Set the random numbers double r1 = prn(seed); double r2 = prn(seed); @@ -26,28 +26,30 @@ double maxwell_spectrum(double T, uint64_t* seed) { return -T * (std::log(r1) + std::log(r2) * c * c); } -double watt_spectrum(double a, double b, uint64_t* seed) { +double watt_spectrum(double a, double b, uint64_t* seed) +{ double w = maxwell_spectrum(a, seed); - return w + 0.25 * a * a * b + (2. * prn(seed) - 1.) * std::sqrt(a * a * b * w); + return w + 0.25 * a * a * b + + uniform_distribution(-1., 1., seed) * std::sqrt(a * a * b * w); } - -double normal_variate(double mean, double standard_deviation, uint64_t* seed) { +double normal_variate(double mean, double standard_deviation, uint64_t* seed) +{ // Sample a normal variate using Marsaglia's polar method double x, y, r2; do { - x = 2 * prn(seed) - 1.; - y = 2 * prn(seed) - 1.; - r2 = x*x + y*y; + x = uniform_distribution(-1., 1., seed); + y = uniform_distribution(-1., 1., seed); + r2 = x * x + y * y; } while (r2 > 1 || r2 == 0); - double z = std::sqrt(-2.0 * std::log(r2)/r2); - return mean + standard_deviation*z*x; + double z = std::sqrt(-2.0 * std::log(r2) / r2); + return mean + standard_deviation * z * x; } - -double muir_spectrum(double e0, double m_rat, double kt, uint64_t* seed) { +double muir_spectrum(double e0, double m_rat, double kt, uint64_t* seed) +{ // https://permalink.lanl.gov/object/tr?what=info:lanl-repo/lareport/LA-05411-MS - double sigma = std::sqrt(4.*e0*kt/m_rat); + double sigma = std::sqrt(4. * e0 * kt / m_rat); return normal_variate(e0, sigma, seed); } diff --git a/src/random_lcg.cpp b/src/random_lcg.cpp index c3bd0d5add..38174be657 100644 --- a/src/random_lcg.cpp +++ b/src/random_lcg.cpp @@ -11,7 +11,7 @@ int64_t master_seed {1}; // LCG parameters constexpr uint64_t prn_mult {6364136223846793005ULL}; // multiplication constexpr uint64_t prn_add {1442695040888963407ULL}; // additive factor, c -constexpr uint64_t prn_stride {152917LL}; // stride between +constexpr uint64_t prn_stride {152917LL}; // stride between particles //============================================================================== // PRN diff --git a/tests/unit_tests/test_lib.py b/tests/unit_tests/test_lib.py index d1c07ce716..87d4843c77 100644 --- a/tests/unit_tests/test_lib.py +++ b/tests/unit_tests/test_lib.py @@ -390,7 +390,7 @@ def test_reset(lib_run): for i in range(20): openmc.lib.next_batch() - # Make sure there are 5 realizations for the 5 active batches. + # Make sure there are 15 realizations for the 15 active batches. assert openmc.lib.num_realizations() == 15 assert openmc.lib.tallies[2].num_realizations == 15 _, keff_sd1 = openmc.lib.keff()