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

Respond to my original comments on #1670

Browse source
This commit is contained in:
Paul Romano 2020-12-22 16:20:19 -06:00
parent 9188786094
commit 04a65235fb
2 changed files with 15 additions and 16 deletions
Download patch file Download diff file Expand all files Collapse all files

17
src/math_functions.cpp
View file

@ -692,16 +692,19 @@ double maxwell_spectrum(double T, uint64_t* seed) {
double normal_variate(double mean, double standard_deviation, uint64_t* seed) {
// now correct method to sample a normal distubted number
double v1 = 2 * prn(seed) - 1.;
double v2 = 2 * prn(seed) - 1.;
double r = std::pow(v1, 2) + std::pow(v2, 2);
if ( r == 0 || r > 1 ) return normal_variate(mean,standard_deviation,seed);
double z = std::sqrt(-2.0 * std::log(r)/r);
return mean + standard_deviation*z*v1;
// 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;
} while (r2 > 1 || r2 == 0);
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) {
// https://permalink.lanl.gov/object/tr?what=info:lanl-repo/lareport/LA-05411-MS
double sigma = std::sqrt(4.*e0*kt/m_rat);
return normal_variate(e0, sigma, seed);
}

14
tests/unit_tests/test_math.py
View file

@ -1,10 +1,11 @@
import numpy as np
import pytest
import scipy as sp
from scipy.stats import shapiro
import openmc
import openmc.lib
import pytest
def test_t_percentile():
# Permutations include 1 DoF, 2 DoF, and > 2 DoF
@ -216,21 +217,16 @@ def test_normal_dist():
# make sigma = 1.0
b = 1.0
# import shapiro test
from scipy.stats import shapiro
samples = []
num_samples = 10000
# sample some numbers
for i in range(num_samples):
# sample the normal distribution from openmc
samples.append(openmc.lib.math.normal_variate(a, b, prn_seed))
prn_seed = prn_seed + 1
prn_seed += 1
stat,p = shapiro(samples)
assert stat > 0.97
# cleanup
del samples
stat, p = shapiro(samples)
assert p > 0.05
def test_broaden_wmp_polynomials():