crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 13:45:36 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Moved math python api to point to the C++ and simplified the fortran methods

Browse source
This commit is contained in:
Adam G Nelson 2018-05-18 19:07:34 -04:00
parent 38c13994bc
commit 890a1b4f55
9 changed files with 67 additions and 211 deletions
Download patch file Download diff file Expand all files Collapse all files

71
openmc/capi/math.py
View file

@ -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

11
src/api.F90
View file

@ -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

2
src/distribution_multivariate.F90
View file

@ -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

168
src/math.F90
View file

@ -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

3
src/math_functions.h
View file

@ -1,3 +1,6 @@
//! \file math_functions.h
//! A collection of elementary math functions.
#ifndef MATH_FUNCTIONS_H
#define MATH_FUNCTIONS_H

5
src/mgxs_header.F90
View file

@ -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

13
src/physics.F90
View file

@ -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)

2
src/physics_mg.F90
View file

@ -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

3
src/scattdata_header.F90
View file

@ -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