diff --git a/openmc/capi/math.py b/openmc/capi/math.py index 3b0b723b32..d1d7abdf64 100644 --- a/openmc/capi/math.py +++ b/openmc/capi/math.py @@ -6,37 +6,33 @@ from numpy.ctypeslib import ndpointer from . import _dll -_dll.t_percentile.restype = c_double -_dll.t_percentile.argtypes = [POINTER(c_double), POINTER(c_int)] +_dll.t_percentile_c.restype = c_double +_dll.t_percentile_c.argtypes = [c_double, c_int] -_dll.calc_pn.restype = None -_dll.calc_pn.argtypes = [POINTER(c_int), POINTER(c_double), - ndpointer(c_double)] +_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 = [POINTER(c_int), POINTER(c_double), - POINTER(c_double)] +_dll.evaluate_legendre_c.restype = c_double +_dll.evaluate_legendre_c.argtypes = [c_int, POINTER(c_double), c_double] -_dll.calc_rn.restype = None -_dll.calc_rn.argtypes = [POINTER(c_int), ndpointer(c_double), - ndpointer(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 = [POINTER(c_int), POINTER(c_double), POINTER(c_double), - ndpointer(c_double)] +_dll.calc_zn_c.restype = None +_dll.calc_zn_c.argtypes = [c_int, c_double, c_double, ndpointer(c_double)] -_dll.rotate_angle.restype = None -_dll.rotate_angle.argtypes = [ndpointer(c_double), POINTER(c_double), - ndpointer(c_double), POINTER(c_double)] -_dll.maxwell_spectrum.restype = c_double -_dll.maxwell_spectrum.argtypes = [POINTER(c_double)] +_dll.rotate_angle_c.restype = None +_dll.rotate_angle_c.argtypes = [ndpointer(c_double), c_double, + POINTER(c_double)] +_dll.maxwell_spectrum_c.restype = c_double +_dll.maxwell_spectrum_c.argtypes = [c_double] -_dll.watt_spectrum.restype = c_double -_dll.watt_spectrum.argtypes = [POINTER(c_double), POINTER(c_double)] +_dll.watt_spectrum_c.restype = c_double +_dll.watt_spectrum_c.argtypes = [c_double, c_double] -_dll.broaden_wmp_polynomials.restype = None -_dll.broaden_wmp_polynomials.argtypes = [POINTER(c_double), POINTER(c_double), - POINTER(c_int), ndpointer(c_double)] +_dll.broaden_wmp_polynomials_c.restype = None +_dll.broaden_wmp_polynomials_c.argtypes = [c_double, c_double, c_int, + ndpointer(c_double)] def t_percentile(p, df): @@ -57,7 +53,7 @@ def t_percentile(p, df): """ - return _dll.t_percentile(c_double(p), c_int(df)) + return _dll.t_percentile_c(p, df) def calc_pn(n, x): @@ -78,7 +74,7 @@ def calc_pn(n, x): """ pnx = np.empty(n + 1, dtype=np.float64) - _dll.calc_pn(c_int(n), c_double(x), pnx) + _dll.calc_pn_c(n, x, pnx) return pnx @@ -101,9 +97,9 @@ def evaluate_legendre(data, x): """ data_arr = np.array(data, dtype=np.float64) - return _dll.evaluate_legendre(c_int(len(data)), - data_arr.ctypes.data_as(POINTER(c_double)), - c_double(x)) + return _dll.evaluate_legendre_c(len(data), + data_arr.ctypes.data_as(POINTER(c_double)), + x) def calc_rn(n, uvw): @@ -127,7 +123,7 @@ def calc_rn(n, uvw): 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_int(n), uvw_arr, rn) + _dll.calc_rn_c(n, uvw_arr, rn) return rn @@ -155,7 +151,7 @@ def calc_zn(n, rho, phi): num_bins = ((n + 1) * (n + 2)) // 2 zn = np.zeros(num_bins, dtype=np.float64) - _dll.calc_zn(c_int(n), c_double(rho), c_double(phi), zn) + _dll.calc_zn_c(n, rho, phi, zn) return zn @@ -179,13 +175,13 @@ def rotate_angle(uvw0, mu, phi=None): """ - uvw = np.zeros(3, dtype=np.float64) uvw0_arr = np.array(uvw0, dtype=np.float64) if phi is None: - _dll.rotate_angle(uvw0_arr, c_double(mu), uvw, None) + _dll.rotate_angle_c(uvw0_arr, mu, None) else: - _dll.rotate_angle(uvw0_arr, c_double(mu), uvw, c_double(phi)) + _dll.rotate_angle_c(uvw0_arr, mu, c_double(phi)) + uvw = uvw0_arr return uvw @@ -206,7 +202,7 @@ def maxwell_spectrum(T): """ - return _dll.maxwell_spectrum(c_double(T)) + return _dll.maxwell_spectrum_c(T) def watt_spectrum(a, b): @@ -226,7 +222,7 @@ def watt_spectrum(a, b): """ - return _dll.watt_spectrum(c_double(a), c_double(b)) + return _dll.watt_spectrum_c(a, b) def broaden_wmp_polynomials(E, dopp, n): @@ -250,6 +246,5 @@ def broaden_wmp_polynomials(E, dopp, n): """ factors = np.zeros(n, dtype=np.float64) - _dll.broaden_wmp_polynomials(c_double(E), c_double(dopp), c_int(n), - factors) + _dll.broaden_wmp_polynomials_c(E, dopp, n, factors) return factors diff --git a/src/api.F90 b/src/api.F90 index a06e2f4cc4..da50007d68 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -95,17 +95,6 @@ module openmc_api public :: openmc_tally_set_nuclides public :: openmc_tally_set_scores public :: openmc_tally_set_type - public :: t_percentile - public :: calc_pn - public :: calc_rn - public :: calc_zn - public :: evaluate_legendre - public :: rotate_angle - public :: maxwell_spectrum - public :: watt_spectrum - public :: faddeeva - public :: w_derivative - public :: broaden_wmp_polynomials contains diff --git a/src/distribution_multivariate.F90 b/src/distribution_multivariate.F90 index 4d7d42cc3d..19db1c297d 100644 --- a/src/distribution_multivariate.F90 +++ b/src/distribution_multivariate.F90 @@ -119,7 +119,7 @@ contains else ! Sample azimuthal angle phi = this % phi % sample() - call rotate_angle(this % reference_uvw, mu, uvw, phi) + uvw = rotate_angle(this % reference_uvw, mu, phi) end if end function polar_azimuthal_sample diff --git a/src/math.F90 b/src/math.F90 index 05942a86e3..33f1ab4d30 100644 --- a/src/math.F90 +++ b/src/math.F90 @@ -26,22 +26,22 @@ module math interface - pure function t_percentile_c_intfc(p, df) bind(C, name='t_percentile_c') & + pure function t_percentile(p, df) bind(C, name='t_percentile_c') & result(t) use ISO_C_BINDING implicit none real(C_DOUBLE), value, intent(in) :: p integer(C_INT), value, intent(in) :: df real(C_DOUBLE) :: t - end function t_percentile_c_intfc + end function t_percentile - pure subroutine calc_pn_c_intfc(n, x, pnx) bind(C, name='calc_pn_c') + pure subroutine calc_pn(n, x, pnx) bind(C, name='calc_pn_c') use ISO_C_BINDING implicit none integer(C_INT), value, intent(in) :: n real(C_DOUBLE), value, intent(in) :: x real(C_DOUBLE), intent(out) :: pnx(n + 1) - end subroutine calc_pn_c_intfc + end subroutine calc_pn pure function evaluate_legendre_c_intfc(n, data, x) & bind(C, name='evaluate_legendre_c') result(val) @@ -53,22 +53,22 @@ module math real(C_DOUBLE) :: val end function evaluate_legendre_c_intfc - pure subroutine calc_rn_c_intfc(n, uvw, rn) bind(C, name='calc_rn_c') + pure subroutine calc_rn(n, uvw, rn) bind(C, name='calc_rn_c') use ISO_C_BINDING implicit none integer(C_INT), value, intent(in) :: n real(C_DOUBLE), intent(in) :: uvw(3) real(C_DOUBLE), intent(out) :: rn(2 * n + 1) - end subroutine calc_rn_c_intfc + end subroutine calc_rn - pure subroutine calc_zn_c_intfc(n, rho, phi, zn) bind(C, name='calc_zn_c') + pure subroutine calc_zn(n, rho, phi, zn) bind(C, name='calc_zn_c') use ISO_C_BINDING implicit none integer(C_INT), value, intent(in) :: n real(C_DOUBLE), value, intent(in) :: rho real(C_DOUBLE), value, intent(in) :: phi real(C_DOUBLE), intent(out) :: zn(((n + 1) * (n + 2)) / 2) - end subroutine calc_zn_c_intfc + end subroutine calc_zn subroutine rotate_angle_c_intfc(uvw, mu, phi) bind(C, name='rotate_angle_c') use ISO_C_BINDING @@ -78,24 +78,24 @@ module math real(C_DOUBLE), optional, intent(in) :: phi end subroutine rotate_angle_c_intfc - function maxwell_spectrum_c_intfc(T) bind(C, name='maxwell_spectrum_c') & + function maxwell_spectrum(T) bind(C, name='maxwell_spectrum_c') & result(E_out) use ISO_C_BINDING implicit none real(C_DOUBLE), value, intent(in) :: T real(C_DOUBLE) :: E_out - end function maxwell_spectrum_c_intfc + end function maxwell_spectrum - function watt_spectrum_c_intfc(a, b) bind(C, name='watt_spectrum_c') & + function watt_spectrum(a, b) bind(C, name='watt_spectrum_c') & result(E_out) use ISO_C_BINDING implicit none real(C_DOUBLE), value, intent(in) :: a real(C_DOUBLE), value, intent(in) :: b real(C_DOUBLE) :: E_out - end function watt_spectrum_c_intfc + end function watt_spectrum - subroutine broaden_wmp_polynomials_c_intfc(E, dopp, n, factors) & + subroutine broaden_wmp_polynomials(E, dopp, n, factors) & bind(C, name='broaden_wmp_polynomials_c') use ISO_C_BINDING implicit none @@ -103,7 +103,7 @@ module math real(C_DOUBLE), value, intent(in) :: dopp integer(C_INT), value, intent(in) :: n real(C_DOUBLE), intent(inout) :: factors(n) - end subroutine broaden_wmp_polynomials_c_intfc + end subroutine broaden_wmp_polynomials function faddeeva_w(z, relerr) bind(C, name='Faddeeva_w') result(w) use ISO_C_BINDING @@ -116,52 +116,13 @@ module math contains -!=============================================================================== -! T_PERCENTILE calculates the percentile of the Student's t distribution with a -! specified probability level and number of degrees of freedom -!=============================================================================== - - pure function t_percentile(p, df) result(t) bind(C) - - real(C_DOUBLE), intent(in) :: p ! probability level - integer(C_INT), intent(in) :: df ! degrees of freedom - real(C_DOUBLE) :: t ! corresponding t-value - - t = t_percentile_c_intfc(p, df) - - end function t_percentile - -!=============================================================================== -! CALC_PN calculates the n-th order Legendre polynomial at the value of x. -! Since this function is called repeatedly during the neutron transport process, -! neither n or x is checked to see if they are in the applicable range. -! This is left to the client developer to use where applicable. x is to be in -! the domain of [-1,1], and 0<=n<=5. If x is outside of the range, the return -! value will be outside the expected range; if n is outside the stated range, -! the return value will be 1.0. -!=============================================================================== - - pure subroutine calc_pn(n, x, pnx) bind(C) - - integer(C_INT), intent(in) :: n ! Legendre order requested - real(C_DOUBLE), intent(in) :: x ! Independent variable the Legendre is to - ! be evaluated at; x must be in the - ! domain [-1,1] - real(C_DOUBLE), intent(out) :: pnx(n + 1) ! The Legendre polys of order n - ! evaluated at x - - call calc_pn_c_intfc(n, x, pnx) - - end subroutine calc_pn - !=============================================================================== ! EVALUATE_LEGENDRE Find the value of f(x) given a set of Legendre coefficients ! and the value of x !=============================================================================== - pure function evaluate_legendre(n, data, x) result(val) bind(C) - integer(C_INT), intent(in) :: n - real(C_DOUBLE), intent(in) :: data(n) + pure function evaluate_legendre(data, x) result(val) bind(C) + real(C_DOUBLE), intent(in) :: data(:) real(C_DOUBLE), intent(in) :: x real(C_DOUBLE) :: val @@ -169,91 +130,23 @@ contains end function evaluate_legendre -!=============================================================================== -! CALC_RN calculates the n-th order real spherical harmonics for a given angle -! (in terms of (u,v,w)). All Rn,m values are provided (where -n<=m<=n) -!=============================================================================== - - pure subroutine calc_rn(n, uvw, rn) bind(C) - - integer(C_INT), intent(in) :: n ! Order requested - real(C_DOUBLE), intent(in) :: uvw(3) ! Direction of travel; - ! assumed to be on unit sphere - real(C_DOUBLE), intent(out) :: rn(2*n + 1) ! The resultant R_n(uvw) - - call calc_rn_c_intfc(n, uvw, rn) - - end subroutine calc_rn - -!=============================================================================== -! CALC_ZN calculates the n-th order modified Zernike polynomial moment for a -! given angle (rho, theta) location in the unit disk. The normlization of the -! polynomials is such that the integral of Z_pq*Z_pq over the unit disk is -! exactly pi -!=============================================================================== - - pure subroutine calc_zn(n, rho, phi, zn) bind(C) - integer(C_INT), intent(in) :: n ! Maximum order - real(C_DOUBLE), intent(in) :: rho ! Radial location in the unit disk - real(C_DOUBLE), intent(in) :: phi ! Theta (radians) location in the unit disk - ! The resulting list of coefficients - real(C_DOUBLE), intent(out) :: zn(((n + 1) * (n + 2)) / 2) - - call calc_zn_c_intfc(n, rho, phi, zn) - end subroutine calc_zn - !=============================================================================== ! ROTATE_ANGLE rotates direction cosines through a polar angle whose cosine is ! mu and through an azimuthal angle sampled uniformly. Note that this is done ! with direct sampling rather than rejection as is done in MCNP and SERPENT. !=============================================================================== - subroutine rotate_angle(uvw0, mu, uvw, phi) bind(C) - real(C_DOUBLE), intent(in) :: uvw0(3) ! directional cosine - real(C_DOUBLE), intent(in) :: mu ! cosine of angle in lab or CM - real(C_DOUBLE), intent(out) :: uvw(3) ! rotated directional cosine - real(C_DOUBLE), target, optional :: phi ! azimuthal angle + function rotate_angle(uvw0, mu, phi) result(uvw) + real(C_DOUBLE), intent(in) :: uvw0(3) ! directional cosine + real(C_DOUBLE), intent(in) :: mu ! cosine of angle in lab or CM + real(C_DOUBLE), intent(in), optional :: phi ! azimuthal angle - real(C_DOUBLE), pointer :: phi_ptr + real(C_DOUBLE) :: uvw(3) ! rotated directional cosine uvw = uvw0 call rotate_angle_c_intfc(uvw, mu, phi) - end subroutine rotate_angle - -!=============================================================================== -! MAXWELL_SPECTRUM samples an energy from the Maxwell fission distribution based -! on a direct sampling scheme. The probability distribution function for a -! Maxwellian is given as p(x) = 2/(T*sqrt(pi))*sqrt(x/T)*exp(-x/T). This PDF can -! be sampled using rule C64 in the Monte Carlo Sampler LA-9721-MS. -!=============================================================================== - - function maxwell_spectrum(T) result(E_out) bind(C) - - real(C_DOUBLE), intent(in) :: T ! tabulated function of incoming E - real(C_DOUBLE) :: E_out ! sampled energy - - E_out = maxwell_spectrum_c_intfc(T) - - end function maxwell_spectrum - -!=============================================================================== -! WATT_SPECTRUM samples the outgoing energy from a Watt energy-dependent fission -! spectrum. Although fitted parameters exist for many nuclides, generally the -! continuous tabular distributions (LAW 4) should be used in lieu of the Watt -! spectrum. This direct sampling scheme is an unpublished scheme based on the -! original Watt spectrum derivation (See F. Brown's MC lectures). -!=============================================================================== - - function watt_spectrum(a, b) result(E_out) bind(C) - - real(C_DOUBLE), intent(in) :: a ! Watt parameter a - real(C_DOUBLE), intent(in) :: b ! Watt parameter b - real(C_DOUBLE) :: E_out ! energy of emitted neutron - - E_out = watt_spectrum_c_intfc(a, b) - - end function watt_spectrum + end function rotate_angle !=============================================================================== ! FADDEEVA the Faddeeva function, using Stephen Johnson's implementation @@ -305,21 +198,4 @@ contains end select end function w_derivative -!=============================================================================== -! BROADEN_WMP_POLYNOMIALS Doppler broadens the windowed multipole curvefit. The -! curvefit is a polynomial of the form -! a/E + b/sqrt(E) + c + d sqrt(E) ... -!=============================================================================== - - subroutine broaden_wmp_polynomials(E, dopp, n, factors) bind(C) - real(C_DOUBLE), intent(in) :: E ! Energy to evaluate at - real(C_DOUBLE), intent(in) :: dopp ! sqrt(atomic weight ratio / kT), - ! kT given in eV. - integer(C_INT), intent(in) :: n ! number of components to polynomial - real(C_DOUBLE), intent(out):: factors(n) ! output leading coefficient - - call broaden_wmp_polynomials_c_intfc(E, dopp, n, factors) - - end subroutine broaden_wmp_polynomials - end module math diff --git a/src/math_functions.h b/src/math_functions.h index 18bd026405..68e89251e3 100644 --- a/src/math_functions.h +++ b/src/math_functions.h @@ -1,3 +1,6 @@ +//! \file math_functions.h +//! A collection of elementary math functions. + #ifndef MATH_FUNCTIONS_H #define MATH_FUNCTIONS_H diff --git a/src/mgxs_header.F90 b/src/mgxs_header.F90 index bcdef6c792..5abeccbcfa 100644 --- a/src/mgxs_header.F90 +++ b/src/mgxs_header.F90 @@ -1103,9 +1103,7 @@ contains end if scatt_coeffs(gin) % data(imu, gout) = & - evaluate_legendre( & - size(input_scatt(gin) % data, dim=1), & - input_scatt(gin) % data(:, gout), mu) + evaluate_legendre(input_scatt(gin) % data(:, gout), mu) ! Ensure positivity of distribution if (scatt_coeffs(gin) % data(imu, gout) < ZERO) & @@ -2097,7 +2095,6 @@ contains scatt_coeffs(gin, iazi, ipol) % data(imu, gout) = & evaluate_legendre(& - size(input_scatt(gin, iazi, ipol) % data, dim=1), & input_scatt(gin, iazi, ipol) % data(:, gout), mu) ! Ensure positivity of distribution diff --git a/src/physics.F90 b/src/physics.F90 index 8183435cb7..b614b81851 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -495,8 +495,7 @@ contains ! Rotate neutron velocity vector to new angle -- note that the speed of the ! neutron in CM does not change in elastic scattering. However, the speed ! will change when we convert back to LAB - call rotate_angle(uvw_cm, mu_cm, v_n) - v_n = vel * v_n + v_n = vel * rotate_angle(uvw_cm, mu_cm) ! Transform back to LAB frame v_n = v_n + v_cm @@ -786,7 +785,7 @@ contains if (abs(mu) > ONE) mu = sign(ONE,mu) ! change direction of particle - call rotate_angle(uvw, mu, uvw) + uvw = rotate_angle(uvw, mu) end subroutine sab_scatter @@ -977,8 +976,7 @@ contains if (abs(mu) < ONE) then ! set and accept target velocity E_t = E_t / awr - call rotate_angle(uvw, mu, v_target) - v_target = sqrt(E_t) * v_target + v_target = sqrt(E_t) * rotate_angle(uvw, mu) exit ARES_REJECT_LOOP end if end do ARES_REJECT_LOOP @@ -1058,8 +1056,7 @@ contains ! Determine velocity vector of target nucleus based on neutron's velocity ! and the sampled angle between them - call rotate_angle(uvw, mu, v_target) - v_target = vt * v_target + v_target = vt * rotate_angle(uvw, mu) end subroutine sample_cxs_target_velocity @@ -1330,7 +1327,7 @@ contains p % mu = mu ! change direction of particle - call rotate_angle(p % coord(1) % uvw, mu, p % coord(1) % uvw) + p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, mu) ! evaluate yield yield = rxn % products(1) % yield % evaluate(E_in) diff --git a/src/physics_mg.F90 b/src/physics_mg.F90 index 05493af057..eb82085941 100644 --- a/src/physics_mg.F90 +++ b/src/physics_mg.F90 @@ -151,7 +151,7 @@ contains p % E = energy_bin_avg(p % g) ! Convert change in angle (mu) to new direction - call rotate_angle(p % coord(1) % uvw, p % mu, p % coord(1) % uvw) + p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, p % mu) ! Set event component p % event = EVENT_SCATTER diff --git a/src/scattdata_header.F90 b/src/scattdata_header.F90 index a3108a2df6..511e2a2376 100644 --- a/src/scattdata_header.F90 +++ b/src/scattdata_header.F90 @@ -445,8 +445,7 @@ contains if (gout < this % gmin(gin) .or. gout > this % gmax(gin)) then f = ZERO else - f = evaluate_legendre(size(this % dist(gin) % data, dim=1), & - this % dist(gin) % data(:, gout), mu) + f = evaluate_legendre(this % dist(gin) % data(:, gout), mu) end if end function scattdatalegendre_calc_f