From 600969ce2bb726b812b21d8f5bba263f7dc67817 Mon Sep 17 00:00:00 2001 From: Adam G Nelson Date: Sun, 29 Apr 2018 06:35:56 -0400 Subject: [PATCH] Exposed math.F90 to C and then created python wrappers around these functions. Seem to work --- openmc/capi/__init__.py | 1 + openmc/capi/math.py | 263 ++++++++++++++++++++++++++ src/api.F90 | 13 ++ src/distribution_multivariate.F90 | 2 +- src/math.F90 | 211 ++++++++++----------- src/physics.F90 | 13 +- src/physics_mg.F90 | 2 +- src/tallies/tally.F90 | 9 +- src/tallies/tally_filter_sph_harm.F90 | 2 +- 9 files changed, 397 insertions(+), 119 deletions(-) create mode 100644 openmc/capi/math.py diff --git a/openmc/capi/__init__.py b/openmc/capi/__init__.py index 217e782a84..672a6e811f 100644 --- a/openmc/capi/__init__.py +++ b/openmc/capi/__init__.py @@ -47,3 +47,4 @@ from .mesh import * from .filter import * from .tally import * from .settings import settings +from .math import * diff --git a/openmc/capi/math.py b/openmc/capi/math.py new file mode 100644 index 0000000000..d51858fd3b --- /dev/null +++ b/openmc/capi/math.py @@ -0,0 +1,263 @@ +from ctypes import (c_int, c_double, POINTER) + +import numpy as np +from numpy.ctypeslib import ndpointer + +from . import _dll + +_dll.normal_percentile.restype = c_double +_dll.normal_percentile.argtypes = [POINTER(c_double)] +_dll.t_percentile.restype = c_double +_dll.t_percentile.argtypes = [POINTER(c_double), POINTER(c_int)] +_dll.calc_pn.restype = c_double +_dll.calc_pn.argtypes = [POINTER(c_int), POINTER(c_double)] +_dll.calc_rn.restype = None +_dll.calc_rn.argtypes = [POINTER(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.evaluate_legendre.restype = c_double +_dll.evaluate_legendre.argtypes = [ndpointer(c_double), POINTER(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.watt_spectrum.restype = c_double +_dll.watt_spectrum.argtypes = [POINTER(c_double), POINTER(c_double)] +# COMPLEX DATA TYPES? +# _dll.faddeeva.restype = c_double +# _dll.faddeeva.argtypes = [POINTER(c_double)] + +# _dll.w_derivative.restype = c_double +# _dll.w_derivative.argtypes = [POINTER(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)] + + +def normal_percentile(p): + """ Calculate the percentile of the standard normal distribution with a + specified probability level + + Parameters + ---------- + p : float + Probability level + + Returns + ------- + float + Corresponding z-value + + """ + + return _dll.normal_percentile(c_double(p)) + + +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(c_double(p), c_int(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 + + """ + + return _dll.calc_pn(c_int(n), 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 (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 = 2 * 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) + 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 + Maximum order + rho : float + 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 + zn = np.zeros(num_bins, dtype=np.float64) + _dll.calc_zn(c_int(n), c_double(rho), c_double(phi), zn) + return zn + + +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(data_arr, c_double(x)) + + +def rotate_angle(uvw0, mu, phi=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, optional + Azimuthal angle; if None, one will be sampled uniformly + + Returns + ------- + numpy.ndarray + Rotated direction cosine + + """ + + uvw = np.zeros(3, dtype=np.float64) + uvw0_arr = np.array(uvw0, dtype=np.float64) + _dll.rotate_angle(uvw0_arr, c_double(mu), uvw, c_double(phi)) + return uvw + + +def maxwell_spectrum(T): + """ Samples an energy from the Maxwell fission distribution based + on a direct sampling scheme. + + Parameters + ---------- + T : float + Spectrum parameter + + Returns + ------- + float + Sampled outgoing energy + + """ + + return _dll.maxwell_spectrum(c_double(T)) + + +def watt_spectrum(a, b): + """ Samples an energy from the Watt energy-dependent fission spectrum. + + Parameters + ---------- + a : float + Spectrum parameter a + b : float + Spectrum parameter b + + Returns + ------- + float + Sampled outgoing energy + + """ + + return _dll.watt_spectrum(c_double(a), c_double(b)) + + +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(c_double(E), c_double(dopp), c_int(n), + factors) + return factors diff --git a/src/api.F90 b/src/api.F90 index df183c292f..30aca68343 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -12,6 +12,7 @@ module openmc_api use geometry_header use hdf5_interface use material_header + use math use mesh_header use message_passing use nuclide_header @@ -97,6 +98,18 @@ module openmc_api public :: openmc_tally_set_nuclides public :: openmc_tally_set_scores public :: openmc_tally_set_type + public :: normal_percentile + 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 650ab26fa2..4d7d42cc3d 100644 --- a/src/distribution_multivariate.F90 +++ b/src/distribution_multivariate.F90 @@ -119,7 +119,7 @@ contains else ! Sample azimuthal angle phi = this % phi % sample() - uvw(:) = rotate_angle(this % reference_uvw, mu, phi) + call rotate_angle(this % reference_uvw, mu, uvw, phi) end if end function polar_azimuthal_sample diff --git a/src/math.F90 b/src/math.F90 index ee8cd0530b..47e6acefb0 100644 --- a/src/math.F90 +++ b/src/math.F90 @@ -6,6 +6,19 @@ module math use random_lcg, only: prn implicit none + private + public :: normal_percentile + 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 !=============================================================================== ! FADDEEVA_W evaluates the scaled complementary error function. This @@ -29,24 +42,24 @@ contains ! distribution with a specified probability level !=============================================================================== - elemental function normal_percentile(p) result(z) + pure function normal_percentile(p) result(z) bind(C) - real(8), intent(in) :: p ! probability level - real(8) :: z ! corresponding z-value + real(C_DOUBLE), intent(in) :: p ! probability level + real(C_DOUBLE) :: z ! corresponding z-value - real(8) :: q - real(8) :: r - real(8), parameter :: p_low = 0.02425_8 - real(8), parameter :: a(6) = (/ & + real(C_DOUBLE) :: q + real(C_DOUBLE) :: r + real(C_DOUBLE), parameter :: p_low = 0.02425_8 + real(C_DOUBLE), parameter :: a(6) = (/ & -3.969683028665376e1_8, 2.209460984245205e2_8, -2.759285104469687e2_8, & 1.383577518672690e2_8, -3.066479806614716e1_8, 2.506628277459239e0_8 /) - real(8), parameter :: b(5) = (/ & + real(C_DOUBLE), parameter :: b(5) = (/ & -5.447609879822406e1_8, 1.615858368580409e2_8, -1.556989798598866e2_8, & 6.680131188771972e1_8, -1.328068155288572e1_8 /) - real(8), parameter :: c(6) = (/ & + real(C_DOUBLE), parameter :: c(6) = (/ & -7.784894002430293e-3_8, -3.223964580411365e-1_8, -2.400758277161838_8, & -2.549732539343734_8, 4.374664141464968_8, 2.938163982698783_8 /) - real(8), parameter :: d(4) = (/ & + real(C_DOUBLE), parameter :: d(4) = (/ & 7.784695709041462e-3_8, 3.224671290700398e-1_8, & 2.445134137142996_8, 3.754408661907416_8 /) @@ -88,16 +101,16 @@ contains ! specified probability level and number of degrees of freedom !=============================================================================== - elemental function t_percentile(p, df) result(t) + pure function t_percentile(p, df) result(t) bind(C) - real(8), intent(in) :: p ! probability level - integer, intent(in) :: df ! degrees of freedom - real(8) :: t ! corresponding t-value + real(C_DOUBLE), intent(in) :: p ! probability level + integer(C_INT), intent(in) :: df ! degrees of freedom + real(C_DOUBLE) :: t ! corresponding t-value - real(8) :: n ! degrees of freedom as a real(8) - real(8) :: k ! n - 2 - real(8) :: z ! percentile of normal distribution - real(8) :: z2 ! z * z + real(C_DOUBLE) :: n ! degrees of freedom as a real(8) + real(C_DOUBLE) :: k ! n - 2 + real(C_DOUBLE) :: z ! percentile of normal distribution + real(C_DOUBLE) :: z2 ! z * z if (df == 1) then ! For one degree of freedom, the t-distribution becomes a Cauchy @@ -140,12 +153,14 @@ contains ! the return value will be 1.0. !=============================================================================== - elemental function calc_pn(n,x) result(pnx) + pure function calc_pn(n,x) result(pnx) bind(C) - integer, intent(in) :: n ! Legendre order requested - real(8), intent(in) :: x ! Independent variable the Legendre is to be - ! evaluated at; x must be in the domain [-1,1] - real(8) :: pnx ! The Legendre poly of order n evaluated at x + 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) :: pnx ! The Legendre poly of order n evaluated + ! at x select case(n) case(1) @@ -185,14 +200,16 @@ contains ! (in terms of (u,v,w)). All Rn,m values are provided (where -n<=m<=n) !=============================================================================== - pure function calc_rn(n,uvw) result(rn) + subroutine calc_rn(n, uvw, rn) bind(C) - integer, intent(in) :: n ! Order requested - real(8), intent(in) :: uvw(3) ! Direction of travel, assumed to be on unit sphere - real(8) :: rn(2*n + 1) ! The resultant R_n(uvw) + 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) :: rn(2*n + 1) ! The resultant R_n(uvw) - real(8) :: phi, w ! Azimuthal and Cosine of Polar angles (from uvw) - real(8) :: w2m1 ! (w^2 - 1), frequently used in these + + real(C_DOUBLE) :: phi, w ! Azimuthal and Cosine of Polar angles (from uvw) + real(C_DOUBLE) :: w2m1 ! (w^2 - 1), frequently used in these w = uvw(3) ! z = cos(polar) if (uvw(1) == ZERO) then @@ -572,7 +589,7 @@ contains rn = ONE end select - end function calc_rn + end subroutine calc_rn !=============================================================================== ! CALC_ZN calculates the n-th order modified Zernike polynomial moment for a @@ -581,31 +598,31 @@ contains ! exactly pi !=============================================================================== - subroutine calc_zn(n, rho, phi, zn) + subroutine calc_zn(n, rho, phi, zn) bind(C) ! This procedure uses the modified Kintner's method for calculating Zernike ! polynomials as outlined in Chong, C. W., Raveendran, P., & Mukundan, ! R. (2003). A comparative analysis of algorithms for fast computation of ! Zernike moments. Pattern Recognition, 36(3), 731-742. - integer, intent(in) :: n ! Maximum order - real(8), intent(in) :: rho ! Radial location in the unit disk - real(8), intent(in) :: phi ! Theta (radians) location in the unit disk - real(8), intent(out) :: zn(:) ! The resulting list of coefficients + 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 + real(C_DOUBLE), intent(out) :: zn(:) ! The resulting list of coefficients - real(8) :: sin_phi, cos_phi ! Sine and Cosine of phi - real(8) :: sin_phi_vec(n+1) ! Contains sin(n*phi) - real(8) :: cos_phi_vec(n+1) ! Contains cos(n*phi) - real(8) :: zn_mat(n+1, n+1) ! Matrix form of the coefficients which is - ! easier to work with - real(8) :: k1, k2, k3, k4 ! Variables for R_m_n calculation - real(8) :: sqrt_norm ! normalization for radial moments - integer :: i,p,q ! Loop counters + real(C_DOUBLE) :: sin_phi, cos_phi ! Sine and Cosine of phi + real(C_DOUBLE) :: sin_phi_vec(n+1) ! Contains sin(n*phi) + real(C_DOUBLE) :: cos_phi_vec(n+1) ! Contains cos(n*phi) + real(C_DOUBLE) :: zn_mat(n+1, n+1) ! Matrix form of the coefficients which is + ! easier to work with + real(C_DOUBLE) :: k1, k2, k3, k4 ! Variables for R_m_n calculation + real(C_DOUBLE) :: sqrt_norm ! normalization for radial moments + integer(C_INT) :: i,p,q ! Loop counters - real(8), parameter :: SQRT_N_1(0:10) = [& + real(C_DOUBLE), parameter :: SQRT_N_1(0:10) = [& sqrt(1.0_8), sqrt(2.0_8), sqrt(3.0_8), sqrt(4.0_8), & sqrt(5.0_8), sqrt(6.0_8), sqrt(7.0_8), sqrt(8.0_8), & sqrt(9.0_8), sqrt(10.0_8), sqrt(11.0_8)] - real(8), parameter :: SQRT_2N_2(0:10) = SQRT_N_1*sqrt(2.0_8) + real(C_DOUBLE), parameter :: SQRT_2N_2(0:10) = SQRT_N_1*sqrt(2.0_8) ! n == radial degree ! m == azimuthal frequency @@ -681,41 +698,17 @@ contains end do end subroutine calc_zn -!=============================================================================== -! EXPAND_HARMONIC expands a given series of real spherical harmonics -!=============================================================================== - - pure function expand_harmonic(data, order, uvw) result(val) - real(8), intent(in) :: data(:) - integer, intent(in) :: order - real(8), intent(in) :: uvw(3) - real(8) :: val - - integer :: l, lm_lo, lm_hi - - val = data(1) - lm_lo = 2 - lm_hi = 4 - do l = 1, order - 1 - val = val + sqrt(TWO * real(l,8) + ONE) * & - dot_product(calc_rn(l,uvw), data(lm_lo:lm_hi)) - lm_lo = lm_hi + 1 - lm_hi = lm_lo + 2 * (l + 1) - end do - - end function expand_harmonic - !=============================================================================== ! EVALUATE_LEGENDRE Find the value of f(x) given a set of Legendre coefficients ! and the value of x !=============================================================================== - pure function evaluate_legendre(data, x) result(val) - real(8), intent(in) :: data(:) - real(8), intent(in) :: x - real(8) :: val + 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 - integer :: l + integer(C_INT) :: l val = HALF * data(1) do l = 1, size(data) - 1 @@ -730,20 +723,20 @@ contains ! with direct sampling rather than rejection as is done in MCNP and SERPENT. !=============================================================================== - function rotate_angle(uvw0, mu, phi) result(uvw) - real(8), intent(in) :: uvw0(3) ! directional cosine - real(8), intent(in) :: mu ! cosine of angle in lab or CM - real(8), optional :: phi ! azimuthal angle - real(8) :: uvw(3) ! rotated directional cosine + 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), optional :: phi ! azimuthal angle - real(8) :: phi_ ! azimuthal angle - real(8) :: sinphi ! sine of azimuthal angle - real(8) :: cosphi ! cosine of azimuthal angle - real(8) :: a ! sqrt(1 - mu^2) - real(8) :: b ! sqrt(1 - w^2) - real(8) :: u0 ! original cosine in x direction - real(8) :: v0 ! original cosine in y direction - real(8) :: w0 ! original cosine in z direction + real(C_DOUBLE) :: phi_ ! azimuthal angle + real(C_DOUBLE) :: sinphi ! sine of azimuthal angle + real(C_DOUBLE) :: cosphi ! cosine of azimuthal angle + real(C_DOUBLE) :: a ! sqrt(1 - mu^2) + real(C_DOUBLE) :: b ! sqrt(1 - w^2) + real(C_DOUBLE) :: u0 ! original cosine in x direction + real(C_DOUBLE) :: v0 ! original cosine in y direction + real(C_DOUBLE) :: w0 ! original cosine in z direction ! Copy original directional cosines u0 = uvw0(1) @@ -776,7 +769,7 @@ contains uvw(3) = mu*w0 + a*(v0*w0*cosphi - u0*sinphi)/b end if - end function rotate_angle + end subroutine rotate_angle !=============================================================================== ! MAXWELL_SPECTRUM samples an energy from the Maxwell fission distribution based @@ -785,13 +778,13 @@ contains ! be sampled using rule C64 in the Monte Carlo Sampler LA-9721-MS. !=============================================================================== - function maxwell_spectrum(T) result(E_out) + function maxwell_spectrum(T) result(E_out) bind(C) - real(8), intent(in) :: T ! tabulated function of incoming E - real(8) :: E_out ! sampled energy + real(C_DOUBLE), intent(in) :: T ! tabulated function of incoming E + real(C_DOUBLE) :: E_out ! sampled energy - real(8) :: r1, r2, r3 ! random numbers - real(8) :: c ! cosine of pi/2*r3 + real(C_DOUBLE) :: r1, r2, r3 ! random numbers + real(C_DOUBLE) :: c ! cosine of pi/2*r3 r1 = prn() r2 = prn() @@ -813,13 +806,13 @@ contains ! original Watt spectrum derivation (See F. Brown's MC lectures). !=============================================================================== - function watt_spectrum(a, b) result(E_out) + function watt_spectrum(a, b) result(E_out) bind(C) - real(8), intent(in) :: a ! Watt parameter a - real(8), intent(in) :: b ! Watt parameter b - real(8) :: E_out ! energy of emitted neutron + 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 - real(8) :: w ! sampled from Maxwellian + 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) @@ -830,9 +823,9 @@ contains ! FADDEEVA the Faddeeva function, using Stephen Johnson's implementation !=============================================================================== - function faddeeva(z) result(wv) - complex(C_DOUBLE_COMPLEX), intent(in) :: z ! The point to evaluate Z at - complex(8) :: wv ! The resulting w(z) value + function faddeeva(z) result(wv) bind(C) + complex(C_DOUBLE_COMPLEX), intent(in) :: z ! The point to evaluate Z at + complex(C_DOUBLE_COMPLEX) :: wv ! The resulting w(z) value real(C_DOUBLE) :: relerr ! Target relative error in inner loop of MIT ! Faddeeva @@ -860,10 +853,10 @@ contains end function faddeeva - recursive function w_derivative(z, order) result(wv) + recursive function w_derivative(z, order) result(wv) bind(C) complex(C_DOUBLE_COMPLEX), intent(in) :: z ! The point to evaluate Z at - integer, intent(in) :: order - complex(8) :: wv ! The resulting w(z) value + integer(C_INT), intent(in) :: order + complex(C_DOUBLE_COMPLEX) :: wv ! The resulting w(z) value select case(order) case (0) @@ -882,12 +875,12 @@ contains ! a/E + b/sqrt(E) + c + d sqrt(E) ... !=============================================================================== - subroutine broaden_wmp_polynomials(E, dopp, n, factors) - real(8), intent(in) :: E ! Energy to evaluate at - real(8), intent(in) :: dopp ! sqrt(atomic weight ratio / kT), + 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, intent(in) :: n ! number of components to polynomial - real(8), intent(out):: factors(n) ! output leading coefficient + integer(C_INT), intent(in) :: n ! number of components to polynomial + real(C_DOUBLE), intent(out):: factors(n) ! output leading coefficient integer :: i diff --git a/src/physics.F90 b/src/physics.F90 index b614b81851..8183435cb7 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -495,7 +495,8 @@ 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 - v_n = vel * rotate_angle(uvw_cm, mu_cm) + call rotate_angle(uvw_cm, mu_cm, v_n) + v_n = vel * v_n ! Transform back to LAB frame v_n = v_n + v_cm @@ -785,7 +786,7 @@ contains if (abs(mu) > ONE) mu = sign(ONE,mu) ! change direction of particle - uvw = rotate_angle(uvw, mu) + call rotate_angle(uvw, mu, uvw) end subroutine sab_scatter @@ -976,7 +977,8 @@ contains if (abs(mu) < ONE) then ! set and accept target velocity E_t = E_t / awr - v_target = sqrt(E_t) * rotate_angle(uvw, mu) + call rotate_angle(uvw, mu, v_target) + v_target = sqrt(E_t) * v_target exit ARES_REJECT_LOOP end if end do ARES_REJECT_LOOP @@ -1056,7 +1058,8 @@ contains ! Determine velocity vector of target nucleus based on neutron's velocity ! and the sampled angle between them - v_target = vt * rotate_angle(uvw, mu) + call rotate_angle(uvw, mu, v_target) + v_target = vt * v_target end subroutine sample_cxs_target_velocity @@ -1327,7 +1330,7 @@ contains p % mu = mu ! change direction of particle - p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, mu) + call rotate_angle(p % coord(1) % uvw, mu, p % coord(1) % uvw) ! evaluate yield yield = rxn % products(1) % yield % evaluate(E_in) diff --git a/src/physics_mg.F90 b/src/physics_mg.F90 index eb82085941..05493af057 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 - p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, p % mu) + call rotate_angle(p % coord(1) % uvw, p % mu, p % coord(1) % uvw) ! Set event component p % event = EVENT_SCATTER diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index f4d12be7e1..f7d3c8a08e 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -46,6 +46,9 @@ module tally end subroutine score_analog_tally_ end interface + real(C_DOUBLE) :: rn(2 * MAX_ANG_ORDER + 1) +!$omp threadprivate(rn) + contains !=============================================================================== @@ -2131,11 +2134,12 @@ contains num_nm = 2 * n + 1 ! multiply score by the angular flux moments and store + call calc_rn(n, p % last_uvw, rn(1:num_nm)) !$omp critical (score_general_scatt_yn) t % results(RESULT_VALUE, score_index: score_index + num_nm - 1, & filter_index) = t % results(RESULT_VALUE, & score_index: score_index + num_nm - 1, filter_index) & - + score * calc_pn(n, p % mu) * calc_rn(n, p % last_uvw) + + score * calc_pn(n, p % mu) * rn(1:num_nm) !$omp end critical (score_general_scatt_yn) end do i = i + (t % moment_order(i) + 1)**2 - 1 @@ -2159,11 +2163,12 @@ contains num_nm = 2 * n + 1 ! multiply score by the angular flux moments and store + call calc_rn(n, uvw, rn(1:num_nm)) !$omp critical (score_general_flux_tot_yn) t % results(RESULT_VALUE, score_index: score_index + num_nm - 1, & filter_index) = t % results(RESULT_VALUE, & score_index: score_index + num_nm - 1, filter_index) & - + score * calc_rn(n, uvw) + + score * rn(1:num_nm) !$omp end critical (score_general_flux_tot_yn) end do i = i + (t % moment_order(i) + 1)**2 - 1 diff --git a/src/tallies/tally_filter_sph_harm.F90 b/src/tallies/tally_filter_sph_harm.F90 index 5c47c2c71f..5faf2b6fdb 100644 --- a/src/tallies/tally_filter_sph_harm.F90 +++ b/src/tallies/tally_filter_sph_harm.F90 @@ -91,7 +91,7 @@ contains ! Calculate n-th order spherical harmonics for (u,v,w) num_nm = 2*n + 1 - rn(1:num_nm) = calc_rn(n, p % last_uvw) + call calc_rn(n, p % last_uvw, rn(1:num_nm)) ! Append matching (bin,weight) for each moment do i = 1, num_nm