From feecf0134eda55b9cac11107c55ef46c7450ac4f Mon Sep 17 00:00:00 2001 From: Adam G Nelson Date: Sun, 6 May 2018 11:28:05 -0400 Subject: [PATCH] Incorporated maxwell and watt_spectrum C++ code, tests pass --- src/math.F90 | 18 ++---------------- src/math_functions.cpp | 12 ++++++------ src/math_functions.h | 4 ++++ tests/unit_tests/test_math.py | 8 ++++---- 4 files changed, 16 insertions(+), 26 deletions(-) diff --git a/src/math.F90 b/src/math.F90 index 6534b66b03..d7e335284c 100644 --- a/src/math.F90 +++ b/src/math.F90 @@ -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 diff --git a/src/math_functions.cpp b/src/math_functions.cpp index 6ddf944a76..1466bde547 100644 --- a/src/math_functions.cpp +++ b/src/math_functions.cpp @@ -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 w_derivative_c(std::complex z, int order){ // std::complex wv; // The resultant w(z) value -// const std::complex twoi_sqrtpi(0.0, 2.0 / std::sqrt(M_PI)); +// const std::complex 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++) { diff --git a/src/math_functions.h b/src/math_functions.h index bea30cf916..aaba082615 100644 --- a/src/math_functions.h +++ b/src/math_functions.h @@ -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 diff --git a/tests/unit_tests/test_math.py b/tests/unit_tests/test_math.py index c8d33e6012..bfa021cebd 100644 --- a/tests/unit_tests/test_math.py +++ b/tests/unit_tests/test_math.py @@ -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():