Fix bug in selecting data temperatures to read with temperature inter…

…polation (#2734)

src/nuclide.cpp

@@ -176,14 +176,14 @@ Nuclide::Nuclide(hid_t group, const vector<double>& temperature)
           if (!contains(temps_to_read, temps_available.front())) {
             temps_to_read.push_back(std::round(temps_available.front()));
           }
-          break;
+          continue;
         }
         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;
+          continue;
         }
         fatal_error(
           "Nuclear data library does not contain cross sections for " + name_ +

tests/unit_tests/test_temp_interp.py

@@ -238,3 +238,38 @@ def test_temperature_interpolation_tolerance(model):
     # All calculated k-effectives should be equal
     assert default_k == pytest.approx(interpolated_k)
     assert interpolated_k == pytest.approx(cell_k)
+
+
+def test_temperature_slightly_above(run_in_tmpdir):
+    """In this test, we have two materials at temperatures close to actual data
+    temperatures. However, one is slightly above the highest temperature which
+    invokes separate logic. The k-effective value should be somewhere between
+    k=2 (if the temperature were only 600 K) and k=1 (if the temperature were
+    only 900 K)."""
+
+    make_fake_cross_section()
+
+    model = openmc.Model()
+    mat1 = openmc.Material()
+    mat1.add_nuclide('U235', 1.0)
+    mat1.temperature = 900.1
+    mat2 = openmc.Material()
+    mat2.add_nuclide('U235', 1.0)
+    mat2.temperature = 600.0
+    model.materials.extend([mat1, mat2])
+    model.materials.cross_sections = str(Path('cross_sections_fake.xml').resolve())
+
+    sph1 = openmc.Sphere(r=1.0)
+    sph2 = openmc.Sphere(r=4.0, boundary_type='reflective')
+    cell1 = openmc.Cell(fill=mat1, region=-sph1)
+    cell2 = openmc.Cell(fill=mat2, region=+sph1 & -sph2)
+    model.geometry = openmc.Geometry([cell1, cell2])
+
+    model.settings.particles = 1000
+    model.settings.inactive = 0
+    model.settings.batches = 10
+    model.settings.temperature = {'method': 'interpolation'}
+
+    sp_filename = model.run()
+    with openmc.StatePoint(sp_filename) as sp:
+        assert 1.1 < sp.keff.n < 1.9