From 4ba9acd459391a1d9806d112480ce47489b49fab Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 1 Aug 2022 18:50:49 -0500 Subject: [PATCH] Adding test results for additional interpolation types --- docs/source/io_formats/statepoint.rst | 7 +- include/openmc/interpolate.h | 65 +++++++++++++++++++ openmc/filter.py | 11 ++-- src/material.cpp | 4 +- src/physics.cpp | 13 ++-- src/tallies/filter_energyfunc.cpp | 24 +++++-- .../filter_energyfun/inputs_true.dat | 20 ++++++ .../filter_energyfun/results_true.dat | 4 ++ .../regression_tests/filter_energyfun/test.py | 35 +++++++--- 9 files changed, 148 insertions(+), 35 deletions(-) create mode 100644 include/openmc/interpolate.h diff --git a/docs/source/io_formats/statepoint.rst b/docs/source/io_formats/statepoint.rst index 71025d1e64..09a531ecd3 100644 --- a/docs/source/io_formats/statepoint.rst +++ b/docs/source/io_formats/statepoint.rst @@ -80,7 +80,7 @@ The current version of the statepoint file format is 17.0. dimension. - **Unstructured Mesh Only:** - **filename** (*char[]*) -- Name of the mesh file. - - **library** (*char[]*) -- Mesh library used to represent the + - **library** (*char[]*) -- Mesh library used to represent the mesh ("moab" or "libmesh"). - **length_multiplier** (*double*) Scaling factor applied to the mesh. - **volumes** (*double[]*) -- Volume of each mesh cell. @@ -108,9 +108,8 @@ The current version of the statepoint file format is 17.0. interpolation. Only used for 'energyfunction' filters. - **y** (*double[]*) -- Interpolant values for energyfunction interpolation. Only used for 'energyfunction' filters. - - **interpolation** (*int*) -- Interpolation type. Either - 1 (linear-linear) or 2 (log-log). Only used for 'energyfunction' - filters. + - **interpolation** (*int*) -- Interpolation type. Only used for + 'energyfunction' filters. **/tallies/derivatives/derivative /** diff --git a/include/openmc/interpolate.h b/include/openmc/interpolate.h new file mode 100644 index 0000000000..b88ab2f2e7 --- /dev/null +++ b/include/openmc/interpolate.h @@ -0,0 +1,65 @@ + +#ifndef OPENMC_INTERPOLATE_H +#define OPENMC_INTERPOLATE_H + +#include +#include + +#include "openmc/search.h" + +namespace openmc { + +inline double interpolate_lin_lin( + double x0, double x1, double y0, double y1, double x) +{ + return y0 + (x - x0) / (x1 - x0) * (y1 - y0); +} + +inline double interpolate_lin_log( + double x0, double x1, double y0, double y1, double x) +{ + return y0 + log(x / x0) / log(x1 / x0) * (y1 - y0); +} + +inline double interpolate_log_lin( + double x0, double x1, double y0, double y1, double x) +{ + return y0 * exp((x - x0) / (x1 - x0) * log(y1 / y0)); +} + +inline double interpolate_log_log( + double x0, double x1, double y0, double y1, double x) +{ + double f = log(x / x0) / log(x1 / x0); + return y0 * exp(f * log(y1 / y0)); +} + +inline double interpolate(const std::vector& xs, + const std::vector& 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& xs, + const std::vector& 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); +} + +} // namespace openmc + +#endif \ No newline at end of file diff --git a/openmc/filter.py b/openmc/filter.py index 989f5d762a..858238bcfc 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -12,6 +12,7 @@ import pandas as pd import openmc import openmc.checkvalue as cv from .cell import Cell +from .data.function import INTERPOLATION_SCHEME from .material import Material from .mixin import IDManagerMixin from .surface import Surface @@ -1877,7 +1878,6 @@ class EnergyFunctionFilter(Filter): A grid of interpolant values in [eV] interpolation : str The type of interpolation to be used. - One of ('linear-linear', 'log-log') filter_id : int Unique identifier for the filter @@ -1889,7 +1889,6 @@ class EnergyFunctionFilter(Filter): A grid of interpolant values in [eV] interpolation : str The type of interpolation to be used. - One of ('linear-linear', 'log-log') id : int Unique identifier for the filter num_bins : Integral @@ -1958,7 +1957,7 @@ class EnergyFunctionFilter(Filter): + group['type'][()].decode() + " instead") energy = group['energy'][()] - y = group['y'][()] + y = group['y'][()] filter_id = int(group.name.split('/')[-1].lstrip('filter ')) out = cls(energy, y, filter_id=filter_id) @@ -2045,7 +2044,7 @@ class EnergyFunctionFilter(Filter): @interpolation.setter def interpolation(self, val): cv.check_type('interpolation', val, str) - cv.check_value('interpolation', val, ('linear-linear', 'log-log')) + cv.check_value('interpolation', val, INTERPOLATION_SCHEME.values()) self._interpolation = val def to_xml_element(self): @@ -2076,11 +2075,11 @@ class EnergyFunctionFilter(Filter): def from_xml_element(cls, elem, **kwargs): filter_id = int(elem.get('id')) energy = [float(x) for x in get_text(elem, 'energy').split()] - y = [float(x) for x in get_text(elem, 'y').split()] + y = [float(x) for x in get_text(elem, 'y').split()] out = cls(energy, y, filter_id=filter_id) if elem.find('interpolation') is not None: out.interpolation = elem.find('interpolation') - return out + return out def can_merge(self, other): return False diff --git a/src/material.cpp b/src/material.cpp index 30dfa5ed50..8c75b079d3 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -734,8 +734,8 @@ void Material::init_bremsstrahlung() // Loop over photon energies double c = 0.0; for (int i = 0; i < j; ++i) { - // Integrate the CDF from the PDF using the trapezoidal rule in log-log - // space + // Integra te the CDF from the PDF using the trapezoidal rule in + // log-log space double w_l = std::log(data::ttb_e_grid(i)); double w_r = std::log(data::ttb_e_grid(i + 1)); double x_l = std::log(ttb->pdf(j, i)); diff --git a/src/physics.cpp b/src/physics.cpp index bcdcf7ecf9..5e9afa0903 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -7,6 +7,7 @@ #include "openmc/eigenvalue.h" #include "openmc/endf.h" #include "openmc/error.h" +#include "openmc/interpolate.h" #include "openmc/material.h" #include "openmc/math_functions.h" #include "openmc/message_passing.h" @@ -883,15 +884,9 @@ Direction sample_target_velocity(const Nuclide& nuc, double E, Direction u, if (sampling_method == ResScatMethod::dbrc) { // interpolate xs since we're not exactly at the energy indices - double xs_low = nuc.elastic_0K_[i_E_low]; - double m = (nuc.elastic_0K_[i_E_low + 1] - xs_low) / - (nuc.energy_0K_[i_E_low + 1] - nuc.energy_0K_[i_E_low]); - xs_low += m * (E_low - nuc.energy_0K_[i_E_low]); - double xs_up = nuc.elastic_0K_[i_E_up]; - m = (nuc.elastic_0K_[i_E_up + 1] - xs_up) / - (nuc.energy_0K_[i_E_up + 1] - nuc.energy_0K_[i_E_up]); - xs_up += m * (E_up - nuc.energy_0K_[i_E_up]); - + double xs_low = + interpolate(nuc.energy_0K_, nuc.elastic_0K_, i_E_low, E_low); + double xs_up = interpolate(nuc.energy_0K_, nuc.elastic_0K_, i_E_up, E_up); // get max 0K xs value over range of practical relative energies double xs_max = *std::max_element( &nuc.elastic_0K_[i_E_low + 1], &nuc.elastic_0K_[i_E_up + 1]); diff --git a/src/tallies/filter_energyfunc.cpp b/src/tallies/filter_energyfunc.cpp index b0a42d099b..f883975a36 100644 --- a/src/tallies/filter_energyfunc.cpp +++ b/src/tallies/filter_energyfunc.cpp @@ -3,6 +3,7 @@ #include #include "openmc/error.h" +#include "openmc/interpolate.h" #include "openmc/search.h" #include "openmc/settings.h" #include "openmc/xml_interface.h" @@ -31,12 +32,15 @@ void EnergyFunctionFilter::from_xml(pugi::xml_node node) std::string interpolation = get_node_value(node, "interpolation"); if (interpolation == "linear-linear") { interpolation_ = Interpolation::lin_lin; + } else if (interpolation == "linear-log") { + interpolation_ = Interpolation::lin_log; + } else if (interpolation == "log-linear") { + interpolation_ = Interpolation::log_lin; } else if (interpolation == "log-log") { interpolation_ = Interpolation::log_log; } else { fatal_error(fmt::format( - "Invalid interpolation type '{}' specified on EnergyFunctionFilter {}. " - "Must be 'linear-linear' or 'log-log'.", + "Found invalid interpolation type '{}' on EnergyFunctionFilter {}.", interpolation, id())); } } @@ -77,12 +81,20 @@ void EnergyFunctionFilter::get_all_bins( double f, w; switch (interpolation_) { case Interpolation::lin_lin: - f = (p.E_last() - energy_[i]) / (energy_[i + 1] - energy_[i]); - w = (1 - f) * y_[i] + f * y_[i + 1]; + 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: - f = log(p.E_last() / energy_[i]) / log(energy_[i + 1] / energy_[i]); - w = y_[i] * exp(f * log(y_[i + 1] / y_[i])); + w = interpolate_log_log( + energy_[i], energy_[i + 1], y_[i], y_[i + 1], p.E_last()); break; default: fatal_error(fmt::format( diff --git a/tests/regression_tests/filter_energyfun/inputs_true.dat b/tests/regression_tests/filter_energyfun/inputs_true.dat index 35d15d61c0..8c3db35186 100644 --- a/tests/regression_tests/filter_energyfun/inputs_true.dat +++ b/tests/regression_tests/filter_energyfun/inputs_true.dat @@ -29,6 +29,16 @@ 0.1 0.1 0.1333 0.158 0.18467 0.25618 0.4297 0.48 0.48 log-log + + 1e-05 0.369 1000.0 100000.0 600000.0 1000000.0 2000000.0 4000000.0 30000000.0 + 0.1 0.1 0.1333 0.158 0.18467 0.25618 0.4297 0.48 0.48 + linear-log + + + 1e-05 0.369 1000.0 100000.0 600000.0 1000000.0 2000000.0 4000000.0 30000000.0 + 0.1 0.1 0.1333 0.158 0.18467 0.25618 0.4297 0.48 0.48 + log-linear + Am241 (n,gamma) @@ -43,4 +53,14 @@ Am241 (n,gamma) + + 4 + Am241 + (n,gamma) + + + 5 + Am241 + (n,gamma) + diff --git a/tests/regression_tests/filter_energyfun/results_true.dat b/tests/regression_tests/filter_energyfun/results_true.dat index 4fedf2a0d7..c0b59680e8 100644 --- a/tests/regression_tests/filter_energyfun/results_true.dat +++ b/tests/regression_tests/filter_energyfun/results_true.dat @@ -2,3 +2,7 @@ 0 448ee8dfd19c4f Am241 ((n,gamma) / (n,gamma)) 1.74e-01 6.83e-03 energyfunction nuclide score mean std. dev. 0 37e006ae6b2e74 Am241 (n,gamma) 8.16e-02 2.24e-03 + energyfunction nuclide score mean std. dev. +0 b4e2ac84068d2d Am241 (n,gamma) 8.19e-02 2.25e-03 + energyfunction nuclide score mean std. dev. +0 dacf88242512ea Am241 (n,gamma) 7.95e-02 2.19e-03 diff --git a/tests/regression_tests/filter_energyfun/test.py b/tests/regression_tests/filter_energyfun/test.py index b09f0e5ec3..c297c094af 100644 --- a/tests/regression_tests/filter_energyfun/test.py +++ b/tests/regression_tests/filter_energyfun/test.py @@ -32,7 +32,7 @@ def model(): with pytest.raises(ValueError): filt1.interpolation = '5th order polynomial' - + # Also make a filter with the .from_tabulated1d constructor. Make sure # the filters are identical. tab1d = openmc.data.Tabulated1D(x, y) @@ -42,13 +42,22 @@ def model(): filt3 = openmc.EnergyFunctionFilter(x, y) filt3.interpolation = 'log-log' + filt4 = openmc.EnergyFunctionFilter(x , y) + filt4.interpolation = 'linear-log' + + filt5 = openmc.EnergyFunctionFilter(x, y) + filt5.interpolation = 'log-linear' + + filters = [filt1, filt3, filt4, filt5] # Make tallies - tallies = [openmc.Tally(), openmc.Tally(), openmc.Tally()] + tallies = [openmc.Tally() for i in range(5)] for t in tallies: t.scores = ['(n,gamma)'] t.nuclides = ['Am241'] - tallies[1].filters = [filt1] - tallies[2].filters = [filt3] + + for t, f in zip(tallies[1:], filters): + t.filters = [f] + model.tallies.extend(tallies) return model @@ -64,9 +73,11 @@ class FilterEnergyFunHarness(PyAPITestHarness): br_tally = sp.tallies[2] / sp.tallies[1] dataframes_string += br_tally.get_pandas_dataframe().to_string() + '\n' - # Write out the log-log interpolation as well - log_tally = sp.tallies[3] - dataframes_string += log_tally.get_pandas_dataframe().to_string() + '\n' + + for t_id in (3, 4, 5): + # Write out the log-log interpolation as well + ef_tally = sp.tallies[t_id] + dataframes_string += ef_tally.get_pandas_dataframe().to_string() + '\n' # Output the tally in a Pandas DataFrame. return dataframes_string @@ -78,7 +89,7 @@ class FilterEnergyFunHarness(PyAPITestHarness): sp = openmc.StatePoint(self._sp_name) # statepoint file round-trip checks - + # linear-linear interpolation tally sp_lin_lin_tally = sp.get_tally(id=2) sp_lin_lin_filt = sp_lin_lin_tally.find_filter(openmc.EnergyFunctionFilter) @@ -105,6 +116,14 @@ class FilterEnergyFunHarness(PyAPITestHarness): # we expect the log-log tally to have a higher value assert all(sp_lin_lin_tally.mean < sp_log_log_tally.mean) + sp_lin_log_tally = self._model.tallies[3] + sp_lin_log_filt = sp_lin_log_tally.find_filter(openmc.EnergyFunctionFilter) + assert sp_lin_log_filt.interpolation == 'linear-log' + + sp_log_lin_tally = self._model.tallies[4] + sp_log_lin_filt = sp_log_lin_tally.find_filter(openmc.EnergyFunctionFilter) + assert sp_log_lin_filt.interpolation == 'log-linear' + def test_filter_energyfun(model): harness = FilterEnergyFunHarness('statepoint.5.h5', model)