From 0587aeeaa4bf834e89eb9309c4729f315d001f03 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 18 Oct 2022 11:58:24 -0500 Subject: [PATCH 1/3] Allow EnergyFunctionFilter interpolation to be set from C API --- include/openmc/capi.h | 3 + include/openmc/tallies/filter_energyfunc.h | 6 +- openmc/lib/filter.py | 26 ++++- src/tallies/filter_energyfunc.cpp | 114 +++++++++++++++------ tests/unit_tests/test_lib.py | 5 + 5 files changed, 119 insertions(+), 35 deletions(-) diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 29e900965..ce58bed12 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -38,6 +38,9 @@ int openmc_energyfunc_filter_get_energy( int openmc_energyfunc_filter_get_y(int32_t index, size_t* n, const double** y); int openmc_energyfunc_filter_set_data( int32_t index, size_t n, const double* energies, const double* y); +int openmc_energyfunc_filter_set_interpolation( + int32_t index, const char* interp); +int openmc_energyfunc_filter_get_interpolation(int32_t index, int* interp); int openmc_extend_cells(int32_t n, int32_t* index_start, int32_t* index_end); int openmc_extend_filters(int32_t n, int32_t* index_start, int32_t* index_end); int openmc_extend_materials( diff --git a/include/openmc/tallies/filter_energyfunc.h b/include/openmc/tallies/filter_energyfunc.h index f650976ef..373fee8dd 100644 --- a/include/openmc/tallies/filter_energyfunc.h +++ b/include/openmc/tallies/filter_energyfunc.h @@ -40,8 +40,9 @@ public: const vector& energy() const { return energy_; } const vector& y() const { return y_; } - Interpolation interpolation_; + Interpolation interpolation() const { return interpolation_; } void set_data(gsl::span energy, gsl::span y); + void set_interpolation(const std::string& interpolation); private: //---------------------------------------------------------------------------- @@ -52,6 +53,9 @@ private: //! Interpolant values. vector y_; + + //! Interpolation scheme + Interpolation interpolation_ {Interpolation::lin_lin}; }; } // namespace openmc diff --git a/openmc/lib/filter.py b/openmc/lib/filter.py index 5d2231cc2..04da92459 100644 --- a/openmc/lib/filter.py +++ b/openmc/lib/filter.py @@ -7,6 +7,7 @@ import numpy as np from numpy.ctypeslib import as_array from openmc.exceptions import AllocationError, InvalidIDError +from openmc.data.function import INTERPOLATION_SCHEME from . import _dll from .core import _FortranObjectWithID from .error import _error_handler @@ -16,9 +17,9 @@ from .mesh import _get_mesh __all__ = [ 'Filter', 'AzimuthalFilter', 'CellFilter', 'CellbornFilter', 'CellfromFilter', - 'CellInstanceFilter', 'CollisionFilter', 'DistribcellFilter', 'DelayedGroupFilter', - 'EnergyFilter', 'EnergyoutFilter', 'EnergyFunctionFilter', 'LegendreFilter', - 'MaterialFilter', 'MeshFilter', 'MeshSurfaceFilter', 'MuFilter', 'ParticleFilter', + 'CellInstanceFilter', 'CollisionFilter', 'DistribcellFilter', 'DelayedGroupFilter', + 'EnergyFilter', 'EnergyoutFilter', 'EnergyFunctionFilter', 'LegendreFilter', + 'MaterialFilter', 'MeshFilter', 'MeshSurfaceFilter', 'MuFilter', 'ParticleFilter', 'PolarFilter', 'SphericalHarmonicsFilter', 'SpatialLegendreFilter', 'SurfaceFilter', 'UniverseFilter', 'ZernikeFilter', 'ZernikeRadialFilter', 'filters' ] @@ -47,6 +48,12 @@ _dll.openmc_energyfunc_filter_get_y.resttpe = c_int _dll.openmc_energyfunc_filter_get_y.errcheck = _error_handler _dll.openmc_energyfunc_filter_get_y.argtypes = [ c_int32, POINTER(c_size_t), POINTER(POINTER(c_double))] +_dll.openmc_energyfunc_filter_get_interpolation.resttpe = c_int +_dll.openmc_energyfunc_filter_get_interpolation.errcheck = _error_handler +_dll.openmc_energyfunc_filter_get_interpolation.argtypes = [c_int32, POINTER(c_int)] +_dll.openmc_energyfunc_filter_set_interpolation.resttpe = c_int +_dll.openmc_energyfunc_filter_set_interpolation.errcheck = _error_handler +_dll.openmc_energyfunc_filter_set_interpolation.argtypes = [c_int32, c_char_p] _dll.openmc_filter_get_id.argtypes = [c_int32, POINTER(c_int32)] _dll.openmc_filter_get_id.restype = c_int _dll.openmc_filter_get_id.errcheck = _error_handler @@ -247,6 +254,8 @@ class EnergyFunctionFilter(Filter): Independent variable for the interpolation y : numpy.ndarray Dependent variable for the interpolation + interpolation : {'histogram', 'linear-linear', 'linear-log', 'log-linear', 'log-log', 'quadratic', 'cubic'} + Interpolation scheme """ energy_array = np.asarray(energy) y_array = np.asarray(y) @@ -264,6 +273,17 @@ class EnergyFunctionFilter(Filter): def y(self): return self._get_attr(_dll.openmc_energyfunc_filter_get_y) + @property + def interpolation(self) -> str: + interp = c_int() + _dll.openmc_energyfunc_filter_get_interpolation(self._index, interp) + return INTERPOLATION_SCHEME[interp.value] + + @interpolation.setter + def interpolation(self, interp: str): + interp_ptr = c_char_p(interp.encode()) + _dll.openmc_energyfunc_filter_set_interpolation(self._index, interp_ptr) + def _get_attr(self, cfunc): array_p = POINTER(c_double)() n = c_size_t() diff --git a/src/tallies/filter_energyfunc.cpp b/src/tallies/filter_energyfunc.cpp index b961105a2..93ae24b2a 100644 --- a/src/tallies/filter_energyfunc.cpp +++ b/src/tallies/filter_energyfunc.cpp @@ -25,42 +25,14 @@ void EnergyFunctionFilter::from_xml(pugi::xml_node node) fatal_error("y values not specified for EnergyFunction filter."); auto y = get_node_array(node, "y"); + this->set_data(energy, 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_interpolation(interpolation); } - - this->set_data(energy, y); } void EnergyFunctionFilter::set_data( @@ -86,6 +58,38 @@ void EnergyFunctionFilter::set_data( } } +void EnergyFunctionFilter::set_interpolation(const std::string& 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.", + this->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.", + this->id())); + interpolation_ = Interpolation::cubic; + } else { + fatal_error(fmt::format( + "Found invalid interpolation type '{}' on EnergyFunctionFilter {}.", + interpolation, this->id())); + } +} + void EnergyFunctionFilter::get_all_bins( const Particle& p, TallyEstimator estimator, FilterMatch& match) const { @@ -105,7 +109,8 @@ void EnergyFunctionFilter::to_statepoint(hid_t filter_group) const 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_)); + write_attribute( + y_dataset, "interpolation", static_cast(interpolation_)); close_dataset(y_dataset); } @@ -189,4 +194,51 @@ extern "C" int openmc_energyfunc_filter_get_y( return 0; } +extern "C" int openmc_energyfunc_filter_set_interpolation( + int32_t index, const char* interp) +{ + // ensure this is a valid index to allocated filter + if (int err = verify_filter(index)) + return err; + + // get a pointer to the filter + const auto& filt_base = model::tally_filters[index].get(); + // downcast to EnergyFunctionFilter + auto* filt = dynamic_cast(filt_base); + + // check if a valid filter was produced + if (!filt) { + set_errmsg( + "Tried to set interpolation data for non-energy function filter."); + return OPENMC_E_INVALID_TYPE; + } + + // Set interpolation + filt->set_interpolation(interp); + return 0; +} + +extern "C" int openmc_energyfunc_filter_get_interpolation( + int32_t index, int* interp) +{ + // ensure this is a valid index to allocated filter + if (int err = verify_filter(index)) + return err; + + // get a pointer to the filter + const auto& filt_base = model::tally_filters[index].get(); + // downcast to EnergyFunctionFilter + auto* filt = dynamic_cast(filt_base); + + // check if a valid filter was produced + if (!filt) { + set_errmsg( + "Tried to set interpolation data for non-energy function filter."); + return OPENMC_E_INVALID_TYPE; + } + + *interp = static_cast(filt->interpolation()); + return 0; +} + } // namespace openmc diff --git a/tests/unit_tests/test_lib.py b/tests/unit_tests/test_lib.py index 5063646f7..0633ceb18 100644 --- a/tests/unit_tests/test_lib.py +++ b/tests/unit_tests/test_lib.py @@ -283,6 +283,11 @@ def test_energy_function_filter(lib_init): assert len(efunc.y) == 2 assert (efunc.y == [0.0, 2.0]).all() + # Default should be lin-lin + assert efunc.interpolation == 'linear-linear' + efunc.interpolation = 'histogram' + assert efunc.interpolation == 'histogram' + def test_tally(lib_init): t = openmc.lib.tallies[1] From 2ae5ccb031c481d4e97c40c843186396c0fd8ffc Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 18 Oct 2022 12:02:02 -0500 Subject: [PATCH 2/3] Fix EnergyFunctionFilter.from_xml_element --- openmc/filter.py | 2 +- tests/unit_tests/test_filters.py | 15 +++++++++++++++ 2 files changed, 16 insertions(+), 1 deletion(-) diff --git a/openmc/filter.py b/openmc/filter.py index d83b9f6ba..d6571cd5a 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -2095,7 +2095,7 @@ class EnergyFunctionFilter(Filter): 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') + out.interpolation = elem.find('interpolation').text return out def can_merge(self, other): diff --git a/tests/unit_tests/test_filters.py b/tests/unit_tests/test_filters.py index a20b129e9..8a6f095ef 100644 --- a/tests/unit_tests/test_filters.py +++ b/tests/unit_tests/test_filters.py @@ -254,3 +254,18 @@ def test_lethargy_bin_width(): energy_bins = openmc.mgxs.GROUP_STRUCTURES['VITAMIN-J-175'] assert f.lethargy_bin_width[0] == np.log10(energy_bins[1]/energy_bins[0]) assert f.lethargy_bin_width[-1] == np.log10(energy_bins[-1]/energy_bins[-2]) + + +def test_energyfunc(): + f = openmc.EnergyFunctionFilter( + [0.0, 10.0, 2.0e3, 1.0e6, 20.0e6], + [1.0, 0.9, 0.8, 0.7, 0.6], + 'histogram' + ) + + # Make sure XML roundtrip works + elem = f.to_xml_element() + new_f = openmc.EnergyFunctionFilter.from_xml_element(elem) + np.testing.assert_allclose(f.energy, new_f.energy) + np.testing.assert_allclose(f.y, new_f.y) + assert f.interpolation == new_f.interpolation From 1814f826c9ec1e621374457886ed7900d82b27ce Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 18 Oct 2022 12:02:22 -0500 Subject: [PATCH 3/3] Remove deprecated cylinder_from_points name in docs --- docs/source/pythonapi/model.rst | 1 - 1 file changed, 1 deletion(-) diff --git a/docs/source/pythonapi/model.rst b/docs/source/pythonapi/model.rst index 5dd3d3293..636935c6b 100644 --- a/docs/source/pythonapi/model.rst +++ b/docs/source/pythonapi/model.rst @@ -11,7 +11,6 @@ Convenience Functions :template: myfunction.rst openmc.model.borated_water - openmc.model.cylinder_from_points openmc.model.hexagonal_prism openmc.model.rectangular_prism openmc.model.subdivide