From 4fc79565150eabd407d7bfea19134bc052d1ba37 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 13 Jan 2021 15:56:19 -0600 Subject: [PATCH] Move broaden_wmp_polynomials to wmp.cpp. Make a few methods const --- include/openmc/math_functions.h | 13 --------- include/openmc/wmp.h | 18 +++++++++++-- src/math_functions.cpp | 45 ------------------------------- src/wmp.cpp | 48 +++++++++++++++++++++++++++++++-- 4 files changed, 62 insertions(+), 62 deletions(-) diff --git a/include/openmc/math_functions.h b/include/openmc/math_functions.h index e2d776dc2..79cd06325 100644 --- a/include/openmc/math_functions.h +++ b/include/openmc/math_functions.h @@ -203,19 +203,6 @@ extern "C" double normal_variate(double mean, double std_dev, uint64_t* seed); extern "C" double muir_spectrum(double e0, double m_rat, double kt, uint64_t* seed); -//============================================================================== -//! Doppler broadens the windowed multipole curvefit. -//! -//! 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 -//============================================================================== - -extern "C" void broaden_wmp_polynomials(double E, double dopp, int n, double factors[]); - //============================================================================== //! Constructs a natural cubic spline. //! diff --git a/include/openmc/wmp.h b/include/openmc/wmp.h index cd82d7a81..014ee7fe7 100644 --- a/include/openmc/wmp.h +++ b/include/openmc/wmp.h @@ -53,7 +53,7 @@ public: //! \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); + std::tuple evaluate(double E, double sqrtkT) const; //! \brief Evaluates the windowed multipole equations for the derivative of //! cross sections in the resolved resonance regions with respect to @@ -63,7 +63,7 @@ public: //! \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); + std::tuple evaluate_deriv(double E, double sqrtkT) const; // Data members std::string name_; //!< Name of nuclide @@ -97,6 +97,20 @@ void check_wmp_version(hid_t file); //! \param[in] i_nuclide Index in global nuclides array void read_multipole_data(int i_nuclide); + +//============================================================================== +//! Doppler broadens the windowed multipole curvefit. +//! +//! 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 +//============================================================================== + +extern "C" void broaden_wmp_polynomials(double E, double dopp, int n, double factors[]); + } // namespace openmc #endif // OPENMC_WMP_H diff --git a/src/math_functions.cpp b/src/math_functions.cpp index a3168fca1..ac60c6ce0 100644 --- a/src/math_functions.cpp +++ b/src/math_functions.cpp @@ -718,51 +718,6 @@ double watt_spectrum(double a, double b, uint64_t* seed) { } -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; - double half_inv_dopp2 = 0.5 / (dopp * dopp); - double quarter_inv_dopp4 = half_inv_dopp2 * half_inv_dopp2; - - double erf_beta; // error function of beta - double exp_m_beta2; // exp(-beta**2) - if (beta > 6.0) { - // Save time, ERF(6) is 1 to machine precision. - // beta/sqrtpi*exp(-beta**2) is also approximately 1 machine epsilon. - erf_beta = 1.; - exp_m_beta2 = 0.; - } else { - erf_beta = std::erf(beta); - exp_m_beta2 = std::exp(-beta * beta); - } - - // Assume that, for sure, we'll use a second order (1/E, 1/V, const) - // fit, and no less. - - factors[0] = erf_beta / E; - factors[1] = 1. / sqrtE; - factors[2] = factors[0] * (half_inv_dopp2 + E) + exp_m_beta2 / - (beta * SQRT_PI); - - // Perform recursive broadening of high order components - for (int i = 0; i < n - 3; i++) { - double ip1_dbl = i + 1; - if (i != 0) { - factors[i + 3] = -factors[i - 1] * (ip1_dbl - 1.) * ip1_dbl * - quarter_inv_dopp4 + factors[i + 1] * - (E + (1. + 2. * ip1_dbl) * half_inv_dopp2); - } else { - // Although it's mathematically identical, factors[0] will contain - // nothing, and we don't want to have to worry about memory. - factors[i + 3] = factors[i + 1] * - (E + (1. + 2. * ip1_dbl) * half_inv_dopp2); - } - } -} - - void spline(int n, const double x[], const double y[], double z[]) { std::vector c_new(n-1); diff --git a/src/wmp.cpp b/src/wmp.cpp index 5d44039d1..993e8f60d 100644 --- a/src/wmp.cpp +++ b/src/wmp.cpp @@ -78,7 +78,7 @@ WindowedMultipole::WindowedMultipole(hid_t group) } std::tuple -WindowedMultipole::evaluate(double E, double sqrtkT) +WindowedMultipole::evaluate(double E, double sqrtkT) const { using namespace std::complex_literals; @@ -161,7 +161,7 @@ WindowedMultipole::evaluate(double E, double sqrtkT) } std::tuple -WindowedMultipole::evaluate_deriv(double E, double sqrtkT) +WindowedMultipole::evaluate_deriv(double E, double sqrtkT) const { // ========================================================================== // Bookkeeping @@ -261,4 +261,48 @@ void read_multipole_data(int i_nuclide) file_close(file); } +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; + double half_inv_dopp2 = 0.5 / (dopp * dopp); + double quarter_inv_dopp4 = half_inv_dopp2 * half_inv_dopp2; + + double erf_beta; // error function of beta + double exp_m_beta2; // exp(-beta**2) + if (beta > 6.0) { + // Save time, ERF(6) is 1 to machine precision. + // beta/sqrtpi*exp(-beta**2) is also approximately 1 machine epsilon. + erf_beta = 1.; + exp_m_beta2 = 0.; + } else { + erf_beta = std::erf(beta); + exp_m_beta2 = std::exp(-beta * beta); + } + + // Assume that, for sure, we'll use a second order (1/E, 1/V, const) + // fit, and no less. + + factors[0] = erf_beta / E; + factors[1] = 1. / sqrtE; + factors[2] = factors[0] * (half_inv_dopp2 + E) + exp_m_beta2 / + (beta * SQRT_PI); + + // Perform recursive broadening of high order components + for (int i = 0; i < n - 3; i++) { + double ip1_dbl = i + 1; + if (i != 0) { + factors[i + 3] = -factors[i - 1] * (ip1_dbl - 1.) * ip1_dbl * + quarter_inv_dopp4 + factors[i + 1] * + (E + (1. + 2. * ip1_dbl) * half_inv_dopp2); + } else { + // Although it's mathematically identical, factors[0] will contain + // nothing, and we don't want to have to worry about memory. + factors[i + 3] = factors[i + 1] * + (E + (1. + 2. * ip1_dbl) * half_inv_dopp2); + } + } +} + } // namespace openmc