From 47828091aa96cf6944317ae5c4ea7f615a25bd60 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 26 Apr 2019 14:32:04 -0500 Subject: [PATCH] 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 6bf14e51b5..cd82d7a811 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 10921fca39..5d44039d16 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); }