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

Merge pull request #1746 from paulromano/wmp-performance-opt

Browse source 
Performance optimizations for windowed multipole
This commit is contained in:
Gavin Ridley 2021-01-20 12:59:30 -05:00 • committed by GitHub
commit c262407af5
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
7 changed files with 122 additions and 111 deletions
Download patch file Download diff file Expand all files Collapse all files

6
docs/source/io_formats/data_wmp.rst
View file

@ -17,13 +17,13 @@ Windowed Multipole Library Format
If 1, Doppler broaden curve fit for window with corresponding index.
If 0, do not.
- **curvefit** (*double[][][]*)
Curve fit coefficients. Indexed by (reaction type, coefficient index,
window index).
Curve fit coefficients. Indexed by (window index, coefficient index,
reaction type).
- **data** (*complex[][]*)
Complex poles and residues. Each pole has a corresponding set of
residues. For example, the :math:`i`-th pole and corresponding residues
are stored as
.. math::
\text{data}[:,i] = [\text{pole},~\text{residue}_1,~\text{residue}_2,
~\ldots]

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

38
include/openmc/wmp.h
View file

@ -35,6 +35,13 @@ constexpr std::array<int, 2> WMP_VERSION {1, 1};
class WindowedMultipole {
public:
// Types
struct WindowInfo {
int index_start; // Index of starting pole
int index_end; // Index of ending pole
bool broaden_poly; // Whether to broaden polynomial curvefit
};
// Constructors, destructors
WindowedMultipole(hid_t group);
@ -46,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
@ -56,20 +63,19 @@ 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
bool fissionable_; //!< Is the nuclide fissionable?
xt::xtensor<std::complex<double>, 2> data_; //!< Poles and residues
double sqrt_awr_; //!< Square root of atomic weight ratio
double E_min_; //!< Minimum energy in [eV]
double E_max_; //!< Maximum energy in [eV]
double spacing_; //!< Spacing in sqrt(E) space
double sqrt_awr_; //!< Square root of atomic weight ratio
double inv_spacing_; //!< 1 / 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, 3> curvefit_; //!< Fitting function (reaction, coeff index, window index)
xt::xtensor<bool, 1> broaden_poly_; //!< Whether to broaden curvefit
bool fissionable_; //!< Is the nuclide fissionable?
std::vector<WindowInfo> window_info_; // Information about a window
xt::xtensor<double, 3> curvefit_; // Curve fit coefficients (window, poly order, reaction)
xt::xtensor<std::complex<double>, 2> data_; //!< Poles and residues
// Constant data
static constexpr int MAX_POLY_COEFFICIENTS =
@ -91,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

2
src/cell.cpp
View file

@ -514,7 +514,7 @@ CSGCell::distance(Position r, Direction u, int32_t on_surface, Particle* p) cons
// Check if this distance is the new minimum.
if (d < min_dist) {
if (std::abs(d - min_dist) / min_dist >= FP_PRECISION) {
if (min_dist - d >= FP_PRECISION*min_dist) {
min_dist = d;
i_surf = -token;
}

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

2
src/physics.cpp
View file

@ -973,7 +973,7 @@ sample_cxs_target_velocity(double awr, double E, Direction u, double kT, uint64_
if (prn(seed) < alpha) {
// With probability alpha, we sample the distribution p(y) =
// y*e^(-y). This can be done with sampling scheme C45 frmo the Monte
// y*e^(-y). This can be done with sampling scheme C45 from the Monte
// Carlo sampler
beta_vt_sq = -std::log(r1*r2);

127
src/wmp.cpp
View file

@ -24,7 +24,8 @@ WindowedMultipole::WindowedMultipole(hid_t group)
name_ = object_name(group).substr(1);
// Read scalar values.
read_dataset(group, "spacing", spacing_);
read_dataset(group, "spacing", inv_spacing_);
inv_spacing_ = 1.0 / inv_spacing_;
read_dataset(group, "sqrtAWR", sqrt_awr_);
read_dataset(group, "E_min", E_min_);
read_dataset(group, "E_max", E_max_);
@ -39,19 +40,22 @@ WindowedMultipole::WindowedMultipole(hid_t group)
// Read the "windows" array and use its shape to figure out the number of
// windows.
read_dataset(group, "windows", windows_);
int n_windows = windows_.shape()[0];
xt::xtensor<int, 2> windows;
read_dataset(group, "windows", windows);
int n_windows = windows.shape()[0];
windows -= 1; // Adjust to 0-based indices
// Read the "broaden_poly" arrays.
read_dataset(group, "broaden_poly", broaden_poly_);
if (n_windows != broaden_poly_.shape()[0]) {
xt::xtensor<bool, 1> broaden_poly;
read_dataset(group, "broaden_poly", broaden_poly);
if (n_windows != broaden_poly.shape()[0]) {
fatal_error("broaden_poly array shape is not consistent with the windows "
"array shape in WMP library for " + name_ + ".");
}
// Read the "curvefit" array.
read_dataset(group, "curvefit", curvefit_);
if (n_windows != broaden_poly_.shape()[0]) {
if (n_windows != curvefit_.shape()[0]) {
fatal_error("curvefit array shape is not consistent with the windows "
"array shape in WMP library for " + name_ + ".");
}
@ -63,10 +67,18 @@ WindowedMultipole::WindowedMultipole(hid_t group)
"Need to compile with WindowedMultipole::MAX_POLY_COEFFICIENTS = {}",
fit_order_ + 1));
}
// Move window information into a vector
window_info_.resize(n_windows);
for (int i = 0; i < n_windows; ++i) {
window_info_[i].index_start = windows(i, 0);
window_info_[i].index_end = windows(i, 1);
window_info_[i].broaden_poly = broaden_poly[i];
}
}
std::tuple<double, double, double>
WindowedMultipole::evaluate(double E, double sqrtkT)
WindowedMultipole::evaluate(double E, double sqrtkT) const
{
using namespace std::complex_literals;
@ -78,10 +90,9 @@ WindowedMultipole::evaluate(double E, double sqrtkT)
double invE = 1.0 / E;
// Locate window containing energy
int i_window = std::min(windows_.shape()[0] - 1,
static_cast<unsigned long>((sqrtE - std::sqrt(E_min_)) / spacing_));
int startw = windows_(i_window, 0) - 1;
int endw = windows_(i_window, 1) - 1;
int i_window = std::min(window_info_.size() - 1,
static_cast<size_t>((sqrtE - std::sqrt(E_min_)) * inv_spacing_));
const auto& window {window_info_[i_window]};
// Initialize the ouptut cross sections
double sig_s = 0.0;
@ -91,7 +102,7 @@ WindowedMultipole::evaluate(double E, double sqrtkT)
// ==========================================================================
// Add the contribution from the curvefit polynomial.
if (sqrtkT > 0.0 && broaden_poly_(i_window)) {
if (sqrtkT > 0.0 && window.broaden_poly) {
// Broaden the curvefit.
double dopp = sqrt_awr_ / sqrtkT;
std::array<double, MAX_POLY_COEFFICIENTS> broadened_polynomials;
@ -121,9 +132,9 @@ WindowedMultipole::evaluate(double E, double sqrtkT)
if (sqrtkT == 0.0) {
// If at 0K, use asymptotic form.
for (int i_pole = startw; i_pole <= endw; ++i_pole) {
for (int i_pole = window.index_start; i_pole <= window.index_end; ++i_pole) {
std::complex<double> psi_chi = -1.0i / (data_(i_pole, MP_EA) - sqrtE);
std::complex<double> c_temp = psi_chi / E;
std::complex<double> c_temp = psi_chi * invE;
sig_s += (data_(i_pole, MP_RS) * c_temp).real();
sig_a += (data_(i_pole, MP_RA) * c_temp).real();
if (fissionable_) {
@ -133,15 +144,13 @@ WindowedMultipole::evaluate(double E, double sqrtkT)
} 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();
}
for (int i_pole = window.index_start; i_pole <= window.index_end; ++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();
}
}
}
@ -150,7 +159,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
@ -166,9 +175,8 @@ WindowedMultipole::evaluate_deriv(double E, double sqrtkT)
}
// 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;
int i_window = (sqrtE - std::sqrt(E_min_)) * inv_spacing_;
const auto& window {window_info_[i_window]};
// Initialize the ouptut cross sections.
double sig_s = 0.0;
@ -184,20 +192,19 @@ WindowedMultipole::evaluate_deriv(double E, double sqrtkT)
// 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();
}
for (int i_pole = window.index_start; i_pole <= window.index_end; ++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);
}
double norm = -0.5*sqrt_awr_ / std::sqrt(K_BOLTZMANN) * std::pow(T, -1.5);
sig_s *= norm;
sig_a *= norm;
sig_f *= norm;
return std::make_tuple(sig_s, sig_a, sig_f);
}
@ -250,4 +257,46 @@ void read_multipole_data(int i_nuclide)
file_close(file);
}
void broaden_wmp_polynomials(double E, double dopp, int n, double factors[])
{
// Broadening of polynomials follows procedure outlined in C. Josey, P. Ducru,
// B. Forget, and K. Smith, "Windowed multipole for cross section Doppler
// broadening," J. Comput. Phys., 307, 715-727 (2016).
// https://doi.org/10.1016/j.jcp.2015.08.013
// 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);
if (n > 3) factors[3] = factors[1] * (E + 3.0 * half_inv_dopp2);
// Perform recursive broadening of high order components (Eq. 16)
for (int i = 1; i < n - 3; i++) {
double ip1_dbl = i + 1;
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);
}
}
} // namespace openmc