Address @jtramm comments

This commit is contained in:
Paul Romano 2021-05-13 07:18:46 -05:00
parent 43804fc09d
commit f1173e82c1
3 changed files with 19 additions and 17 deletions

View file

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

View file

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

View file

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