diff --git a/openmc/lib/math.py b/openmc/lib/math.py index 029f1c4c9ff..8c62f241624 100644 --- a/openmc/lib/math.py +++ b/openmc/lib/math.py @@ -1,5 +1,4 @@ -from ctypes import c_int, c_double, POINTER, c_uint64 -from random import getrandbits +from ctypes import c_int, c_double import numpy as np from numpy.ctypeslib import ndpointer @@ -7,137 +6,18 @@ from . import _dll -_dll.t_percentile.restype = c_double -_dll.t_percentile.argtypes = [c_double, c_int] - -_dll.calc_pn_c.restype = None -_dll.calc_pn_c.argtypes = [c_int, c_double, ndpointer(c_double)] - -_dll.evaluate_legendre.restype = c_double -_dll.evaluate_legendre.argtypes = [c_int, POINTER(c_double), c_double] - -_dll.calc_rn_c.restype = None -_dll.calc_rn_c.argtypes = [c_int, ndpointer(c_double), ndpointer(c_double)] - _dll.calc_zn.restype = None _dll.calc_zn.argtypes = [c_int, c_double, c_double, ndpointer(c_double)] _dll.calc_zn_rad.restype = None _dll.calc_zn_rad.argtypes = [c_int, c_double, ndpointer(c_double)] -_dll.rotate_angle_c.restype = None -_dll.rotate_angle_c.argtypes = [ndpointer(c_double), c_double, - POINTER(c_double), POINTER(c_uint64)] -_dll.maxwell_spectrum.restype = c_double -_dll.maxwell_spectrum.argtypes = [c_double, POINTER(c_uint64)] - -_dll.watt_spectrum.restype = c_double -_dll.watt_spectrum.argtypes = [c_double, c_double, POINTER(c_uint64)] - -_dll.broaden_wmp_polynomials.restype = None -_dll.broaden_wmp_polynomials.argtypes = [c_double, c_double, c_int, - ndpointer(c_double)] - -_dll.normal_variate.restype = c_double -_dll.normal_variate.argtypes = [c_double, c_double, POINTER(c_uint64)] - -def t_percentile(p, df): - """ Calculate the percentile of the Student's t distribution with a - specified probability level and number of degrees of freedom - - Parameters - ---------- - p : float - Probability level - df : int - Degrees of freedom - - Returns - ------- - float - Corresponding t-value - - """ - - return _dll.t_percentile(p, df) - - -def calc_pn(n, x): - """ Calculate the n-th order Legendre polynomial at the value of x. - - Parameters - ---------- - n : int - Legendre order - x : float - Independent variable to evaluate the Legendre at - - Returns - ------- - float - Corresponding Legendre polynomial result - - """ - - pnx = np.empty(n + 1, dtype=np.float64) - _dll.calc_pn_c(n, x, pnx) - return pnx - - -def evaluate_legendre(data, x): - """ Finds the value of f(x) given a set of Legendre coefficients - and the value of x. - - Parameters - ---------- - data : iterable of float - Legendre coefficients - x : float - Independent variable to evaluate the Legendre at - - Returns - ------- - float - Corresponding Legendre expansion result - - """ - - data_arr = np.array(data, dtype=np.float64) - return _dll.evaluate_legendre(len(data)-1, - data_arr.ctypes.data_as(POINTER(c_double)), x) - - -def calc_rn(n, uvw): - """ Calculate the n-th order real Spherical Harmonics for a given angle; - all Rn,m values are provided for all n (where -n <= m <= n). - - Parameters - ---------- - n : int - Harmonics order - uvw : iterable of float - Independent variable to evaluate the Legendre at - - Returns - ------- - numpy.ndarray - Corresponding real harmonics value - - """ - - num_nm = (n + 1) * (n + 1) - rn = np.empty(num_nm, dtype=np.float64) - uvw_arr = np.array(uvw, dtype=np.float64) - _dll.calc_rn_c(n, uvw_arr, rn) - return rn - def calc_zn(n, rho, phi): """ Calculate the n-th order modified Zernike polynomial moment for a given angle (rho, theta) location in the unit disk. The normalization of the polynomials is such that the integral of Z_pq*Z_pq over the unit disk is exactly pi - Parameters ---------- n : int @@ -146,12 +26,10 @@ def calc_zn(n, rho, phi): Radial location in the unit disk phi : float Theta (radians) location in the unit disk - Returns ------- numpy.ndarray Corresponding resulting list of coefficients - """ num_bins = ((n + 1) * (n + 2)) // 2 @@ -165,161 +43,19 @@ def calc_zn_rad(n, rho): moment with no azimuthal dependency (m=0) for a given radial location in the unit disk. The normalization of the polynomials is such that the integral of Z_pq*Z_pq over the unit disk is exactly pi. - Parameters ---------- n : int Maximum order rho : float Radial location in the unit disk - Returns ------- numpy.ndarray Corresponding resulting list of coefficients - """ num_bins = n // 2 + 1 zn_rad = np.zeros(num_bins, dtype=np.float64) _dll.calc_zn_rad(n, rho, zn_rad) return zn_rad - - -def rotate_angle(uvw0, mu, phi, prn_seed=None): - """ Rotates direction cosines through a polar angle whose cosine is - mu and through an azimuthal angle sampled uniformly. - - Parameters - ---------- - uvw0 : iterable of float - Original direction cosine - mu : float - Polar angle cosine to rotate - phi : float - Azimuthal angle; if None, one will be sampled uniformly - prn_seed : int - Pseudorandom number generator (PRNG) seed; if None, one will be - generated randomly. - - Returns - ------- - numpy.ndarray - Rotated direction cosine - - """ - - if prn_seed is None: - prn_seed = getrandbits(63) - - uvw0_arr = np.array(uvw0, dtype=np.float64) - if phi is None: - _dll.rotate_angle_c(uvw0_arr, mu, None, c_uint64(prn_seed)) - else: - _dll.rotate_angle_c(uvw0_arr, mu, c_double(phi), c_uint64(prn_seed)) - - uvw = uvw0_arr - - return uvw - - -def maxwell_spectrum(T, prn_seed=None): - """ Samples an energy from the Maxwell fission distribution based - on a direct sampling scheme. - - Parameters - ---------- - T : float - Spectrum parameter - prn_seed : int - Pseudorandom number generator (PRNG) seed; if None, one will be - generated randomly. - - Returns - ------- - float - Sampled outgoing energy - - """ - - if prn_seed is None: - prn_seed = getrandbits(63) - - return _dll.maxwell_spectrum(T, c_uint64(prn_seed)) - - -def watt_spectrum(a, b, prn_seed=None): - """ Samples an energy from the Watt energy-dependent fission spectrum. - - Parameters - ---------- - a : float - Spectrum parameter a - b : float - Spectrum parameter b - prn_seed : int - Pseudorandom number generator (PRNG) seed; if None, one will be - generated randomly. - - Returns - ------- - float - Sampled outgoing energy - - """ - - if prn_seed is None: - prn_seed = getrandbits(63) - - return _dll.watt_spectrum(a, b, c_uint64(prn_seed)) - - -def normal_variate(mean_value, std_dev, prn_seed=None): - """ Samples an energy from the Normal distribution. - - Parameters - ---------- - mean_value : float - Mean of the Normal distribution - std_dev : float - Standard deviation of the normal distribution - prn_seed : int - Pseudorandom number generator (PRNG) seed; if None, one will be - generated randomly. - - Returns - ------- - float - Sampled outgoing normally distributed value - - """ - - if prn_seed is None: - prn_seed = getrandbits(63) - - return _dll.normal_variate(mean_value, std_dev, c_uint64(prn_seed)) - - -def broaden_wmp_polynomials(E, dopp, n): - """ Doppler broadens the windowed multipole curvefit. The curvefit is a - polynomial of the form a/E + b/sqrt(E) + c + d sqrt(E) ... - - Parameters - ---------- - E : float - Energy to evaluate at - dopp : float - sqrt(atomic weight ratio / kT), with kT given in eV - n : int - Number of components to the polynomial - - Returns - ------- - numpy.ndarray - Resultant leading coefficients - - """ - - factors = np.zeros(n, dtype=np.float64) - _dll.broaden_wmp_polynomials(E, dopp, n, factors) - return factors diff --git a/tests/cpp_unit_tests/CMakeLists.txt b/tests/cpp_unit_tests/CMakeLists.txt index 24ec1b7a90e..f0f5f2853ad 100644 --- a/tests/cpp_unit_tests/CMakeLists.txt +++ b/tests/cpp_unit_tests/CMakeLists.txt @@ -3,6 +3,7 @@ set(TEST_NAMES test_file_utils test_tally test_interpolate + test_math # Add additional unit test files here ) diff --git a/tests/cpp_unit_tests/test_math.cpp b/tests/cpp_unit_tests/test_math.cpp new file mode 100644 index 00000000000..1ad7c4b7097 --- /dev/null +++ b/tests/cpp_unit_tests/test_math.cpp @@ -0,0 +1,357 @@ +#include +#include + +#include +#include +#include + +#include "openmc/math_functions.h" +#include "openmc/random_dist.h" +#include "openmc/random_lcg.h" +#include "openmc/wmp.h" + +TEST_CASE("Test t_percentile") +{ + // The reference solutions come from scipy.stats.t.ppf + std::vector> ref_ts { + {-15.894544844102773, -0.32491969623407446, 0.000000000000000, + 0.32491969623407446, 15.894544844102759}, + {-4.848732214442601, -0.2886751346880066, 0.000000000000000, + 0.2886751346880066, 4.848732214442598}, + {-2.756508521909475, -0.2671808657039658, 0.000000000000000, + 0.2671808657039658, 2.7565085219094745}}; + + // Permutations include 1 DoF, 2 DoF, and > 2 DoF + // We will test 5 p-values at 3-DoF values + std::vector test_ps {0.02, 0.4, 0.5, 0.6, 0.98}; + std::vector test_dfs {1, 2, 5}; + + for (int i = 0; i < test_dfs.size(); i++) { + int df = test_dfs[i]; + + std::vector test_ts; + + for (double p : test_ps) { + double test_t = openmc::t_percentile(p, df); + test_ts.push_back(test_t); + } + + // The 5 DoF approximation in openmc.lib.math.t_percentile is off by up to + // 8e-3 from the scipy solution, so test that one separately with looser + // tolerance + double tolerance = (df > 2) ? 1e-2 : 1e-6; + + REQUIRE_THAT( + ref_ts[i], Catch::Matchers::Approx(test_ts).epsilon(tolerance)); + } +} + +TEST_CASE("Test calc_pn") +{ + // The reference solutions come from scipy.special.eval_legendre + std::vector> ref_vals { + {1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1}, + {1, -0.5, -0.125, 0.4375, -0.289062, -0.0898438, 0.323242, -0.223145, + -0.0736389, 0.267899, -0.188229}, + {1, 0, -0.5, -0, 0.375, 0, -0.3125, -0, 0.273438, 0, -0.246094}, + {1, 0.5, -0.125, -0.4375, -0.289062, 0.0898438, 0.323242, 0.223145, + -0.0736389, -0.267899, -0.188229}, + {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}}; + + int max_order = 10; + std::vector test_xs = {-1.0, -0.5, 0.0, 0.5, 1.0}; + + std::vector> test_vals; + for (double x : test_xs) { + std::vector test_val(max_order + 1); + openmc::calc_pn_c(max_order, x, test_val.data()); + test_vals.push_back(test_val); + } + + for (int i = 0; i < ref_vals.size(); i++) { + REQUIRE_THAT(ref_vals[i], Catch::Matchers::Approx(test_vals[i])); + } +} + +TEST_CASE("Test evaluate_legendre") +{ + // The reference solutions come from numpy.polynomial.legendre.legval + std::vector ref_vals { + 5.5, -0.45597649, -1.35351562, -2.7730999, 60.5}; + + int max_order = 10; + std::vector test_xs = {-1.0, -0.5, 0.0, 0.5, 1.0}; + + // Set the coefficients back to 1s for the test values since + // evaluate legendre incorporates the (2l+1)/2 term on its own + std::vector test_coeffs(max_order + 1, 1.0); + + std::vector test_vals; + for (double x : test_xs) { + test_vals.push_back( + openmc::evaluate_legendre(test_coeffs.size() - 1, test_coeffs.data(), x)); + } + + REQUIRE_THAT(ref_vals, Catch::Matchers::Approx(test_vals)); +} + +TEST_CASE("Test calc_rn") +{ + std::vector ref_vals {1.000000000000000, -0.019833838076210, + 0.980066577841242, -0.197676811654084, 0.006790834062088, + -0.033668438114859, 0.940795745502164, -0.335561350977312, + 0.033500236162691, -0.001831975566765, 0.014882082223994, + -0.046185860057145, 0.883359726009014, -0.460318044571973, + 0.073415616482180, -0.005922278973373, 0.000448625292461, + -0.004750335422039, 0.025089695062177, -0.057224052171859, + 0.809468042300133, -0.570331780454957, 0.123771351522967, + -0.015356543011155, 0.001061098599927, -0.000104097571795, + 0.001319047965347, -0.009263463267120, 0.037043163155191, + -0.066518621473934, 0.721310852552881, -0.662967447756079, + 0.182739660926192, -0.029946258412359, 0.003119841820746, + -0.000190549327031, 0.000023320052630, -0.000338370521658, + 0.002878809439524, -0.015562587450914, 0.050271226423217, + -0.073829294593737, 0.621486505922182, -0.735830327235834, + 0.247995745731425, -0.050309614442385, 0.006809024629381, + -0.000619383085285, 0.000034086826414, -0.000005093626712, + 0.000082405610567, -0.000809532012556, 0.005358034016708, + -0.023740240859138, 0.064242405926477, -0.078969918083157, + 0.512915839160049, -0.787065093668736, 0.316917738015632, + -0.076745744765114, 0.012672942183651, -0.001481838409317, + 0.000120451946983, -0.000006047366709, 0.000001091052697, + -0.000019334294214, 0.000213051604838, -0.001640234119608, + 0.008982263900105, -0.033788039035668, 0.078388909900756, + -0.081820779415058, 0.398746190829636, -0.815478614863816, + 0.386704633068855, -0.109227544713261, 0.021245051959237, + -0.003002428416676, 0.000311416667310, -0.000022954482885, + 0.000001059646310, -0.000000230023931, 0.000004408854505, + -0.000053457925526, 0.000464152759861, -0.002976305522860, + 0.013958017448970, -0.045594791382625, 0.092128969315914, + -0.082334538374971, 0.282248459574595, -0.820599067736528, + 0.454486474163594, -0.147395565311743, 0.033013815809602, + -0.005448090715661, 0.000678450207914, -0.000063467485444, + 0.000004281943868, -0.000000182535754, 0.000000047847775, + -0.000000982664801, 0.000012933320414, -0.000124076425457, + 0.000901739739837, -0.004982323311961, 0.020457776068931, + -0.058948376674391, 0.104888993733747, -0.080538298991650, + 0.166710763818175, -0.802696588503912, 0.517433650833039, + -0.190564076304612, 0.048387190622376, -0.009120081648146, + 0.001318069323039, -0.000147308722683, 0.000012561029621, + -0.000000779794781, 0.000000030722703}; + + int max_order = 10; + + double azi = 0.1; // Longitude + double pol = 0.2; // Latitude + double mu = std::cos(pol); + + std::vector test_uvw {std::sin(pol) * std::cos(azi), + std::sin(pol) * std::sin(azi), std::cos(pol)}; + + std::vector test_vals((max_order + 1) * (max_order + 1), 0); + openmc::calc_rn_c(max_order, test_uvw.data(), test_vals.data()); + + REQUIRE_THAT(ref_vals, Catch::Matchers::Approx(test_vals)); +} + +TEST_CASE("Test calc_zn") +{ + std::vector ref_vals {1.00000000e+00, 2.39712769e-01, 4.38791281e-01, + 2.10367746e-01, -5.00000000e-01, 1.35075576e-01, 1.24686873e-01, + -2.99640962e-01, -5.48489101e-01, 8.84215021e-03, 5.68310892e-02, + -4.20735492e-01, -1.25000000e-01, -2.70151153e-01, -2.60091773e-02, + 1.87022545e-02, -3.42888902e-01, 1.49820481e-01, 2.74244551e-01, + -2.43159131e-02, -2.50357380e-02, 2.20500013e-03, -1.98908812e-01, + 4.07587508e-01, 4.37500000e-01, 2.61708929e-01, 9.10321205e-02, + -1.54686328e-02, -2.74049397e-03, -7.94845816e-02, 4.75368705e-01, + 7.11647284e-02, 1.30266162e-01, 3.37106977e-02, 1.06401886e-01, + -7.31606787e-03, -2.95625975e-03, -1.10250006e-02, 3.55194307e-01, + -1.44627826e-01, -2.89062500e-01, -9.28644588e-02, -1.62557358e-01, + 7.73431638e-02, -2.55329539e-03, -1.90923851e-03, 1.57578403e-02, + 1.72995854e-01, -3.66267690e-01, -1.81657333e-01, -3.32521518e-01, + -2.59738162e-02, -2.31580576e-01, 4.20673902e-02, -4.11710546e-04, + -9.36449487e-04, 1.92156884e-02, 2.82515641e-02, -3.90713738e-01, + -1.69280296e-01, -8.98437500e-02, -1.08693628e-01, 1.78813094e-01, + -1.98191857e-01, 1.65964201e-02, 2.77013853e-04}; + + int n = 10; + double rho = 0.5; + double phi = 0.5; + + int nums = ((n + 1) * (n + 2)) / 2; + + std::vector test_vals(nums, 0); + openmc::calc_zn(n, rho, phi, test_vals.data()); + + REQUIRE_THAT(ref_vals, Catch::Matchers::Approx(test_vals)); +} + +TEST_CASE("Test calc_zn_rad") +{ + std::vector ref_vals {1.00000000e+00, -5.00000000e-01, + -1.25000000e-01, 4.37500000e-01, -2.89062500e-01, -8.98437500e-02}; + + int n = 10; + double rho = 0.5; + + int nums = n / 2 + 1; + std::vector test_vals(nums, 0); + openmc::calc_zn_rad(n, rho, test_vals.data()); + + REQUIRE_THAT(ref_vals, Catch::Matchers::Approx(test_vals)); +} + +TEST_CASE("Test rotate_angle") +{ + std::vector uvw0 {1.0, 0.0, 0.0}; + double phi = 0.0; + + uint64_t prn_seed = 1; + openmc::prn(&prn_seed); + + SECTION("Test rotate_angle mu is 0") + { + std::vector ref_uvw {0.0, 0.0, -1.0}; + + double mu = 0.0; + + std::vector test_uvw(uvw0); + openmc::rotate_angle_c(test_uvw.data(), mu, &phi, &prn_seed); + + REQUIRE_THAT(ref_uvw, Catch::Matchers::Approx(test_uvw)); + } + + SECTION("Test rotate_angle mu is 1") + { + std::vector ref_uvw = {1.0, 0.0, 0.0}; + + double mu = 1.0; + + std::vector test_uvw(uvw0); + openmc::rotate_angle_c(test_uvw.data(), mu, &phi, &prn_seed); + + REQUIRE_THAT(ref_uvw, Catch::Matchers::Approx(test_uvw)); + } + + // Now to test phi is None + SECTION("Test rotate_angle no phi") + { + // When seed = 1, phi will be sampled as 1.9116495709698769 + // The resultant reference is from hand-calculations given the above + std::vector ref_uvw = { + 0.9, -0.422746750548505, 0.10623175090659095}; + + double mu = 0.9; + prn_seed = 1; + + std::vector test_uvw(uvw0); + openmc::rotate_angle_c(test_uvw.data(), mu, NULL, &prn_seed); + + REQUIRE_THAT(ref_uvw, Catch::Matchers::Approx(test_uvw)); + } +} + +TEST_CASE("Test maxwell_spectrum") +{ + double ref_val = 0.27767406743161277; + + double T = 0.5; + uint64_t prn_seed = 1; + + double test_val = openmc::maxwell_spectrum(T, &prn_seed); + + REQUIRE(ref_val == test_val); +} + +TEST_CASE("Test watt_spectrum") +{ + double ref_val = 0.30957476387766697; + + double a = 0.5; + double b = 0.75; + uint64_t prn_seed = 1; + + double test_val = openmc::watt_spectrum(a, b, &prn_seed); + + REQUIRE(ref_val == test_val); +} + +TEST_CASE("Test normal_variate") +{ + + // Generate a series of normally distributed random numbers and test + // whether their mean and standard deviation are close to the expected value + SECTION("Test with non-zero standard deviation") + { + uint64_t seed = 1; + + double mean = 0.0; + double standard_deviation = 1.0; + + int num_samples = 10000; + double sum = 0.0; + double sum_squared_difference = 0.0; + + for (int i = 0; i < num_samples; ++i) { + double sample = openmc::normal_variate(mean, standard_deviation, &seed); + sum += sample; + sum_squared_difference += (sample - mean) * (sample - mean); + } + + double actual_mean = sum / num_samples; + double actual_standard_deviation = + std::sqrt(sum_squared_difference / num_samples); + + REQUIRE_THAT(mean, Catch::Matchers::WithinAbs(actual_mean, 0.1)); + REQUIRE_THAT(standard_deviation, + Catch::Matchers::WithinAbs(actual_standard_deviation, 0.1)); + } + + // When the standard deviation is zero + // the generated random number should always be equal to the mean + SECTION("Test with zero standard deviation") + { + uint64_t seed = 1; + double mean = 5.0; + double standard_deviation = 0.0; + + for (int i = 0; i < 10; ++i) { + double sample = openmc::normal_variate(mean, standard_deviation, &seed); + REQUIRE(sample == mean); + } + } +} + +TEST_CASE("Test broaden_wmp_polynomials") +{ + double test_E = 0.5; + int n = 6; + + // Two branches of the code to worry about, beta > 6 and otherwise + // beta = sqrtE * dopp + SECTION("Test broaden_wmp_polynomials beta > 6") + { + std::vector ref_val { + 2., 1.41421356, 1.0001, 0.70731891, 0.50030001, 0.353907}; + + double test_dopp = 100.0; // approximately U235 at room temperature + + std::vector test_val(n, 0); + openmc::broaden_wmp_polynomials(test_E, test_dopp, n, test_val.data()); + + REQUIRE_THAT(ref_val, Catch::Matchers::Approx(test_val)); + } + + SECTION("Test broaden_wmp_polynomials beta < 6") + { + std::vector ref_val = { + 1.99999885, 1.41421356, 1.04, 0.79195959, 0.6224, 0.50346003}; + + double test_dopp = 5.0; + + std::vector test_val(n, 0); + openmc::broaden_wmp_polynomials(test_E, test_dopp, n, test_val.data()); + + REQUIRE_THAT(ref_val, Catch::Matchers::Approx(test_val)); + } +} diff --git a/tests/unit_tests/test_math.py b/tests/unit_tests/test_math.py deleted file mode 100644 index c50057fcd6f..00000000000 --- a/tests/unit_tests/test_math.py +++ /dev/null @@ -1,247 +0,0 @@ -import numpy as np -import pytest -import scipy as sp -from scipy.stats import shapiro - -import openmc -import openmc.lib - - -def test_t_percentile(): - # Permutations include 1 DoF, 2 DoF, and > 2 DoF - # We will test 5 p-values at 3-DoF values - test_ps = [0.02, 0.4, 0.5, 0.6, 0.98] - test_dfs = [1, 2, 5] - - # The reference solutions come from Scipy - ref_ts = [[sp.stats.t.ppf(p, df) for p in test_ps] for df in test_dfs] - - test_ts = [[openmc.lib.math.t_percentile(p, df) for p in test_ps] - for df in test_dfs] - - # The 5 DoF approximation in openmc.lib.math.t_percentile is off by up to - # 8e-3 from the scipy solution, so test that one separately with looser - # tolerance - assert np.allclose(ref_ts[:-1], test_ts[:-1]) - assert np.allclose(ref_ts[-1], test_ts[-1], atol=1e-2) - - -def test_calc_pn(): - max_order = 10 - test_xs = np.linspace(-1., 1., num=5, endpoint=True) - - # Reference solutions from scipy - ref_vals = np.array([sp.special.eval_legendre(n, test_xs) - for n in range(0, max_order + 1)]) - - test_vals = [] - for x in test_xs: - test_vals.append(openmc.lib.math.calc_pn(max_order, x).tolist()) - - test_vals = np.swapaxes(np.array(test_vals), 0, 1) - - assert np.allclose(ref_vals, test_vals) - - -def test_evaluate_legendre(): - max_order = 10 - # Coefficients are set to 1, but will incorporate the (2l+1)/2 norm factor - # for the reference solution - test_coeffs = [0.5 * (2. * l + 1.) for l in range(max_order + 1)] - test_xs = np.linspace(-1., 1., num=5, endpoint=True) - - ref_vals = np.polynomial.legendre.legval(test_xs, test_coeffs) - - # Set the coefficients back to 1s for the test values since - # evaluate legendre incorporates the (2l+1)/2 term on its own - test_coeffs = [1. for l in range(max_order + 1)] - - test_vals = np.array([openmc.lib.math.evaluate_legendre(test_coeffs, x) - for x in test_xs]) - - assert np.allclose(ref_vals, test_vals) - - -def test_calc_rn(): - max_order = 10 - test_ns = np.array([i for i in range(0, max_order + 1)]) - azi = 0.1 # Longitude - pol = 0.2 # Latitude - test_uvw = np.array([np.sin(pol) * np.cos(azi), - np.sin(pol) * np.sin(azi), - np.cos(pol)]) - - # Reference solutions from the equations - ref_vals = [] - - def coeff(n, m): - return np.sqrt((2. * n + 1) * sp.special.factorial(n - m) / - (sp.special.factorial(n + m))) - - def pnm_bar(n, m, mu): - val = coeff(n, m) - if m != 0: - val *= np.sqrt(2.) - val *= sp.special.lpmv([m], [n], [mu]) - return val[0] - - ref_vals = [] - for n in test_ns: - for m in range(-n, n + 1): - if m < 0: - ylm = pnm_bar(n, np.abs(m), np.cos(pol)) * \ - np.sin(np.abs(m) * azi) - else: - ylm = pnm_bar(n, m, np.cos(pol)) * np.cos(m * azi) - - # Un-normalize for comparison - ylm /= np.sqrt(2. * n + 1.) - ref_vals.append(ylm) - - test_vals = [] - test_vals = openmc.lib.math.calc_rn(max_order, test_uvw) - - assert np.allclose(ref_vals, test_vals) - - -def test_calc_zn(): - n = 10 - rho = 0.5 - phi = 0.5 - - # Reference solution from running the C++ implementation - ref_vals = np.array([ - 1.00000000e+00, 2.39712769e-01, 4.38791281e-01, - 2.10367746e-01, -5.00000000e-01, 1.35075576e-01, - 1.24686873e-01, -2.99640962e-01, -5.48489101e-01, - 8.84215021e-03, 5.68310892e-02, -4.20735492e-01, - -1.25000000e-01, -2.70151153e-01, -2.60091773e-02, - 1.87022545e-02, -3.42888902e-01, 1.49820481e-01, - 2.74244551e-01, -2.43159131e-02, -2.50357380e-02, - 2.20500013e-03, -1.98908812e-01, 4.07587508e-01, - 4.37500000e-01, 2.61708929e-01, 9.10321205e-02, - -1.54686328e-02, -2.74049397e-03, -7.94845816e-02, - 4.75368705e-01, 7.11647284e-02, 1.30266162e-01, - 3.37106977e-02, 1.06401886e-01, -7.31606787e-03, - -2.95625975e-03, -1.10250006e-02, 3.55194307e-01, - -1.44627826e-01, -2.89062500e-01, -9.28644588e-02, - -1.62557358e-01, 7.73431638e-02, -2.55329539e-03, - -1.90923851e-03, 1.57578403e-02, 1.72995854e-01, - -3.66267690e-01, -1.81657333e-01, -3.32521518e-01, - -2.59738162e-02, -2.31580576e-01, 4.20673902e-02, - -4.11710546e-04, -9.36449487e-04, 1.92156884e-02, - 2.82515641e-02, -3.90713738e-01, -1.69280296e-01, - -8.98437500e-02, -1.08693628e-01, 1.78813094e-01, - -1.98191857e-01, 1.65964201e-02, 2.77013853e-04]) - - test_vals = openmc.lib.math.calc_zn(n, rho, phi) - - assert np.allclose(ref_vals, test_vals) - - -def test_calc_zn_rad(): - n = 10 - rho = 0.5 - - # Reference solution from running the C++ implementation - ref_vals = np.array([ - 1.00000000e+00, -5.00000000e-01, -1.25000000e-01, - 4.37500000e-01, -2.89062500e-01,-8.98437500e-02]) - - test_vals = openmc.lib.math.calc_zn_rad(n, rho) - - assert np.allclose(ref_vals, test_vals) - - -def test_rotate_angle(): - uvw0 = np.array([1., 0., 0.]) - phi = 0. - mu = 0. - - # reference: mu of 0 pulls the vector the bottom, so: - ref_uvw = np.array([0., 0., -1.]) - - test_uvw = openmc.lib.math.rotate_angle(uvw0, mu, phi) - - assert np.array_equal(ref_uvw, test_uvw) - - # Repeat for mu = 1 (no change) - mu = 1. - ref_uvw = np.array([1., 0., 0.]) - - test_uvw = openmc.lib.math.rotate_angle(uvw0, mu, phi) - - assert np.array_equal(ref_uvw, test_uvw) - - # Now to test phi is None - mu = 0.9 - phi = None - prn_seed = 1 - - # When seed = 1, phi will be sampled as 1.9116495709698769 - # The resultant reference is from hand-calculations given the above - ref_uvw = [0.9, -0.422746750548505, 0.10623175090659095] - test_uvw = openmc.lib.math.rotate_angle(uvw0, mu, phi, prn_seed) - - assert np.allclose(ref_uvw, test_uvw) - - -def test_maxwell_spectrum(): - prn_seed = 1 - T = 0.5 - ref_val = 0.27767406743161277 - test_val = openmc.lib.math.maxwell_spectrum(T, prn_seed) - - assert ref_val == test_val - - -def test_watt_spectrum(): - prn_seed = 1 - a = 0.5 - b = 0.75 - ref_val = 0.30957476387766697 - test_val = openmc.lib.math.watt_spectrum(a, b, prn_seed) - - assert ref_val == test_val - - -def test_normal_dist(): - # When standard deviation is zero, sampled value should be mean - prn_seed = 1 - mean = 14.08 - stdev = 0.0 - ref_val = 14.08 - test_val = openmc.lib.math.normal_variate(mean, stdev, prn_seed) - assert ref_val == pytest.approx(test_val) - - # Use Shapiro-Wilk test to ensure normality of sampled vairates - stdev = 1.0 - samples = [] - num_samples = 10000 - for _ in range(num_samples): - # sample the normal distribution from openmc - samples.append(openmc.lib.math.normal_variate(mean, stdev, prn_seed)) - prn_seed += 1 - stat, p = shapiro(samples) - assert p > 0.05 - - -def test_broaden_wmp_polynomials(): - # Two branches of the code to worry about, beta > 6 and otherwise - # beta = sqrtE * dopp - # First lets do beta > 6 - test_E = 0.5 - test_dopp = 100. # approximately U235 at room temperature - n = 6 - - ref_val = [2., 1.41421356, 1.0001, 0.70731891, 0.50030001, 0.353907] - test_val = openmc.lib.math.broaden_wmp_polynomials(test_E, test_dopp, n) - - assert np.allclose(ref_val, test_val) - - # now beta < 6 - test_dopp = 5. - ref_val = [1.99999885, 1.41421356, 1.04, 0.79195959, 0.6224, 0.50346003] - test_val = openmc.lib.math.broaden_wmp_polynomials(test_E, test_dopp, n) - - assert np.allclose(ref_val, test_val)