Adding test results for additional interpolation types

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 <id>/**
 

include/openmc/interpolate.h

@@ -0,0 +1,65 @@
+
+#ifndef OPENMC_INTERPOLATE_H
+#define OPENMC_INTERPOLATE_H
+
+#include <cmath>
+#include <vector>
+
+#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<double>& xs,
+  const std::vector<double>& 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<double>& xs,
+  const std::vector<double>& 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

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

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

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]);

src/tallies/filter_energyfunc.cpp

@@ -3,6 +3,7 @@
 #include <fmt/core.h>
 
 #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(

tests/regression_tests/filter_energyfun/inputs_true.dat

@@ -29,6 +29,16 @@
     <y>0.1 0.1 0.1333 0.158 0.18467 0.25618 0.4297 0.48 0.48</y>
     <interpolation>log-log</interpolation>
   </filter>
+  <filter id="4" type="energyfunction">
+    <energy>1e-05 0.369 1000.0 100000.0 600000.0 1000000.0 2000000.0 4000000.0 30000000.0</energy>
+    <y>0.1 0.1 0.1333 0.158 0.18467 0.25618 0.4297 0.48 0.48</y>
+    <interpolation>linear-log</interpolation>
+  </filter>
+  <filter id="5" type="energyfunction">
+    <energy>1e-05 0.369 1000.0 100000.0 600000.0 1000000.0 2000000.0 4000000.0 30000000.0</energy>
+    <y>0.1 0.1 0.1333 0.158 0.18467 0.25618 0.4297 0.48 0.48</y>
+    <interpolation>log-linear</interpolation>
+  </filter>
   <tally id="1">
     <nuclides>Am241</nuclides>
     <scores>(n,gamma)</scores>
@@ -43,4 +53,14 @@
     <nuclides>Am241</nuclides>
     <scores>(n,gamma)</scores>
   </tally>
+  <tally id="4">
+    <filters>4</filters>
+    <nuclides>Am241</nuclides>
+    <scores>(n,gamma)</scores>
+  </tally>
+  <tally id="5">
+    <filters>5</filters>
+    <nuclides>Am241</nuclides>
+    <scores>(n,gamma)</scores>
+  </tally>
 </tallies>

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

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)