Account for temperature in reaction rate calculation

include/openmc/nuclide.h

@@ -35,6 +35,11 @@ public:
     std::vector<double> 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<double>& temperature);
   ~Nuclide();
@@ -61,10 +66,12 @@ public:
   //! \brief Calculate reaction rate based on group-wise flux distribution
   //
   //! \param[in] MT ENDF MT value for desired reaction
+  //! \param[in] temperature Temperature in [K]
   //! \param[in] energy Energy group boundaries in [eV]
   //! \param[in] flux Flux in each energy group (not normalized per eV)
   //! \return Reaction rate
-  double collapse_rate(int MT, gsl::span<const double> energy, gsl::span<const double> flux) const;
+  double collapse_rate(int MT, double temperature, gsl::span<const double> energy,
+    gsl::span<const double> flux) const;
 
   // Data members
   std::string name_; //!< Name of nuclide, e.g. "U235"
@@ -113,6 +120,8 @@ public:
 private:
   void create_derived(const Function1D* prompt_photons, const Function1D* delayed_photons);
 
+  InterpResult find_temperature(double T) const;
+
   static int XS_TOTAL;
   static int XS_ABSORPTION;
   static int XS_FISSION;

include/openmc/reaction.h

@@ -29,12 +29,13 @@ public:
 
   //! \brief Calculate reaction rate based on group-wise flux distribution
   //
+  //! \param[in] i_temp Temperature index
   //! \param[in] energy Energy group boundaries in [eV]
   //! \param[in] flux Flux in each energy group (not normalized per eV)
   //! \param[in] grid Nuclide energy grid
   //! \return Reaction rate
-  double collapse_rate(gsl::span<const double> energy, gsl::span<const double> flux,
-    const std::vector<double>& grid) const;
+  double collapse_rate(gsl::index i_temp, gsl::span<const double> energy,
+    gsl::span<const double> flux, const std::vector<double>& grid) const;
 
   //! Cross section at a single temperature
   struct TemperatureXS {

openmc/deplete/helpers.py

@@ -194,7 +194,7 @@ class FluxCollapseHelper(ReactionRateHelper):
             for mt, i_react in zip(self._mts, react_index):
                 # Use flux to collapse reaction rate (per N)
                 nuc = openmc.lib.nuclides[name]
-                rate_per_nuc = nuc.collapse_rate(mt, self._energies, flux)
+                rate_per_nuc = nuc.collapse_rate(mt, mat.temperature, self._energies, flux)
 
                 # Multiply by density to get absolute reaction rate
                 self._results_cache[i_nuc, i_react] = rate_per_nuc * density

openmc/lib/nuclide.py

@@ -25,8 +25,10 @@ _dll.openmc_load_nuclide.errcheck = _error_handler
 _dll.openmc_nuclide_name.argtypes = [c_int, POINTER(c_char_p)]
 _dll.openmc_nuclide_name.restype = c_int
 _dll.openmc_nuclide_name.errcheck = _error_handler
-_dll.openmc_nuclide_collapse_rate.argtypes = [c_int, c_int, _array_1d_dble,
-    _array_1d_dble, c_int, POINTER(c_double)]
+_dll.openmc_nuclide_collapse_rate.argtypes = [c_int, c_int, c_double,
+    _array_1d_dble, _array_1d_dble, c_int, POINTER(c_double)]
+_dll.openmc_nuclide_collapse_rate.restype = c_int
+_dll.openmc_nuclide_collapse_rate.errcheck = _error_handler
 _dll.nuclides_size.restype = c_size_t
 
 
@@ -77,13 +79,15 @@ class Nuclide(_FortranObject):
         _dll.openmc_nuclide_name(self._index, name)
         return name.value.decode()
 
-    def collapse_rate(self, MT, energy, flux):
+    def collapse_rate(self, MT, temperature, energy, flux):
         """Calculate reaction rate based on group-wise flux distribution
 
         Parameters
         ----------
         MT : int
             ENDF MT value of the desired reaction
+        temperature : float
+            Temperature in [K] at which to evaluate cross sections
         energy : iterable of float
             Energy group boundaries in [eV]
         flux : iterable of float
@@ -98,7 +102,8 @@ class Nuclide(_FortranObject):
         energy = np.asarray(energy, dtype=float)
         flux = np.asarray(flux, dtype=float)
         xs = c_double()
-        _dll.openmc_nuclide_collapse_rate(self._index, MT, energy, flux, len(flux), xs)
+        _dll.openmc_nuclide_collapse_rate(self._index, MT, temperature, energy,
+                                          flux, len(flux), xs)
         return xs.value
 
 

src/dagmc.cpp

@@ -263,10 +263,8 @@ void load_dagmc_geometry()
       MB_CHK_ERR_CONT(rval);
       double temp = std::stod(temp_value);
       c->sqrtkT_.push_back(std::sqrt(K_BOLTZMANN * temp));
-    } else if (mat->temperature_ > 0.0) {
-      c->sqrtkT_.push_back(std::sqrt(K_BOLTZMANN * mat->temperature_));
     } else {
-      c->sqrtkT_.push_back(std::sqrt(K_BOLTZMANN * settings::temperature_default));
+      c->sqrtkT_.push_back(std::sqrt(K_BOLTZMANN * mat->temperature()));
     }
 
   }

src/nuclide.cpp

@@ -15,6 +15,8 @@
 #include "openmc/string_utils.h"
 #include "openmc/thermal.h"
 
+#include <fmt/core.h>
+
 #include "xtensor/xbuilder.hpp"
 #include "xtensor/xview.hpp"
 
@@ -908,7 +910,53 @@ void Nuclide::calculate_urr_xs(int i_temp, Particle& p) const
 
 }
 
-double Nuclide::collapse_rate(int MT, gsl::span<const double> energy, gsl::span<const double> flux) const
+Nuclide::InterpResult 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;
+  double kT = K_BOLTZMANN * T;
+  switch (settings::temperature_method) {
+  case TemperatureMethod::NEAREST:
+    {
+      double max_diff = INFTY;
+      for (gsl::index t = 0; t < kTs_.size(); ++t) {
+        double diff = std::abs(kTs_[t] - kT);
+        if (diff < max_diff) {
+          res.i = 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]);
+  }
+
+  Ensures(res.i >= 0 && res.i < kTs_.size());
+  Ensures(res.f >= 0.0 && res.f <= 1.0);
+
+  return res;
+}
+
+double Nuclide::collapse_rate(int MT, double temperature, gsl::span<const double> energy,
+  gsl::span<const double> flux) const
 {
   Expects(MT > 0);
   Expects(energy.size() > 0);
@@ -916,10 +964,24 @@ double Nuclide::collapse_rate(int MT, gsl::span<const double> energy, gsl::span<
 
   int i_rx = reaction_index_[MT];
   if (i_rx < 0) return 0.0;
-
   const auto& rx = reactions_[i_rx];
-  const auto& grid = grid_[0].energy;
-  return rx->collapse_rate(energy, flux, grid);
+
+  // Determine temperature index
+  auto res = 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);
+
+  if (res.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);
+  } else {
+    // If interpolation factor is zero, return reaction rate at lower temperature
+    return rr_low;
+  }
 }
 
 //==============================================================================
@@ -1044,14 +1106,22 @@ openmc_nuclide_name(int index, const char** name)
 }
 
 extern "C" int
-openmc_nuclide_collapse_rate(int index, int MT, const double* energy, const double* flux, int n, double* xs)
+openmc_nuclide_collapse_rate(int index, int MT, double temperature,
+  const double* energy, const double* flux, int n, double* xs)
 {
   if (index < 0 || index >= data::nuclides.size()) {
     set_errmsg("Index in nuclides vector is out of bounds.");
     return OPENMC_E_OUT_OF_BOUNDS;
   }
 
-  *xs = data::nuclides[index]->collapse_rate(MT, {energy, energy + n + 1}, {flux, flux + n});
+  try {
+    *xs = data::nuclides[index]->collapse_rate(MT, temperature,
+      {energy, energy + n + 1}, {flux, flux + n});
+  } catch (const std::out_of_range& e) {
+    fmt::print("Cuaght error\n");
+    set_errmsg(e.what());
+    return OPENMC_E_OUT_OF_BOUNDS;
+  }
   return 0;
 }
 

src/reaction.cpp

@@ -85,11 +85,9 @@ Reaction::Reaction(hid_t group, const std::vector<int>& temperatures)
 }
 
 double
-Reaction::collapse_rate(gsl::span<const double> energy, gsl::span<const double> flux, const std::vector<double>& grid) const
+Reaction::collapse_rate(gsl::index i_temp, gsl::span<const double> energy,
+  gsl::span<const double> flux, const std::vector<double>& grid) const
 {
-  // TODO: Figure out how to deal with temperature
-  int i_temp = 0;
-
   // Find index corresponding to first energy
   const auto& xs = xs_[i_temp].value;
   int i_low = lower_bound_index(grid.cbegin(), grid.cend(), energy.front());