Add rtol, atol to MonotonicHermiteSpline

This allows the user to control the clipping, allowing some overshoot

Also restore the FCI MMS test with this interpolator

Author: ZedThree

include/bout/interpolation_xz.hxx

@@ -24,6 +24,7 @@
 #ifndef BOUT_INTERP_XZ_H
 #define BOUT_INTERP_XZ_H
 
+#include "bout/bout_types.hxx"
 #include "bout/mask.hxx"
 
 #define USE_NEW_WEIGHTS 1
@@ -166,7 +167,8 @@ protected:
 #endif
 
 public:
-  XZHermiteSpline(Mesh* mesh = nullptr) : XZHermiteSpline(0, mesh) {}
+  XZHermiteSpline(Mesh* mesh = nullptr, [[maybe_unused]] Options* options = nullptr)
+      : XZHermiteSpline(0, mesh) {}
   XZHermiteSpline(int y_offset = 0, Mesh* mesh = nullptr);
   XZHermiteSpline(const BoutMask& mask, int y_offset = 0, Mesh* mesh = nullptr)
       : XZHermiteSpline(y_offset, mesh) {
@@ -210,9 +212,29 @@ public:
 /// but also degrades accuracy near maxima and minima.
 /// Perhaps should only impose near boundaries, since that is where
 /// problems most obviously occur.
+///
+/// You can control how tight the clipping to the range of the neighbouring cell
+/// values through ``rtol`` and ``atol``:
+///
+///     diff = (max_of_neighours - min_of_neighours) * rtol + atol
+///
+/// and the interpolated value is instead clipped to the range
+/// ``[min_of_neighours - diff, max_of_neighours + diff]``
 class XZMonotonicHermiteSpline : public XZHermiteSpline {
+  /// Absolute tolerance for clipping
+  BoutReal atol = 0.0;
+  /// Relative tolerance for clipping
+  BoutReal rtol = 1.0;
+
 public:
-  XZMonotonicHermiteSpline(Mesh* mesh = nullptr) : XZHermiteSpline(0, mesh) {
+  XZMonotonicHermiteSpline(Mesh* mesh = nullptr, Options* options = nullptr)
+      : XZHermiteSpline(0, mesh),
+        atol{(*options)["atol"]
+                 .doc("Absolute tolerance for clipping overshoot")
+                 .withDefault(0.0)},
+        rtol{(*options)["rtol"]
+                 .doc("Relative tolerance for clipping overshoot")
+                 .withDefault(1.0)} {
     if (localmesh->getNXPE() > 1) {
       throw BoutException("Do not support MPI splitting in X");
     }
@@ -245,7 +267,8 @@ class XZLagrange4pt : public XZInterpolation {
   Field3D t_x, t_z;
 
 public:
-  XZLagrange4pt(Mesh* mesh = nullptr) : XZLagrange4pt(0, mesh) {}
+  XZLagrange4pt(Mesh* mesh = nullptr, [[maybe_unused]] Options* options = nullptr)
+      : XZLagrange4pt(0, mesh) {}
   XZLagrange4pt(int y_offset = 0, Mesh* mesh = nullptr);
   XZLagrange4pt(const BoutMask& mask, int y_offset = 0, Mesh* mesh = nullptr)
       : XZLagrange4pt(y_offset, mesh) {
@@ -278,7 +301,8 @@ class XZBilinear : public XZInterpolation {
   Field3D w0, w1, w2, w3;
 
 public:
-  XZBilinear(Mesh* mesh = nullptr) : XZBilinear(0, mesh) {}
+  XZBilinear(Mesh* mesh = nullptr, [[maybe_unused]] Options* options = nullptr)
+      : XZBilinear(0, mesh) {}
   XZBilinear(int y_offset = 0, Mesh* mesh = nullptr);
   XZBilinear(const BoutMask& mask, int y_offset = 0, Mesh* mesh = nullptr)
       : XZBilinear(y_offset, mesh) {
@@ -302,18 +326,18 @@ public:
 };
 
 class XZInterpolationFactory
-    : public Factory<XZInterpolation, XZInterpolationFactory, Mesh*> {
+    : public Factory<XZInterpolation, XZInterpolationFactory, Mesh*, Options*> {
 public:
   static constexpr auto type_name = "XZInterpolation";
   static constexpr auto section_name = "xzinterpolation";
   static constexpr auto option_name = "type";
   static constexpr auto default_type = "hermitespline";
 
   ReturnType create(Options* options = nullptr, Mesh* mesh = nullptr) const {
-    return Factory::create(getType(options), mesh);
+    return Factory::create(getType(options), mesh, options);
   }
-  ReturnType create(const std::string& type, [[maybe_unused]] Options* options) const {
-    return Factory::create(type, nullptr);
+  ReturnType create(const std::string& type, Options* options) const {
+    return Factory::create(type, nullptr, options);
   }
 
   static void ensureRegistered();

src/mesh/interpolation/monotonic_hermite_spline_xz.cxx

@@ -25,7 +25,7 @@
 #include "bout/interpolation_xz.hxx"
 #include "bout/mesh.hxx"
 
-#include <vector>
+#include <algorithm>
 
 Field3D XZMonotonicHermiteSpline::interpolate(const Field3D& f,
                                               const std::string& region) const {
@@ -96,20 +96,17 @@ Field3D XZMonotonicHermiteSpline::interpolate(const Field3D& f,
     // Perhaps should only impose near boundaries, since that is where
     // problems most obviously occur.
     const BoutReal localmax = BOUTMAX(f[ic], f[icxp], f[iczp], f[icxpzp]);
-
     const BoutReal localmin = BOUTMIN(f[ic], f[icxp], f[iczp], f[icxpzp]);
 
     ASSERT2(std::isfinite(localmax) || i.x() < localmesh->xstart
             || i.x() > localmesh->xend);
     ASSERT2(std::isfinite(localmin) || i.x() < localmesh->xstart
             || i.x() > localmesh->xend);
 
-    if (result > localmax) {
-      result = localmax;
-    }
-    if (result < localmin) {
-      result = localmin;
-    }
+    const auto diff = ((localmax - localmin) * rtol) + atol;
+
+    result = std::min(result, localmax + diff);
+    result = std::max(result, localmin - diff);
 
     f_interp[iyp] = result;
   }

tests/MMS/spatial/fci/runtest

@@ -236,6 +236,17 @@ if __name__ == "__main__":
             "yperiodic": True,
             "args": "mesh:ddy:first=C2 mesh:paralleltransform:xzinterpolation:type=lagrange4pt",
         },
+        "nslice=1 monotonichermitespline": {
+            "nslice": 1,
+            "order": 2,
+            "yperiodic": True,
+            "args": (
+                "mesh:ddy:first=C2 "
+                "mesh:paralleltransform:xzinterpolation:type=monotonichermitespline "
+                "mesh:paralleltransform:xzinterpolation:rtol=1e-3 "
+                "mesh:paralleltransform:xzinterpolation:atol=5e-3"
+            ),
+        },
     }
 
     for name, case in cases.items():

tests/integrated/test-interpolate/runtest

@@ -25,7 +25,6 @@ methods = {
     "hermitespline": 3,
     "lagrange4pt": 3,
     "bilinear": 2,
-    "monotonichermitespline": 2,
 }