From c704ae5ac18ebd6d79cf2f1c62e12d271cd994f6 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 13 Jan 2021 14:50:30 -0600 Subject: [PATCH 1/6] Avoid calling std::abs and divsd instruction in CSGCell::distance --- src/cell.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cell.cpp b/src/cell.cpp index 6ca39749f..418258d8b 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -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; } From 425de8943fe29e3ebbfa96fcd9288e08823f636c Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 26 Apr 2019 11:04:59 -0500 Subject: [PATCH 2/6] Minor optimizations for WindowedMultipole --- docs/source/io_formats/data_wmp.rst | 6 +++--- include/openmc/wmp.h | 2 +- src/physics.cpp | 2 +- src/wmp.cpp | 20 +++++++++++--------- 4 files changed, 16 insertions(+), 14 deletions(-) diff --git a/docs/source/io_formats/data_wmp.rst b/docs/source/io_formats/data_wmp.rst index 8d7cf6ba0..e7261f3fb 100644 --- a/docs/source/io_formats/data_wmp.rst +++ b/docs/source/io_formats/data_wmp.rst @@ -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] diff --git a/include/openmc/wmp.h b/include/openmc/wmp.h index 2378f89da..6bf14e51b 100644 --- a/include/openmc/wmp.h +++ b/include/openmc/wmp.h @@ -65,7 +65,7 @@ public: 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 inv_spacing_; //!< 1 / 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) diff --git a/src/physics.cpp b/src/physics.cpp index 2636207fb..201a893c6 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -972,7 +972,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); diff --git a/src/wmp.cpp b/src/wmp.cpp index 708cf923f..10921fca3 100644 --- a/src/wmp.cpp +++ b/src/wmp.cpp @@ -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_); @@ -41,6 +42,7 @@ WindowedMultipole::WindowedMultipole(hid_t group) // 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_); @@ -51,7 +53,7 @@ WindowedMultipole::WindowedMultipole(hid_t group) // 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_ + "."); } @@ -79,9 +81,9 @@ WindowedMultipole::evaluate(double E, double sqrtkT) // Locate window containing energy int i_window = std::min(windows_.shape()[0] - 1, - static_cast((sqrtE - std::sqrt(E_min_)) / spacing_)); - int startw = windows_(i_window, 0) - 1; - int endw = windows_(i_window, 1) - 1; + static_cast((sqrtE - std::sqrt(E_min_)) * inv_spacing_)); + int startw = windows_(i_window, 0); + int endw = windows_(i_window, 1); // Initialize the ouptut cross sections double sig_s = 0.0; @@ -123,7 +125,7 @@ WindowedMultipole::evaluate(double E, double sqrtkT) // 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; + std::complex 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_) { @@ -166,9 +168,9 @@ 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_; + int startw = windows_(i_window, 0); + int endw = windows_(i_window, 1); // Initialize the ouptut cross sections. double sig_s = 0.0; From 47828091aa96cf6944317ae5c4ea7f615a25bd60 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 26 Apr 2019 14:32:04 -0500 Subject: [PATCH 3/6] Improve memory layout slightly for windowed multipole evaluations --- include/openmc/wmp.h | 18 +++++++---- src/wmp.cpp | 75 +++++++++++++++++++++++++------------------- 2 files changed, 54 insertions(+), 39 deletions(-) diff --git a/include/openmc/wmp.h b/include/openmc/wmp.h index 6bf14e51b..cd82d7a81 100644 --- a/include/openmc/wmp.h +++ b/include/openmc/wmp.h @@ -35,6 +35,13 @@ constexpr std::array 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); @@ -60,16 +67,15 @@ public: // Data members std::string name_; //!< Name of nuclide - bool fissionable_; //!< Is the nuclide fissionable? - xt::xtensor, 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 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 windows_; //!< Indices of pole at start/end of window - xt::xtensor curvefit_; //!< Fitting function (reaction, coeff index, window index) - xt::xtensor broaden_poly_; //!< Whether to broaden curvefit + bool fissionable_; //!< Is the nuclide fissionable? + std::vector window_info_; // Information about a window + xt::xtensor curvefit_; // Curve fit coefficients (window, poly order, reaction) + xt::xtensor, 2> data_; //!< Poles and residues // Constant data static constexpr int MAX_POLY_COEFFICIENTS = diff --git a/src/wmp.cpp b/src/wmp.cpp index 10921fca3..5d44039d1 100644 --- a/src/wmp.cpp +++ b/src/wmp.cpp @@ -40,13 +40,15 @@ 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]; - windows_ -= 1; // Adjust to 0-based indices + xt::xtensor 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 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_ + "."); } @@ -65,6 +67,14 @@ 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 @@ -80,10 +90,11 @@ 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((sqrtE - std::sqrt(E_min_)) * inv_spacing_)); - int startw = windows_(i_window, 0); - int endw = windows_(i_window, 1); + int i_window = std::min(window_info_.size() - 1, + static_cast((sqrtE - std::sqrt(E_min_)) * inv_spacing_)); + const auto& window {window_info_[i_window]}; + int startw = window.index_start; + int endw = window.index_end; // Initialize the ouptut cross sections double sig_s = 0.0; @@ -93,7 +104,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 broadened_polynomials; @@ -135,15 +146,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 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(); - } + 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(); } } } @@ -169,8 +178,9 @@ WindowedMultipole::evaluate_deriv(double E, double sqrtkT) // Locate us int i_window = (sqrtE - std::sqrt(E_min_)) * inv_spacing_; - int startw = windows_(i_window, 0); - int endw = windows_(i_window, 1); + const auto& window {window_info_[i_window]}; + int startw = window.index_start; + int endw = window.index_end; // Initialize the ouptut cross sections. double sig_s = 0.0; @@ -186,20 +196,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 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(); - } + 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); } + 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); } From 4fc79565150eabd407d7bfea19134bc052d1ba37 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 13 Jan 2021 15:56:19 -0600 Subject: [PATCH 4/6] 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 From d26bdefd1f1d0f75e39180305d68c21cb4cdd41b Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 18 Jan 2021 21:46:37 -0600 Subject: [PATCH 5/6] Improvement in broaden_wmp_polynomials suggested by @gridley --- src/wmp.cpp | 24 +++++++++++------------- 1 file changed, 11 insertions(+), 13 deletions(-) diff --git a/src/wmp.cpp b/src/wmp.cpp index 993e8f60d..7dd24294d 100644 --- a/src/wmp.cpp +++ b/src/wmp.cpp @@ -263,6 +263,11 @@ void read_multipole_data(int i_nuclide) 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; @@ -283,25 +288,18 @@ void broaden_wmp_polynomials(double E, double dopp, int n, double factors[]) // 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 - for (int i = 0; i < n - 3; i++) { + // Perform recursive broadening of high order components (Eq. 16) + for (int i = 1; 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); - } + 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); } } From dc743da596914e16ebc18af1aa44dc5871ad3a24 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 19 Jan 2021 13:40:38 -0600 Subject: [PATCH 6/6] Remove startw/endw per @gridley suggestion --- src/wmp.cpp | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/src/wmp.cpp b/src/wmp.cpp index 7dd24294d..6a56ddd0f 100644 --- a/src/wmp.cpp +++ b/src/wmp.cpp @@ -93,8 +93,6 @@ WindowedMultipole::evaluate(double E, double sqrtkT) const int i_window = std::min(window_info_.size() - 1, static_cast((sqrtE - std::sqrt(E_min_)) * inv_spacing_)); const auto& window {window_info_[i_window]}; - int startw = window.index_start; - int endw = window.index_end; // Initialize the ouptut cross sections double sig_s = 0.0; @@ -134,7 +132,7 @@ WindowedMultipole::evaluate(double E, double sqrtkT) const 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 psi_chi = -1.0i / (data_(i_pole, MP_EA) - sqrtE); std::complex c_temp = psi_chi * invE; sig_s += (data_(i_pole, MP_RS) * c_temp).real(); @@ -146,7 +144,7 @@ WindowedMultipole::evaluate(double E, double sqrtkT) const } else { // At temperature, use Faddeeva function-based form. double dopp = sqrt_awr_ / sqrtkT; - 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 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(); @@ -179,8 +177,6 @@ WindowedMultipole::evaluate_deriv(double E, double sqrtkT) const // Locate us int i_window = (sqrtE - std::sqrt(E_min_)) * inv_spacing_; const auto& window {window_info_[i_window]}; - int startw = window.index_start; - int endw = window.index_end; // Initialize the ouptut cross sections. double sig_s = 0.0; @@ -196,7 +192,7 @@ WindowedMultipole::evaluate_deriv(double E, double sqrtkT) const // Add the contribution from the poles in this window. double dopp = sqrt_awr_ / sqrtkT; - 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 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();