From 3e2ddceaac8484503c57fef4248af50399067b96 Mon Sep 17 00:00:00 2001 From: josh Date: Tue, 11 Oct 2022 03:03:50 +0000 Subject: [PATCH] Add tolerance for interp temps outside of bounds --- docs/source/io_formats/settings.rst | 11 +++++++-- openmc/settings.py | 16 +++++++------ src/cell.cpp | 4 ++-- src/nuclide.cpp | 35 +++++++++++++++++++++++++++++ src/thermal.cpp | 35 +++++++++++++++++++++-------- 5 files changed, 81 insertions(+), 20 deletions(-) diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 1a8f63716f..7a794bf003 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -832,7 +832,9 @@ cell, the nearest temperature at which cross sections are given is to be applied, within a given tolerance (see :ref:`temperature_tolerance`). A value of "interpolation" indicates that cross sections are to be linear-linear interpolated between temperatures at which nuclear data are present (see -:ref:`temperature_treatment`). +:ref:`temperature_treatment`). With the "interpolation" method, temperatures +outside of the bounds of the nuclear data may be accepted, provided they still +fall within the tolerance (see :ref:`temperature_tolerance`). *Default*: "nearest" @@ -871,7 +873,12 @@ The ```` element specifies a tolerance in Kelvin that is to be applied when the "nearest" temperature method is used. For example, if a cell temperature is 340 K and the tolerance is 15 K, then the closest temperature in the range of 325 K to 355 K will be used to evaluate cross -sections. +sections. If the ```` is "interpolation", the tolerance +specified applies to cell temperatures outside of the data bounds. For example, +If a cell is specified at 695K, a tolerance of 15K and data only available at +700K and 1000K, the cell's cross sections will be evaluated at 700K, since +desired temperature of 695K is within the tolerance of the actual data despite +not being bounded on both sides. *Default*: 10 K diff --git a/openmc/settings.py b/openmc/settings.py index ad5a9b28ae..9d9b421ca3 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -186,13 +186,15 @@ class Settings: 'default', 'method', 'range', 'tolerance', and 'multipole'. The value for 'default' should be a float representing the default temperature in Kelvin. The value for 'method' should be 'nearest' or 'interpolation'. - If the method is 'nearest', 'tolerance' indicates a range of temperature - within which cross sections may be used. The value for 'range' should be - a pair a minimum and maximum temperatures which are used to indicate - that cross sections be loaded at all temperatures within the - range. 'multipole' is a boolean indicating whether or not the windowed - multipole method should be used to evaluate resolved resonance cross - sections. + If the method is 'nearest', 'tolerance' indicates a range of + temperature within which cross sections may be used. If the method is + 'interpolation', 'tolerance' indicates the range of temperatures outside + of the available cross section temperatures where cross sections will + evaluate to the nearer bound. The value for 'range' should be a pair a + minimum and maximum temperatures which are used to indicate that cross + sections be loaded at all temperatures within the range. 'multipole' is + a boolean indicating whether or not the windowed multipole method should + be used to evaluate resolved resonance cross sections. trace : tuple or list Show detailed information about a single particle, indicated by three integers: the batch number, generation number, and particle number diff --git a/src/cell.cpp b/src/cell.cpp index 4ff35498c6..d69f79a3e7 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -334,10 +334,10 @@ double Cell::temperature(int32_t instance) const void Cell::set_temperature(double T, int32_t instance, bool set_contained) { if (settings::temperature_method == TemperatureMethod::INTERPOLATION) { - if (T < data::temperature_min) { + if (T < data::temperature_min - settings::temperature_tolerance) { throw std::runtime_error {"Temperature is below minimum temperature at " "which data is available."}; - } else if (T > data::temperature_max) { + } else if (T > data::temperature_max + settings::temperature_tolerance) { throw std::runtime_error {"Temperature is above maximum temperature at " "which data is available."}; } diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 564e27d446..4c1706509b 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -169,6 +169,22 @@ Nuclide::Nuclide(hid_t group, const vector& temperature) } if (!found_pair) { + // If no pairs found, check if the desired temperature falls just + // outside of data + if (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())); + } + break; + } + if (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())); + } + break; + } fatal_error( "Nuclear data library does not contain cross sections for " + name_ + " at temperatures that bound " + std::to_string(T_desired) + " K."); @@ -646,6 +662,16 @@ void Nuclide::calculate_xs( } break; case TemperatureMethod::INTERPOLATION: + // If current kT outside of the bounds of available, snap to the bound + if (kT < kTs_.front()) { + i_temp = 0; + break; + } + if (kT > kTs_.back()) { + i_temp = kTs_.size() - 1; + break; + } + // Find temperatures that bound the actual temperature for (i_temp = 0; i_temp < kTs_.size() - 1; ++i_temp) { if (kTs_[i_temp] <= kT && kT < kTs_[i_temp + 1]) @@ -969,6 +995,15 @@ std::pair Nuclide::find_temperature(double T) const } break; case TemperatureMethod::INTERPOLATION: + // If current kT outside of the bounds of available, snap to the bound + if (kT < kTs_.front()) { + i_temp = 0; + break; + } + if (kT > kTs_.back()) { + i_temp = kTs_.size() - 1; + break; + } // Find temperatures that bound the actual temperature while (kTs_[i_temp + 1] < kT && i_temp + 1 < n - 1) ++i_temp; diff --git a/src/thermal.cpp b/src/thermal.cpp index 1a9592d0af..66e1fb0095 100644 --- a/src/thermal.cpp +++ b/src/thermal.cpp @@ -117,6 +117,16 @@ ThermalScattering::ThermalScattering( } } if (!found) { + // If no pairs found, check if the desired temperature falls within + // bounds' tolerance + if (T - temps_available[0] <= -settings::temperature_tolerance) { + temps_to_read.push_back(std::round(temps_available[0])); + break; + } + if (T - temps_available[n - 1] <= settings::temperature_tolerance) { + temps_to_read.push_back(std::round(temps_available[n - 1])); + break; + } fatal_error( fmt::format("Nuclear data library does not contain cross " "sections for {} at temperatures that bound {} K.", @@ -159,20 +169,27 @@ void ThermalScattering::calculate_xs(double E, double sqrtkT, int* i_temp, auto n = kTs_.size(); if (n > 1) { - // Find temperatures that bound the actual temperature - while (kTs_[i + 1] < kT && i + 1 < n - 1) - ++i; - if (settings::temperature_method == TemperatureMethod::NEAREST) { + while (kTs_[i + 1] < kT && i + 1 < n - 1) + ++i; // Pick closer of two bounding temperatures if (kT - kTs_[i] > kTs_[i + 1] - kT) ++i; - } else { - // Randomly sample between temperature i and i+1 - double f = (kT - kTs_[i]) / (kTs_[i + 1] - kTs_[i]); - if (f > prn(seed)) - ++i; + // If current kT outside of the bounds of available, snap to the bound + if (kT < kTs_.front()) { + i = 0; + } else if (kT > kTs_.back()) { + i = kTs_.size() - 1; + } else { + // Find temperatures that bound the actual temperature + while (kTs_[i + 1] < kT && i + 1 < n - 1) + ++i; + // Randomly sample between temperature i and i+1 + double f = (kT - kTs_[i]) / (kTs_[i + 1] - kTs_[i]); + if (f > prn(seed)) + ++i; + } } }