Use temperature bounds in multigroup mode temperature

include/openmc/settings.h

@@ -167,6 +167,11 @@ void read_settings_xml();
 //! \param[in] root XML node for <settings>
 void read_settings_xml(pugi::xml_node root);
 
+//! Select temperatures to read based on what is needed and available
+void select_temperatures(std::string& object_name,
+  const vector<double>& available_temperatures,
+  const vector<double>& requested_temperatures, vector<int>& temps_to_read);
+
 void free_memory_settings();
 
 } // namespace openmc

src/mgxs.cpp

@@ -71,8 +71,7 @@ void Mgxs::metadata_from_hdf5(hid_t xs_id, const vector<double>& temperature,
     dset_names[i] = new char[151];
   }
   get_datasets(kT_group, dset_names);
-  vector<size_t> shape = {num_temps};
-  xt::xarray<double> temps_available(shape);
+  vector<double> temps_available(num_temps);
   for (int i = 0; i < num_temps; i++) {
     read_double(kT_group, dset_names[i], &temps_available[i], true);
 
@@ -85,63 +84,9 @@ void Mgxs::metadata_from_hdf5(hid_t xs_id, const vector<double>& temperature,
   delete[] dset_names;
   std::sort(temps_available.begin(), temps_available.end());
 
-  // If only one temperature is available, lets just use nearest temperature
-  // interpolation
-  if ((num_temps == 1) &&
-      (settings::temperature_method == TemperatureMethod::INTERPOLATION)) {
-    warning("Cross sections for " + strtrim(name) + " are only available " +
-            "at one temperature.  Reverting to the nearest temperature " +
-            "method.");
-    settings::temperature_method = TemperatureMethod::NEAREST;
-  }
-
-  switch (settings::temperature_method) {
-  case TemperatureMethod::NEAREST:
-    // Determine actual temperatures to read
-    for (const auto& T : temperature) {
-      auto i_closest = xt::argmin(xt::abs(temps_available - T))[0];
-      double temp_actual = temps_available[i_closest];
-      if (std::fabs(temp_actual - T) < settings::temperature_tolerance) {
-        if (std::find(temps_to_read.begin(), temps_to_read.end(),
-              std::round(temp_actual)) == temps_to_read.end()) {
-          temps_to_read.push_back(std::round(temp_actual));
-        }
-      } else {
-        fatal_error(fmt::format(
-          "MGXS library does not contain cross sections "
-          "for {} at or near {} K. Available temperatures "
-          "are {} K. Consider making use of openmc.Settings.temperature "
-          "to specify how intermediate temperatures are treated.",
-          in_name, std::round(T), concatenate(temps_available)));
-      }
-    }
-    break;
-
-  case TemperatureMethod::INTERPOLATION:
-    for (int i = 0; i < temperature.size(); i++) {
-      for (int j = 0; j < num_temps; j++) {
-        if (j == (num_temps - 1)) {
-          fatal_error("MGXS Library does not contain cross sections for " +
-                      in_name + " at temperatures that bound " +
-                      std::to_string(std::round(temperature[i])));
-        }
-        if ((temps_available[j] <= temperature[i]) &&
-            (temperature[i] < temps_available[j + 1])) {
-          if (std::find(temps_to_read.begin(), temps_to_read.end(),
-                temps_available[j]) == temps_to_read.end()) {
-            temps_to_read.push_back(temps_available[j]);
-          }
-
-          if (std::find(temps_to_read.begin(), temps_to_read.end(),
-                temps_available[j + 1]) == temps_to_read.end()) {
-            temps_to_read.push_back(temps_available[j + 1]);
-          }
-          break;
-        }
-      }
-    }
-  }
-  std::sort(temps_to_read.begin(), temps_to_read.end());
+  // Code shared with nuclide, placed in settings.cpp
+  openmc::select_temperatures(
+    in_name, temps_available, temperature, temps_to_read);
 
   // Get the library's temperatures
   int n_temperature = temps_to_read.size();

src/nuclide.cpp

@@ -90,131 +90,10 @@ Nuclide::Nuclide(hid_t group, const vector<double>& temperature)
   }
   std::sort(temps_available.begin(), temps_available.end());
 
-  // If only one temperature is available, revert to nearest temperature
-  if (temps_available.size() == 1 &&
-      settings::temperature_method == TemperatureMethod::INTERPOLATION) {
-    if (mpi::master) {
-      warning("Cross sections for " + name_ +
-              " are only available at one "
-              "temperature. Reverting to nearest temperature method.");
-    }
-    settings::temperature_method = TemperatureMethod::NEAREST;
-  }
-
-  // Determine actual temperatures to read -- start by checking whether a
-  // temperature range was given (indicated by T_max > 0), in which case all
-  // temperatures in the range are loaded irrespective of what temperatures
-  // actually appear in the model
+  // Code shared with mgxs, placed in settings.cpp
   vector<int> temps_to_read;
-  int n = temperature.size();
-  double T_min = n > 0 ? settings::temperature_range[0] : 0.0;
-  double T_max = n > 0 ? settings::temperature_range[1] : INFTY;
-  if (T_max > 0.0) {
-    // Determine first available temperature below or equal to T_min
-    auto T_min_it =
-      std::upper_bound(temps_available.begin(), temps_available.end(), T_min);
-    if (T_min_it != temps_available.begin())
-      --T_min_it;
-
-    // Determine first available temperature above or equal to T_max
-    auto T_max_it =
-      std::lower_bound(temps_available.begin(), temps_available.end(), T_max);
-    if (T_max_it != temps_available.end())
-      ++T_max_it;
-
-    // Add corresponding temperatures to vector
-    for (auto it = T_min_it; it != T_max_it; ++it) {
-      temps_to_read.push_back(std::round(*it));
-    }
-  }
-
-  switch (settings::temperature_method) {
-  case TemperatureMethod::NEAREST:
-    // Find nearest temperatures
-    for (double T_desired : temperature) {
-
-      // Determine closest temperature
-      double min_delta_T = INFTY;
-      double T_actual = 0.0;
-      for (auto T : temps_available) {
-        double delta_T = std::abs(T - T_desired);
-        if (delta_T < min_delta_T) {
-          T_actual = T;
-          min_delta_T = delta_T;
-        }
-      }
-
-      if (std::abs(T_actual - T_desired) < settings::temperature_tolerance) {
-        if (!contains(temps_to_read, std::round(T_actual))) {
-          temps_to_read.push_back(std::round(T_actual));
-
-          // Write warning for resonance scattering data if 0K is not available
-          if (std::abs(T_actual - T_desired) > 0 && T_desired == 0 &&
-              mpi::master) {
-            warning(name_ +
-                    " does not contain 0K data needed for resonance "
-                    "scattering options selected. Using data at " +
-                    std::to_string(T_actual) + " K instead.");
-          }
-        }
-      } else {
-        fatal_error(fmt::format(
-          "Nuclear data library does not contain cross sections "
-          "for {}  at or near {} K. Available temperatures "
-          "are {} K. Consider making use of openmc.Settings.temperature "
-          "to specify how intermediate temperatures are treated.",
-          name_, std::to_string(T_desired), concatenate(temps_available)));
-      }
-    }
-    break;
-
-  case TemperatureMethod::INTERPOLATION:
-    // If temperature interpolation or multipole is selected, get a list of
-    // bounding temperatures for each actual temperature present in the model
-    for (double T_desired : temperature) {
-      bool found_pair = false;
-      for (int j = 0; j < temps_available.size() - 1; ++j) {
-        if (temps_available[j] <= T_desired &&
-            T_desired < temps_available[j + 1]) {
-          int T_j = temps_available[j];
-          int T_j1 = temps_available[j + 1];
-          if (!contains(temps_to_read, T_j)) {
-            temps_to_read.push_back(T_j);
-          }
-          if (!contains(temps_to_read, T_j1)) {
-            temps_to_read.push_back(T_j1);
-          }
-          found_pair = true;
-        }
-      }
-
-      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(temps_available.front());
-          }
-          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(temps_available.back());
-          }
-          continue;
-        }
-        fatal_error(
-          "Nuclear data library does not contain cross sections for " + name_ +
-          " at temperatures that bound " + std::to_string(T_desired) + " K.");
-      }
-    }
-    break;
-  }
-
-  // Sort temperatures to read
-  std::sort(temps_to_read.begin(), temps_to_read.end());
+  openmc::select_temperatures(
+    name_, temps_available, temperature, temps_to_read);
 
   data::temperature_min =
     std::min(data::temperature_min, static_cast<double>(temps_to_read.front()));

src/settings.cpp

@@ -1090,6 +1090,138 @@ void read_settings_xml(pugi::xml_node root)
   }
 }
 
+void select_temperatures(std::string& object_name,
+  const vector<double>& available_temperatures,
+  const vector<double>& requested_temperatures, vector<int>& temps_to_read)
+{
+  // If only one temperature is available, revert to nearest temperature
+  if (available_temperatures.size() == 1 &&
+      settings::temperature_method == TemperatureMethod::INTERPOLATION) {
+    if (mpi::master) {
+      warning("Cross sections for " + object_name +
+              " are only available at one "
+              "temperature. Reverting to nearest temperature method.");
+    }
+    settings::temperature_method = TemperatureMethod::NEAREST;
+  }
+
+  // Determine actual temperatures to read -- start by checking whether a
+  // temperature range was given (indicated by T_max > 0), in which case all
+  // temperatures in the range are loaded irrespective of what temperatures
+  // actually appear in the model
+  int n = requested_temperatures.size();
+  double T_min = n > 0 ? settings::temperature_range[0] : 0.0;
+  double T_max = n > 0 ? settings::temperature_range[1] : INFTY;
+  if (T_max > 0.0) {
+    // Determine first available temperature below or equal to T_min
+    auto T_min_it = std::upper_bound(
+      available_temperatures.begin(), available_temperatures.end(), T_min);
+    if (T_min_it != available_temperatures.begin())
+      --T_min_it;
+
+    // Determine first available temperature above or equal to T_max
+    auto T_max_it = std::lower_bound(
+      available_temperatures.begin(), available_temperatures.end(), T_max);
+    if (T_max_it != available_temperatures.end())
+      ++T_max_it;
+
+    // Add corresponding temperatures to vector
+    for (auto it = T_min_it; it != T_max_it; ++it) {
+      temps_to_read.push_back(std::round(*it));
+    }
+  }
+
+  switch (settings::temperature_method) {
+  case TemperatureMethod::NEAREST:
+    // Find nearest temperatures
+    for (double T_desired : requested_temperatures) {
+
+      // Determine closest temperature
+      double min_delta_T = INFTY;
+      double T_actual = 0.0;
+      for (auto T : available_temperatures) {
+        double delta_T = std::abs(T - T_desired);
+        if (delta_T < min_delta_T) {
+          T_actual = T;
+          min_delta_T = delta_T;
+        }
+      }
+
+      if (std::abs(T_actual - T_desired) < settings::temperature_tolerance) {
+        if (!contains(temps_to_read, std::round(T_actual))) {
+          temps_to_read.push_back(std::round(T_actual));
+
+          // Write warning for resonance scattering data if 0K is not available
+          if (std::abs(T_actual - T_desired) > 0 && T_desired == 0 &&
+              mpi::master) {
+            warning(object_name +
+                    " does not contain 0K data needed for resonance "
+                    "scattering options selected. Using data at " +
+                    std::to_string(T_actual) + " K instead.");
+          }
+        }
+      } else {
+        fatal_error(fmt::format(
+          "Nuclear data library does not contain cross sections "
+          "for {}  at or near {} K. Available temperatures "
+          "are {} K. Consider making use of openmc.Settings.temperature "
+          "to specify how intermediate temperatures are treated.",
+          object_name, std::to_string(T_desired),
+          concatenate(available_temperatures)));
+      }
+    }
+    break;
+
+  case TemperatureMethod::INTERPOLATION:
+    // If temperature interpolation or multipole is selected, get a list of
+    // bounding temperatures for each actual temperature present in the model
+    for (double T_desired : requested_temperatures) {
+      bool found_pair = false;
+      for (int j = 0; j < available_temperatures.size() - 1; ++j) {
+        if (available_temperatures[j] <= T_desired &&
+            T_desired < available_temperatures[j + 1]) {
+          int T_j = available_temperatures[j];
+          int T_j1 = available_temperatures[j + 1];
+          if (!contains(temps_to_read, T_j)) {
+            temps_to_read.push_back(T_j);
+          }
+          if (!contains(temps_to_read, T_j1)) {
+            temps_to_read.push_back(T_j1);
+          }
+          found_pair = true;
+        }
+      }
+
+      if (!found_pair) {
+        // If no pairs found, check if the desired temperature falls just
+        // outside of data
+        if (std::abs(T_desired - available_temperatures.front()) <=
+            settings::temperature_tolerance) {
+          if (!contains(temps_to_read, available_temperatures.front())) {
+            temps_to_read.push_back(available_temperatures.front());
+          }
+          continue;
+        }
+        if (std::abs(T_desired - available_temperatures.back()) <=
+            settings::temperature_tolerance) {
+          if (!contains(temps_to_read, available_temperatures.back())) {
+            temps_to_read.push_back(available_temperatures.back());
+          }
+          continue;
+        }
+        fatal_error(
+          "Nuclear data library does not contain cross sections for " +
+          object_name + " at temperatures that bound " +
+          std::to_string(T_desired) + " K.");
+      }
+    }
+    break;
+  }
+
+  // Sort temperatures to read
+  std::sort(temps_to_read.begin(), temps_to_read.end());
+}
+
 void free_memory_settings()
 {
   settings::statepoint_batch.clear();