Skip to content

Commit

Permalink
update math function unit test with catch2 (#2955)
Browse files Browse the repository at this point in the history 
Co-authored-by: Paul Romano <[email protected]>
  • Loading branch information
@lhchg @paulromano
lhchg and paulromano authored Jun 19, 2024
1 parent ddc9526 commit 3bedd04
Show file tree
Hide file tree
Showing 4 changed files with 359 additions and 512 deletions.
266 changes: 1 addition & 265 deletions openmc/lib/math.py
View file
Original file line number Diff line number Diff line change
@@ -1,143 +1,23 @@
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

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
Expand All @@ -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
Expand All @@ -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
1 change: 1 addition & 0 deletions tests/cpp_unit_tests/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@ set(TEST_NAMES
test_file_utils
test_tally
test_interpolate
test_math
# Add additional unit test files here
)

Expand Down
Loading

0 comments on commit 3bedd04

Please sign in to comment.