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

Finish converting windowed multipole functions to C++

Browse source
This commit is contained in:
Paul Romano 2018-12-28 07:19:20 -06:00
parent 149588b74e
commit 38ff673cf3
18 changed files with 417 additions and 619 deletions
Download patch file Download diff file Expand all files Collapse all files

8
CMakeLists.txt
View file

@ -294,11 +294,9 @@ endif()
# faddeeva library
#===============================================================================
add_library(faddeeva STATIC vendor/faddeeva/Faddeeva.c)
target_compile_options(faddeeva PRIVATE ${cflags})
set_target_properties(faddeeva PROPERTIES
C_STANDARD 99
C_STANDARD_REQUIRED ON)
add_library(faddeeva STATIC vendor/faddeeva/Faddeeva.cc)
target_include_directories(faddeeva PUBLIC vendor/faddeeva/)
target_compile_options(faddeeva PRIVATE ${cxxflags})
#===============================================================================
# libopenmc

4
include/openmc/endf.h
View file

@ -99,6 +99,10 @@ private:
std::vector<double> factors_; //!< Partial sums of structure factors [eV-b]
};
//! Read 1D function from HDF5 dataset
//! \param[in] group HDF5 group containing dataset
//! \param[in] name Name of dataset
//! \return Unique pointer to 1D function
std::unique_ptr<Function1D> read_function(hid_t group, const char* name);
} // namespace openmc

122
include/openmc/math_functions.h
View file

@ -5,6 +5,7 @@
#define OPENMC_MATH_FUNCTIONS_H
#include <cmath>
#include <complex>
#include <cstdlib>
#include "openmc/constants.h"
@ -18,8 +19,8 @@ namespace openmc {
//! Calculate the percentile of the standard normal distribution with a
//! specified probability level.
//!
//! @param p The probability level
//! @return The requested percentile
//! \param p The probability level
//! \return The requested percentile
//==============================================================================
extern "C" double normal_percentile(double p);
@ -28,9 +29,9 @@ extern "C" double normal_percentile(double p);
//! Calculate the percentile of the Student's t distribution with a specified
//! probability level and number of degrees of freedom.
//!
//! @param p The probability level
//! @param df The degrees of freedom
//! @return The requested percentile
//! \param p The probability level
//! \param df The degrees of freedom
//! \return The requested percentile
//==============================================================================
extern "C" double t_percentile_c(double p, int df);
@ -38,9 +39,9 @@ extern "C" double t_percentile_c(double p, int df);
//==============================================================================
//! Calculate the n-th order Legendre polynomials at the value of x.
//!
//! @param n The maximum order requested
//! @param x The value to evaluate at; x is expected to be within [-1,1]
//! @param pnx The requested Legendre polynomials of order 0 to n (inclusive)
//! \param n The maximum order requested
//! \param x The value to evaluate at; x is expected to be within [-1,1]
//! \param pnx The requested Legendre polynomials of order 0 to n (inclusive)
//! evaluated at x.
//==============================================================================
@ -50,11 +51,11 @@ extern "C" void calc_pn_c(int n, double x, double pnx[]);
//! Find the value of f(x) given a set of Legendre coefficients and the value
//! of x.
//!
//! @param n The maximum order of the expansion
//! @param data The polynomial expansion coefficient data; without the (2l+1)/2
//! \param n The maximum order of the expansion
//! \param data The polynomial expansion coefficient data; without the (2l+1)/2
//! factor.
//! @param x The value to evaluate at; x is expected to be within [-1,1]
//! @return The requested Legendre polynomials of order 0 to n (inclusive)
//! \param x The value to evaluate at; x is expected to be within [-1,1]
//! \return The requested Legendre polynomials of order 0 to n (inclusive)
//! evaluated at x
//==============================================================================
@ -64,9 +65,9 @@ extern "C" double evaluate_legendre(int n, const double data[], double x);
//! Calculate the n-th order real spherical harmonics for a given angle (in
//! terms of (u,v,w)) for all 0<=n and -m<=n<=n.
//!
//! @param n The maximum order requested
//! @param uvw[3] The direction the harmonics are requested at
//! @param rn The requested harmonics of order 0 to n (inclusive)
//! \param n The maximum order requested
//! \param uvw[3] The direction the harmonics are requested at
//! \param rn The requested harmonics of order 0 to n (inclusive)
//! evaluated at uvw.
//==============================================================================
@ -83,10 +84,10 @@ extern "C" void calc_rn_c(int n, const double uvw[3], double rn[]);
//! The normalization of the polynomials is such that the integral of Z_pq^2
//! over the unit disk is exactly pi.
//!
//! @param n The maximum order requested
//! @param rho The radial parameter to specify location on the unit disk
//! @param phi The angle parameter to specify location on the unit disk
//! @param zn The requested moments of order 0 to n (inclusive)
//! \param n The maximum order requested
//! \param rho The radial parameter to specify location on the unit disk
//! \param phi The angle parameter to specify location on the unit disk
//! \param zn The requested moments of order 0 to n (inclusive)
//! evaluated at rho and phi.
//==============================================================================
@ -106,10 +107,10 @@ extern "C" void calc_zn(int n, double rho, double phi, double zn[]);
//! The normalization of the polynomials is such that the integral of Z_pq^2
//! over the unit disk is exactly pi.
//!
//! @param n The maximum order requested
//! @param rho The radial parameter to specify location on the unit disk
//! @param phi The angle parameter to specify location on the unit disk
//! @param zn_rad The requested moments of order 0 to n (inclusive)
//! \param n The maximum order requested
//! \param rho The radial parameter to specify location on the unit disk
//! \param phi The angle parameter to specify location on the unit disk
//! \param zn_rad The requested moments of order 0 to n (inclusive)
//! evaluated at rho and phi when m = 0.
//==============================================================================
@ -122,9 +123,9 @@ extern "C" void calc_zn_rad(int n, double rho, double zn_rad[]);
//! This is done with direct sampling rather than rejection sampling as is done
//! in MCNP and Serpent.
//!
//! @param uvw[3] The initial, and final, direction vector
//! @param mu The cosine of angle in lab or CM
//! @param phi The azimuthal angle; will randomly chosen angle if a nullptr
//! \param uvw[3] The initial, and final, direction vector
//! \param mu The cosine of angle in lab or CM
//! \param phi The azimuthal angle; will randomly chosen angle if a nullptr
//! is passed
//==============================================================================
@ -140,8 +141,8 @@ Direction rotate_angle(Direction u, double mu, double* phi);
//! 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.
//!
//! @param T The tabulated function of the incoming energy
//! @result The sampled outgoing energy
//! \param T The tabulated function of the incoming energy
//! \return The sampled outgoing energy
//==============================================================================
extern "C" double maxwell_spectrum(double T);
@ -154,9 +155,9 @@ extern "C" double maxwell_spectrum(double T);
//! spectrum. This direct sampling scheme is an unpublished scheme based on the
//! original Watt spectrum derivation (See F. Brown's MC lectures).
//!
//! @param a Watt parameter a
//! @param b Watt parameter b
//! @result The sampled outgoing energy
//! \param a Watt parameter a
//! \param b Watt parameter b
//! \return The sampled outgoing energy
//==============================================================================
extern "C" double watt_spectrum(double a, double b);
@ -166,14 +167,13 @@ extern "C" double watt_spectrum(double a, double b);
//!
//! The curvefit is a polynomial of the form a/E + b/sqrt(E) + c + d sqrt(E)...
//!
//! @param E The energy to evaluate the broadening at
//! @param dopp sqrt(atomic weight ratio / kT) with kT given in eV
//! @param n The number of components to the polynomial
//! @param factors The output leading coefficient
//! \param E The energy to evaluate the broadening at
//! \param dopp sqrt(atomic weight ratio / kT) with kT given in eV
//! \param n The number of components to the polynomial
//! \param factors The output leading coefficient
//==============================================================================
extern "C" void broaden_wmp_polynomials_c(double E, double dopp, int n,
double factors[]);
extern "C" void broaden_wmp_polynomials(double E, double dopp, int n, double factors[]);
//==============================================================================
//! Constructs a natural cubic spline.
@ -183,11 +183,11 @@ extern "C" void broaden_wmp_polynomials_c(double E, double dopp, int n,
//! used in any subsequent calls to spline_interpolate or spline_integrate for
//! the same set of x and y values.
//!
//! @param n Number of points
//! @param x Values of the independent variable, which must be strictly
//! \param n Number of points
//! \param x Values of the independent variable, which must be strictly
//! increasing.
//! @param y Values of the dependent variable.
//! @param[out] z The second derivative of the interpolating function at each
//! \param y Values of the dependent variable.
//! \param[out] z The second derivative of the interpolating function at each
//! value of x.
//==============================================================================
@ -196,14 +196,14 @@ extern "C" void spline_c(int n, const double x[], const double y[], double z[]);
//==============================================================================
//! Determine the cubic spline interpolated y-value for a given x-value.
//!
//! @param n Number of points
//! @param x Values of the independent variable, which must be strictly
//! \param n Number of points
//! \param x Values of the independent variable, which must be strictly
//! increasing.
//! @param y Values of the dependent variable.
//! @param z The second derivative of the interpolating function at each
//! \param y Values of the dependent variable.
//! \param z The second derivative of the interpolating function at each
//! value of x.
//! @param xint Point at which to evaluate the cubic spline polynomial
//! @result Interpolated value
//! \param xint Point at which to evaluate the cubic spline polynomial
//! \return Interpolated value
//==============================================================================
extern "C" double spline_interpolate_c(int n, const double x[], const double y[],
@ -213,19 +213,33 @@ extern "C" double spline_interpolate_c(int n, const double x[], const double y[]
//! Evaluate the definite integral of the interpolating cubic spline between
//! the given endpoints.
//!
//! @param n Number of points
//! @param x Values of the independent variable, which must be strictly
//! \param n Number of points
//! \param x Values of the independent variable, which must be strictly
//! increasing.
//! @param y Values of the dependent variable.
//! @param z The second derivative of the interpolating function at each
//! \param y Values of the dependent variable.
//! \param z The second derivative of the interpolating function at each
//! value of x.
//! @param xa Lower limit of integration
//! @param xb Upper limit of integration
//! @result Integral
//! \param xa Lower limit of integration
//! \param xb Upper limit of integration
//! \return Integral
//==============================================================================
extern "C" double spline_integrate_c(int n, const double x[], const double y[],
const double z[], double xa, double xb);
//! Evaluate the Faddeeva function
//!
//! \param z Complex argument
//! \return Faddeeva function evaluated at z
std::complex<double> faddeeva(std::complex<double> z);
//! Evaluate derivative of the Faddeeva function
//!
//! \param z Complex argument
//! \param order Order of the derivative
//! \return Derivative of Faddeeva function evaluated at z
std::complex<double> w_derivative(std::complex<double> z, int order);
} // namespace openmc
#endif // OPENMC_MATH_FUNCTIONS_H

11
include/openmc/nuclide.h
View file

@ -14,7 +14,8 @@
#include "openmc/endf.h"
#include "openmc/reaction.h"
#include "openmc/reaction_product.h"
#include "urr.h"
#include "openmc/urr.h"
#include "openmc/wmp.h"
namespace openmc {
@ -130,6 +131,9 @@ public:
std::vector<EnergyGrid> grid_; //!< Energy grid at each temperature
std::vector<xt::xtensor<double, 2>> xs_; //!< Cross sections at each temperature
// Multipole data
std::unique_ptr<WindowedMultipole> multipole_;
// Fission data
bool fissionable_ {false}; //!< Whether nuclide is fissionable
bool has_partial_fission_ {false}; //!< has partial fission reactions?
@ -193,11 +197,6 @@ extern "C" MaterialMacroXS material_xs;
//==============================================================================
extern "C" void set_micro_xs();
extern "C" void nuclide_multipole_eval(int i_nuclide, double E, double sqrtkT,
double* sig_s, double* sig_a, double* sig_f);
extern "C" bool nuclide_wmp_present(int i_nuclide);
extern "C" double nuclide_wmp_emin(int i_nuclide);
extern "C" double nuclide_wmp_emax(int i_nuclide);
extern "C" void nuclide_calculate_urr_xs(bool use_mp, int i_nuclide,
int i_temp, double E);

42
include/openmc/wmp.h
View file

@ -4,10 +4,31 @@
#include "hdf5.h"
#include "xtensor/xtensor.hpp"
#include <complex>
#include <string>
#include <tuple>
namespace openmc {
//========================================================================
// Constants
//========================================================================
// Constants that determine which value to access
constexpr int MP_EA {0}; // Pole
constexpr int MP_RS {1}; // Residue scattering
constexpr int MP_RA {2}; // Residue absorption
constexpr int MP_RF {3}; // Residue fission
// Polynomial fit indices
constexpr int FIT_S {0}; // Scattering
constexpr int FIT_A {1}; // Absorption
constexpr int FIT_F {2}; // Fission
//========================================================================
// Windowed multipole data
//========================================================================
class WindowedMultipole {
public:
// Types, aliases
@ -16,6 +37,25 @@ public:
// Constructors, destructors
WindowedMultipole(hid_t group);
// Methods
//! Evaluate the windowed multipole equations for cross sections in the
//! resolved resonance regions
//!
//! \param E Incident neutron energy in [eV]
//! \param sqrtkT Square root of temperature times Boltzmann constant
//! \return Tuple of elastic scattering, absorption, and fission cross sections in [b]
std::tuple<double, double, double> evaluate(double E, double sqrtkT);
//! Evaluates the windowed multipole equations for the derivative of cross
//! sections in the resolved resonance regions with respect to temperature.
//!
//! \param E Incident neutron energy in [eV]
//! \param sqrtkT Square root of temperature times Boltzmann constant
//! \return Tuple of derivatives of elastic scattering, absorption, and
//! fission cross sections in [b/K]
std::tuple<double, double, double> evaluate_deriv(double E, double sqrtkT);
// Data members
std::string name_; //!< Name of nuclide
bool fissionable_; //!< Is the nuclide fissionable?
@ -26,7 +66,7 @@ public:
double spacing_; //!< Spacing in sqrt(E) space
int fit_order_; //!< Order of the fit
xt::xtensor<int, 2> windows_; //!< Indices of pole at start/end of window
xt::xtensor<double, 2> curvefit_; //!< Fitting function (reaction, coeff index, window index)
xt::xtensor<double, 3> curvefit_; //!< Fitting function (reaction, coeff index, window index)
xt::xtensor<bool, 1> broaden_poly_; //!< Whether to broaden curvefit
};

8
openmc/capi/math.py
View file

@ -33,9 +33,9 @@ _dll.maxwell_spectrum.argtypes = [c_double]
_dll.watt_spectrum.restype = c_double
_dll.watt_spectrum.argtypes = [c_double, c_double]
_dll.broaden_wmp_polynomials_c.restype = None
_dll.broaden_wmp_polynomials_c.argtypes = [c_double, c_double, c_int,
ndpointer(c_double)]
_dll.broaden_wmp_polynomials.restype = None
_dll.broaden_wmp_polynomials.argtypes = [c_double, c_double, c_int,
ndpointer(c_double)]
def t_percentile(p, df):
@ -274,5 +274,5 @@ def broaden_wmp_polynomials(E, dopp, n):
"""
factors = np.zeros(n, dtype=np.float64)
_dll.broaden_wmp_polynomials_c(E, dopp, n, factors)
_dll.broaden_wmp_polynomials(E, dopp, n, factors)
return factors

15
src/input_xml.F90
View file

@ -2256,6 +2256,14 @@ contains
integer(HID_T) :: file_id
integer(HID_T) :: group_id
interface
subroutine nuclide_load_multipole(ptr, group) bind(C)
import C_PTR, HID_T
type(C_PTR), value :: ptr
integer(HID_T), value :: group
end subroutine
end interface
associate (nuc => nuclides(i_table))
! Look for WMP data in cross_sections.xml
@ -2286,12 +2294,11 @@ contains
! Read nuclide data from HDF5
group_id = open_group(file_id, nuc % name)
allocate(nuc % multipole)
call nuc % multipole % from_hdf5(group_id)
nuc % mp_present = .true.
! Close the group and file.
call nuclide_load_multipole(nuc % ptr, group_id)
call close_group(group_id)
! Close the file
call file_close(file_id)
end associate

71
src/math.F90
View file

@ -11,9 +11,6 @@ module math
public :: calc_pn
public :: calc_rn
public :: rotate_angle
public :: faddeeva
public :: w_derivative
public :: broaden_wmp_polynomials
public :: spline
public :: spline_interpolate
public :: spline_integrate
@ -53,16 +50,6 @@ module math
real(C_DOUBLE), optional, intent(in) :: phi
end subroutine rotate_angle_c_intfc
subroutine broaden_wmp_polynomials(E, dopp, n, factors) &
bind(C, name='broaden_wmp_polynomials_c')
use ISO_C_BINDING
implicit none
real(C_DOUBLE), value, intent(in) :: E
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
subroutine spline(n, x, y, z) bind(C, name='spline_c')
use ISO_C_BINDING
implicit none
@ -96,14 +83,6 @@ module math
real(C_DOUBLE), value, intent(in) :: xb
real(C_DOUBLE) :: s
end function spline_integrate
function faddeeva_w(z, relerr) bind(C, name='Faddeeva_w') result(w)
use ISO_C_BINDING
implicit none
complex(C_DOUBLE_COMPLEX), value :: z
real(C_DOUBLE), value :: relerr
complex(C_DOUBLE_COMPLEX) :: w
end function faddeeva_w
end interface
contains
@ -126,54 +105,4 @@ contains
end function rotate_angle
!===============================================================================
! FADDEEVA the Faddeeva function, using Stephen Johnson's implementation
!===============================================================================
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
! Technically, the value we want is given by the equation:
! w(z) = I/Pi * Integrate[Exp[-t^2]/(z-t), {t, -Infinity, Infinity}]
! as shown in Equation 63 from Hwang, R. N. "A rigorous pole
! representation of multilevel cross sections and its practical
! applications." Nuclear Science and Engineering 96.3 (1987): 192-209.
!
! The MIT Faddeeva function evaluates w(z) = exp(-z^2)erfc(-iz). These
! two forms of the Faddeeva function are related by a transformation.
!
! If we call the integral form w_int, and the function form w_fun:
! For imag(z) > 0, w_int(z) = w_fun(z)
! For imag(z) < 0, w_int(z) = -conjg(w_fun(conjg(z)))
! Note that faddeeva_w will interpret zero as machine epsilon
relerr = ZERO
if (aimag(z) > ZERO) then
wv = faddeeva_w(z, relerr)
else
wv = -conjg(faddeeva_w(conjg(z), relerr))
end if
end function faddeeva
recursive function w_derivative(z, order) result(wv) bind(C)
complex(C_DOUBLE_COMPLEX), intent(in) :: z ! The point to evaluate Z at
integer(C_INT), intent(in) :: order
complex(C_DOUBLE_COMPLEX) :: wv ! The resulting w(z) value
select case(order)
case (0)
wv = faddeeva(z)
case (1)
wv = -TWO * z * faddeeva(z) + TWO * ONEI / SQRT_PI
case default
wv = -TWO * z * w_derivative(z, order-1) &
- TWO * (order-1) * w_derivative(z, order-2)
end select
end function w_derivative
end module math

39
src/math_functions.cpp
View file

@ -1,5 +1,7 @@
#include "openmc/math_functions.h"
#include "Faddeeva.hh"
namespace openmc {
//==============================================================================
@ -682,7 +684,8 @@ double watt_spectrum(double a, double b) {
}
void broaden_wmp_polynomials_c(double E, double dopp, int n, double factors[]) {
void broaden_wmp_polynomials(double E, double dopp, int n, double factors[])
{
// Factors is already pre-allocated
double sqrtE = std::sqrt(E);
double beta = sqrtE * dopp;
@ -817,4 +820,38 @@ double spline_integrate_c(int n, const double x[], const double y[],
return s;
}
std::complex<double> faddeeva(std::complex<double> z)
{
// Technically, the value we want is given by the equation:
// w(z) = I/pi * Integrate[Exp[-t^2]/(z-t), {t, -Infinity, Infinity}]
// as shown in Equation 63 from Hwang, R. N. "A rigorous pole
// representation of multilevel cross sections and its practical
// applications." Nucl. Sci. Eng. 96.3 (1987): 192-209.
//
// The MIT Faddeeva function evaluates w(z) = exp(-z^2)erfc(-iz). These
// two forms of the Faddeeva function are related by a transformation.
//
// If we call the integral form w_int, and the function form w_fun:
// For imag(z) > 0, w_int(z) = w_fun(z)
// For imag(z) < 0, w_int(z) = -conjg(w_fun(conjg(z)))
// Note that Faddeeva::w will interpret zero as machine epsilon
return z.imag() > 0.0 ? Faddeeva::w(z) :
-std::conj(Faddeeva::w(std::conj(z)));
}
std::complex<double> w_derivative(std::complex<double> z, int order)
{
using namespace std::complex_literals;
switch (order) {
case 0:
return faddeeva(z);
case 1:
return -2.0*z*faddeeva(z) + 2.0i / SQRT_PI;
default:
return -2.0*z*w_derivative(z, order-1)
- 2.0*(order-1)*w_derivative(z, order-2);
}
}
} // namespace openmc

121
src/multipole_header.F90
View file

@ -7,129 +7,8 @@ module multipole_header
implicit none
!========================================================================
! Multipole related constants
! Constants that determine which value to access
integer, parameter :: MP_EA = 1 ! Pole
! Residue indices
integer, parameter :: MP_RS = 2, & ! Residue scattering
MP_RA = 3, & ! Residue absorption
MP_RF = 4 ! Residue fission
! Polynomial fit indices
integer, parameter :: FIT_S = 1, & ! Scattering
FIT_A = 2, & ! Absorption
FIT_F = 3 ! Fission
!===============================================================================
! MULTIPOLE contains all the components needed for the windowed multipole
! temperature dependent cross section libraries for the resolved resonance
! region.
!===============================================================================
type MultipoleArray
character(20) :: name ! name of nuclide, e.g. U235
! Isotope Properties
logical :: fissionable ! Is this isotope fissionable?
complex(8), allocatable :: data(:,:) ! Poles and residues
real(8) :: sqrtAWR ! Square root of the atomic
! weight ratio
! Windows
integer :: fit_order ! Order of the fit. 1 linear,
! 2 quadratic, etc.
real(8) :: E_min ! Start energy for the windows
real(8) :: E_max ! End energy for the windows
real(8) :: spacing ! The actual spacing in sqrt(E)
! space.
integer, allocatable :: windows(:, :) ! Contains the indexes of the poles
! at the start and end of the
! window
real(8), allocatable :: curvefit(:,:,:) ! Contains the fitting function.
! (reaction type, coeff index,
! window index)
integer, allocatable :: broaden_poly(:) ! if 1, broaden, if 0, don't.
contains
procedure :: from_hdf5 => multipole_from_hdf5
end type MultipoleArray
contains
!===============================================================================
! FROM_HDF5 loads multipole data from a HDF5 file.
!===============================================================================
subroutine multipole_from_hdf5(this, group_id)
class(MultipoleArray), intent(inout) :: this
integer(HID_T), intent(in) :: group_id
integer :: n_poles, n_residues, n_windows
integer(HSIZE_T) :: dims_1d(1), dims_2d(2), dims_3d(3)
integer(HID_T) :: dset
! Get name of nuclide from group
this % name = get_name(group_id)
! Get rid of leading '/'
this % name = trim(this % name(2:))
! Read scalar values.
call read_dataset(this % spacing, group_id, "spacing")
call read_dataset(this % sqrtAWR, group_id, "sqrtAWR")
call read_dataset(this % E_min, group_id, "E_min")
call read_dataset(this % E_max, group_id, "E_max")
! Read the "data" array. Use its shape to figure out the number of poles
! and residue types in this data.
dset = open_dataset(group_id, "data")
call get_shape(dset, dims_2d)
n_residues = int(dims_2d(1), 4) - 1
n_poles = int(dims_2d(2), 4)
allocate(this % data(n_residues+1, n_poles))
if (n_poles > 0) call read_dataset(this % data, dset)
call close_dataset(dset)
! Check to see if this data includes fission residues.
this % fissionable = (n_residues == 3)
! Read the "windows" array and use its shape to figure out the number of
! windows.
dset = open_dataset(group_id, "windows")
call get_shape(dset, dims_2d)
n_windows = int(dims_2d(2), 4)
allocate(this % windows(dims_2d(1), n_windows))
call read_dataset(this % windows, dset)
call close_dataset(dset)
! Read the "broaden_poly" arrays.
dset = open_dataset(group_id, "broaden_poly")
call get_shape(dset, dims_1d)
if (dims_1d(1) /= n_windows) call fatal_error("broaden_poly array shape is&
&not consistent with the windows array shape in WMP library for"&
// trim(this % name) // ".")
allocate(this % broaden_poly(n_windows))
call read_dataset(this % broaden_poly, dset)
call close_dataset(dset)
! Read the "curvefit" array.
dset = open_dataset(group_id, "curvefit")
call get_shape(dset, dims_3d)
if (dims_3d(3) /= n_windows) call fatal_error("curvefit array shape is not&
&consistent with the windows array shape in WMP library for"&
// trim(this % name) // ".")
allocate(this % curvefit(dims_3d(1), dims_3d(2), dims_3d(3)))
call read_dataset(this % curvefit, dset)
call close_dataset(dset)
this % fit_order = int(dims_3d(2), 4) - 1
end subroutine multipole_from_hdf5
!===============================================================================
! CHECK_WMP_VERSION checks for the right version of WMP data within HDF5
! files

24
src/nuclide.cpp
View file

@ -509,9 +509,8 @@ void Nuclide::calculate_xs(int i_sab, double E, int i_log_union,
// Check to see if there is multipole data present at this energy
bool use_mp = false;
if (nuclide_wmp_present(i_nuclide_)) {
use_mp = (E >= nuclide_wmp_emin(i_nuclide_) &&
E <= nuclide_wmp_emax(i_nuclide_));
if (multipole_) {
use_mp = (E >= multipole_->E_min_ && E <= multipole_->E_max_);
}
int i_temp = -1;
@ -520,7 +519,7 @@ void Nuclide::calculate_xs(int i_sab, double E, int i_log_union,
if (use_mp) {
// Call multipole kernel
double sig_s, sig_a, sig_f;
nuclide_multipole_eval(i_nuclide_, E, sqrtkT, &sig_s, &sig_a, &sig_f);
std::tie(sig_s, sig_a, sig_f) = multipole_->evaluate(E, sqrtkT);
micro_xs.total = sig_s + sig_a;
micro_xs.elastic = sig_s;
@ -900,6 +899,23 @@ extern "C" void nuclide_calculate_elastic_xs_c(Nuclide* nuc)
nuc->calculate_elastic_xs();
}
extern "C" void nuclide_load_multipole(Nuclide* nuc, hid_t group)
{
nuc->multipole_ = std::make_unique<WindowedMultipole>(group);
}
extern "C" void multipole_deriv_eval(Nuclide* nuc, double E, double sqrtkT,
double* sig_s, double* sig_a, double* sig_f)
{
std::tie(*sig_s, *sig_a, *sig_f) = nuc->multipole_->evaluate_deriv(E, sqrtkT);
}
extern "C" bool multipole_in_range(Nuclide* nuc, double E)
{
return nuc->multipole_ && E >= nuc->multipole_->E_min_&&
E <= nuc->multipole_->E_max_;
}
extern "C" void nuclides_clear() { data::nuclides.clear(); }
extern "C" NuclideMicroXS* micro_xs_ptr();

257
src/nuclide_header.F90
View file

@ -11,10 +11,6 @@ module nuclide_header
use endf_header, only: Function1D, Polynomial, Tabulated1D
use error
use hdf5_interface
use math, only: faddeeva, w_derivative, &
broaden_wmp_polynomials
use multipole_header, only: MP_EA, MP_RS, MP_RA, MP_RF, &
FIT_S, FIT_A, FIT_F, MultipoleArray
use message_passing
use random_lcg, only: prn, future_prn, prn_set_stream
use reaction_header, only: Reaction
@ -73,8 +69,7 @@ module nuclide_header
class(Function1D), allocatable :: total_nu
! Multipole data
logical :: mp_present = .false.
type(MultipoleArray), pointer :: multipole => null()
logical :: mp_present = .false.
! Reactions
type(Reaction), allocatable :: reactions(:)
@ -90,7 +85,6 @@ module nuclide_header
type(C_PTR) :: ptr
contains
procedure :: clear => nuclide_clear
procedure :: from_hdf5 => nuclide_from_hdf5
procedure :: init_grid => nuclide_init_grid
procedure :: nu => nuclide_nu
@ -189,6 +183,23 @@ module nuclide_header
character(kind=C_CHAR), intent(in) :: name(*)
type(C_PTR) :: path
end function
subroutine multipole_deriv_eval(ptr, E, sqrtkT, sig_s, sig_a, sig_f) bind(C)
import C_PTR, C_DOUBLE
type(C_PTR), value :: ptr
real(C_DOUBLE), value :: E
real(C_DOUBLE), value :: sqrtkT
real(C_DOUBLE), intent(out) :: sig_s
real(C_DOUBLE), intent(out) :: sig_a
real(C_DOUBLE), intent(out) :: sig_f
end subroutine
function multipole_in_range(ptr, E) result(b) bind(C)
import C_PTR, C_DOUBLE, C_BOOL
type(C_PTR), value :: ptr
real(C_DOUBLE), value :: E
logical(C_BOOL) :: b
end function
end interface
contains
@ -219,17 +230,6 @@ contains
ptr = C_LOC(micro_xs(1))
end function
!===============================================================================
! NUCLIDE_CLEAR resets and deallocates data in Nuclide
!===============================================================================
subroutine nuclide_clear(this)
class(Nuclide), intent(inout) :: this ! The Nuclide object to clear
if (associated(this % multipole)) deallocate(this % multipole)
end subroutine nuclide_clear
subroutine nuclide_from_hdf5(this, group_id, temperature, method, tolerance, &
minmax, master, i_nuclide)
class(Nuclide), intent(inout) :: this
@ -731,192 +731,6 @@ contains
call nuclide_calculate_elastic_xs_c(this % ptr)
end subroutine nuclide_calculate_elastic_xs
!===============================================================================
! MULTIPOLE_EVAL evaluates the windowed multipole equations for cross
! sections in the resolved resonance regions
!===============================================================================
subroutine multipole_eval(multipole, E, sqrtkT, sig_s, sig_a, sig_f)
type(MultipoleArray), intent(in) :: multipole ! The windowed multipole
! object to process.
real(8), intent(in) :: E ! The energy at which to
! evaluate the cross section
real(8), intent(in) :: sqrtkT ! The temperature in the form
! sqrt(kT), at which
! to evaluate the XS.
real(8), intent(out) :: sig_s ! Scattering cross section
real(8), intent(out) :: sig_a ! Absorption cross section
real(8), intent(out) :: sig_f ! Fission cross section
complex(8) :: psi_chi ! The value of the psi-chi function for the
! asymptotic form
complex(8) :: c_temp ! complex temporary variable
complex(8) :: w_val ! The faddeeva function evaluated at Z
complex(8) :: Z ! sqrt(atomic weight ratio / kT) * (sqrt(E) - pole)
real(8) :: broadened_polynomials(multipole % fit_order + 1)
real(8) :: sqrtE ! sqrt(E), eV
real(8) :: invE ! 1/E, eV
real(8) :: dopp ! sqrt(atomic weight ratio / kT) = 1 / (2 sqrt(xi))
real(8) :: temp ! real temporary value
integer :: i_pole ! index of pole
integer :: i_poly ! index of curvefit
integer :: i_window ! index of window
integer :: startw ! window start pointer (for poles)
integer :: endw ! window end pointer
! ==========================================================================
! Bookkeeping
! Define some frequently used variables.
sqrtE = sqrt(E)
invE = ONE / E
! Locate us.
i_window = floor((sqrtE - sqrt(multipole % E_min)) / multipole % spacing &
+ ONE)
startw = multipole % windows(1, i_window)
endw = multipole % windows(2, i_window)
! Initialize the ouptut cross sections.
sig_s = ZERO
sig_a = ZERO
sig_f = ZERO
! ==========================================================================
! Add the contribution from the curvefit polynomial.
if (sqrtkT /= ZERO .and. multipole % broaden_poly(i_window) == 1) then
! Broaden the curvefit.
dopp = multipole % sqrtAWR / sqrtkT
call broaden_wmp_polynomials(E, dopp, multipole % fit_order + 1, &
broadened_polynomials)
do i_poly = 1, multipole % fit_order+1
sig_s = sig_s + multipole % curvefit(FIT_S, i_poly, i_window) &
* broadened_polynomials(i_poly)
sig_a = sig_a + multipole % curvefit(FIT_A, i_poly, i_window) &
* broadened_polynomials(i_poly)
if (multipole % fissionable) then
sig_f = sig_f + multipole % curvefit(FIT_F, i_poly, i_window) &
* broadened_polynomials(i_poly)
end if
end do
else ! Evaluate as if it were a polynomial
temp = invE
do i_poly = 1, multipole % fit_order+1
sig_s = sig_s + multipole % curvefit(FIT_S, i_poly, i_window) * temp
sig_a = sig_a + multipole % curvefit(FIT_A, i_poly, i_window) * temp
if (multipole % fissionable) then
sig_f = sig_f + multipole % curvefit(FIT_F, i_poly, i_window) * temp
end if
temp = temp * sqrtE
end do
end if
! ==========================================================================
! Add the contribution from the poles in this window.
if (sqrtkT == ZERO) then
! If at 0K, use asymptotic form.
do i_pole = startw, endw
psi_chi = -ONEI / (multipole % data(MP_EA, i_pole) - sqrtE)
c_temp = psi_chi / E
sig_s = sig_s + real(multipole % data(MP_RS, i_pole) * c_temp)
sig_a = sig_a + real(multipole % data(MP_RA, i_pole) * c_temp)
if (multipole % fissionable) then
sig_f = sig_f + real(multipole % data(MP_RF, i_pole) * c_temp)
end if
end do
else
! At temperature, use Faddeeva function-based form.
dopp = multipole % sqrtAWR / sqrtkT
if (endw >= startw) then
do i_pole = startw, endw
Z = (sqrtE - multipole % data(MP_EA, i_pole)) * dopp
w_val = faddeeva(Z) * dopp * invE * SQRT_PI
sig_s = sig_s + real(multipole % data(MP_RS, i_pole) * w_val)
sig_a = sig_a + real(multipole % data(MP_RA, i_pole) * w_val)
if (multipole % fissionable) then
sig_f = sig_f + real(multipole % data(MP_RF, i_pole) * w_val)
end if
end do
end if
end if
end subroutine multipole_eval
!===============================================================================
! MULTIPOLE_DERIV_EVAL evaluates the windowed multipole equations for the
! derivative of cross sections in the resolved resonance regions with respect to
! temperature.
!===============================================================================
subroutine multipole_deriv_eval(multipole, E, sqrtkT, sig_s, sig_a, sig_f)
type(MultipoleArray), intent(in) :: multipole ! The windowed multipole
! object to process.
real(8), intent(in) :: E ! The energy at which to
! evaluate the cross section
real(8), intent(in) :: sqrtkT ! The temperature in the form
! sqrt(kT), at which to
! evaluate the XS.
real(8), intent(out) :: sig_s ! Scattering cross section
real(8), intent(out) :: sig_a ! Absorption cross section
real(8), intent(out) :: sig_f ! Fission cross section
complex(8) :: w_val ! The faddeeva function evaluated at Z
complex(8) :: Z ! sqrt(atomic weight ratio / kT) * (sqrt(E) - pole)
real(8) :: sqrtE ! sqrt(E), eV
real(8) :: invE ! 1/E, eV
real(8) :: dopp ! sqrt(atomic weight ratio / kT)
integer :: i_pole ! index of pole
integer :: i_window ! index of window
integer :: startw ! window start pointer (for poles)
integer :: endw ! window end pointer
real(8) :: T
! ==========================================================================
! Bookkeeping
! Define some frequently used variables.
sqrtE = sqrt(E)
invE = ONE / E
T = sqrtkT**2 / K_BOLTZMANN
if (sqrtkT == ZERO) call fatal_error("Windowed multipole temperature &
&derivatives are not implemented for 0 Kelvin cross sections.")
! Locate us
i_window = floor((sqrtE - sqrt(multipole % E_min)) / multipole % spacing &
+ ONE)
startw = multipole % windows(1, i_window)
endw = multipole % windows(2, i_window)
! Initialize the ouptut cross sections.
sig_s = ZERO
sig_a = ZERO
sig_f = ZERO
! TODO Polynomials: Some of the curvefit polynomials Doppler broaden so
! rigorously we should be computing the derivative of those. But in
! practice, those derivatives are only large at very low energy and they
! have no effect on reactor calculations.
! ==========================================================================
! Add the contribution from the poles in this window.
dopp = multipole % sqrtAWR / sqrtkT
if (endw >= startw) then
do i_pole = startw, endw
Z = (sqrtE - multipole % data(MP_EA, i_pole)) * dopp
w_val = -invE * SQRT_PI * HALF * w_derivative(Z, 2)
sig_s = sig_s + real(multipole % data(MP_RS, i_pole) * w_val)
sig_a = sig_a + real(multipole % data(MP_RA, i_pole) * w_val)
if (multipole % fissionable) then
sig_f = sig_f + real(multipole % data(MP_RF, i_pole) * w_val)
end if
end do
sig_s = -HALF*multipole % sqrtAWR / sqrt(K_BOLTZMANN) * T**(-1.5) * sig_s
sig_a = -HALF*multipole % sqrtAWR / sqrt(K_BOLTZMANN) * T**(-1.5) * sig_a
sig_f = -HALF*multipole % sqrtAWR / sqrt(K_BOLTZMANN) * T**(-1.5) * sig_f
end if
end subroutine multipole_deriv_eval
!===============================================================================
! CHECK_DATA_VERSION checks for the right version of nuclear data within HDF5
! files
@ -947,8 +761,6 @@ contains
!===============================================================================
subroutine free_memory_nuclide()
integer :: i
interface
subroutine library_clear() bind(C)
end subroutine
@ -960,9 +772,6 @@ contains
! Deallocate cross section data, listings, and cache
if (allocated(nuclides)) then
! First call the clear routines
do i = 1, size(nuclides)
call nuclides(i) % clear()
end do
deallocate(nuclides)
call nuclides_clear()
end if
@ -1085,34 +894,4 @@ contains
end if
end function openmc_nuclide_name
function nuclide_wmp_present(i_nuclide) result(b) bind(C)
integer(C_INT), value :: i_nuclide
logical(C_BOOL) :: b
b = nuclides(i_nuclide + 1) % mp_present
end function
function nuclide_wmp_emin(i_nuclide) result(E) bind(C)
integer(C_INT), value :: i_nuclide
real(C_DOUBLE) :: E
E = nuclides(i_nuclide + 1) % multipole % E_min
end function
function nuclide_wmp_emax(i_nuclide) result(E) bind(C)
integer(C_INT), value :: i_nuclide
real(C_DOUBLE) :: E
E = nuclides(i_nuclide + 1) % multipole % E_max
end function
subroutine nuclide_multipole_eval(i_nuclide, E, sqrtkT, sig_s, sig_a, sig_f) bind(C)
integer(C_INT), value :: i_nuclide
real(C_DOUBLE), value :: E
real(C_DOUBLE), value :: sqrtkT
real(C_DOUBLE), intent(out) :: sig_s
real(C_DOUBLE), intent(out) :: sig_a
real(C_DOUBLE), intent(out) :: sig_f
call multipole_eval(nuclides(i_nuclide + 1) % multipole, E, sqrtkT, &
sig_s, sig_a, sig_f)
end subroutine
end module nuclide_header

7
src/physics.cpp
View file

@ -511,8 +511,8 @@ Reaction* sample_fission(int i_nuclide, double E)
// Check to see if we are in a windowed multipole range. WMP only supports
// the first fission reaction.
if (nuclide_wmp_present(i_nuclide)) {
if (E >= nuclide_wmp_emin(i_nuclide) && E <= nuclide_wmp_emax(i_nuclide)) {
if (nuc->multipole_) {
if (E >= nuc->multipole_->E_min_ && E <= nuc->multipole_->E_max_) {
return nuc->fission_rx_[0];
}
}
@ -635,8 +635,7 @@ void scatter(Particle* p, int i_nuclide, int i_nuc_mat)
// NON-S(A,B) ELASTIC SCATTERING
// Determine temperature
double kT = nuclide_wmp_present(i_nuclide) ?
p->sqrtkT*p->sqrtkT : nuc->kTs_[i_temp];
double kT = nuc->multipole_ ? p->sqrtkT*p->sqrtkT : nuc->kTs_[i_temp];
// Perform collision physics for elastic scattering
elastic_scatter(i_nuclide, nuc->reactions_[0].get(), kT,

102
src/tallies/tally.F90
View file

@ -3248,10 +3248,8 @@ contains
dsig_s = ZERO
dsig_a = ZERO
associate (nuc => nuclides(p % event_nuclide))
if (nuc % mp_present .and. &
p % last_E >= nuc % multipole % E_min .and. &
p % last_E <= nuc % multipole % E_max) then
call multipole_deriv_eval(nuc % multipole, p % last_E, &
if (multipole_in_range(nuc % ptr, p % last_E)) then
call multipole_deriv_eval(nuc % ptr, p % last_E, &
p % sqrtkT, dsig_s, dsig_a, dsig_f)
end if
end associate
@ -3275,10 +3273,8 @@ contains
dsig_s = ZERO
associate (nuc => nuclides(p % event_nuclide))
if (nuc % mp_present .and. &
p % last_E >= nuc % multipole % E_min .and. &
p % last_E <= nuc % multipole % E_max) then
call multipole_deriv_eval(nuc % multipole, p % last_E, &
if (multipole_in_range(nuc % ptr, p % last_E)) then
call multipole_deriv_eval(nuc % ptr, p % last_E, &
p % sqrtkT, dsig_s, dsig_a, dsig_f)
end if
end associate
@ -3300,10 +3296,8 @@ contains
dsig_a = ZERO
associate (nuc => nuclides(p % event_nuclide))
if (nuc % mp_present .and. &
p % last_E >= nuc % multipole % E_min .and. &
p % last_E <= nuc % multipole % E_max) then
call multipole_deriv_eval(nuc % multipole, p % last_E, &
if (multipole_in_range(nuc % ptr, p % last_E)) then
call multipole_deriv_eval(nuc % ptr, p % last_E, &
p % sqrtkT, dsig_s, dsig_a, dsig_f)
end if
end associate
@ -3325,10 +3319,8 @@ contains
dsig_f = ZERO
associate (nuc => nuclides(p % event_nuclide))
if (nuc % mp_present .and. &
p % last_E >= nuc % multipole % E_min .and. &
p % last_E <= nuc % multipole % E_max) then
call multipole_deriv_eval(nuc % multipole, p % last_E, &
if (multipole_in_range(nuc % ptr, p % last_E)) then
call multipole_deriv_eval(nuc % ptr, p % last_E, &
p % sqrtkT, dsig_s, dsig_a, dsig_f)
end if
end associate
@ -3350,10 +3342,8 @@ contains
dsig_f = ZERO
associate (nuc => nuclides(p % event_nuclide))
if (nuc % mp_present .and. &
p % last_E >= nuc % multipole % E_min .and. &
p % last_E <= nuc % multipole % E_max) then
call multipole_deriv_eval(nuc % multipole, p % last_E, &
if (multipole_in_range(nuc % ptr, p % last_E)) then
call multipole_deriv_eval(nuc % ptr, p % last_E, &
p % sqrtkT, dsig_s, dsig_a, dsig_f)
end if
end associate
@ -3386,11 +3376,9 @@ contains
associate(mat => materials(p % material))
do l = 1, mat % n_nuclides
associate (nuc => nuclides(mat % nuclide(l)))
if (nuc % mp_present .and. &
p % last_E >= nuc % multipole % E_min .and. &
p % last_E <= nuc % multipole % E_max .and. &
if (multipole_in_range(nuc % ptr, p % last_E) .and. &
micro_xs(mat % nuclide(l)) % total > ZERO) then
call multipole_deriv_eval(nuc % multipole, p % last_E, &
call multipole_deriv_eval(nuc % ptr, p % last_E, &
p % sqrtkT, dsig_s, dsig_a, dsig_f)
cum_dsig = cum_dsig + (dsig_s + dsig_a) &
* mat % atom_density(l)
@ -3405,10 +3393,8 @@ contains
dsig_s = ZERO
dsig_a = ZERO
associate (nuc => nuclides(i_nuclide))
if (nuc % mp_present .and. &
p % last_E >= nuc % multipole % E_min .and. &
p % last_E <= nuc % multipole % E_max) then
call multipole_deriv_eval(nuc % multipole, p % last_E, &
if (multipole_in_range(nuc % ptr, p % last_E)) then
call multipole_deriv_eval(nuc % ptr, p % last_E, &
p % sqrtkT, dsig_s, dsig_a, dsig_f)
end if
end associate
@ -3426,12 +3412,10 @@ contains
associate(mat => materials(p % material))
do l = 1, mat % n_nuclides
associate (nuc => nuclides(mat % nuclide(l)))
if (nuc % mp_present .and. &
p % last_E >= nuc % multipole % E_min .and. &
p % last_E <= nuc % multipole % E_max .and. &
if (multipole_in_range(nuc % ptr, p % last_E) .and. &
(micro_xs(mat % nuclide(l)) % total &
- micro_xs(mat % nuclide(l)) % absorption) > ZERO) then
call multipole_deriv_eval(nuc % multipole, p % last_E, &
call multipole_deriv_eval(nuc % ptr, p % last_E, &
p % sqrtkT, dsig_s, dsig_a, dsig_f)
cum_dsig = cum_dsig + dsig_s * mat % atom_density(l)
end if
@ -3445,10 +3429,8 @@ contains
then
dsig_s = ZERO
associate (nuc => nuclides(i_nuclide))
if (nuc % mp_present .and. &
p % last_E >= nuc % multipole % E_min .and. &
p % last_E <= nuc % multipole % E_max) then
call multipole_deriv_eval(nuc % multipole, p % last_E, &
if (multipole_in_range(nuc % ptr, p % last_E)) then
call multipole_deriv_eval(nuc % ptr, p % last_E, &
p % sqrtkT, dsig_s, dsig_a, dsig_f)
end if
end associate
@ -3467,11 +3449,9 @@ contains
associate(mat => materials(p % material))
do l = 1, mat % n_nuclides
associate (nuc => nuclides(mat % nuclide(l)))
if (nuc % mp_present .and. &
p % last_E >= nuc % multipole % E_min .and. &
p % last_E <= nuc % multipole % E_max .and. &
if (multipole_in_range(nuc % ptr, p % last_E) .and. &
micro_xs(mat % nuclide(l)) % absorption > ZERO) then
call multipole_deriv_eval(nuc % multipole, p % last_E, &
call multipole_deriv_eval(nuc % ptr, p % last_E, &
p % sqrtkT, dsig_s, dsig_a, dsig_f)
cum_dsig = cum_dsig + dsig_a * mat % atom_density(l)
end if
@ -3484,10 +3464,8 @@ contains
.and. material_xs % absorption > ZERO) then
dsig_a = ZERO
associate (nuc => nuclides(i_nuclide))
if (nuc % mp_present .and. &
p % last_E >= nuc % multipole % E_min .and. &
p % last_E <= nuc % multipole % E_max) then
call multipole_deriv_eval(nuc % multipole, p % last_E, &
if (multipole_in_range(nuc % ptr, p % last_E)) then
call multipole_deriv_eval(nuc % ptr, p % last_E, &
p % sqrtkT, dsig_s, dsig_a, dsig_f)
end if
end associate
@ -3505,11 +3483,9 @@ contains
associate(mat => materials(p % material))
do l = 1, mat % n_nuclides
associate (nuc => nuclides(mat % nuclide(l)))
if (nuc % mp_present .and. &
p % last_E >= nuc % multipole % E_min .and. &
p % last_E <= nuc % multipole % E_max .and. &
if (multipole_in_range(nuc % ptr, p % last_E) .and. &
micro_xs(mat % nuclide(l)) % fission > ZERO) then
call multipole_deriv_eval(nuc % multipole, p % last_E, &
call multipole_deriv_eval(nuc % ptr, p % last_E, &
p % sqrtkT, dsig_s, dsig_a, dsig_f)
cum_dsig = cum_dsig + dsig_f * mat % atom_density(l)
end if
@ -3522,10 +3498,8 @@ contains
.and. material_xs % fission > ZERO) then
dsig_f = ZERO
associate (nuc => nuclides(i_nuclide))
if (nuc % mp_present .and. &
p % last_E >= nuc % multipole % E_min .and. &
p % last_E <= nuc % multipole % E_max) then
call multipole_deriv_eval(nuc % multipole, p % last_E, &
if (multipole_in_range(nuc % ptr, p % last_E)) then
call multipole_deriv_eval(nuc % ptr, p % last_E, &
p % sqrtkT, dsig_s, dsig_a, dsig_f)
end if
end associate
@ -3543,11 +3517,9 @@ contains
associate(mat => materials(p % material))
do l = 1, mat % n_nuclides
associate (nuc => nuclides(mat % nuclide(l)))
if (nuc % mp_present .and. &
p % last_E >= nuc % multipole % E_min .and. &
p % last_E <= nuc % multipole % E_max .and. &
if (multipole_in_range(nuc % ptr, p % last_E) .and. &
micro_xs(mat % nuclide(l)) % nu_fission > ZERO) then
call multipole_deriv_eval(nuc % multipole, p % last_E, &
call multipole_deriv_eval(nuc % ptr, p % last_E, &
p % sqrtkT, dsig_s, dsig_a, dsig_f)
cum_dsig = cum_dsig + dsig_f * mat % atom_density(l) &
* micro_xs(mat % nuclide(l)) % nu_fission &
@ -3562,10 +3534,8 @@ contains
.and. material_xs % nu_fission > ZERO) then
dsig_f = ZERO
associate (nuc => nuclides(i_nuclide))
if (nuc % mp_present .and. &
p % last_E >= nuc % multipole % E_min .and. &
p % last_E <= nuc % multipole % E_max) then
call multipole_deriv_eval(nuc % multipole, p % last_E, &
if (multipole_in_range(nuc % ptr, p % last_E)) then
call multipole_deriv_eval(nuc % ptr, p % last_E, &
p % sqrtkT, dsig_s, dsig_a, dsig_f)
end if
end associate
@ -3634,13 +3604,11 @@ contains
if (mat % id() == deriv % diff_material) then
do l=1, mat % n_nuclides
associate (nuc => nuclides(mat % nuclide(l)))
if (nuc % mp_present .and. &
p % E >= nuc % multipole % E_min .and. &
p % E <= nuc % multipole % E_max) then
if (multipole_in_range(nuc % ptr, p % E)) then
! phi is proportional to e^(-Sigma_tot * dist)
! (1 / phi) * (d_phi / d_T) = - (d_Sigma_tot / d_T) * dist
! (1 / phi) * (d_phi / d_T) = - N (d_sigma_tot / d_T) * dist
call multipole_deriv_eval(nuc % multipole, p % E, &
call multipole_deriv_eval(nuc % ptr, p % E, &
p % sqrtkT, dsig_s, dsig_a, dsig_f)
deriv % flux_deriv = deriv % flux_deriv &
- distance * (dsig_s + dsig_a) * mat % atom_density(l)
@ -3721,13 +3689,11 @@ contains
do l=1, mat % n_nuclides
associate (nuc => nuclides(mat % nuclide(l)))
if (mat % nuclide(l) == p % event_nuclide .and. &
nuc % mp_present .and. &
p % last_E >= nuc % multipole % E_min .and. &
p % last_E <= nuc % multipole % E_max) then
multipole_in_range(nuc % ptr, p % last_E)) then
! phi is proportional to Sigma_s
! (1 / phi) * (d_phi / d_T) = (d_Sigma_s / d_T) / Sigma_s
! (1 / phi) * (d_phi / d_T) = (d_sigma_s / d_T) / sigma_s
call multipole_deriv_eval(nuc % multipole, p % last_E, &
call multipole_deriv_eval(nuc % ptr, p % last_E, &
p % sqrtkT, dsig_s, dsig_a, dsig_f)
deriv % flux_deriv = deriv % flux_deriv + dsig_s&
/ (micro_xs(mat % nuclide(l)) % total &

140
src/wmp.cpp
View file

@ -1,6 +1,10 @@
#include "openmc/wmp.h"
#include "openmc/constants.h"
#include "openmc/hdf5_interface.h"
#include "openmc/math_functions.h"
#include <cmath>
namespace openmc {
@ -44,4 +48,140 @@ WindowedMultipole::WindowedMultipole(hid_t group)
fit_order_ = curvefit_.shape()[1] - 1;
}
std::tuple<double, double, double>
WindowedMultipole::evaluate(double E, double sqrtkT)
{
using namespace std::complex_literals;
// ==========================================================================
// Bookkeeping
// Define some frequently used variables.
double sqrtE = std::sqrt(E);
double invE = 1.0 / E;
// Locate window containing energy
int i_window = (sqrtE - std::sqrt(E_min_)) / spacing_;
int startw = windows_(i_window, 0) - 1;
int endw = windows_(i_window, 1) - 1;
// Initialize the ouptut cross sections
double sig_s = 0.0;
double sig_a = 0.0;
double sig_f = 0.0;
// ==========================================================================
// Add the contribution from the curvefit polynomial.
if (sqrtkT > 0.0 && broaden_poly_(i_window)) {
// Broaden the curvefit.
double dopp = sqrt_awr_ / sqrtkT;
std::vector<double> broadened_polynomials(fit_order_ + 1);
broaden_wmp_polynomials(E, dopp, fit_order_ + 1, broadened_polynomials.data());
for (int i_poly = 0; i_poly < fit_order_ + 1; ++i_poly) {
sig_s += curvefit_(i_window, i_poly, FIT_S) * broadened_polynomials[i_poly];
sig_a += curvefit_(i_window, i_poly, FIT_A) * broadened_polynomials[i_poly];
if (fissionable_) {
sig_f += curvefit_(i_window, i_poly, FIT_F) * broadened_polynomials[i_poly];
}
}
} else {
// Evaluate as if it were a polynomial
double temp = invE;
for (int i_poly = 0; i_poly < fit_order_ + 1; ++i_poly) {
sig_s += curvefit_(i_window, i_poly, FIT_S) * temp;
sig_a += curvefit_(i_window, i_poly, FIT_A) * temp;
if (fissionable_) {
sig_f += curvefit_(i_window, i_poly, FIT_F) * temp;
}
temp *= sqrtE;
}
}
// ==========================================================================
// Add the contribution from the poles in this window.
if (sqrtkT == 0.0) {
// If at 0K, use asymptotic form.
for (int i_pole = startw; i_pole <= endw; ++i_pole) {
std::complex<double> psi_chi = -1.0i / (data_(i_pole, MP_EA) - sqrtE);
std::complex<double> c_temp = psi_chi / E;
sig_s += (data_(i_pole, MP_RS) * c_temp).real();
sig_a += (data_(i_pole, MP_RA) * c_temp).real();
if (fissionable_) {
sig_f += (data_(i_pole, MP_RF) * c_temp).real();
}
}
} else {
// At temperature, use Faddeeva function-based form.
double dopp = sqrt_awr_ / sqrtkT;
if (endw >= startw) {
for (int i_pole = startw; i_pole <= endw; ++i_pole) {
std::complex<double> z = (sqrtE - data_(i_pole, MP_EA)) * dopp;
std::complex<double> w_val = faddeeva(z) * dopp * invE * SQRT_PI;
sig_s += (data_(i_pole, MP_RS) * w_val).real();
sig_a += (data_(i_pole, MP_RA) * w_val).real();
if (fissionable_) {
sig_f += (data_(i_pole, MP_RF) * w_val).real();
}
}
}
}
return std::make_tuple(sig_s, sig_a, sig_f);
}
std::tuple<double, double, double>
WindowedMultipole::evaluate_deriv(double E, double sqrtkT)
{
// ==========================================================================
// Bookkeeping
// Define some frequently used variables.
double sqrtE = std::sqrt(E);
double invE = 1.0 / E;
double T = sqrtkT*sqrtkT / K_BOLTZMANN;
if (sqrtkT == 0.0) {
fatal_error("Windowed multipole temperature derivatives are not implemented"
"for 0 Kelvin cross sections.");
}
// Locate us
int i_window = (sqrtE - std::sqrt(E_min_)) / spacing_;
int startw = windows_(i_window, 0) - 1;
int endw = windows_(i_window, 1) - 1;
// Initialize the ouptut cross sections.
double sig_s = 0.0;
double sig_a = 0.0;
double sig_f = 0.0;
// TODO Polynomials: Some of the curvefit polynomials Doppler broaden so
// rigorously we should be computing the derivative of those. But in
// practice, those derivatives are only large at very low energy and they
// have no effect on reactor calculations.
// ==========================================================================
// Add the contribution from the poles in this window.
double dopp = sqrt_awr_ / sqrtkT;
if (endw >= startw) {
for (int i_pole = startw; i_pole <= endw; ++i_pole) {
std::complex<double> z = (sqrtE - data_(i_pole, MP_EA)) * dopp;
std::complex<double> w_val = -invE * SQRT_PI * 0.5 * w_derivative(z, 2);
sig_s += (data_(i_pole, MP_RS) * w_val).real();
sig_a += (data_(i_pole, MP_RA) * w_val).real();
if (fissionable_) {
sig_f += (data_(i_pole, MP_RF) * w_val).real();
}
}
sig_s *= -0.5*sqrt_awr_ / std::sqrt(K_BOLTZMANN) * std::pow(T, -1.5);
sig_a *= -0.5*sqrt_awr_ / std::sqrt(K_BOLTZMANN) * std::pow(T, -1.5);
sig_f *= -0.5*sqrt_awr_ / std::sqrt(K_BOLTZMANN) * std::pow(T, -1.5);
}
return std::make_tuple(sig_s, sig_a, sig_f);
}
} // namespace openmc

10
tests/regression_tests/diff_tally/results_true.dat
View file

@ -92,16 +92,16 @@ d_material,d_nuclide,d_variable,score,mean,std. dev.
1,,temperature,total,5.7228698e-05,1.7465295e-04
1,,temperature,absorption,2.9495471e-05,3.1637786e-05
1,,temperature,scatter,2.7733227e-05,1.4323068e-04
1,,temperature,fission,-5.6710689e-06,1.2800905e-05
1,,temperature,fission,-5.6710690e-06,1.2800905e-05
1,,temperature,nu-fission,-1.3819116e-05,3.1189136e-05
1,,temperature,total,-5.8315815e-06,1.8705738e-05
1,,temperature,absorption,-5.3204426e-06,1.6688104e-05
1,,temperature,scatter,-5.1113883e-07,2.0213875e-06
1,,temperature,fission,-5.6723577e-06,1.2800214e-05
1,,temperature,absorption,-5.3204427e-06,1.6688104e-05
1,,temperature,scatter,-5.1113882e-07,2.0213875e-06
1,,temperature,fission,-5.6723578e-06,1.2800214e-05
1,,temperature,nu-fission,-1.3822264e-05,3.1187458e-05
1,,temperature,total,-5.5283703e-04,7.6408422e-04
1,,temperature,absorption,-6.2179157e-06,1.0365194e-05
1,,temperature,scatter,-5.4661911e-04,7.5373133e-04
1,,temperature,scatter,-5.4661912e-04,7.5373133e-04
1,,temperature,fission,0.0000000e+00,0.0000000e+00
1,,temperature,nu-fission,0.0000000e+00,0.0000000e+00
1,,temperature,total,0.0000000e+00,0.0000000e+00

3
vendor/faddeeva/Faddeeva.c vendored
View file

@ -1,3 +0,0 @@
/* The Faddeeva.cc file contains macros to let it compile as C code
(assuming C99 complex-number support), so just #include it. */
#include "Faddeeva.cc"

52
vendor/faddeeva/Faddeeva.h → vendor/faddeeva/Faddeeva.hh vendored
View file

@ -1,5 +1,5 @@
/* Copyright (c) 2012 Massachusetts Institute of Technology
*
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
@ -7,62 +7,56 @@
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject to
* the following conditions:
*
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
*
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
* LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
* OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
/* Available at: http://ab-initio.mit.edu/Faddeeva
Header file for Faddeeva.c; see Faddeeva.cc for more information. */
Header file for Faddeeva.cc; see that file for more information. */
#ifndef FADDEEVA_H
#define FADDEEVA_H 1
#ifndef FADDEEVA_HH
#define FADDEEVA_HH 1
// Require C99 complex-number support
#include <complex.h>
#include <complex>
#ifdef __cplusplus
extern "C"
{
#endif /* __cplusplus */
namespace Faddeeva {
// compute w(z) = exp(-z^2) erfc(-iz) [ Faddeeva / scaled complex error func ]
extern double complex Faddeeva_w(double complex z,double relerr);
extern double Faddeeva_w_im(double x); // special-case code for Im[w(x)] of real x
extern std::complex<double> w(std::complex<double> z,double relerr=0);
extern double w_im(double x); // special-case code for Im[w(x)] of real x
// Various functions that we can compute with the help of w(z)
// compute erfcx(z) = exp(z^2) erfc(z)
extern double complex Faddeeva_erfcx(double complex z, double relerr);
extern double Faddeeva_erfcx_re(double x); // special case for real x
extern std::complex<double> erfcx(std::complex<double> z, double relerr=0);
extern double erfcx(double x); // special case for real x
// compute erf(z), the error function of complex arguments
extern double complex Faddeeva_erf(double complex z, double relerr);
extern double Faddeeva_erf_re(double x); // special case for real x
extern std::complex<double> erf(std::complex<double> z, double relerr=0);
extern double erf(double x); // special case for real x
// compute erfi(z) = -i erf(iz), the imaginary error function
extern double complex Faddeeva_erfi(double complex z, double relerr);
extern double Faddeeva_erfi_re(double x); // special case for real x
extern std::complex<double> erfi(std::complex<double> z, double relerr=0);
extern double erfi(double x); // special case for real x
// compute erfc(z) = 1 - erf(z), the complementary error function
extern double complex Faddeeva_erfc(double complex z, double relerr);
extern double Faddeeva_erfc_re(double x); // special case for real x
extern std::complex<double> erfc(std::complex<double> z, double relerr=0);
extern double erfc(double x); // special case for real x
// compute Dawson(z) = sqrt(pi)/2 * exp(-z^2) * erfi(z)
extern double complex Faddeeva_Dawson(double complex z, double relerr);
extern double Faddeeva_Dawson_re(double x); // special case for real x
extern std::complex<double> Dawson(std::complex<double> z, double relerr=0);
extern double Dawson(double x); // special case for real x
#ifdef __cplusplus
}
#endif /* __cplusplus */
} // namespace Faddeeva
#endif // FADDEEVA_H
#endif // FADDEEVA_HH