Add tolerance for interp temps outside of bounds

docs/source/io_formats/settings.rst

@@ -832,7 +832,9 @@ cell, the nearest temperature at which cross sections are given is to be
 applied, within a given tolerance (see :ref:`temperature_tolerance`). A value of
 "interpolation" indicates that cross sections are to be linear-linear
 interpolated between temperatures at which nuclear data are present (see
-:ref:`temperature_treatment`).
+:ref:`temperature_treatment`). With the "interpolation" method, temperatures
+outside of the bounds of the nuclear data may be accepted, provided they still
+fall within the tolerance (see :ref:`temperature_tolerance`).
 
   *Default*: "nearest"
 
@@ -871,7 +873,12 @@ The ``<temperature_tolerance>`` element specifies a tolerance in Kelvin that is
 to be applied when the "nearest" temperature method is used. For example, if a
 cell temperature is 340 K and the tolerance is 15 K, then the closest
 temperature in the range of 325 K to 355 K will be used to evaluate cross
-sections.
+sections. If the ``<temperature_method>`` is "interpolation", the tolerance
+specified applies to cell temperatures outside of the data bounds. For example,
+If a cell is specified at 695K, a tolerance of 15K and data only available at
+700K and 1000K, the cell's cross sections will be evaluated at 700K, since
+desired temperature of 695K is within the tolerance of the actual data despite
+not being bounded on both sides.
 
   *Default*: 10 K
 

openmc/settings.py

@@ -186,13 +186,15 @@ class Settings:
         'default', 'method', 'range', 'tolerance', and 'multipole'. The value
         for 'default' should be a float representing the default temperature in
         Kelvin. The value for 'method' should be 'nearest' or 'interpolation'.
-        If the method is 'nearest', 'tolerance' indicates a range of temperature
-        within which cross sections may be used. The value for 'range' should be
-        a pair a minimum and maximum temperatures which are used to indicate
-        that cross sections be loaded at all temperatures within the
-        range. 'multipole' is a boolean indicating whether or not the windowed
-        multipole method should be used to evaluate resolved resonance cross
-        sections.
+        If the method is 'nearest', 'tolerance' indicates a range of
+        temperature within which cross sections may be used. If the method is
+        'interpolation', 'tolerance' indicates the range of temperatures outside
+        of the available cross section temperatures where cross sections will
+        evaluate to the nearer bound. The value for 'range' should be a pair a
+        minimum and maximum temperatures which are used to indicate that cross
+        sections be loaded at all temperatures within the range. 'multipole' is
+        a boolean indicating whether or not the windowed multipole method should
+        be used to evaluate resolved resonance cross sections.
     trace : tuple or list
         Show detailed information about a single particle, indicated by three
         integers: the batch number, generation number, and particle number

src/cell.cpp

@@ -334,10 +334,10 @@ double Cell::temperature(int32_t instance) const
 void Cell::set_temperature(double T, int32_t instance, bool set_contained)
 {
   if (settings::temperature_method == TemperatureMethod::INTERPOLATION) {
-    if (T < data::temperature_min) {
+    if (T < data::temperature_min - settings::temperature_tolerance) {
       throw std::runtime_error {"Temperature is below minimum temperature at "
                                 "which data is available."};
-    } else if (T > data::temperature_max) {
+    } else if (T > data::temperature_max + settings::temperature_tolerance) {
       throw std::runtime_error {"Temperature is above maximum temperature at "
                                 "which data is available."};
     }

src/nuclide.cpp

@@ -169,6 +169,22 @@ Nuclide::Nuclide(hid_t group, const vector<double>& temperature)
       }
 
       if (!found_pair) {
+        // If no pairs found, check if the desired temperature falls just
+        // outside of data
+        if (T_desired - temps_available.front() <=
+            -settings::temperature_tolerance) {
+          if (!contains(temps_to_read, temps_available.front())) {
+            temps_to_read.push_back(std::round(temps_available.front()));
+          }
+          break;
+        }
+        if (T_desired - temps_available.back() <=
+            settings::temperature_tolerance) {
+          if (!contains(temps_to_read, temps_available.back())) {
+            temps_to_read.push_back(std::round(temps_available.back()));
+          }
+          break;
+        }
         fatal_error(
           "Nuclear data library does not contain cross sections for " + name_ +
           " at temperatures that bound " + std::to_string(T_desired) + " K.");
@@ -646,6 +662,16 @@ void Nuclide::calculate_xs(
     } break;
 
     case TemperatureMethod::INTERPOLATION:
+      // If current kT outside of the bounds of available, snap to the bound
+      if (kT < kTs_.front()) {
+        i_temp = 0;
+        break;
+      }
+      if (kT > kTs_.back()) {
+        i_temp = kTs_.size() - 1;
+        break;
+      }
+
       // Find temperatures that bound the actual temperature
       for (i_temp = 0; i_temp < kTs_.size() - 1; ++i_temp) {
         if (kTs_[i_temp] <= kT && kT < kTs_[i_temp + 1])
@@ -969,6 +995,15 @@ std::pair<gsl::index, double> Nuclide::find_temperature(double T) const
   } break;
 
   case TemperatureMethod::INTERPOLATION:
+    // If current kT outside of the bounds of available, snap to the bound
+    if (kT < kTs_.front()) {
+      i_temp = 0;
+      break;
+    }
+    if (kT > kTs_.back()) {
+      i_temp = kTs_.size() - 1;
+      break;
+    }
     // Find temperatures that bound the actual temperature
     while (kTs_[i_temp + 1] < kT && i_temp + 1 < n - 1)
       ++i_temp;

src/thermal.cpp

@@ -117,6 +117,16 @@ ThermalScattering::ThermalScattering(
         }
       }
       if (!found) {
+        // If no pairs found, check if the desired temperature falls within
+        // bounds' tolerance
+        if (T - temps_available[0] <= -settings::temperature_tolerance) {
+          temps_to_read.push_back(std::round(temps_available[0]));
+          break;
+        }
+        if (T - temps_available[n - 1] <= settings::temperature_tolerance) {
+          temps_to_read.push_back(std::round(temps_available[n - 1]));
+          break;
+        }
         fatal_error(
           fmt::format("Nuclear data library does not contain cross "
                       "sections for {} at temperatures that bound {} K.",
@@ -159,20 +169,27 @@ void ThermalScattering::calculate_xs(double E, double sqrtkT, int* i_temp,
 
   auto n = kTs_.size();
   if (n > 1) {
-    // Find temperatures that bound the actual temperature
-    while (kTs_[i + 1] < kT && i + 1 < n - 1)
-      ++i;
-
     if (settings::temperature_method == TemperatureMethod::NEAREST) {
+      while (kTs_[i + 1] < kT && i + 1 < n - 1)
+        ++i;
       // Pick closer of two bounding temperatures
       if (kT - kTs_[i] > kTs_[i + 1] - kT)
         ++i;
-
     } else {
-      // Randomly sample between temperature i and i+1
-      double f = (kT - kTs_[i]) / (kTs_[i + 1] - kTs_[i]);
-      if (f > prn(seed))
-        ++i;
+      // If current kT outside of the bounds of available, snap to the bound
+      if (kT < kTs_.front()) {
+        i = 0;
+      } else if (kT > kTs_.back()) {
+        i = kTs_.size() - 1;
+      } else {
+        // Find temperatures that bound the actual temperature
+        while (kTs_[i + 1] < kT && i + 1 < n - 1)
+          ++i;
+        // Randomly sample between temperature i and i+1
+        double f = (kT - kTs_[i]) / (kTs_[i + 1] - kTs_[i]);
+        if (f > prn(seed))
+          ++i;
+      }
     }
   }