crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 14:15:42 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Updating interpolation function signature.

Browse source
This commit is contained in:
Patrick Shriwise 2022-09-06 22:49:00 -05:00
parent 5168b3392a
commit 35aa433b9a
2 changed files with 33 additions and 55 deletions
Download patch file Download diff file Expand all files Collapse all files

63
include/openmc/interpolate.h
View file

@ -34,41 +34,11 @@ inline double interpolate_log_log(
return y0 * exp(f * log(y1 / y0));
}
inline double interpolate(const std::vector<double>& xs,
const std::vector<double>& ys, int idx, double x,
Interpolation i = Interpolation::lin_lin)
{
switch (i) {
case Interpolation::lin_lin:
return interpolate_lin_lin(xs[idx], xs[idx + 1], ys[idx], ys[idx + 1], x);
case Interpolation::log_log:
return interpolate_log_log(xs[idx], xs[idx + 1], ys[idx], ys[idx + 1], x);
case Interpolation::lin_log:
return interpolate_lin_log(xs[idx], xs[idx + 1], ys[idx], ys[idx + 1], x);
case Interpolation::log_lin:
return interpolate_log_lin(xs[idx], xs[idx + 1], ys[idx], ys[idx + 1], x);
default:
fatal_error("Unsupported interpolation");
}
}
inline double interpolate(const std::vector<double>& xs,
const std::vector<double>& ys, double x,
Interpolation i = Interpolation::lin_lin)
{
int idx = lower_bound_index(xs.begin(), xs.end(), x);
return interpolate(xs, ys, idx, x, i);
}
inline double interpolate_lagrangian(
const std::vector<double>& xs, const std::vector<double>& ys, double x)
const std::vector<double>& xs, const std::vector<double>& ys, int idx, double x, int order)
{
int idx = lower_bound_index(xs.begin(), xs.end(), x);
std::vector<double> coeffs;
int order = 3;
for (int i = 0; i < order + 1; i++) {
double numerator {1.0};
double denominator {1.0};
@ -85,6 +55,37 @@ inline double interpolate_lagrangian(
coeffs.begin(), coeffs.end(), ys.begin() + idx, 0.0);
}
inline double interpolate(const std::vector<double>& xs,
const std::vector<double>& ys, double x,
Interpolation i = Interpolation::lin_lin)
{
int idx = lower_bound_index(xs.begin(), xs.end(), x);
switch (i) {
case Interpolation::lin_lin:
return interpolate_lin_lin(xs[idx], xs[idx + 1], ys[idx], ys[idx + 1], x);
case Interpolation::log_log:
return interpolate_log_log(xs[idx], xs[idx + 1], ys[idx], ys[idx + 1], x);
case Interpolation::lin_log:
return interpolate_lin_log(xs[idx], xs[idx + 1], ys[idx], ys[idx + 1], x);
case Interpolation::log_lin:
return interpolate_log_lin(xs[idx], xs[idx + 1], ys[idx], ys[idx + 1], x);
case Interpolation::quadratic:
// move back one point if x is in the last interval of the x-grid
if (idx == xs.size() - 2 && idx > 0) idx--;
return interpolate_lagrangian(xs, ys, x, idx, 2);
case Interpolation::cubic:
// if x is not in the first interval of the x-grid, move back one
if (idx > 0) idx--;
// if the index was the last interval of the x-grid, move it back one more
if (idx == xs.size() - 3) idx--;
return interpolate_lagrangian(xs, ys, x, idx, 3);
default:
fatal_error("Unsupported interpolation");
}
}
} // namespace openmc
#endif

25
src/tallies/filter_energyfunc.cpp
View file

@ -75,31 +75,8 @@ void EnergyFunctionFilter::get_all_bins(
const Particle& p, TallyEstimator estimator, FilterMatch& match) const
{
if (p.E_last() >= energy_.front() && p.E_last() <= energy_.back()) {
// Search for the incoming energy bin.
auto i = lower_bound_index(energy_.begin(), energy_.end(), p.E_last());
double f, w;
switch (interpolation_) {
case Interpolation::lin_lin:
w = interpolate_lin_lin(
energy_[i], energy_[i + 1], y_[i], y_[i + 1], p.E_last());
break;
case Interpolation::lin_log:
w = interpolate_lin_log(
energy_[i], energy_[i + 1], y_[i], y_[i + 1], p.E_last());
break;
case Interpolation::log_lin:
w = interpolate_log_lin(
energy_[i], energy_[i + 1], y_[i], y_[i + 1], p.E_last());
break;
case Interpolation::log_log:
w = interpolate_log_log(
energy_[i], energy_[i + 1], y_[i], y_[i + 1], p.E_last());
break;
default:
fatal_error(fmt::format(
"Invalid interpolation scheme found on EnergyFunctionFilter {}", id()));
}
double w = interpolate(energy_, y_, p.E_last(), interpolation_);
// Interpolate on the lin-lin grid.
match.bins_.push_back(0);