Merge pull request #2265 from joshmay1/temp_interp_tol

Adds a tolerance for temperatures slightly out of bounds when interpolating

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 is only available
+at 700K and 1000K, the cell's cross sections will be evaluated at 700K, since
+the 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

@@ -187,13 +187,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 of 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 of
+        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

@@ -99,12 +99,12 @@ 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) {
-      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 (T < (data::temperature_min - settings::temperature_tolerance)) {
+      throw std::runtime_error {fmt::format("Temperature of {} K is below minimum temperature at "
+                                "which data is available of {} K.", T, data::temperature_min)};
+    } else if (T > (data::temperature_max + settings::temperature_tolerance)) {
+      throw std::runtime_error {fmt::format("Temperature of {} K is above maximum temperature at "
+                                "which data is available of {} K.", T, data::temperature_max)};
     }
   }
 

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 (std::abs(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 (std::abs(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,10 +117,24 @@ ThermalScattering::ThermalScattering(
         }
       }
       if (!found) {
-        fatal_error(
+        // If no pairs found, check if the desired temperature falls within
+        // bounds' tolerance
+        if (std::abs(T - temps_available[0]) <= settings::temperature_tolerance){
+          if (std::find(temps_to_read.begin(), temps_to_read.end(), std::round(temps_available[0])) ==
+              temps_to_read.end()) {
+          temps_to_read.push_back(std::round(temps_available[0]));
+        }}
+        else if (std::abs(T - temps_available[n - 1]) <= settings::temperature_tolerance){
+          if (std::find(temps_to_read.begin(), temps_to_read.end(), std::round(temps_available[n - 1])) ==
+              temps_to_read.end()){
+          temps_to_read.push_back(std::round(temps_available[n - 1]));
+        }}
+        else {
+          fatal_error(
           fmt::format("Nuclear data library does not contain cross "
                       "sections for {} at temperatures that bound {} K.",
             name_, std::round(T)));
+            }
       }
     }
   }
@@ -159,20 +173,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;
+      }
     }
   }
 

tests/unit_tests/test_temp_interp.py

@@ -10,7 +10,7 @@ import pytest
 
 
 def make_fake_cross_section():
-    """Create fake U235 nuclide
+    """Create fake U235 nuclide with a fake thermal scattering library attached
 
     This nuclide is designed to have k_inf=1 at 300 K, k_inf=2 at 600 K, and
     k_inf=1 at 900 K. The absorption cross section is also constant with
@@ -81,8 +81,62 @@ def make_fake_cross_section():
     # Export HDF5 file
     u235_fake.export_to_hdf5('U235_fake.h5', 'w')
 
+    # Create a fake thermal scattering library attached to the fake U235 data
+    c_U_fake = openmc.data.ThermalScattering("c_U_fake", 1.9968, 4.9, [0.0253])
+    c_U_fake.nuclides = ['U235']
+
+    # Create elastic reaction
+    bragg_edges = [0.00370672, 0.00494229]
+    factors = [0.00375735, 0.01386287]
+    coherent_xs = openmc.data.CoherentElastic(bragg_edges, factors)
+    incoherent_xs_294 = openmc.data.Tabulated1D([0.00370672, 0.00370672], [0.00370672, 0.00370672])
+    elastic_xs_base = openmc.data.Sum((coherent_xs, incoherent_xs_294))
+    elastic_xs = {'294K': elastic_xs_base, '600K': elastic_xs_base}
+    coherent_dist = openmc.data.CoherentElasticAE(coherent_xs)
+    incoherent_dist_294 = openmc.data.IncoherentElasticAEDiscrete([
+        [-0.6, -0.18, 0.18, 0.6], [-0.6, -0.18, 0.18, 0.6]
+    ])
+    incoherent_dist_600 = openmc.data.IncoherentElasticAEDiscrete([
+        [-0.1, -0.2, 0.2, 0.1], [-0.1, -0.2, 0.2, 0.1]
+    ])
+    elastic_dist = {
+        '294K': openmc.data.MixedElasticAE(coherent_dist, incoherent_dist_294),
+        '600K': openmc.data.MixedElasticAE(coherent_dist, incoherent_dist_600)
+        }
+    c_U_fake.elastic = openmc.data.ThermalScatteringReaction(elastic_xs, elastic_dist)
+
+    # Create inelastic reaction
+    inelastic_xs = {
+        '294K': openmc.data.Tabulated1D([1.0e-5, 4.9], [13.4, 3.35]),
+        '600K': openmc.data.Tabulated1D([1.0e-2, 10], [1.4, 5])
+        }
+    breakpoints = [3]
+    interpolation = [2]
+    energy = [1.0e-5, 4.3e-2, 4.9]
+    energy_out = [
+        openmc.data.Tabular([0.0002, 0.067, 0.146, 0.366], [0.25, 0.25, 0.25, 0.25]),
+        openmc.data.Tabular([0.0001, 0.009, 0.137, 0.277], [0.25, 0.25, 0.25, 0.25]),
+        openmc.data.Tabular([0.0579, 4.555, 4.803, 4.874], [0.25, 0.25, 0.25, 0.25]),
+    ]
+    for eout in energy_out:
+        eout.normalize()
+        eout.c = eout.cdf()
+    discrete = openmc.stats.Discrete([-0.9, -0.6, -0.3, -0.1, 0.1, 0.3, 0.6, 0.9], [1/8]*8)
+    discrete.c = discrete.cdf()[1:]
+    mu = [[discrete]*4]*3
+    dist = openmc.data.IncoherentInelasticAE(
+        breakpoints, interpolation, energy, energy_out, mu)
+    inelastic_dist = {'294K': dist, '600K': dist}
+    inelastic = openmc.data.ThermalScatteringReaction(inelastic_xs, inelastic_dist)
+    c_U_fake.inelastic = inelastic
+
+    # Export HDF5 file
+    c_U_fake.export_to_hdf5("c_U_fake.h5")
+
+    # Create a data library of the fake nuclide and its thermal scattering data
     lib = openmc.data.DataLibrary()
     lib.register_file('U235_fake.h5')
+    lib.register_file("c_U_fake.h5")
     lib.export_to_xml('cross_sections_fake.xml')
 
 
@@ -119,21 +173,23 @@ def model(tmp_path_factory):
 
 
 @pytest.mark.parametrize(
-    ["method", "temperature", "fission_expected"],
+    ["method", "temperature", "fission_expected", "tolerance"],
     [
-        ("nearest", 300.0, 0.5),
-        ("nearest", 600.0, 1.0),
-        ("nearest", 900.0, 0.5),
-        ("interpolation", 360.0, 0.6),
-        ("interpolation", 450.0, 0.75),
-        ("interpolation", 540.0, 0.9),
-        ("interpolation", 660.0, 0.9),
-        ("interpolation", 750.0, 0.75),
-        ("interpolation", 840.0, 0.6),
+        ("nearest", 300.0, 0.5, 10),
+        ("nearest", 600.0, 1.0, 10),
+        ("nearest", 900.0, 0.5, 10),
+        ("interpolation", 360.0, 0.6, 10),
+        ("interpolation", 450.0, 0.75, 10),
+        ("interpolation", 540.0, 0.9, 10),
+        ("interpolation", 660.0, 0.9, 10),
+        ("interpolation", 750.0, 0.75, 10),
+        ("interpolation", 840.0, 0.6, 10),
+        ("interpolation", 295.0, 0.5, 10),
+        ("interpolation", 990.0, 0.5, 100),
     ]
 )
-def test_interpolation(model, method, temperature, fission_expected):
-    model.settings.temperature = {'method': method, 'default': temperature}
+def test_interpolation(model, method, temperature, fission_expected, tolerance):
+    model.settings.temperature = {'method': method, 'default': temperature, "tolerance": tolerance}
     sp_filename = model.run()
     with openmc.StatePoint(sp_filename) as sp:
         t = sp.tallies[model.tallies[0].id]
@@ -152,3 +208,33 @@ def test_interpolation(model, method, temperature, fission_expected):
             assert k.n == pytest.approx(nu*fission_expected)
         else:
             assert abs(k.n - nu*fission_expected) <= 3*k.s
+
+
+def test_temperature_interpolation_tolerance(model):
+    """Test applying global and cell temperatures with thermal scattering libraries
+    """
+    model.materials[0].add_s_alpha_beta("c_U_fake")
+
+    # Default k-effective, using the thermal scattering data's minimum available temperature
+    model.settings.temperature = {'method': "nearest", 'default': 294, "tolerance": 50}
+    sp_filename = model.run()
+    with openmc.StatePoint(sp_filename) as sp:
+        default_k = sp.keff.n
+
+    # Get k-effective with temperature below the minimum but in interpolation mode
+    model.settings.temperature = {'method': "interpolation", 'default': 255, "tolerance": 50}
+    sp_filename = model.run()
+    with openmc.StatePoint(sp_filename) as sp:
+        interpolated_k = sp.keff.n
+
+    # Get the k-effective with the temperature applied to the cell, instead of globally
+    model.settings.temperature = {'method': "interpolation", 'default': 500, "tolerance": 50}
+    for cell in model.geometry.get_all_cells().values():
+        cell.temperature = 275
+    sp_filename = model.run()
+    with openmc.StatePoint(sp_filename) as sp:
+        cell_k = sp.keff.n
+
+    # All calculated k-effectives should be equal
+    assert default_k == pytest.approx(interpolated_k)
+    assert interpolated_k == pytest.approx(cell_k)