diff --git a/docs/source/io_formats/statepoint.rst b/docs/source/io_formats/statepoint.rst index c1598f2c3..a05034a2d 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. @@ -109,6 +109,9 @@ The current version of the statepoint file format is 17.0. - **y** (*double[]*) -- Interpolant values for energyfunction interpolation. Only used for 'energyfunction' filters. + :Attributes: - **interpolation** (*int*) -- Interpolation type. Only used for + 'energyfunction' filters. + **/tallies/derivatives/derivative /** :Datasets: - **independent variable** (*char[]*) -- Independent variable of diff --git a/include/openmc/constants.h b/include/openmc/constants.h index a7362ea9e..fa2251b84 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -318,7 +318,11 @@ enum class Interpolation { lin_lin = 2, lin_log = 3, log_lin = 4, - log_log = 5 + log_log = 5, + // 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 new file mode 100644 index 000000000..05353f130 --- /dev/null +++ b/include/openmc/interpolate.h @@ -0,0 +1,95 @@ +#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 + 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 * 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 = std::log(x / x0) / std::log(x1 / x0); + return y0 * std::exp(f * std::log(y1 / y0)); +} + +inline double interpolate_lagrangian(gsl::span xs, + gsl::span ys, int idx, double x, int order) +{ + double output {0.0}; + + for (int i = 0; i < order + 1; i++) { + double numerator {1.0}; + double denominator {1.0}; + for (int j = 0; j < order + 1; j++) { + if (i == j) + continue; + numerator *= (x - xs[idx + j]); + denominator *= (xs[idx + i] - xs[idx + j]); + } + output += (numerator / denominator) * ys[i]; + } + + return output; +} + +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); + + if (idx == xs.size()) + 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: + 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, 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, idx, x, 3); + default: + fatal_error("Unsupported interpolation"); + } +} + +} // namespace openmc + +#endif \ No newline at end of file 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 7874a775b..d83b9f6ba 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 + Interpolation scheme: {'histogram', 'linear-linear', 'linear-log', + 'log-linear', 'log-log', 'quadratic', 'cubic'} filter_id : int Unique identifier for the filter @@ -1884,6 +1887,9 @@ class EnergyFunctionFilter(Filter): A grid of energy values in [eV] y : iterable of Real A grid of interpolant values in [eV] + interpolation : str + Interpolation scheme: {'histogram', 'linear-linear', 'linear-log', + 'log-linear', 'log-log', 'quadratic', 'cubic'} id : int Unique identifier for the filter num_bins : Integral @@ -1891,14 +1897,25 @@ class EnergyFunctionFilter(Filter): """ - def __init__(self, energy, y, filter_id=None): + # keys selected to match those in function.py where possible + # 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 self.y = y self.id = filter_id + self.interpolation = interpolation 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: @@ -1932,12 +1949,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 @@ -1949,10 +1968,16 @@ 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 ')) - return cls(energy, y, filter_id=filter_id) + out = cls(energy, y, filter_id=filter_id) + if 'interpolation' in y_grp.attrs: + out.interpolation = \ + cls.INTERPOLATION_SCHEMES[y_grp.attrs['interpolation'][()]] + + return out @classmethod def from_tabulated1d(cls, tab1d): @@ -1974,10 +1999,11 @@ 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') - - return cls(tab1d.x, tab1d.y) + 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, cls.INTERPOLATION_SCHEMES[interpolation_val]) @property def energy(self): @@ -1987,6 +2013,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 +2051,19 @@ 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, 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): """Return XML Element representing the Filter. @@ -2040,6 +2083,9 @@ 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 @@ -2047,7 +2093,10 @@ class EnergyFunctionFilter(Filter): 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) + 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 diff --git a/src/tallies/filter_energyfunc.cpp b/src/tallies/filter_energyfunc.cpp index fd595e332..b961105a2 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" @@ -25,6 +26,40 @@ void EnergyFunctionFilter::from_xml(pugi::xml_node node) auto y = get_node_array(node, "y"); + // default to linear-linear interpolation + interpolation_ = Interpolation::lin_lin; + if (check_for_node(node, "interpolation")) { + std::string interpolation = get_node_value(node, "interpolation"); + if (interpolation == "histogram") { + interpolation_ = Interpolation::histogram; + } else 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 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( + "Found invalid interpolation type '{}' on EnergyFunctionFilter {}.", + interpolation, id())); + } + } + this->set_data(energy, y); } @@ -55,15 +90,12 @@ 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()); - // Compute the interpolation factor between the nearest bins. - double f = (p.E_last() - energy_[i]) / (energy_[i + 1] - energy_[i]); + double w = interpolate(energy_, y_, p.E_last(), interpolation_); // 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 +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_); + 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 diff --git a/tests/regression_tests/filter_energyfun/inputs_true.dat b/tests/regression_tests/filter_energyfun/inputs_true.dat index 0f506f3f5..db7f91ed0 100644 --- a/tests/regression_tests/filter_energyfun/inputs_true.dat +++ b/tests/regression_tests/filter_energyfun/inputs_true.dat @@ -22,6 +22,42 @@ 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 + + + 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 + + + 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 + + + 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 + + + 0.0 5000000.0 10000000.0 15000000.0 + 0.2 0.7 0.7 0.2 + histogram Am241 @@ -32,4 +68,39 @@ Am241 (n,gamma) + + 3 + Am241 + (n,gamma) + + + 4 + Am241 + (n,gamma) + + + 5 + Am241 + (n,gamma) + + + 6 + Am241 + (n,gamma) + + + 7 + Am241 + (n,gamma) + + + 8 + Am241 + (n,gamma) + + + 9 + 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..0aadac82e 100644 --- a/tests/regression_tests/filter_energyfun/results_true.dat +++ b/tests/regression_tests/filter_energyfun/results_true.dat @@ -1,2 +1,10 @@ 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 + 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 + 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 f61f05d78..295b8ebd8 100644 --- a/tests/regression_tests/filter_energyfun/test.py +++ b/tests/regression_tests/filter_energyfun/test.py @@ -28,20 +28,66 @@ 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): + filt1.interpolation = '🥏' + # Also make a filter with the .from_tabulated1d constructor. Make sure # the filters are identical. tab1d = openmc.data.Tabulated1D(x, y) filt2 = openmc.EnergyFunctionFilter.from_tabulated1d(tab1d) assert filt1 == filt2, 'Error with the .from_tabulated1d constructor' + 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' + + # 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' + + filt9 = openmc.EnergyFunctionFilter(x, y) + filt9.interpolation = 'histogram' + + filters = [filt1, filt3, filt4, filt5, filt6, filt7, filt8, filt9] # Make tallies - tallies = [openmc.Tally(), openmc.Tally()] + tallies = [openmc.Tally() for _ in range(len(filters) + 1)] for t in tallies: t.scores = ['(n,gamma)'] t.nuclides = ['Am241'] - tallies[1].filters = [filt1] + + for t, f in zip(tallies[1:], filters): + t.filters = [f] + 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 @@ -50,13 +96,74 @@ 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' + + for t_id in (3, 4, 5, 6): + ef_tally = sp.tallies[t_id] + dataframes_string += ef_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() + + # Read the statepoint file. + 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) + + model_lin_lin_tally = self._model.tallies[1] + model_lin_lin_filt = model_lin_lin_tally.find_filter(openmc.EnergyFunctionFilter) + + 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) + + 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' + + # 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) + + # 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) - harness.main() + harness.main() \ No newline at end of file