From 1e436577ea6d202e2124dee318f212fe043b62e9 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 7 Aug 2020 10:51:28 -0500 Subject: [PATCH] Use pair instead of struct as return type of find_temperature --- include/openmc/nuclide.h | 12 +++++----- src/nuclide.cpp | 49 +++++++++++++++++----------------------- 2 files changed, 27 insertions(+), 34 deletions(-) diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h index 0155c318f4..69bb50cc8f 100644 --- a/include/openmc/nuclide.h +++ b/include/openmc/nuclide.h @@ -7,6 +7,7 @@ #include #include // for unique_ptr #include +#include // for pair #include #include @@ -35,11 +36,6 @@ public: std::vector energy; }; - struct InterpResult { - gsl::index i; //!< Index in tabulated data - double f; //!< Interpolation factor between i and i+1 - }; - // Constructors/destructors Nuclide(hid_t group, const std::vector& temperature); ~Nuclide(); @@ -120,7 +116,11 @@ public: private: void create_derived(const Function1D* prompt_photons, const Function1D* delayed_photons); - InterpResult find_temperature(double T) const; + //! Determine temperature index and interpolation factor + // + //! \param[in] T Temperature in [K] + //! \return Temperature index and interpolation factor + std::pair find_temperature(double T) const; static int XS_TOTAL; static int XS_ABSORPTION; diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 0d69adb8ec..02b0ed3653 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -910,49 +910,40 @@ void Nuclide::calculate_urr_xs(int i_temp, Particle& p) const } -Nuclide::InterpResult Nuclide::find_temperature(double T) const +std::pair Nuclide::find_temperature(double T) const { Expects(T >= 0.0); - // Make sure value is within bounds - double T_min = kTs_.front() / K_BOLTZMANN; - double T_max = kTs_.back() / K_BOLTZMANN; - if (T < T_min || T > T_max) { - throw std::out_of_range{fmt::format( - "Temperature out of range (min={:.1f}, max={:.1f}).", T_min, T_max - )}; - } - // Determine temperature index - InterpResult res; + gsl::index i_temp = 0; + double f = 0.0; double kT = K_BOLTZMANN * T; + gsl::index n = kTs_.size(); switch (settings::temperature_method) { case TemperatureMethod::NEAREST: { double max_diff = INFTY; - for (gsl::index t = 0; t < kTs_.size(); ++t) { + for (gsl::index t = 0; t < n; ++t) { double diff = std::abs(kTs_[t] - kT); if (diff < max_diff) { - res.i = t; + i_temp = t; max_diff = diff; } } } - res.f = 0.0; break; case TemperatureMethod::INTERPOLATION: // Find temperatures that bound the actual temperature - for (res.i = 0; res.i < kTs_.size() - 1; ++res.i) { - if (kTs_[res.i] <= kT && kT < kTs_[res.i + 1]) break; - } - res.f = (kT - kTs_[res.i]) / (kTs_[res.i + 1] - kTs_[res.i]); + while (kTs_[i_temp + 1] < kT && i_temp + 1 < n - 1) ++i_temp; + + // Determine interpolation factor + f = (kT - kTs_[i_temp]) / (kTs_[i_temp + 1] - kTs_[i_temp]); } - Ensures(res.i >= 0 && res.i < kTs_.size()); - Ensures(res.f >= 0.0 && res.f <= 1.0); + Ensures(i_temp >= 0 && i_temp < n); - return res; + return {i_temp, f}; } double Nuclide::collapse_rate(int MT, double temperature, gsl::span energy, @@ -967,17 +958,19 @@ double Nuclide::collapse_rate(int MT, double temperature, gsl::spanfind_temperature(temperature); + gsl::index i_temp; + double f; + std::tie(i_temp, f) = this->find_temperature(temperature); // Get reaction rate at lower temperature - const auto& grid_low = grid_[res.i].energy; - double rr_low = rx->collapse_rate(res.i, energy, flux, grid_low); + const auto& grid_low = grid_[i_temp].energy; + double rr_low = rx->collapse_rate(i_temp, energy, flux, grid_low); - if (res.f > 0.0) { + if (f > 0.0) { // Interpolate between reaction rate at lower and higher temperature - const auto& grid_high = grid_[res.i + 1].energy; - double rr_high = rx->collapse_rate(res.i + 1, energy, flux, grid_high); - return rr_low + res.f*(rr_high - rr_low); + const auto& grid_high = grid_[i_temp + 1].energy; + double rr_high = rx->collapse_rate(i_temp + 1, energy, flux, grid_high); + return rr_low + f*(rr_high - rr_low); } else { // If interpolation factor is zero, return reaction rate at lower temperature return rr_low;