crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 21:55:41 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Use pair instead of struct as return type of find_temperature

Browse source
This commit is contained in:
Paul Romano 2020-08-07 10:51:28 -05:00
parent 0a61bc0cbe
commit 1e436577ea
2 changed files with 27 additions and 34 deletions
Download patch file Download diff file Expand all files Collapse all files

12
include/openmc/nuclide.h
View file

@ -7,6 +7,7 @@
#include <array>
#include <memory> // for unique_ptr
#include <unordered_map>
#include <utility> // for pair
#include <vector>
#include <gsl/gsl>
@ -35,11 +36,6 @@ public:
std::vector<double> energy;
};
struct InterpResult {
gsl::index i; //!< Index in tabulated data
double f; //!< Interpolation factor between i and i+1
};
// Constructors/destructors
Nuclide(hid_t group, const std::vector<double>& temperature);
~Nuclide();
@ -120,7 +116,11 @@ public:
private:
void create_derived(const Function1D* prompt_photons, const Function1D* delayed_photons);
InterpResult find_temperature(double T) const;
//! Determine temperature index and interpolation factor
//
//! \param[in] T Temperature in [K]
//! \return Temperature index and interpolation factor
std::pair<gsl::index, double> find_temperature(double T) const;
static int XS_TOTAL;
static int XS_ABSORPTION;

49
src/nuclide.cpp
View file

@ -910,49 +910,40 @@ void Nuclide::calculate_urr_xs(int i_temp, Particle& p) const
}
Nuclide::InterpResult Nuclide::find_temperature(double T) const
std::pair<gsl::index, double> Nuclide::find_temperature(double T) const
{
Expects(T >= 0.0);
// Make sure value is within bounds
double T_min = kTs_.front() / K_BOLTZMANN;
double T_max = kTs_.back() / K_BOLTZMANN;
if (T < T_min || T > T_max) {
throw std::out_of_range{fmt::format(
"Temperature out of range (min={:.1f}, max={:.1f}).", T_min, T_max
)};
}
// Determine temperature index
InterpResult res;
gsl::index i_temp = 0;
double f = 0.0;
double kT = K_BOLTZMANN * T;
gsl::index n = kTs_.size();
switch (settings::temperature_method) {
case TemperatureMethod::NEAREST:
{
double max_diff = INFTY;
for (gsl::index t = 0; t < kTs_.size(); ++t) {
for (gsl::index t = 0; t < n; ++t) {
double diff = std::abs(kTs_[t] - kT);
if (diff < max_diff) {
res.i = t;
i_temp = t;
max_diff = diff;
}
}
}
res.f = 0.0;
break;
case TemperatureMethod::INTERPOLATION:
// Find temperatures that bound the actual temperature
for (res.i = 0; res.i < kTs_.size() - 1; ++res.i) {
if (kTs_[res.i] <= kT && kT < kTs_[res.i + 1]) break;
}
res.f = (kT - kTs_[res.i]) / (kTs_[res.i + 1] - kTs_[res.i]);
while (kTs_[i_temp + 1] < kT && i_temp + 1 < n - 1) ++i_temp;
// Determine interpolation factor
f = (kT - kTs_[i_temp]) / (kTs_[i_temp + 1] - kTs_[i_temp]);
}
Ensures(res.i >= 0 && res.i < kTs_.size());
Ensures(res.f >= 0.0 && res.f <= 1.0);
Ensures(i_temp >= 0 && i_temp < n);
return res;
return {i_temp, f};
}
double Nuclide::collapse_rate(int MT, double temperature, gsl::span<const double> energy,
@ -967,17 +958,19 @@ double Nuclide::collapse_rate(int MT, double temperature, gsl::span<const double
const auto& rx = reactions_[i_rx];
// Determine temperature index
auto res = this->find_temperature(temperature);
gsl::index i_temp;
double f;
std::tie(i_temp, f) = this->find_temperature(temperature);
// Get reaction rate at lower temperature
const auto& grid_low = grid_[res.i].energy;
double rr_low = rx->collapse_rate(res.i, energy, flux, grid_low);
const auto& grid_low = grid_[i_temp].energy;
double rr_low = rx->collapse_rate(i_temp, energy, flux, grid_low);
if (res.f > 0.0) {
if (f > 0.0) {
// Interpolate between reaction rate at lower and higher temperature
const auto& grid_high = grid_[res.i + 1].energy;
double rr_high = rx->collapse_rate(res.i + 1, energy, flux, grid_high);
return rr_low + res.f*(rr_high - rr_low);
const auto& grid_high = grid_[i_temp + 1].energy;
double rr_high = rx->collapse_rate(i_temp + 1, energy, flux, grid_high);
return rr_low + f*(rr_high - rr_low);
} else {
// If interpolation factor is zero, return reaction rate at lower temperature
return rr_low;