From 026d172255e64e45fc7e6ca37d23557ac455f469 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 28 Jul 2022 10:04:23 -0500 Subject: [PATCH 01/37] Adding interpolation property to EnergyFunctionFilter. --- openmc/filter.py | 41 ++++++++++++++++++++++++++++++++++------- 1 file changed, 34 insertions(+), 7 deletions(-) diff --git a/openmc/filter.py b/openmc/filter.py index 7874a775b..0d8458be5 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -1884,6 +1884,9 @@ class EnergyFunctionFilter(Filter): A grid of energy values in [eV] y : iterable of Real A grid of interpolant values in [eV] + interpolation : str + Used to indicate the type of interpolation used. + One of ('linear-linear', 'log-log') id : int Unique identifier for the filter num_bins : Integral @@ -1895,6 +1898,7 @@ class EnergyFunctionFilter(Filter): self.energy = energy self.y = y self.id = filter_id + self.interpolation = 'linear-linear' def __eq__(self, other): if type(self) is not type(other): @@ -1949,10 +1953,14 @@ 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 ')) - return cls(energy, y, filter_id=filter_id) + out = cls(energy, y, filter_id=filter_id) + if 'interpolation' in group: + out.interpolation = group['interpolation'][()] + + return out @classmethod def from_tabulated1d(cls, tab1d): @@ -1974,10 +1982,13 @@ class EnergyFunctionFilter(Filter): if tab1d.n_regions > 1: raise ValueError('Only Tabulated1Ds with a single interpolation ' 'region are supported') - if tab1d.interpolation[0] != 2: - raise ValueError('Only linear-linear Tabulated1Ds are supported') + if tab1d.interpolation[0] not in (2, 5): + raise ValueError('Only linear-linear or log-log Tabulated1Ds are supported') + out = cls(tab1d.x, tab1d.y) + if tab1d.interpolation[0] == 5: + out.interpolation = 'log-log' - return cls(tab1d.x, tab1d.y) + return out @property def energy(self): @@ -1987,6 +1998,10 @@ class EnergyFunctionFilter(Filter): def y(self): return self._y + @property + def interpolation(self): + return self._interpolation + @property def bins(self): raise AttributeError('EnergyFunctionFilters have no bins.') @@ -2021,6 +2036,12 @@ class EnergyFunctionFilter(Filter): def bins(self, bins): raise RuntimeError('EnergyFunctionFilters have no bins.') + @interpolation.setter + def interpolation(self, val): + cv.check_type('interpolation', val, str) + cv.check_value('interpolation', val, ('linear-linear', 'log-log')) + self._interpolation = val + def to_xml_element(self): """Return XML Element representing the Filter. @@ -2040,14 +2061,20 @@ class EnergyFunctionFilter(Filter): subelement = ET.SubElement(element, 'y') subelement.text = ' '.join(str(y) for y in self.y) + subelement = ET.SubElement(element, 'interpolation') + subelement.text = self.interpolation + return element @classmethod 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()] - return cls(energy, y, filter_id=filter_id) + 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 def can_merge(self, other): return False From 8f26df7f5a91545903331916981c9d51a45880e4 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 28 Jul 2022 10:04:57 -0500 Subject: [PATCH 02/37] Adding interpolation value check and updating test inputs --- tests/regression_tests/filter_energyfun/inputs_true.dat | 1 + tests/regression_tests/filter_energyfun/test.py | 2 ++ 2 files changed, 3 insertions(+) diff --git a/tests/regression_tests/filter_energyfun/inputs_true.dat b/tests/regression_tests/filter_energyfun/inputs_true.dat index 0f506f3f5..15de0f027 100644 --- a/tests/regression_tests/filter_energyfun/inputs_true.dat +++ b/tests/regression_tests/filter_energyfun/inputs_true.dat @@ -22,6 +22,7 @@ 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-linear Am241 diff --git a/tests/regression_tests/filter_energyfun/test.py b/tests/regression_tests/filter_energyfun/test.py index f61f05d78..c76b525d7 100644 --- a/tests/regression_tests/filter_energyfun/test.py +++ b/tests/regression_tests/filter_energyfun/test.py @@ -28,6 +28,8 @@ def model(): # Make an EnergyFunctionFilter directly from the x and y lists. filt1 = openmc.EnergyFunctionFilter(x, y) + assert filt1.interpolation == 'linear-linear' + # Also make a filter with the .from_tabulated1d constructor. Make sure # the filters are identical. tab1d = openmc.data.Tabulated1D(x, y) From 0d551cbb21f76894ea70f9e712f81d09f4b3c7f2 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 28 Jul 2022 11:14:30 -0500 Subject: [PATCH 03/37] Adding support for log-log interpolation on the C++ side --- include/openmc/tallies/filter_energyfunc.h | 2 ++ openmc/filter.py | 5 +++-- src/tallies/filter_energyfunc.cpp | 26 ++++++++++++++++++---- 3 files changed, 27 insertions(+), 6 deletions(-) diff --git a/include/openmc/tallies/filter_energyfunc.h b/include/openmc/tallies/filter_energyfunc.h index c066dfa17..f650976ef 100644 --- a/include/openmc/tallies/filter_energyfunc.h +++ b/include/openmc/tallies/filter_energyfunc.h @@ -1,6 +1,7 @@ #ifndef OPENMC_TALLIES_FILTER_ENERGYFUNC_H #define OPENMC_TALLIES_FILTER_ENERGYFUNC_H +#include "openmc/constants.h" #include "openmc/tallies/filter.h" #include "openmc/vector.h" @@ -39,6 +40,7 @@ public: const vector& energy() const { return energy_; } const vector& y() const { return y_; } + Interpolation interpolation_; void set_data(gsl::span energy, gsl::span y); private: diff --git a/openmc/filter.py b/openmc/filter.py index 0d8458be5..f44ac99f8 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -1957,8 +1957,9 @@ class EnergyFunctionFilter(Filter): filter_id = int(group.name.split('/')[-1].lstrip('filter ')) out = cls(energy, y, filter_id=filter_id) - if 'interpolation' in group: - out.interpolation = group['interpolation'][()] + if 'interpolation' in group.attrs: + out.interpolation = \ + openmc.data.INTERPOLATION_SCHEME[group.attrs['interpolation']] return out diff --git a/src/tallies/filter_energyfunc.cpp b/src/tallies/filter_energyfunc.cpp index fd595e332..2fcca7110 100644 --- a/src/tallies/filter_energyfunc.cpp +++ b/src/tallies/filter_energyfunc.cpp @@ -22,9 +22,15 @@ void EnergyFunctionFilter::from_xml(pugi::xml_node node) if (!check_for_node(node, "y")) fatal_error("y values not specified for EnergyFunction filter."); - + auto y = get_node_array(node, "y"); + // use linear-linear interpolation by default + interpolation_ = Interpolation::lin_lin; + if (check_for_node(node, "interpolation")) { + std::string interpolation = get_node_value(node, "interpolation"); + } + this->set_data(energy, y); } @@ -58,12 +64,23 @@ void EnergyFunctionFilter::get_all_bins( // Search for the incoming energy bin. auto i = lower_bound_index(energy_.begin(), energy_.end(), p.E_last()); - // Compute the interpolation factor between the nearest bins. - double f = (p.E_last() - energy_[i]) / (energy_[i + 1] - energy_[i]); + 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]; + 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])); + break; + default: + fatal_error(fmt::format("Invalid interpolation scheme found on EnergyFunctionFilter {}", id())); + } // Interpolate on the lin-lin grid. match.bins_.push_back(0); - match.weights_.push_back((1 - f) * y_[i] + f * y_[i + 1]); + match.weights_.push_back(w); } } @@ -72,6 +89,7 @@ void EnergyFunctionFilter::to_statepoint(hid_t filter_group) const Filter::to_statepoint(filter_group); write_dataset(filter_group, "energy", energy_); write_dataset(filter_group, "y", y_); + write_attr_int(filter_group, "interpolation", interpolation_); } std::string EnergyFunctionFilter::text_label(int bin) const From 0c1b1875b7df89f4bb9ae8eee53d116b11816d3b Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 28 Jul 2022 11:31:12 -0500 Subject: [PATCH 04/37] Adding a quick test for the interpolation value checking --- src/tallies/filter_energyfunc.cpp | 23 ++++++++++--------- .../regression_tests/filter_energyfun/test.py | 3 +++ 2 files changed, 15 insertions(+), 11 deletions(-) diff --git a/src/tallies/filter_energyfunc.cpp b/src/tallies/filter_energyfunc.cpp index 2fcca7110..d783366bc 100644 --- a/src/tallies/filter_energyfunc.cpp +++ b/src/tallies/filter_energyfunc.cpp @@ -22,7 +22,7 @@ void EnergyFunctionFilter::from_xml(pugi::xml_node node) if (!check_for_node(node, "y")) fatal_error("y values not specified for EnergyFunction filter."); - + auto y = get_node_array(node, "y"); // use linear-linear interpolation by default @@ -66,16 +66,17 @@ 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]; - 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])); - break; - default: - fatal_error(fmt::format("Invalid interpolation scheme found on EnergyFunctionFilter {}", id())); + case Interpolation::lin_lin: + f = (p.E_last() - energy_[i]) / (energy_[i + 1] - energy_[i]); + w = (1 - f) * y_[i] + f * y_[i + 1]; + 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])); + break; + default: + fatal_error(fmt::format( + "Invalid interpolation scheme found on EnergyFunctionFilter {}", id())); } // Interpolate on the lin-lin grid. diff --git a/tests/regression_tests/filter_energyfun/test.py b/tests/regression_tests/filter_energyfun/test.py index c76b525d7..ba5a9b16f 100644 --- a/tests/regression_tests/filter_energyfun/test.py +++ b/tests/regression_tests/filter_energyfun/test.py @@ -30,6 +30,9 @@ def model(): assert filt1.interpolation == 'linear-linear' + 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) From a0de07bf707957f615cd163ce1ea96de433e052b Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 28 Jul 2022 18:23:51 -0500 Subject: [PATCH 05/37] Correcting hdf5 function signature --- src/tallies/filter_energyfunc.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tallies/filter_energyfunc.cpp b/src/tallies/filter_energyfunc.cpp index d783366bc..37c1c3178 100644 --- a/src/tallies/filter_energyfunc.cpp +++ b/src/tallies/filter_energyfunc.cpp @@ -90,7 +90,7 @@ void EnergyFunctionFilter::to_statepoint(hid_t filter_group) const Filter::to_statepoint(filter_group); write_dataset(filter_group, "energy", energy_); write_dataset(filter_group, "y", y_); - write_attr_int(filter_group, "interpolation", interpolation_); + write_attribute(filter_group, "interpolation", static_cast(interpolation_)); } std::string EnergyFunctionFilter::text_label(int bin) const From b4cf42b092eb456a7ad7b6480939a4c0bcb21b94 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Fri, 29 Jul 2022 23:37:34 -0500 Subject: [PATCH 06/37] Adding entry to statepoint documentation. --- docs/source/io_formats/statepoint.rst | 3 +++ src/tallies/filter_energyfunc.cpp | 3 ++- 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/docs/source/io_formats/statepoint.rst b/docs/source/io_formats/statepoint.rst index c1598f2c3..2332ebf21 100644 --- a/docs/source/io_formats/statepoint.rst +++ b/docs/source/io_formats/statepoint.rst @@ -108,6 +108,9 @@ 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. One of + (1 - linear-linear or 2 - log-log). Only used for 'energyfunction' + filters. **/tallies/derivatives/derivative /** diff --git a/src/tallies/filter_energyfunc.cpp b/src/tallies/filter_energyfunc.cpp index 37c1c3178..819b5d13e 100644 --- a/src/tallies/filter_energyfunc.cpp +++ b/src/tallies/filter_energyfunc.cpp @@ -90,7 +90,8 @@ void EnergyFunctionFilter::to_statepoint(hid_t filter_group) const Filter::to_statepoint(filter_group); write_dataset(filter_group, "energy", energy_); write_dataset(filter_group, "y", y_); - write_attribute(filter_group, "interpolation", static_cast(interpolation_)); + write_attribute( + filter_group, "interpolation", static_cast(interpolation_)); } std::string EnergyFunctionFilter::text_label(int bin) const From d88648634e0d9f84ab2642fc77ee20053ef33f65 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Fri, 29 Jul 2022 23:49:58 -0500 Subject: [PATCH 07/37] Adding check for round-trip of spec for eff --- .../regression_tests/filter_energyfun/test.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/tests/regression_tests/filter_energyfun/test.py b/tests/regression_tests/filter_energyfun/test.py index ba5a9b16f..b5ab1387c 100644 --- a/tests/regression_tests/filter_energyfun/test.py +++ b/tests/regression_tests/filter_energyfun/test.py @@ -61,6 +61,24 @@ class FilterEnergyFunHarness(PyAPITestHarness): # Output the tally in a Pandas DataFrame. return br_tally.get_pandas_dataframe().to_string() + '\n' + def _compare_results(self): + super()._compare_results() + + # Read the statepoint file. + sp = openmc.StatePoint(self._sp_name) + + # check that values are round-tripped correctly from + # the statepoint file + sp_tally = sp.get_tally(id=2) + sp_filt = sp_tally.find_filter(openmc.EnergyFunctionFilter) + + model_tally = self._model.tallies[1] + model_filt = model_tally.find_filter(openmc.EnergyFunctionFilter) + + assert sp_filt.interpolation == model_filt.interpolation + assert all(sp_filt.energy == model_filt.energy) + assert all(sp_filt.y == model_filt.y) + def test_filter_energyfun(model): harness = FilterEnergyFunHarness('statepoint.5.h5', model) From 2ef49a9c9b66736cc69eb4a7a7ff31d7cd9b4407 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sat, 30 Jul 2022 00:18:10 -0500 Subject: [PATCH 08/37] Updating tests for eff --- openmc/filter.py | 6 ++++-- .../regression_tests/filter_energyfun/inputs_true.dat | 10 ++++++++++ .../filter_energyfun/results_true.dat | 2 +- tests/regression_tests/filter_energyfun/test.py | 11 +++++++++-- 4 files changed, 24 insertions(+), 5 deletions(-) diff --git a/openmc/filter.py b/openmc/filter.py index f44ac99f8..f5268ad81 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -1894,11 +1894,11 @@ class EnergyFunctionFilter(Filter): """ - def __init__(self, energy, y, filter_id=None): + def __init__(self, energy, y, interpolation='linear-linear', filter_id=None): self.energy = energy self.y = y self.id = filter_id - self.interpolation = 'linear-linear' + self.interpolation = interpolation def __eq__(self, other): if type(self) is not type(other): @@ -1936,12 +1936,14 @@ class EnergyFunctionFilter(Filter): string = type(self).__name__ + '\n' string += '{: <16}=\t{}\n'.format('\tEnergy', self.energy) string += '{: <16}=\t{}\n'.format('\tInterpolant', self.y) + string += '{: <16}=\t{}\n'.format('\tInterpolation', self.interpolation) return hash(string) def __repr__(self): string = type(self).__name__ + '\n' string += '{: <16}=\t{}\n'.format('\tEnergy', self.energy) string += '{: <16}=\t{}\n'.format('\tInterpolant', self.y) + string += '{: <16}=\t{}\n'.format('\tInterpolation', self.interpolation) string += '{: <16}=\t{}\n'.format('\tID', self.id) return string diff --git a/tests/regression_tests/filter_energyfun/inputs_true.dat b/tests/regression_tests/filter_energyfun/inputs_true.dat index 15de0f027..35d15d61c 100644 --- a/tests/regression_tests/filter_energyfun/inputs_true.dat +++ b/tests/regression_tests/filter_energyfun/inputs_true.dat @@ -24,6 +24,11 @@ 0.1 0.1 0.1333 0.158 0.18467 0.25618 0.4297 0.48 0.48 linear-linear + + 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-log + Am241 (n,gamma) @@ -33,4 +38,9 @@ Am241 (n,gamma) + + 3 + Am241 + (n,gamma) + diff --git a/tests/regression_tests/filter_energyfun/results_true.dat b/tests/regression_tests/filter_energyfun/results_true.dat index a1fda93a8..a75b3fd4f 100644 --- a/tests/regression_tests/filter_energyfun/results_true.dat +++ b/tests/regression_tests/filter_energyfun/results_true.dat @@ -1,2 +1,2 @@ energyfunction nuclide score mean std. dev. -0 d2effa26cb3cf2 Am241 ((n,gamma) / (n,gamma)) 1.74e-01 6.83e-03 +0 448ee8dfd19c4f Am241 ((n,gamma) / (n,gamma)) 1.74e-01 6.83e-03 diff --git a/tests/regression_tests/filter_energyfun/test.py b/tests/regression_tests/filter_energyfun/test.py index b5ab1387c..cb658b902 100644 --- a/tests/regression_tests/filter_energyfun/test.py +++ b/tests/regression_tests/filter_energyfun/test.py @@ -39,14 +39,21 @@ def model(): filt2 = openmc.EnergyFunctionFilter.from_tabulated1d(tab1d) assert filt1 == filt2, 'Error with the .from_tabulated1d constructor' + filt3 = openmc.EnergyFunctionFilter(x, y) + filt3.interpolation = 'log-log' + print(filt3) + # Make tallies - tallies = [openmc.Tally(), openmc.Tally()] + tallies = [openmc.Tally(), openmc.Tally(), openmc.Tally()] for t in tallies: t.scores = ['(n,gamma)'] t.nuclides = ['Am241'] tallies[1].filters = [filt1] + tallies[2].filters = [filt3] model.tallies.extend(tallies) + + return model @@ -75,7 +82,7 @@ class FilterEnergyFunHarness(PyAPITestHarness): model_tally = self._model.tallies[1] model_filt = model_tally.find_filter(openmc.EnergyFunctionFilter) - assert sp_filt.interpolation == model_filt.interpolation + assert sp_filt.interpolation == 'linear-linear' assert all(sp_filt.energy == model_filt.energy) assert all(sp_filt.y == model_filt.y) From c09f539cd5bd87ab3201c495a9c96552a0e5fe68 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sat, 30 Jul 2022 00:34:37 -0500 Subject: [PATCH 09/37] Updating check to log-log --- openmc/filter.py | 4 ++-- src/tallies/filter_energyfunc.cpp | 12 +++++++++++- tests/regression_tests/filter_energyfun/test.py | 10 ++++------ 3 files changed, 17 insertions(+), 9 deletions(-) diff --git a/openmc/filter.py b/openmc/filter.py index f5268ad81..432ea1599 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -1959,9 +1959,9 @@ class EnergyFunctionFilter(Filter): filter_id = int(group.name.split('/')[-1].lstrip('filter ')) out = cls(energy, y, filter_id=filter_id) - if 'interpolation' in group.attrs: + if 'interpolation' in group: out.interpolation = \ - openmc.data.INTERPOLATION_SCHEME[group.attrs['interpolation']] + openmc.data.INTERPOLATION_SCHEME[group['interpolation'][()]] return out diff --git a/src/tallies/filter_energyfunc.cpp b/src/tallies/filter_energyfunc.cpp index 819b5d13e..fb76fb215 100644 --- a/src/tallies/filter_energyfunc.cpp +++ b/src/tallies/filter_energyfunc.cpp @@ -29,6 +29,16 @@ void EnergyFunctionFilter::from_xml(pugi::xml_node node) interpolation_ = Interpolation::lin_lin; if (check_for_node(node, "interpolation")) { std::string interpolation = get_node_value(node, "interpolation"); + if (interpolation == "linear-linear") { + interpolation_ = Interpolation::lin_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'.", + interpolation, id())); + } } this->set_data(energy, y); @@ -90,7 +100,7 @@ void EnergyFunctionFilter::to_statepoint(hid_t filter_group) const Filter::to_statepoint(filter_group); write_dataset(filter_group, "energy", energy_); write_dataset(filter_group, "y", y_); - write_attribute( + write_dataset( filter_group, "interpolation", static_cast(interpolation_)); } diff --git a/tests/regression_tests/filter_energyfun/test.py b/tests/regression_tests/filter_energyfun/test.py index cb658b902..e0a191ce8 100644 --- a/tests/regression_tests/filter_energyfun/test.py +++ b/tests/regression_tests/filter_energyfun/test.py @@ -52,8 +52,6 @@ def model(): tallies[2].filters = [filt3] model.tallies.extend(tallies) - - return model @@ -74,15 +72,15 @@ class FilterEnergyFunHarness(PyAPITestHarness): # Read the statepoint file. sp = openmc.StatePoint(self._sp_name) - # check that values are round-tripped correctly from + # check that information is round-tripped correctly from # the statepoint file - sp_tally = sp.get_tally(id=2) + sp_tally = sp.get_tally(id=3) sp_filt = sp_tally.find_filter(openmc.EnergyFunctionFilter) - model_tally = self._model.tallies[1] + model_tally = self._model.tallies[2] model_filt = model_tally.find_filter(openmc.EnergyFunctionFilter) - assert sp_filt.interpolation == 'linear-linear' + assert sp_filt.interpolation == 'log-log' assert all(sp_filt.energy == model_filt.energy) assert all(sp_filt.y == model_filt.y) From aeca1173e7f9b4a51c774ecfbd52f13050f41983 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sat, 30 Jul 2022 00:43:20 -0500 Subject: [PATCH 10/37] Expanding on eff checks for log-log case. --- .../regression_tests/filter_energyfun/test.py | 35 ++++++++++++++----- 1 file changed, 26 insertions(+), 9 deletions(-) diff --git a/tests/regression_tests/filter_energyfun/test.py b/tests/regression_tests/filter_energyfun/test.py index e0a191ce8..9e8189f21 100644 --- a/tests/regression_tests/filter_energyfun/test.py +++ b/tests/regression_tests/filter_energyfun/test.py @@ -72,17 +72,34 @@ class FilterEnergyFunHarness(PyAPITestHarness): # Read the statepoint file. sp = openmc.StatePoint(self._sp_name) - # check that information is round-tripped correctly from - # the statepoint file - sp_tally = sp.get_tally(id=3) - sp_filt = sp_tally.find_filter(openmc.EnergyFunctionFilter) + # 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) - model_tally = self._model.tallies[2] - model_filt = model_tally.find_filter(openmc.EnergyFunctionFilter) + model_lin_lin_tally = self._model.tallies[1] + model_lin_lin_filt = model_lin_lin_tally.find_filter(openmc.EnergyFunctionFilter) - assert sp_filt.interpolation == 'log-log' - assert all(sp_filt.energy == model_filt.energy) - assert all(sp_filt.y == model_filt.y) + assert sp_lin_lin_filt.interpolation == 'linear-linear' + assert all(sp_lin_lin_filt.energy == model_lin_lin_filt.energy) + assert all(sp_lin_lin_filt.y == model_lin_lin_filt.y) + + # log-log interpolation tally + sp_log_log_tally = sp.get_tally(id=3) + sp_log_log_filt = sp_log_log_tally.find_filter(openmc.EnergyFunctionFilter) + + model_log_log_tally = self._model.tallies[2] + model_log_log_filt = model_log_log_tally.find_filter(openmc.EnergyFunctionFilter) + + assert sp_log_log_filt.interpolation == 'log-log' + assert all(sp_log_log_filt.energy == model_log_log_filt.energy) + assert all(sp_log_log_filt.y == model_log_log_filt.y) + + + # because the values of y are monotonically increasing, + # we expect the log-log tally to have a higher value + assert all(sp_lin_lin_tally.mean < sp_log_log_tally.mean) def test_filter_energyfun(model): From 7ea116c89fce2a63c853b82021483f14cdb228e4 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sat, 30 Jul 2022 00:51:52 -0500 Subject: [PATCH 11/37] Adding log-lg tally output and updating results. --- tests/regression_tests/filter_energyfun/results_true.dat | 2 ++ tests/regression_tests/filter_energyfun/test.py | 9 +++++++-- 2 files changed, 9 insertions(+), 2 deletions(-) diff --git a/tests/regression_tests/filter_energyfun/results_true.dat b/tests/regression_tests/filter_energyfun/results_true.dat index a75b3fd4f..4fedf2a0d 100644 --- a/tests/regression_tests/filter_energyfun/results_true.dat +++ b/tests/regression_tests/filter_energyfun/results_true.dat @@ -1,2 +1,4 @@ energyfunction nuclide score mean std. dev. 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 diff --git a/tests/regression_tests/filter_energyfun/test.py b/tests/regression_tests/filter_energyfun/test.py index 9e8189f21..93eed3327 100644 --- a/tests/regression_tests/filter_energyfun/test.py +++ b/tests/regression_tests/filter_energyfun/test.py @@ -60,11 +60,17 @@ class FilterEnergyFunHarness(PyAPITestHarness): # Read the statepoint file. sp = openmc.StatePoint(self._sp_name) + dataframes_string = "" # Use tally arithmetic to compute the branching ratio. 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' # Output the tally in a Pandas DataFrame. - return br_tally.get_pandas_dataframe().to_string() + '\n' + return dataframes_string def _compare_results(self): super()._compare_results() @@ -96,7 +102,6 @@ class FilterEnergyFunHarness(PyAPITestHarness): assert all(sp_log_log_filt.energy == model_log_log_filt.energy) assert all(sp_log_log_filt.y == model_log_log_filt.y) - # because the values of y are monotonically increasing, # we expect the log-log tally to have a higher value assert all(sp_lin_lin_tally.mean < sp_log_log_tally.mean) From 66ff5b873de7095ef9d9704aca0f46dd8ad7d823 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sun, 31 Jul 2022 00:51:34 -0500 Subject: [PATCH 12/37] Adding missing docstring --- openmc/filter.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/openmc/filter.py b/openmc/filter.py index 432ea1599..989f5d762 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -1875,6 +1875,9 @@ class EnergyFunctionFilter(Filter): A grid of energy values in [eV] y : iterable of Real 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 @@ -1885,7 +1888,7 @@ class EnergyFunctionFilter(Filter): y : iterable of Real A grid of interpolant values in [eV] interpolation : str - Used to indicate the type of interpolation used. + The type of interpolation to be used. One of ('linear-linear', 'log-log') id : int Unique identifier for the filter From bdf73abb06fc211ee10d5fc49946f7502f47ddb0 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sun, 31 Jul 2022 01:05:42 -0500 Subject: [PATCH 13/37] Some cleanup --- docs/source/io_formats/statepoint.rst | 4 ++-- src/tallies/filter_energyfunc.cpp | 2 +- tests/regression_tests/filter_energyfun/test.py | 1 - 3 files changed, 3 insertions(+), 4 deletions(-) diff --git a/docs/source/io_formats/statepoint.rst b/docs/source/io_formats/statepoint.rst index 2332ebf21..71025d1e6 100644 --- a/docs/source/io_formats/statepoint.rst +++ b/docs/source/io_formats/statepoint.rst @@ -108,8 +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. One of - (1 - linear-linear or 2 - log-log). Only used for 'energyfunction' + - **interpolation** (*int*) -- Interpolation type. Either + 1 (linear-linear) or 2 (log-log). Only used for 'energyfunction' filters. **/tallies/derivatives/derivative /** diff --git a/src/tallies/filter_energyfunc.cpp b/src/tallies/filter_energyfunc.cpp index fb76fb215..b0a42d099 100644 --- a/src/tallies/filter_energyfunc.cpp +++ b/src/tallies/filter_energyfunc.cpp @@ -25,7 +25,7 @@ void EnergyFunctionFilter::from_xml(pugi::xml_node node) auto y = get_node_array(node, "y"); - // use linear-linear interpolation by default + // default to linear-linear interpolation interpolation_ = Interpolation::lin_lin; if (check_for_node(node, "interpolation")) { std::string interpolation = get_node_value(node, "interpolation"); diff --git a/tests/regression_tests/filter_energyfun/test.py b/tests/regression_tests/filter_energyfun/test.py index 93eed3327..b09f0e5ec 100644 --- a/tests/regression_tests/filter_energyfun/test.py +++ b/tests/regression_tests/filter_energyfun/test.py @@ -41,7 +41,6 @@ def model(): filt3 = openmc.EnergyFunctionFilter(x, y) filt3.interpolation = 'log-log' - print(filt3) # Make tallies tallies = [openmc.Tally(), openmc.Tally(), openmc.Tally()] From d78a9b188594a262bb84cb758fd8f13f75887966 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 1 Aug 2022 17:34:39 -0500 Subject: [PATCH 14/37] Adding new interpolation entry. --- include/openmc/constants.h | 3 ++- openmc/data/function.py | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/include/openmc/constants.h b/include/openmc/constants.h index a7362ea9e..9c9a5f9b3 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -318,7 +318,8 @@ enum class Interpolation { lin_lin = 2, lin_log = 3, log_lin = 4, - log_log = 5 + log_log = 5, + cubic = 6 }; enum class RunMode { diff --git a/openmc/data/function.py b/openmc/data/function.py index b0390d19c..f03319eed 100644 --- a/openmc/data/function.py +++ b/openmc/data/function.py @@ -13,7 +13,7 @@ from openmc.mixin import EqualityMixin from .data import EV_PER_MEV INTERPOLATION_SCHEME = {1: 'histogram', 2: 'linear-linear', 3: 'linear-log', - 4: 'log-linear', 5: 'log-log'} + 4: 'log-linear', 5: 'log-log', 6 : 'cubic'} def sum_functions(funcs): From 4ba9acd459391a1d9806d112480ce47489b49fab Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 1 Aug 2022 18:50:49 -0500 Subject: [PATCH 15/37] 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 71025d1e6..09a531ecd 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 000000000..b88ab2f2e --- /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 989f5d762..858238bcf 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 30dfa5ed5..8c75b079d 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 bcdcf7ecf..5e9afa090 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 b0a42d099..f883975a3 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 35d15d61c..8c3db3518 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 4fedf2a0d..c0b59680e 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 b09f0e5ec..c297c094a 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) From 48a9ac9a2cb9df43d391f3a0dbe6948f26536afb Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sat, 6 Aug 2022 06:48:47 -0500 Subject: [PATCH 16/37] Adding initial lagrangian interpolation. --- include/openmc/interpolate.h | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/include/openmc/interpolate.h b/include/openmc/interpolate.h index b88ab2f2e..068314033 100644 --- a/include/openmc/interpolate.h +++ b/include/openmc/interpolate.h @@ -60,6 +60,31 @@ inline double interpolate(const std::vector& xs, return interpolate(xs, ys, idx, x, i); } +inline double interpolate_lagrangian( + const std::vector& xs, const std::vector& ys, double x) +{ + int idx = lower_bound_index(xs.begin(), xs.end(), x); + + std::vector coeffs; + + int order = 3; + + for (int i = 0; i < order + 1; i++) { + double numerator {1.0}; + double denominator {1.0}; + for (int j = 0; j < order; j++) { + if (i == j) + continue; + numerator *= (x - xs[idx + j]); + denominator *= (xs[idx + i] - xs[idx + j]); + } + coeffs.push_back(numerator / denominator); + } + + return std::inner_product( + coeffs.begin(), coeffs.end(), ys.begin() + idx, 0.0); +} + } // namespace openmc #endif \ No newline at end of file From 445fb0af3627b728a35de8fffa33f17501b85ead Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 8 Aug 2022 07:38:37 -0500 Subject: [PATCH 17/37] Adjusting interpolation enumeration values. --- include/openmc/constants.h | 3 ++- openmc/data/function.py | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/include/openmc/constants.h b/include/openmc/constants.h index 9c9a5f9b3..32cd8186a 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -319,7 +319,8 @@ enum class Interpolation { lin_log = 3, log_lin = 4, log_log = 5, - cubic = 6 + quadratic = 6, + cubic = 7 }; enum class RunMode { diff --git a/openmc/data/function.py b/openmc/data/function.py index f03319eed..59f549d6f 100644 --- a/openmc/data/function.py +++ b/openmc/data/function.py @@ -13,7 +13,7 @@ from openmc.mixin import EqualityMixin from .data import EV_PER_MEV INTERPOLATION_SCHEME = {1: 'histogram', 2: 'linear-linear', 3: 'linear-log', - 4: 'log-linear', 5: 'log-log', 6 : 'cubic'} + 4: 'log-linear', 5: 'log-log', 6 : 'quadratic', 7 : 'cubic'} def sum_functions(funcs): From d300bac7340f473dc8d4cf6b586ea92a69ecac55 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 6 Sep 2022 20:30:43 -0500 Subject: [PATCH 18/37] Reverting some changes outside of the EFF. --- src/material.cpp | 4 ++-- src/physics.cpp | 13 +++++++++---- 2 files changed, 11 insertions(+), 6 deletions(-) diff --git a/src/material.cpp b/src/material.cpp index 8c75b079d..30dfa5ed5 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) { - // Integra te the CDF from the PDF using the trapezoidal rule in - // log-log space + // Integrate 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 5e9afa090..bcdcf7ecf 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -7,7 +7,6 @@ #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" @@ -884,9 +883,15 @@ 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 = - 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); + 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]); + // 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]); From 5168b3392a2a5b4e3aa54531dc1ba0a7bedf0487 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 6 Sep 2022 20:31:00 -0500 Subject: [PATCH 19/37] Accepting other interpolation types for converted Tabulated1D distributions --- openmc/filter.py | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/openmc/filter.py b/openmc/filter.py index 858238bcf..9815572a3 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -1987,10 +1987,18 @@ class EnergyFunctionFilter(Filter): if tab1d.n_regions > 1: raise ValueError('Only Tabulated1Ds with a single interpolation ' 'region are supported') - if tab1d.interpolation[0] not in (2, 5): - raise ValueError('Only linear-linear or log-log Tabulated1Ds are supported') + if tab1d.interpolation[0] not in (2, 3, 4, 5): + raise ValueError('Only linear-linear, linear-log, log-linear, and ' + 'log-log Tabulated1Ds are supported') out = cls(tab1d.x, tab1d.y) - if tab1d.interpolation[0] == 5: + + if tab1d.interpolation[0] == 2: + out.interpolation = 'linear-lienar' + elif tab1d.interpolation[0] == 3: + out.interpolation = 'linear-log' + elif tab1d.interpolation[0] == 4: + out.interpolation = 'log-linear' + elif tab1d.interpolation[0] == 5: out.interpolation = 'log-log' return out From 35aa433b9a460c6cf5093d3864fb50e8c56d140c Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 6 Sep 2022 22:49:00 -0500 Subject: [PATCH 20/37] Updating interpolation function signature. --- include/openmc/interpolate.h | 63 ++++++++++++++++--------------- src/tallies/filter_energyfunc.cpp | 25 +----------- 2 files changed, 33 insertions(+), 55 deletions(-) diff --git a/include/openmc/interpolate.h b/include/openmc/interpolate.h index 068314033..2764b7ffc 100644 --- a/include/openmc/interpolate.h +++ b/include/openmc/interpolate.h @@ -34,41 +34,11 @@ inline double interpolate_log_log( 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); -} - inline double interpolate_lagrangian( - const std::vector& xs, const std::vector& ys, double x) + const std::vector& xs, const std::vector& ys, int idx, double x, int order) { - int idx = lower_bound_index(xs.begin(), xs.end(), x); - std::vector 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& xs, + const std::vector& 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 \ No newline at end of file diff --git a/src/tallies/filter_energyfunc.cpp b/src/tallies/filter_energyfunc.cpp index f883975a3..a09854e84 100644 --- a/src/tallies/filter_energyfunc.cpp +++ b/src/tallies/filter_energyfunc.cpp @@ -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); From fd7683fe31db6791670c6d512a0739445d9397e6 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 6 Sep 2022 22:50:27 -0500 Subject: [PATCH 21/37] Correction to str literal 'linear-linear' --- openmc/filter.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/openmc/filter.py b/openmc/filter.py index 9815572a3..8b0c84c7e 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -1992,8 +1992,9 @@ class EnergyFunctionFilter(Filter): 'log-log Tabulated1Ds are supported') out = cls(tab1d.x, tab1d.y) + # set interpolation type if tab1d.interpolation[0] == 2: - out.interpolation = 'linear-lienar' + out.interpolation = 'linear-linear' elif tab1d.interpolation[0] == 3: out.interpolation = 'linear-log' elif tab1d.interpolation[0] == 4: From fbce363784e2b9d4b95c246317bccc2230f18703 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Wed, 7 Sep 2022 23:05:15 -0500 Subject: [PATCH 22/37] Using a separate dictionary for interpolation values --- include/openmc/interpolate.h | 3 +++ openmc/data/function.py | 2 +- openmc/filter.py | 10 +++++++--- tests/regression_tests/filter_energyfun/test.py | 2 -- 4 files changed, 11 insertions(+), 6 deletions(-) diff --git a/include/openmc/interpolate.h b/include/openmc/interpolate.h index 2764b7ffc..d02d87ff2 100644 --- a/include/openmc/interpolate.h +++ b/include/openmc/interpolate.h @@ -61,6 +61,9 @@ inline double interpolate(const std::vector& xs, { int idx = lower_bound_index(xs.begin(), xs.end(), x); + if (idx == xs.size()) + idx--; + switch (i) { case Interpolation::lin_lin: return interpolate_lin_lin(xs[idx], xs[idx + 1], ys[idx], ys[idx + 1], x); diff --git a/openmc/data/function.py b/openmc/data/function.py index 59f549d6f..b0390d19c 100644 --- a/openmc/data/function.py +++ b/openmc/data/function.py @@ -13,7 +13,7 @@ from openmc.mixin import EqualityMixin from .data import EV_PER_MEV INTERPOLATION_SCHEME = {1: 'histogram', 2: 'linear-linear', 3: 'linear-log', - 4: 'log-linear', 5: 'log-log', 6 : 'quadratic', 7 : 'cubic'} + 4: 'log-linear', 5: 'log-log'} def sum_functions(funcs): diff --git a/openmc/filter.py b/openmc/filter.py index 8b0c84c7e..c1ce573c9 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -12,7 +12,6 @@ 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 @@ -1896,6 +1895,11 @@ class EnergyFunctionFilter(Filter): """ + # keys selected to match those in function.py where possible + INTERPOLATION_SCHEMES = {2: 'linear-linear', 3: 'linear-log', + 4: 'log-linear', 5: 'log-log', + 6 : 'quadratic', 7 : 'cubic'} + def __init__(self, energy, y, interpolation='linear-linear', filter_id=None): self.energy = energy self.y = y @@ -1963,7 +1967,7 @@ class EnergyFunctionFilter(Filter): out = cls(energy, y, filter_id=filter_id) if 'interpolation' in group: out.interpolation = \ - openmc.data.INTERPOLATION_SCHEME[group['interpolation'][()]] + cls.INTERPOLATION_SCHEMES[group['interpolation'][()]] return out @@ -2053,7 +2057,7 @@ class EnergyFunctionFilter(Filter): @interpolation.setter def interpolation(self, val): cv.check_type('interpolation', val, str) - cv.check_value('interpolation', val, INTERPOLATION_SCHEME.values()) + cv.check_value('interpolation', val, self.INTERPOLATION_SCHEMES.values()) self._interpolation = val def to_xml_element(self): diff --git a/tests/regression_tests/filter_energyfun/test.py b/tests/regression_tests/filter_energyfun/test.py index c297c094a..9c2be7be3 100644 --- a/tests/regression_tests/filter_energyfun/test.py +++ b/tests/regression_tests/filter_energyfun/test.py @@ -73,9 +73,7 @@ class FilterEnergyFunHarness(PyAPITestHarness): br_tally = sp.tallies[2] / sp.tallies[1] dataframes_string += br_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' From 3c3ce6e2e7bafd0dc8e6b28cc3dd10a852d150c8 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 8 Sep 2022 22:28:54 -0500 Subject: [PATCH 23/37] A couple of corrections to the lagrangian interpolation function --- include/openmc/interpolate.h | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/include/openmc/interpolate.h b/include/openmc/interpolate.h index d02d87ff2..3f5d2d361 100644 --- a/include/openmc/interpolate.h +++ b/include/openmc/interpolate.h @@ -42,7 +42,7 @@ inline double interpolate_lagrangian( for (int i = 0; i < order + 1; i++) { double numerator {1.0}; double denominator {1.0}; - for (int j = 0; j < order; j++) { + for (int j = 0; j < order + 1; j++) { if (i == j) continue; numerator *= (x - xs[idx + j]); @@ -55,7 +55,7 @@ inline double interpolate_lagrangian( coeffs.begin(), coeffs.end(), ys.begin() + idx, 0.0); } -inline double interpolate(const std::vector& xs, +double interpolate(const std::vector& xs, const std::vector& ys, double x, Interpolation i = Interpolation::lin_lin) { @@ -76,13 +76,13 @@ inline double interpolate(const std::vector& xs, 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); + return interpolate_lagrangian(xs, ys, idx, x, 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); + return interpolate_lagrangian(xs, ys, idx, x, 3); default: fatal_error("Unsupported interpolation"); } From 84276c54bfc1b076591f274f64cd9890acac36d7 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 8 Sep 2022 22:29:25 -0500 Subject: [PATCH 24/37] Adding additional interpolation types to the constructor --- src/tallies/filter_energyfunc.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/tallies/filter_energyfunc.cpp b/src/tallies/filter_energyfunc.cpp index a09854e84..0abebf885 100644 --- a/src/tallies/filter_energyfunc.cpp +++ b/src/tallies/filter_energyfunc.cpp @@ -38,6 +38,10 @@ void EnergyFunctionFilter::from_xml(pugi::xml_node node) interpolation_ = Interpolation::log_lin; } else if (interpolation == "log-log") { interpolation_ = Interpolation::log_log; + } else if (interpolation == "quadratic") { + interpolation_ = Interpolation::quadratic; + } else if (interpolation == "cubic") { + interpolation_ = Interpolation::cubic; } else { fatal_error(fmt::format( "Found invalid interpolation type '{}' on EnergyFunctionFilter {}.", From ac043cb67a3da93d16761d44d23451f15588b7ba Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 8 Sep 2022 22:37:03 -0500 Subject: [PATCH 25/37] Adding checks for minimum number of data points for quadratic/cubic interpolation --- src/tallies/filter_energyfunc.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/tallies/filter_energyfunc.cpp b/src/tallies/filter_energyfunc.cpp index 0abebf885..170c05262 100644 --- a/src/tallies/filter_energyfunc.cpp +++ b/src/tallies/filter_energyfunc.cpp @@ -39,8 +39,12 @@ void EnergyFunctionFilter::from_xml(pugi::xml_node node) } else if (interpolation == "log-log") { interpolation_ = Interpolation::log_log; } else if (interpolation == "quadratic") { + if (energy.size() < 3) + fatal_error(fmt::format("Quadratic interpolation on EnergyFunctionFilter {} requires at least 3 data points.", id())); interpolation_ = Interpolation::quadratic; } else if (interpolation == "cubic") { + if (energy.size() < 4) + fatal_error(fmt::format("Cubic interpolation on EnergyFunctionFilter {} requires at least 4 data points.", id())); interpolation_ = Interpolation::cubic; } else { fatal_error(fmt::format( From 3b00ed3864d2fabff601f31d67994f4828c0af28 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 8 Sep 2022 22:37:57 -0500 Subject: [PATCH 26/37] Adding checks for data size based on interpolation type --- openmc/filter.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/openmc/filter.py b/openmc/filter.py index c1ce573c9..7a85e2d45 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -2058,6 +2058,13 @@ class EnergyFunctionFilter(Filter): def interpolation(self, val): cv.check_type('interpolation', val, str) cv.check_value('interpolation', val, self.INTERPOLATION_SCHEMES.values()) + + if val == 'quadratic' and len(self.energy) < 3: + raise ValueError('Quadratic interpolation requires 3 or more values.') + + if val == 'cubic' and len(self.energy) < 4: + raise ValueError('Cubic interpolation requires 3 or more values.') + self._interpolation = val def to_xml_element(self): From b68e6e96e614546d2cad91e8b023f994bdbd5d74 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 8 Sep 2022 22:38:22 -0500 Subject: [PATCH 27/37] Updating tests to include checks for quadratic/cubic interpolation --- .../filter_energyfun/inputs_true.dat | 30 +++++++++++++++++++ .../filter_energyfun/results_true.dat | 2 ++ .../regression_tests/filter_energyfun/test.py | 29 +++++++++++++++--- 3 files changed, 57 insertions(+), 4 deletions(-) diff --git a/tests/regression_tests/filter_energyfun/inputs_true.dat b/tests/regression_tests/filter_energyfun/inputs_true.dat index 8c3db3518..133b3167f 100644 --- a/tests/regression_tests/filter_energyfun/inputs_true.dat +++ b/tests/regression_tests/filter_energyfun/inputs_true.dat @@ -39,6 +39,21 @@ 0.1 0.1 0.1333 0.158 0.18467 0.25618 0.4297 0.48 0.48 log-linear + + 0.0 5000000.0 10000000.0 15000000.0 + 0.2 0.7 0.7 0.2 + linear-linear + + + 0.0 5000000.0 10000000.0 15000000.0 + 0.2 0.7 0.7 0.2 + quadratic + + + 0.0 5000000.0 10000000.0 15000000.0 + 0.2 0.7 0.7 0.2 + cubic + Am241 (n,gamma) @@ -63,4 +78,19 @@ Am241 (n,gamma) + + 6 + Am241 + (n,gamma) + + + 7 + Am241 + (n,gamma) + + + 8 + Am241 + (n,gamma) + diff --git a/tests/regression_tests/filter_energyfun/results_true.dat b/tests/regression_tests/filter_energyfun/results_true.dat index c0b59680e..0aadac82e 100644 --- a/tests/regression_tests/filter_energyfun/results_true.dat +++ b/tests/regression_tests/filter_energyfun/results_true.dat @@ -6,3 +6,5 @@ 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 + energyfunction nuclide score mean std. dev. +0 fe168c70d9e078 Am241 (n,gamma) 1.06e-01 2.96e-03 diff --git a/tests/regression_tests/filter_energyfun/test.py b/tests/regression_tests/filter_energyfun/test.py index 9c2be7be3..12144c1fd 100644 --- a/tests/regression_tests/filter_energyfun/test.py +++ b/tests/regression_tests/filter_energyfun/test.py @@ -48,9 +48,21 @@ def model(): filt5 = openmc.EnergyFunctionFilter(x, y) filt5.interpolation = 'log-linear' - filters = [filt1, filt3, filt4, filt5] + # define a trapezoidal function for comparison + x = [0.0, 5e6, 1e7, 1.5e7] + y = [0.2, 0.7, 0.7, 0.2] + + filt6 = openmc.EnergyFunctionFilter(x, y) + + filt7 = openmc.EnergyFunctionFilter(x, y) + filt7.interpolation = 'quadratic' + + filt8 = openmc.EnergyFunctionFilter(x, y) + filt8.interpolation = 'cubic' + + filters = [filt1, filt3, filt4, filt5, filt6, filt7, filt8] # Make tallies - tallies = [openmc.Tally() for i in range(5)] + tallies = [openmc.Tally() for _ in range(len(filters) + 1)] for t in tallies: t.scores = ['(n,gamma)'] t.nuclides = ['Am241'] @@ -73,7 +85,7 @@ class FilterEnergyFunHarness(PyAPITestHarness): br_tally = sp.tallies[2] / sp.tallies[1] dataframes_string += br_tally.get_pandas_dataframe().to_string() + '\n' - for t_id in (3, 4, 5): + for t_id in (3, 4, 5, 6): ef_tally = sp.tallies[t_id] dataframes_string += ef_tally.get_pandas_dataframe().to_string() + '\n' @@ -122,7 +134,16 @@ class FilterEnergyFunHarness(PyAPITestHarness): sp_log_lin_filt = sp_log_lin_tally.find_filter(openmc.EnergyFunctionFilter) assert sp_log_lin_filt.interpolation == 'log-linear' + # check that the cubic interpolation provides a higher value + # than linear-linear + contrived_lin_lin_tally = sp.get_tally(id=6) + contrived_quadratic_tally = sp.get_tally(id=7) + contrived_cubic_tally = sp.get_tally(id=8) + + assert all(contrived_lin_lin_tally.mean < contrived_quadratic_tally.mean) + assert all(contrived_lin_lin_tally.mean < contrived_cubic_tally.mean) + def test_filter_energyfun(model): harness = FilterEnergyFunHarness('statepoint.5.h5', model) - harness.main() + harness.main() \ No newline at end of file From 2fa63f33020da122e52460f64f4f4963ebc89209 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 8 Sep 2022 22:41:54 -0500 Subject: [PATCH 28/37] Applying clang format --- include/openmc/interpolate.h | 19 ++++++++++--------- src/tallies/filter_energyfunc.cpp | 9 +++++++-- 2 files changed, 17 insertions(+), 11 deletions(-) diff --git a/include/openmc/interpolate.h b/include/openmc/interpolate.h index 3f5d2d361..8ccdc5bd7 100644 --- a/include/openmc/interpolate.h +++ b/include/openmc/interpolate.h @@ -34,8 +34,8 @@ inline double interpolate_log_log( return y0 * exp(f * log(y1 / y0)); } -inline double interpolate_lagrangian( - const std::vector& xs, const std::vector& ys, int idx, double x, int order) +inline double interpolate_lagrangian(const std::vector& xs, + const std::vector& ys, int idx, double x, int order) { std::vector coeffs; @@ -55,9 +55,8 @@ inline double interpolate_lagrangian( coeffs.begin(), coeffs.end(), ys.begin() + idx, 0.0); } -double interpolate(const std::vector& xs, - const std::vector& ys, double x, - Interpolation i = Interpolation::lin_lin) +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); @@ -75,20 +74,22 @@ double interpolate(const std::vector& xs, 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--; + if (idx == xs.size() - 2 && idx > 0) + idx--; return interpolate_lagrangian(xs, ys, idx, x, 2); case Interpolation::cubic: // if x is not in the first interval of the x-grid, move back one - if (idx > 0) idx--; + 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--; + if (idx == xs.size() - 3) + idx--; return interpolate_lagrangian(xs, ys, idx, x, 3); default: fatal_error("Unsupported interpolation"); } } - } // namespace openmc #endif \ No newline at end of file diff --git a/src/tallies/filter_energyfunc.cpp b/src/tallies/filter_energyfunc.cpp index 170c05262..bf2b2e6cb 100644 --- a/src/tallies/filter_energyfunc.cpp +++ b/src/tallies/filter_energyfunc.cpp @@ -40,11 +40,16 @@ void EnergyFunctionFilter::from_xml(pugi::xml_node node) interpolation_ = Interpolation::log_log; } else if (interpolation == "quadratic") { if (energy.size() < 3) - fatal_error(fmt::format("Quadratic interpolation on EnergyFunctionFilter {} requires at least 3 data points.", id())); + fatal_error( + fmt::format("Quadratic interpolation on EnergyFunctionFilter {} " + "requires at least 3 data points.", + id())); interpolation_ = Interpolation::quadratic; } else if (interpolation == "cubic") { if (energy.size() < 4) - fatal_error(fmt::format("Cubic interpolation on EnergyFunctionFilter {} requires at least 4 data points.", id())); + fatal_error(fmt::format("Cubic interpolation on EnergyFunctionFilter " + "{} requires at least 4 data points.", + id())); interpolation_ = Interpolation::cubic; } else { fatal_error(fmt::format( From 53d3cd005f2267793b9489f0f39a5b8843468f67 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sat, 17 Sep 2022 09:46:33 -0500 Subject: [PATCH 29/37] Apply suggestions from @paulromano Co-authored-by: Paul Romano --- include/openmc/interpolate.h | 9 ++++----- tests/regression_tests/filter_energyfun/test.py | 2 +- 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/include/openmc/interpolate.h b/include/openmc/interpolate.h index 8ccdc5bd7..fd95e96a8 100644 --- a/include/openmc/interpolate.h +++ b/include/openmc/interpolate.h @@ -1,4 +1,3 @@ - #ifndef OPENMC_INTERPOLATE_H #define OPENMC_INTERPOLATE_H @@ -18,20 +17,20 @@ inline double interpolate_lin_lin( 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); + return y0 + std::log(x / x0) / std::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)); + return y0 * std::exp((x - x0) / (x1 - x0) * std::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)); + double f = std::log(x / x0) / std::log(x1 / x0); + return y0 * std::exp(f * std::log(y1 / y0)); } inline double interpolate_lagrangian(const std::vector& xs, diff --git a/tests/regression_tests/filter_energyfun/test.py b/tests/regression_tests/filter_energyfun/test.py index 12144c1fd..44a2c44b0 100644 --- a/tests/regression_tests/filter_energyfun/test.py +++ b/tests/regression_tests/filter_energyfun/test.py @@ -42,7 +42,7 @@ def model(): filt3 = openmc.EnergyFunctionFilter(x, y) filt3.interpolation = 'log-log' - filt4 = openmc.EnergyFunctionFilter(x , y) + filt4 = openmc.EnergyFunctionFilter(x, y) filt4.interpolation = 'linear-log' filt5 = openmc.EnergyFunctionFilter(x, y) From 3c9e359a575741d98278dfd0e5a2c6e12741c870 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sat, 17 Sep 2022 09:46:00 -0500 Subject: [PATCH 30/37] Addressing some comments from @paulromano --- include/openmc/constants.h | 6 ++-- include/openmc/interpolate.h | 14 +++++--- openmc/filter.py | 32 ++++++++----------- src/tallies/filter_energyfunc.cpp | 4 ++- .../regression_tests/filter_energyfun/test.py | 2 +- 5 files changed, 30 insertions(+), 28 deletions(-) diff --git a/include/openmc/constants.h b/include/openmc/constants.h index 32cd8186a..fa2251b84 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -319,8 +319,10 @@ enum class Interpolation { lin_log = 3, log_lin = 4, log_log = 5, - quadratic = 6, - cubic = 7 + // skip 6 b/c ENDF-6 reserves this value for + // "special one-dimensional interpolation law" + quadratic = 7, + cubic = 8 }; enum class RunMode { diff --git a/include/openmc/interpolate.h b/include/openmc/interpolate.h index fd95e96a8..33167fa2a 100644 --- a/include/openmc/interpolate.h +++ b/include/openmc/interpolate.h @@ -33,10 +33,12 @@ inline double interpolate_log_log( return y0 * std::exp(f * std::log(y1 / y0)); } -inline double interpolate_lagrangian(const std::vector& xs, - const std::vector& ys, int idx, double x, int order) +inline double interpolate_lagrangian(gsl::span xs, + gsl::span ys, int idx, double x, int order) { - std::vector coeffs; + Expects(order <= 3); + std::array coeffs; + coeffs.fill(0.0); for (int i = 0; i < order + 1; i++) { double numerator {1.0}; @@ -47,14 +49,14 @@ inline double interpolate_lagrangian(const std::vector& xs, numerator *= (x - xs[idx + j]); denominator *= (xs[idx + i] - xs[idx + j]); } - coeffs.push_back(numerator / denominator); + coeffs[i] = numerator / denominator; } return std::inner_product( coeffs.begin(), coeffs.end(), ys.begin() + idx, 0.0); } -double interpolate(const std::vector& xs, const std::vector& ys, +double interpolate(gsl::span xs, gsl::span ys, double x, Interpolation i = Interpolation::lin_lin) { int idx = lower_bound_index(xs.begin(), xs.end(), x); @@ -63,6 +65,8 @@ double interpolate(const std::vector& xs, const std::vector& ys, idx--; switch (i) { + case Interpolation::histogram: + return ys[idx]; case Interpolation::lin_lin: return interpolate_lin_lin(xs[idx], xs[idx + 1], ys[idx], ys[idx + 1], x); case Interpolation::log_log: diff --git a/openmc/filter.py b/openmc/filter.py index 7a85e2d45..254bf5245 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -1887,7 +1887,9 @@ class EnergyFunctionFilter(Filter): y : iterable of Real A grid of interpolant values in [eV] interpolation : str - The type of interpolation to be used. + The type of interpolation to be used. One of + ('histogram', 'linear-linear', 'linear-log', 'log-linear', + 'log-log', 'quadratic', 'cubic') id : int Unique identifier for the filter num_bins : Integral @@ -1896,9 +1898,12 @@ class EnergyFunctionFilter(Filter): """ # keys selected to match those in function.py where possible - INTERPOLATION_SCHEMES = {2: 'linear-linear', 3: 'linear-log', - 4: 'log-linear', 5: 'log-log', - 6 : 'quadratic', 7 : 'cubic'} + # skip 6 b/c ENDF-6 reserves this value for + # "special one-dimensional interpolation law" + INTERPOLATION_SCHEMES = {1: 'histogram', 2: 'linear-linear', + 3: 'linear-log', 4: 'log-linear', + 5: 'log-log', 7: 'quadratic', + 8: 'cubic'} def __init__(self, energy, y, interpolation='linear-linear', filter_id=None): self.energy = energy @@ -1991,22 +1996,11 @@ class EnergyFunctionFilter(Filter): if tab1d.n_regions > 1: raise ValueError('Only Tabulated1Ds with a single interpolation ' 'region are supported') - if tab1d.interpolation[0] not in (2, 3, 4, 5): - raise ValueError('Only linear-linear, linear-log, log-linear, and ' + interpolation = tab1d.interpolation[0] + if interpolation not in cls.INTERPOLATION_SCHEMES.keys(): + raise ValueError('Only histogram, linear-linear, linear-log, log-linear, and ' 'log-log Tabulated1Ds are supported') - out = cls(tab1d.x, tab1d.y) - - # set interpolation type - if tab1d.interpolation[0] == 2: - out.interpolation = 'linear-linear' - elif tab1d.interpolation[0] == 3: - out.interpolation = 'linear-log' - elif tab1d.interpolation[0] == 4: - out.interpolation = 'log-linear' - elif tab1d.interpolation[0] == 5: - out.interpolation = 'log-log' - - return out + return cls(tab1d.x, tab1d.y, interpolation) @property def energy(self): diff --git a/src/tallies/filter_energyfunc.cpp b/src/tallies/filter_energyfunc.cpp index bf2b2e6cb..e5a6bc9b7 100644 --- a/src/tallies/filter_energyfunc.cpp +++ b/src/tallies/filter_energyfunc.cpp @@ -30,7 +30,9 @@ void EnergyFunctionFilter::from_xml(pugi::xml_node node) interpolation_ = Interpolation::lin_lin; if (check_for_node(node, "interpolation")) { std::string interpolation = get_node_value(node, "interpolation"); - if (interpolation == "linear-linear") { + if (interpolation == "histogram") { + interpolation_ = Interpolation::histogram; + } else if (interpolation == "linear-linear") { interpolation_ = Interpolation::lin_lin; } else if (interpolation == "linear-log") { interpolation_ = Interpolation::lin_log; diff --git a/tests/regression_tests/filter_energyfun/test.py b/tests/regression_tests/filter_energyfun/test.py index 44a2c44b0..4a42ce0cb 100644 --- a/tests/regression_tests/filter_energyfun/test.py +++ b/tests/regression_tests/filter_energyfun/test.py @@ -31,7 +31,7 @@ def model(): assert filt1.interpolation == 'linear-linear' with pytest.raises(ValueError): - filt1.interpolation = '5th order polynomial' + filt1.interpolation = '🥏' # Also make a filter with the .from_tabulated1d constructor. Make sure # the filters are identical. From 848c0f3e0b408e01f14287be2a70283debb76a67 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sat, 17 Sep 2022 09:50:31 -0500 Subject: [PATCH 31/37] Update to eff docstring --- openmc/filter.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/openmc/filter.py b/openmc/filter.py index 254bf5245..326572c4a 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -1887,9 +1887,8 @@ class EnergyFunctionFilter(Filter): y : iterable of Real A grid of interpolant values in [eV] interpolation : str - The type of interpolation to be used. One of - ('histogram', 'linear-linear', 'linear-log', 'log-linear', - 'log-log', 'quadratic', 'cubic') + Interpolation scheme: {'histogram', 'linear-linear', 'linear-log', + 'log-linear', 'log-log', 'quadratic', 'cubic'} id : int Unique identifier for the filter num_bins : Integral From f33ea966b21802e206bb46c04c7aa8e45943e01b Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sat, 17 Sep 2022 11:10:27 -0500 Subject: [PATCH 32/37] Compute interpolation value in loop. Correction to from_tab1D. --- include/openmc/interpolate.h | 9 +++------ openmc/filter.py | 6 +++--- 2 files changed, 6 insertions(+), 9 deletions(-) diff --git a/include/openmc/interpolate.h b/include/openmc/interpolate.h index 33167fa2a..05353f130 100644 --- a/include/openmc/interpolate.h +++ b/include/openmc/interpolate.h @@ -36,9 +36,7 @@ inline double interpolate_log_log( inline double interpolate_lagrangian(gsl::span xs, gsl::span ys, int idx, double x, int order) { - Expects(order <= 3); - std::array coeffs; - coeffs.fill(0.0); + double output {0.0}; for (int i = 0; i < order + 1; i++) { double numerator {1.0}; @@ -49,11 +47,10 @@ inline double interpolate_lagrangian(gsl::span xs, numerator *= (x - xs[idx + j]); denominator *= (xs[idx + i] - xs[idx + j]); } - coeffs[i] = numerator / denominator; + output += (numerator / denominator) * ys[i]; } - return std::inner_product( - coeffs.begin(), coeffs.end(), ys.begin() + idx, 0.0); + return output; } double interpolate(gsl::span xs, gsl::span ys, diff --git a/openmc/filter.py b/openmc/filter.py index 326572c4a..2af5fe5b7 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -1995,11 +1995,11 @@ class EnergyFunctionFilter(Filter): if tab1d.n_regions > 1: raise ValueError('Only Tabulated1Ds with a single interpolation ' 'region are supported') - interpolation = tab1d.interpolation[0] - if interpolation not in cls.INTERPOLATION_SCHEMES.keys(): + interpolation_val = tab1d.interpolation[0] + if interpolation_val not in cls.INTERPOLATION_SCHEMES.keys(): raise ValueError('Only histogram, linear-linear, linear-log, log-linear, and ' 'log-log Tabulated1Ds are supported') - return cls(tab1d.x, tab1d.y, interpolation) + return cls(tab1d.x, tab1d.y, cls.INTERPOLATION_SCHEMES[interpolation_val]) @property def energy(self): From 834450066e9926d0441ca027dda8c23eb2798174 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sat, 17 Sep 2022 11:31:46 -0500 Subject: [PATCH 33/37] Adding check of interpolation parameter in eff __eq__ --- openmc/filter.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/openmc/filter.py b/openmc/filter.py index 2af5fe5b7..35a65f4fe 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -1913,6 +1913,8 @@ class EnergyFunctionFilter(Filter): def __eq__(self, other): if type(self) is not type(other): return False + elif not self.interpolation == other.interpolation: + return False elif not all(self.energy == other.energy): return False else: From 3351918394cad106ce6ca35fa3fc4b53a2389e87 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sat, 17 Sep 2022 11:32:24 -0500 Subject: [PATCH 34/37] Adding a histogram filter test. Testing conversion from Tabulated1D for many interpolation types --- .../filter_energyfun/inputs_true.dat | 10 ++++++++++ .../regression_tests/filter_energyfun/test.py | 18 +++++++++++++++++- 2 files changed, 27 insertions(+), 1 deletion(-) diff --git a/tests/regression_tests/filter_energyfun/inputs_true.dat b/tests/regression_tests/filter_energyfun/inputs_true.dat index 133b3167f..db7f91ed0 100644 --- a/tests/regression_tests/filter_energyfun/inputs_true.dat +++ b/tests/regression_tests/filter_energyfun/inputs_true.dat @@ -54,6 +54,11 @@ 0.2 0.7 0.7 0.2 cubic + + 0.0 5000000.0 10000000.0 15000000.0 + 0.2 0.7 0.7 0.2 + histogram + Am241 (n,gamma) @@ -93,4 +98,9 @@ Am241 (n,gamma) + + 9 + Am241 + (n,gamma) + diff --git a/tests/regression_tests/filter_energyfun/test.py b/tests/regression_tests/filter_energyfun/test.py index 4a42ce0cb..489eec8c0 100644 --- a/tests/regression_tests/filter_energyfun/test.py +++ b/tests/regression_tests/filter_energyfun/test.py @@ -28,6 +28,7 @@ def model(): # Make an EnergyFunctionFilter directly from the x and y lists. filt1 = openmc.EnergyFunctionFilter(x, y) + # check interpolatoin property setter assert filt1.interpolation == 'linear-linear' with pytest.raises(ValueError): @@ -60,7 +61,10 @@ def model(): filt8 = openmc.EnergyFunctionFilter(x, y) filt8.interpolation = 'cubic' - filters = [filt1, filt3, filt4, filt5, filt6, filt7, filt8] + filt9 = openmc.EnergyFunctionFilter(x, y) + filt9.interpolation = 'histogram' + + filters = [filt1, filt3, filt4, filt5, filt6, filt7, filt8, filt9] # Make tallies tallies = [openmc.Tally() for _ in range(len(filters) + 1)] for t in tallies: @@ -72,6 +76,14 @@ def model(): model.tallies.extend(tallies) + interpolation_vals = \ + list(openmc.EnergyFunctionFilter.INTERPOLATION_SCHEMES.keys()) + for i_val in interpolation_vals: + # breakpoint here is fake and unused + t1d = openmc.data.Tabulated1D(x, y, breakpoints=[1], interpolation=[i_val]) + f = openmc.EnergyFunctionFilter.from_tabulated1d(t1d) + assert f.interpolation == openmc.EnergyFunctionFilter.INTERPOLATION_SCHEMES[i_val] + return model @@ -143,6 +155,10 @@ class FilterEnergyFunHarness(PyAPITestHarness): assert all(contrived_lin_lin_tally.mean < contrived_quadratic_tally.mean) assert all(contrived_lin_lin_tally.mean < contrived_cubic_tally.mean) + # check that the histogram tally is less than the quadratic/cubic interpolations + histogram_tally = sp.get_tally(id=9) + assert all(histogram_tally.mean < contrived_quadratic_tally.mean) + assert all(histogram_tally.mean < contrived_cubic_tally.mean) def test_filter_energyfun(model): harness = FilterEnergyFunHarness('statepoint.5.h5', model) From 8228db6d47d22367eaafbad7ba32e417341c320f Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sat, 17 Sep 2022 12:05:48 -0500 Subject: [PATCH 35/37] Writing interpolation attribute as an attribute of the 'y' dataset. --- docs/source/io_formats/statepoint.rst | 5 +++-- openmc/filter.py | 7 ++++--- src/tallies/filter_energyfunc.cpp | 5 +++-- 3 files changed, 10 insertions(+), 7 deletions(-) diff --git a/docs/source/io_formats/statepoint.rst b/docs/source/io_formats/statepoint.rst index 09a531ecd..a05034a2d 100644 --- a/docs/source/io_formats/statepoint.rst +++ b/docs/source/io_formats/statepoint.rst @@ -108,8 +108,9 @@ 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. Only used for - 'energyfunction' filters. + + :Attributes: - **interpolation** (*int*) -- Interpolation type. Only used for + 'energyfunction' filters. **/tallies/derivatives/derivative /** diff --git a/openmc/filter.py b/openmc/filter.py index 35a65f4fe..1928f5ddb 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -1967,13 +1967,14 @@ class EnergyFunctionFilter(Filter): + group['type'][()].decode() + " instead") energy = group['energy'][()] - y = group['y'][()] + y_grp = group['y'] + y = y_grp[()] filter_id = int(group.name.split('/')[-1].lstrip('filter ')) out = cls(energy, y, filter_id=filter_id) - if 'interpolation' in group: + if 'interpolation' in y_grp.attrs: out.interpolation = \ - cls.INTERPOLATION_SCHEMES[group['interpolation'][()]] + cls.INTERPOLATION_SCHEMES[y_grp.attrs['interpolation'][()]] return out diff --git a/src/tallies/filter_energyfunc.cpp b/src/tallies/filter_energyfunc.cpp index e5a6bc9b7..b961105a2 100644 --- a/src/tallies/filter_energyfunc.cpp +++ b/src/tallies/filter_energyfunc.cpp @@ -104,8 +104,9 @@ void EnergyFunctionFilter::to_statepoint(hid_t filter_group) const Filter::to_statepoint(filter_group); write_dataset(filter_group, "energy", energy_); write_dataset(filter_group, "y", y_); - write_dataset( - filter_group, "interpolation", static_cast(interpolation_)); + hid_t y_dataset = open_dataset(filter_group, "y"); + write_attribute(y_dataset, "interpolation", static_cast(interpolation_)); + close_dataset(y_dataset); } std::string EnergyFunctionFilter::text_label(int bin) const From a6c5dd4d8acccb21c5b5eaafc0a25279f1c7c890 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sat, 17 Sep 2022 12:06:12 -0500 Subject: [PATCH 36/37] Formatting --- tests/regression_tests/filter_energyfun/test.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/tests/regression_tests/filter_energyfun/test.py b/tests/regression_tests/filter_energyfun/test.py index 489eec8c0..295b8ebd8 100644 --- a/tests/regression_tests/filter_energyfun/test.py +++ b/tests/regression_tests/filter_energyfun/test.py @@ -80,9 +80,13 @@ def model(): list(openmc.EnergyFunctionFilter.INTERPOLATION_SCHEMES.keys()) for i_val in interpolation_vals: # breakpoint here is fake and unused - t1d = openmc.data.Tabulated1D(x, y, breakpoints=[1], interpolation=[i_val]) + t1d = openmc.data.Tabulated1D(x, + y, + breakpoints=[1], + interpolation=[i_val]) f = openmc.EnergyFunctionFilter.from_tabulated1d(t1d) - assert f.interpolation == openmc.EnergyFunctionFilter.INTERPOLATION_SCHEMES[i_val] + assert f.interpolation == \ + openmc.EnergyFunctionFilter.INTERPOLATION_SCHEMES[i_val] return model From 61a0458d4c27d7e9bdb54abb012c6c0fb54c51a9 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 22 Sep 2022 10:46:55 -0500 Subject: [PATCH 37/37] Include interpolation types in params docstr --- openmc/filter.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/openmc/filter.py b/openmc/filter.py index 1928f5ddb..d83b9f6ba 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -1876,7 +1876,8 @@ class EnergyFunctionFilter(Filter): y : iterable of Real A grid of interpolant values in [eV] interpolation : str - The type of interpolation to be used. + Interpolation scheme: {'histogram', 'linear-linear', 'linear-log', + 'log-linear', 'log-log', 'quadratic', 'cubic'} filter_id : int Unique identifier for the filter