Skip to content

Commit

Permalink
Address review: consolidate code between mgxs.cpp and nuclide.cpp for…
Browse files Browse the repository at this point in the history 
… reading temperatures
  • Loading branch information
@GiudGiud
GiudGiud committed Aug 29, 2024
1 parent f3b6afe commit 3d707b1
Show file tree
Hide file tree
Showing 4 changed files with 140 additions and 202 deletions.
6 changes: 6 additions & 0 deletions include/openmc/settings.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -167,6 +167,12 @@ void read_settings_xml();
//! \param[in] root XML node for <settings>
void read_settings_xml(pugi::xml_node root);

//! Read temperatures from XML
void read_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
Expand Down
84 changes: 3 additions & 81 deletions src/mgxs.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -72,7 +72,7 @@ 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);
vector<double> available_temps(shape.size());
for (int i = 0; i < num_temps; i++) {
read_double(kT_group, dset_names[i], &available_temps[i], true);

Expand All @@ -85,86 +85,8 @@ void Mgxs::metadata_from_hdf5(hid_t xs_id, const vector<double>& 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();
Expand Down
123 changes: 2 additions & 121 deletions src/nuclide.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -90,128 +90,9 @@ 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(
"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<double>(temps_to_read.front()));
Expand Down
129 changes: 129 additions & 0 deletions src/settings.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1052,6 +1052,135 @@ 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(
"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();
Expand Down

0 comments on commit 3d707b1

Please sign in to comment.