From 3d707b12a6efeb8fa6129084093c3c7d56125612 Mon Sep 17 00:00:00 2001 From: Guillaume Giudicelli Date: Thu, 29 Aug 2024 15:14:11 -0600 Subject: [PATCH] Address review: consolidate code between mgxs.cpp and nuclide.cpp for reading temperatures --- include/openmc/settings.h | 6 ++ src/mgxs.cpp | 84 +------------------------ src/nuclide.cpp | 123 +----------------------------------- src/settings.cpp | 129 ++++++++++++++++++++++++++++++++++++++ 4 files changed, 140 insertions(+), 202 deletions(-) diff --git a/include/openmc/settings.h b/include/openmc/settings.h index a95f1ced9f1..c8bf48027d7 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -167,6 +167,12 @@ void read_settings_xml(); //! \param[in] root XML node for void read_settings_xml(pugi::xml_node root); +//! Read temperatures from XML +void read_temperatures(std::string & object_name, + const vector& available_temperatures, + const vector& requested_temperatures, + vector & temps_to_read); + void free_memory_settings(); } // namespace openmc diff --git a/src/mgxs.cpp b/src/mgxs.cpp index 7646c1834bf..7d83ecbc617 100644 --- a/src/mgxs.cpp +++ b/src/mgxs.cpp @@ -72,7 +72,7 @@ void Mgxs::metadata_from_hdf5(hid_t xs_id, const vector& temperature, } get_datasets(kT_group, dset_names); vector shape = {num_temps}; - xt::xarray available_temps(shape); + vector available_temps(shape.size()); for (int i = 0; i < num_temps; i++) { read_double(kT_group, dset_names[i], &available_temps[i], true); @@ -85,86 +85,8 @@ void Mgxs::metadata_from_hdf5(hid_t xs_id, const vector& temperature, delete[] dset_names; std::sort(available_temps.begin(), available_temps.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; - } - - // 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 = 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(available_temps.begin(), available_temps.end(), T_min); - if (T_min_it != available_temps.begin()) - --T_min_it; - - // Determine first available temperature above or equal to T_max - auto T_max_it = - std::lower_bound(available_temps.begin(), available_temps.end(), T_max); - if (T_max_it != available_temps.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: - // Determine actual temperatures to read - for (const auto& T : temperature) { - auto i_closest = xt::argmin(xt::abs(available_temps - T))[0]; - double temp_actual = available_temps[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.", - in_name, std::round(T))); - } - } - 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 ((available_temps[j] <= temperature[i]) && - (temperature[i] < available_temps[j + 1])) { - if (std::find(temps_to_read.begin(), temps_to_read.end(), - std::round(available_temps[j])) == temps_to_read.end()) { - temps_to_read.push_back(std::round((int)available_temps[j])); - } - - if (std::find(temps_to_read.begin(), temps_to_read.end(), - std::round(available_temps[j + 1])) == temps_to_read.end()) { - temps_to_read.push_back(std::round((int)available_temps[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, available_temps, temperature, temps_to_read); // Get the library's temperatures int n_temperature = temps_to_read.size(); diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 04ec110738c..87e24bcad65 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -90,128 +90,9 @@ Nuclide::Nuclide(hid_t group, const vector& 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 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( - "Nuclear data library does not contain cross sections for " + name_ + - " at or near " + std::to_string(T_desired) + " K."); - } - } - 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 = std::round(temps_available[j]); - int T_j1 = std::round(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(std::round(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(std::round(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(temps_to_read.front())); diff --git a/src/settings.cpp b/src/settings.cpp index 6c3226b04a6..affe3c2c526 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -1052,6 +1052,135 @@ void read_settings_xml(pugi::xml_node root) } } +void select_temperatures(std::string & object_name, + const vector& available_temperatures, + const vector& requested_temperatures, + vector & 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( + "Nuclear data library does not contain cross sections for " + object_name + + " at or near " + std::to_string(T_desired) + " K."); + } + } + 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 = std::round(available_temperatures[j]); + int T_j1 = std::round(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(std::round(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(std::round(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();