From 38ff673cf35768e59bf0648d25a1f90a8e7e3e51 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 28 Dec 2018 07:19:20 -0600 Subject: [PATCH] Finish converting windowed multipole functions to C++ --- CMakeLists.txt | 8 +- include/openmc/endf.h | 4 + include/openmc/math_functions.h | 122 +++++---- include/openmc/nuclide.h | 11 +- include/openmc/wmp.h | 42 ++- openmc/capi/math.py | 8 +- src/input_xml.F90 | 15 +- src/math.F90 | 71 ----- src/math_functions.cpp | 39 ++- src/multipole_header.F90 | 121 --------- src/nuclide.cpp | 24 +- src/nuclide_header.F90 | 257 ++---------------- src/physics.cpp | 7 +- src/tallies/tally.F90 | 102 +++---- src/wmp.cpp | 140 ++++++++++ .../diff_tally/results_true.dat | 10 +- vendor/faddeeva/Faddeeva.c | 3 - vendor/faddeeva/{Faddeeva.h => Faddeeva.hh} | 52 ++-- 18 files changed, 417 insertions(+), 619 deletions(-) delete mode 100644 vendor/faddeeva/Faddeeva.c rename vendor/faddeeva/{Faddeeva.h => Faddeeva.hh} (58%) diff --git a/CMakeLists.txt b/CMakeLists.txt index 73c2ebcf8d..7e9e24161c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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 diff --git a/include/openmc/endf.h b/include/openmc/endf.h index 170a647d50..ae480f8cfc 100644 --- a/include/openmc/endf.h +++ b/include/openmc/endf.h @@ -99,6 +99,10 @@ private: std::vector 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 read_function(hid_t group, const char* name); } // namespace openmc diff --git a/include/openmc/math_functions.h b/include/openmc/math_functions.h index ca0646961b..d143c81f1e 100644 --- a/include/openmc/math_functions.h +++ b/include/openmc/math_functions.h @@ -5,6 +5,7 @@ #define OPENMC_MATH_FUNCTIONS_H #include +#include #include #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 faddeeva(std::complex 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 w_derivative(std::complex z, int order); + } // namespace openmc #endif // OPENMC_MATH_FUNCTIONS_H diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h index a8178b11c3..e551aba613 100644 --- a/include/openmc/nuclide.h +++ b/include/openmc/nuclide.h @@ -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 grid_; //!< Energy grid at each temperature std::vector> xs_; //!< Cross sections at each temperature + // Multipole data + std::unique_ptr 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); diff --git a/include/openmc/wmp.h b/include/openmc/wmp.h index c8058d3694..cf537a61d3 100644 --- a/include/openmc/wmp.h +++ b/include/openmc/wmp.h @@ -4,10 +4,31 @@ #include "hdf5.h" #include "xtensor/xtensor.hpp" +#include #include +#include 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 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 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 windows_; //!< Indices of pole at start/end of window - xt::xtensor curvefit_; //!< Fitting function (reaction, coeff index, window index) + xt::xtensor curvefit_; //!< Fitting function (reaction, coeff index, window index) xt::xtensor broaden_poly_; //!< Whether to broaden curvefit }; diff --git a/openmc/capi/math.py b/openmc/capi/math.py index 34beb585fd..47c7cafcdd 100644 --- a/openmc/capi/math.py +++ b/openmc/capi/math.py @@ -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 diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 40d5ca8f9c..be93a164e0 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -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 diff --git a/src/math.F90 b/src/math.F90 index fe4156b9dc..5a1076d9d7 100644 --- a/src/math.F90 +++ b/src/math.F90 @@ -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 diff --git a/src/math_functions.cpp b/src/math_functions.cpp index a8646f09d9..05b6501499 100644 --- a/src/math_functions.cpp +++ b/src/math_functions.cpp @@ -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 faddeeva(std::complex 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 w_derivative(std::complex 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 diff --git a/src/multipole_header.F90 b/src/multipole_header.F90 index 359bb3ac14..f621b2e1c3 100644 --- a/src/multipole_header.F90 +++ b/src/multipole_header.F90 @@ -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& - ¬ 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 diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 53579cc857..95d5b2b745 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -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(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(); diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index cc6f1761c2..11ba75aa18 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -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 diff --git a/src/physics.cpp b/src/physics.cpp index 5c8d2ccf3c..4a303c3ce7 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -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, diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index e014e7402d..1398ec2537 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -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 & diff --git a/src/wmp.cpp b/src/wmp.cpp index 5ad4218c80..77454dd170 100644 --- a/src/wmp.cpp +++ b/src/wmp.cpp @@ -1,6 +1,10 @@ #include "openmc/wmp.h" +#include "openmc/constants.h" #include "openmc/hdf5_interface.h" +#include "openmc/math_functions.h" + +#include namespace openmc { @@ -44,4 +48,140 @@ WindowedMultipole::WindowedMultipole(hid_t group) fit_order_ = curvefit_.shape()[1] - 1; } +std::tuple +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 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 psi_chi = -1.0i / (data_(i_pole, MP_EA) - sqrtE); + std::complex 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 z = (sqrtE - data_(i_pole, MP_EA)) * dopp; + std::complex 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 +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 z = (sqrtE - data_(i_pole, MP_EA)) * dopp; + std::complex 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 diff --git a/tests/regression_tests/diff_tally/results_true.dat b/tests/regression_tests/diff_tally/results_true.dat index 41f3ac4f77..f4c0b3d221 100644 --- a/tests/regression_tests/diff_tally/results_true.dat +++ b/tests/regression_tests/diff_tally/results_true.dat @@ -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 diff --git a/vendor/faddeeva/Faddeeva.c b/vendor/faddeeva/Faddeeva.c deleted file mode 100644 index 78f5714efa..0000000000 --- a/vendor/faddeeva/Faddeeva.c +++ /dev/null @@ -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" diff --git a/vendor/faddeeva/Faddeeva.h b/vendor/faddeeva/Faddeeva.hh similarity index 58% rename from vendor/faddeeva/Faddeeva.h rename to vendor/faddeeva/Faddeeva.hh index 9e26bc1ed9..c4a2e9717f 100644 --- a/vendor/faddeeva/Faddeeva.h +++ b/vendor/faddeeva/Faddeeva.hh @@ -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 +#include -#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 w(std::complex 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 erfcx(std::complex 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 erf(std::complex 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 erfi(std::complex 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 erfc(std::complex 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 Dawson(std::complex 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