Store min/max data temperatures. Use in Cell::set_temperature to check bounds

include/openmc/nuclide.h

@@ -126,6 +126,12 @@ namespace data {
 extern std::array<double, 2> energy_min;
 extern std::array<double, 2> energy_max;
 
+//! Minimum temperature in [K] that nuclide data is available at
+extern double temperature_min;
+
+//! Maximum temperature in [K] that nuclide data is available at
+extern double temperature_max;
+
 extern std::vector<std::unique_ptr<Nuclide>> nuclides;
 extern std::unordered_map<std::string, int> nuclide_map;
 

src/cell.cpp

@@ -16,6 +16,7 @@
 #include "openmc/hdf5_interface.h"
 #include "openmc/lattice.h"
 #include "openmc/material.h"
+#include "openmc/nuclide.h"
 #include "openmc/settings.h"
 #include "openmc/surface.h"
 #include "openmc/xml_interface.h"
@@ -244,6 +245,16 @@ Cell::temperature(int32_t instance) const
 void
 Cell::set_temperature(double T, int32_t instance)
 {
+  if (settings::temperature_method == TEMPERATURE_INTERPOLATION) {
+    if (T < data::temperature_min) {
+      throw std::runtime_error{"Temperature is below minimum temperature at "
+        "which data is available."};
+    } else if (T > data::temperature_max) {
+      throw std::runtime_error{"Temperature is above maximum temperature at "
+        "which data is available."};
+    }
+  }
+
   if (instance >= 0) {
     sqrtkT_.at(instance) = std::sqrt(K_BOLTZMANN * T);
   } else {
@@ -1057,7 +1068,6 @@ openmc_cell_set_fill(int32_t index, int type, int32_t n,
   return 0;
 }
 
-//TODO: make sure data is loaded for this temperature
 extern "C" int
 openmc_cell_set_temperature(int32_t index, double T, const int32_t* instance)
 {

src/cross_sections.cpp

@@ -357,6 +357,12 @@ read_ce_cross_sections(const std::vector<std::vector<double>>& nuc_temps,
     }
   }
 
+  // Show minimum/maximum temperature
+  write_message("Minimum neutron data temperature: " +
+    std::to_string(data::temperature_min) + " K");
+  write_message("Maximum neutron data temperature: " +
+    std::to_string(data::temperature_max) + " K");
+
   // If the user wants multipole, make sure we found a multipole library.
   if (settings::temperature_multipole) {
     bool mp_found = false;

src/finalize.cpp

@@ -117,6 +117,8 @@ int openmc_finalize()
 
   data::energy_max = {INFTY, INFTY};
   data::energy_min = {0.0, 0.0};
+  data::temperature_min = 0.0;
+  data::temperature_max = INFTY;
   model::root_universe = -1;
   openmc::openmc_set_seed(DEFAULT_SEED);
 

src/nuclide.cpp

@@ -18,7 +18,7 @@
 #include "xtensor/xbuilder.hpp"
 #include "xtensor/xview.hpp"
 
-#include <algorithm> // for sort
+#include <algorithm> // for sort, min_element
 #include <string> // for to_string, stoi
 
 namespace openmc {
@@ -30,6 +30,8 @@ namespace openmc {
 namespace data {
 std::array<double, 2> energy_min {0.0, 0.0};
 std::array<double, 2> energy_max {INFTY, INFTY};
+double temperature_min {0.0};
+double temperature_max {INFTY};
 std::vector<std::unique_ptr<Nuclide>> nuclides;
 std::unordered_map<std::string, int> nuclide_map;
 } // namespace data
@@ -154,6 +156,12 @@ Nuclide::Nuclide(hid_t group, const std::vector<double>& temperature, int i_nucl
   // Sort temperatures to read
   std::sort(temps_to_read.begin(), temps_to_read.end());
 
+  double T_min_read = *std::min_element(temps_to_read.cbegin(), temps_to_read.cend());
+  double T_max_read = *std::max_element(temps_to_read.cbegin(), temps_to_read.cend());
+
+  data::temperature_min = std::max(data::temperature_min, T_min_read);
+  data::temperature_max = std::min(data::temperature_max, T_max_read);
+
   hid_t energy_group = open_group(group, "energy");
   for (const auto& T : temps_to_read) {
     std::string dset {std::to_string(T) + "K"};