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

Incorporated maxwell and watt_spectrum C++ code, tests pass

Browse source
This commit is contained in:
Adam G Nelson 2018-05-06 11:28:05 -04:00
parent 2000b4f998
commit feecf0134e
4 changed files with 16 additions and 26 deletions
Download patch file Download diff file Expand all files Collapse all files

18
src/math.F90
View file

@ -364,18 +364,7 @@ contains
real(C_DOUBLE), intent(in) :: T ! tabulated function of incoming E
real(C_DOUBLE) :: E_out ! sampled energy
real(C_DOUBLE) :: r1, r2, r3 ! random numbers
real(C_DOUBLE) :: c ! cosine of pi/2*r3
r1 = prn()
r2 = prn()
r3 = prn()
! determine cosine of pi/2*r
c = cos(PI/TWO*r3)
! determine outgoing energy
E_out = -T*(log(r1) + log(r2)*c*c)
E_out = maxwell_spectrum_cc(T)
end function maxwell_spectrum
@ -393,10 +382,7 @@ contains
real(C_DOUBLE), intent(in) :: b ! Watt parameter b
real(C_DOUBLE) :: E_out ! energy of emitted neutron
real(C_DOUBLE) :: w ! sampled from Maxwellian
w = maxwell_spectrum(a)
E_out = w + a*a*b/4. + (TWO*prn() - ONE)*sqrt(a*a*b*w)
E_out = watt_spectrum_cc(a, b)
end function watt_spectrum

12
src/math_functions.cpp
View file

@ -56,7 +56,7 @@ double __attribute__ ((const)) normal_percentile_c(double p) {
// Refinement based on Newton's method
z = z - (0.5 * std::erfc(-z / std::sqrt(2.0)) - p) * std::sqrt(2.0 * M_PI) *
z = z - (0.5 * std::erfc(-z / std::sqrt(2.0)) - p) * std::sqrt(2.0 * PI) *
std::exp(0.5 * z * z);
return z;
@ -82,7 +82,7 @@ double __attribute__ ((const)) t_percentile_c(double p, int df){
// For one degree of freedom, the t-distribution becomes a Cauchy
// distribution whose cdf we can invert directly
t = std::tan(M_PI*(p - 0.5));
t = std::tan(PI*(p - 0.5));
} else if (df == 2) {
// For two degrees of freedom, the cdf is given by 1/2 + x/(2*sqrt(x^2 +
// 2)). This can be directly inverted to yield the solution below
@ -604,7 +604,7 @@ void rotate_angle_c(double uvw[3], double mu, double phi) {
if (phi != -10.) {
phi_ = phi;
} else {
phi_ = 2. * M_PI * prn();
phi_ = 2. * PI * prn();
}
// Precompute factors to save flops
@ -647,7 +647,7 @@ double maxwell_spectrum_c(double T) {
r3 = prn();
// determine cosine of pi/2*r
c = std::cos(M_PI / 2. * r3);
c = std::cos(PI / 2. * r3);
// determine outgoing energy
E_out = -T * (std::log(r1) + std::log(r2) * c * c);
@ -710,7 +710,7 @@ double watt_spectrum_c(double a, double b) {
// std::complex<double> w_derivative_c(std::complex<double> z, int order){
// std::complex<double> wv; // The resultant w(z) value
// const std::complex<double> twoi_sqrtpi(0.0, 2.0 / std::sqrt(M_PI));
// const std::complex<double> twoi_sqrtpi(0.0, 2.0 / std::sqrt(PI));
// switch(order) {
// case 0:
@ -764,7 +764,7 @@ void broaden_wmp_polynomials_c(double E, double dopp, int n, double factors[]) {
factors[0] = erf_beta / E;
factors[1] = 1. / sqrtE;
factors[2] = factors[0] * (half_inv_dopp2 + E) + exp_m_beta2 /
(beta * std::sqrt(M_PI));
(beta * std::sqrt(PI));
// Perform recursive broadening of high order components
for (i = 0; i < n - 3; i++) {

4
src/math_functions.h
View file

@ -11,6 +11,10 @@
namespace openmc {
// TODO: cmath::M_PI has 3 more digits precision than the Fortran constant we
// use so for now we will reuse the Fortran constant until we are OK with
// modifying test results
const double PI = 3.1415926535898;
//==============================================================================
// NORMAL_PERCENTILE calculates the percentile of the standard normal

8
tests/unit_tests/test_math.py
View file

@ -147,8 +147,8 @@ def test_maxwell_spectrum():
T = 0.5
ref_val = 0.6129982175261098
test_val = openmc.capi.math.maxwell_spectrum(T)
print(test_val)
assert np.isclose(ref_val, test_val)
assert ref_val == test_val
def test_watt_spectrum():
@ -158,8 +158,8 @@ def test_watt_spectrum():
b = 0.75
ref_val = 0.6247242713640233
test_val = openmc.capi.math.watt_spectrum(a, b)
print(test_val)
assert np.isclose(ref_val, test_val)
assert ref_val == test_val
def test_broaden_wmp_polynomials():