Nuclide temperatures - solution to issue #3102 (#3110)

This commit is contained in:
Zoe Prieto 2024-08-15 12:10:26 -03:00 committed by GitHub
parent 9483cce0bc
commit 10c511a0e2
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
5 changed files with 61 additions and 36 deletions

View file

@ -1,6 +1,7 @@
#ifndef OPENMC_STRING_UTILS_H
#define OPENMC_STRING_UTILS_H
#include <sstream>
#include <string>
#include "openmc/vector.h"
@ -15,7 +16,7 @@ std::string to_element(const std::string& name);
void to_lower(std::string& str);
int word_count(std::string const& str);
int word_count(const std::string& str);
vector<std::string> split(const std::string& in);
@ -23,5 +24,20 @@ bool ends_with(const std::string& value, const std::string& ending);
bool starts_with(const std::string& value, const std::string& beginning);
template<typename T>
inline std::string concatenate(const T& values, const std::string& del = ", ")
{
if (values.size() == 0)
return "";
std::stringstream oss;
auto it = values.begin();
oss << *it++;
while (it != values.end()) {
oss << del << *it++;
}
return oss.str();
}
} // namespace openmc
#endif // OPENMC_STRING_UTILS_H

View file

@ -72,18 +72,18 @@ void Mgxs::metadata_from_hdf5(hid_t xs_id, const vector<double>& temperature,
}
get_datasets(kT_group, dset_names);
vector<size_t> shape = {num_temps};
xt::xarray<double> available_temps(shape);
xt::xarray<double> temps_available(shape);
for (int i = 0; i < num_temps; i++) {
read_double(kT_group, dset_names[i], &available_temps[i], true);
read_double(kT_group, dset_names[i], &temps_available[i], true);
// convert eV to Kelvin
available_temps[i] /= K_BOLTZMANN;
temps_available[i] = std::round(temps_available[i] / K_BOLTZMANN);
// Done with dset_names, so delete it
delete[] dset_names[i];
}
delete[] dset_names;
std::sort(available_temps.begin(), available_temps.end());
std::sort(temps_available.begin(), temps_available.end());
// If only one temperature is available, lets just use nearest temperature
// interpolation
@ -99,17 +99,20 @@ void Mgxs::metadata_from_hdf5(hid_t xs_id, const vector<double>& temperature,
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];
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.",
in_name, std::round(T)));
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;
@ -122,16 +125,16 @@ void Mgxs::metadata_from_hdf5(hid_t xs_id, const vector<double>& temperature,
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 ((temps_available[j] <= temperature[i]) &&
(temperature[i] < temps_available[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]));
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(),
std::round(available_temps[j + 1])) == temps_to_read.end()) {
temps_to_read.push_back(std::round((int)available_temps[j + 1]));
temps_available[j + 1]) == temps_to_read.end()) {
temps_to_read.push_back(temps_available[j + 1]);
}
break;
}

View file

@ -86,7 +86,7 @@ Nuclide::Nuclide(hid_t group, const vector<double>& temperature)
for (const auto& name : dset_names) {
double T;
read_dataset(kT_group, name.c_str(), T);
temps_available.push_back(T / K_BOLTZMANN);
temps_available.push_back(std::round(T / K_BOLTZMANN));
}
std::sort(temps_available.begin(), temps_available.end());
@ -158,9 +158,12 @@ Nuclide::Nuclide(hid_t group, const vector<double>& temperature)
}
}
} else {
fatal_error(
"Nuclear data library does not contain cross sections for " + name_ +
" at or near " + std::to_string(T_desired) + " K.");
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;
@ -173,8 +176,8 @@ Nuclide::Nuclide(hid_t group, const vector<double>& temperature)
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]);
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);
}
@ -191,14 +194,14 @@ Nuclide::Nuclide(hid_t group, const vector<double>& temperature)
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()));
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(std::round(temps_available.back()));
temps_to_read.push_back(temps_available.back());
}
continue;
}

View file

@ -2,7 +2,6 @@
#include <algorithm> // for equal
#include <cctype> // for tolower, isspace
#include <sstream>
namespace openmc {
@ -36,7 +35,7 @@ void to_lower(std::string& str)
str[i] = std::tolower(str[i]);
}
int word_count(std::string const& str)
int word_count(const std::string& str)
{
std::stringstream stream(str);
std::string dum;

View file

@ -19,6 +19,7 @@
#include "openmc/secondary_correlated.h"
#include "openmc/secondary_thermal.h"
#include "openmc/settings.h"
#include "openmc/string_utils.h"
namespace openmc {
@ -59,7 +60,7 @@ ThermalScattering::ThermalScattering(
// Read temperature value
double T;
read_dataset(kT_group, dset_names[i].data(), T);
temps_available[i] = T / K_BOLTZMANN;
temps_available[i] = std::round(T / K_BOLTZMANN);
}
std::sort(temps_available.begin(), temps_available.end());
@ -89,9 +90,12 @@ ThermalScattering::ThermalScattering(
temps_to_read.push_back(std::round(temp_actual));
}
} else {
fatal_error(fmt::format("Nuclear data library does not contain cross "
"sections for {} at or near {} K.",
name_, std::round(T)));
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::round(T), concatenate(temps_available)));
}
}
break;
@ -103,8 +107,8 @@ ThermalScattering::ThermalScattering(
bool found = false;
for (int j = 0; j < temps_available.size() - 1; ++j) {
if (temps_available[j] <= T && T < temps_available[j + 1]) {
int T_j = std::round(temps_available[j]);
int T_j1 = std::round(temps_available[j + 1]);
int T_j = temps_available[j];
int T_j1 = temps_available[j + 1];
if (std::find(temps_to_read.begin(), temps_to_read.end(), T_j) ==
temps_to_read.end()) {
temps_to_read.push_back(T_j);
@ -122,14 +126,14 @@ ThermalScattering::ThermalScattering(
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]));
temps_available[0]) == temps_to_read.end()) {
temps_to_read.push_back(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]));
temps_available[n - 1]) == temps_to_read.end()) {
temps_to_read.push_back(temps_available[n - 1]);
}
} else {
fatal_error(