mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 22:26:08 -04:00
Exposed math.F90 to C and then created python wrappers around these functions. Seem to work
This commit is contained in:
Adam G Nelson
parent
a63a62f38d
commit
600969ce2b
9 changed files with 397 additions and 119 deletions
|
|
@ -47,3 +47,4 @@ from .mesh import *
|
|||
from .filter import *
|
||||
from .tally import *
|
||||
from .settings import settings
|
||||
from .math import *
|
||||
|
|
|
|||
263
openmc/capi/math.py
Normal file
263
openmc/capi/math.py
Normal file
|
|
@ -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
|
||||
13
src/api.F90
13
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
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
||||
|
|
|
|||
211
src/math.F90
211
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
|
||||
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue