Merge pull request #2268 from paulromano/energyfunc-fix

Fixes for EnergyFunctionFilter

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

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(

include/openmc/tallies/filter_energyfunc.h

@@ -40,8 +40,9 @@ public:
 
   const vector<double>& energy() const { return energy_; }
   const vector<double>& y() const { return y_; }
-  Interpolation interpolation_;
+  Interpolation interpolation() const { return interpolation_; }
   void set_data(gsl::span<const double> energy, gsl::span<const double> y);
+  void set_interpolation(const std::string& interpolation);
 
 private:
   //----------------------------------------------------------------------------
@@ -52,6 +53,9 @@ private:
 
   //! Interpolant values.
   vector<double> y_;
+
+  //! Interpolation scheme
+  Interpolation interpolation_ {Interpolation::lin_lin};
 };
 
 } // namespace openmc

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):

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()

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<double>(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<int>(y_dataset, "interpolation", static_cast<int>(interpolation_));
+  write_attribute<int>(
+    y_dataset, "interpolation", static_cast<int>(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<EnergyFunctionFilter*>(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<EnergyFunctionFilter*>(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<int>(filt->interpolation());
+  return 0;
+}
+
 } // namespace openmc

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

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]