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

Move broaden_wmp_polynomials to wmp.cpp. Make a few methods const

Browse source
This commit is contained in:
Paul Romano 2021-01-13 15:56:19 -06:00
parent 47828091aa
commit 4fc7956515
4 changed files with 62 additions and 62 deletions
Download patch file Download diff file Expand all files Collapse all files

13
include/openmc/math_functions.h
View file

@ -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.
//!

18
include/openmc/wmp.h
View file

@ -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<double, double, double> evaluate(double E, double sqrtkT);
std::tuple<double, double, double> 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<double, double, double> evaluate_deriv(double E, double sqrtkT);
std::tuple<double, double, double> 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

45
src/math_functions.cpp
View file

@ -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<double> c_new(n-1);

48
src/wmp.cpp
View file

@ -78,7 +78,7 @@ WindowedMultipole::WindowedMultipole(hid_t group)
}
std::tuple<double, double, double>
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<double, double, double>
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