From ef4947f451bf1512f27258d1d126203b38c66c07 Mon Sep 17 00:00:00 2001 From: Guillaume Giudicelli Date: Thu, 29 Aug 2024 15:31:00 -0600 Subject: [PATCH] Address review: consolidate code between mgxs.cpp and nuclide.cpp for reading temperatures --- include/openmc/settings.h | 5 ++ src/mgxs.cpp | 89 ++----------------------- src/nuclide.cpp | 127 +----------------------------------- src/settings.cpp | 132 ++++++++++++++++++++++++++++++++++++++ 4 files changed, 144 insertions(+), 209 deletions(-) diff --git a/include/openmc/settings.h b/include/openmc/settings.h index a95f1ced9f1..fecc737eefc 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -167,6 +167,11 @@ void read_settings_xml(); //! \param[in] root XML node for 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& 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 a6012225348..25b1bd377fc 100644 --- a/src/mgxs.cpp +++ b/src/mgxs.cpp @@ -71,8 +71,7 @@ void Mgxs::metadata_from_hdf5(hid_t xs_id, const vector& temperature, dset_names[i] = new char[151]; } get_datasets(kT_group, dset_names); - vector shape = {num_temps}; - xt::xarray temps_available(shape); + vector temps_available(num_temps); for (int i = 0; i < num_temps; i++) { read_double(kT_group, dset_names[i], &temps_available[i], true); @@ -85,89 +84,9 @@ void Mgxs::metadata_from_hdf5(hid_t xs_id, const vector& 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; - } - - // 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(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(); diff --git a/src/nuclide.cpp b/src/nuclide.cpp index f720f848bc2..42b7d7919d2 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -90,131 +90,10 @@ 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(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(temps_to_read.front())); diff --git a/src/settings.cpp b/src/settings.cpp index ba451e219c5..6bece40f0be 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -1090,6 +1090,138 @@ 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(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();