Add more tests for FCI operators

CMakeLists.txt

@@ -376,6 +376,7 @@ if (zoidberg_FOUND EQUAL 0)
 else()
   set(zoidberg_FOUND OFF)
 endif()
+message(STATUS "Found Zoidberg for FCI tests: ${zoidberg_FOUND}")
 
 option(BOUT_GENERATE_FIELDOPS "Automatically re-generate the Field arithmetic operators from the Python templates. \
 Requires Python3, clang-format, and Jinja2. Turn this OFF to skip generating them if, for example, \

include/bout/difops.hxx

@@ -40,7 +40,9 @@
 #include "bout/field3d.hxx"
 
 #include "bout/bout_types.hxx"
-#include "bout/solver.hxx"
+#include "bout/coordinates.hxx"
+
+class Solver;
 
 /*!
  * Parallel derivative (central differencing) in Y

include/bout/interpolation_xz.hxx

@@ -24,7 +24,9 @@
 #ifndef BOUT_INTERP_XZ_H
 #define BOUT_INTERP_XZ_H
 
-#include "bout/mask.hxx"
+#include <bout/bout_types.hxx>
+#include <bout/generic_factory.hxx>
+#include <bout/mask.hxx>
 
 #define USE_NEW_WEIGHTS 1
 #if BOUT_HAS_PETSC
@@ -166,7 +168,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 +213,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");
     }
@@ -248,7 +271,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) {
@@ -284,7 +308,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) {
@@ -308,18 +333,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/coordinates.cxx

@@ -1577,7 +1577,7 @@ Field3D Coordinates::Div_par(const Field3D& f, CELL_LOC outloc,
 
   // Need Bxy at location of f, which might be different from location of this
   // Coordinates object
-  auto Bxy_floc = f.getCoordinates()->Bxy;
+  const auto& Bxy_floc = f.getCoordinates()->Bxy;
 
   if (!f.hasParallelSlices()) {
     // No yup/ydown fields. The Grad_par operator will

src/mesh/interpolation/hermite_spline_xz.cxx

@@ -171,7 +171,19 @@ void XZHermiteSpline::calcWeights(const Field3D& delta_x, const Field3D& delta_z
     BoutReal t_x = delta_x(x, y, z) - static_cast<BoutReal>(i_corn);
     BoutReal t_z = delta_z(x, y, z) - static_cast<BoutReal>(k_corner(x, y, z));
 
-    // NOTE: A (small) hack to avoid one-sided differences
+    // NOTE: A (small) hack to avoid one-sided differences. We need at
+    // least 2 interior points due to an awkwardness with the
+    // boundaries. The splines need derivatives in x, but we don't
+    // know the value in the boundaries, so _any_ interpolation in the
+    // last interior cell can't be done. Instead, we fudge the
+    // interpolation in the last cell to be at the extreme right-hand
+    // edge of the previous cell (that is, exactly on the last
+    // interior point). However, this doesn't work with only one
+    // interior point, because we have to do something similar to the
+    // _first_ cell, and these two fudges cancel out and we end up
+    // indexing into the boundary anyway.
+    // TODO(peter): Can we remove this if we apply (dirichlet?) BCs to
+    // the X derivatives? Note that we need at least _2_
     if (i_corn >= xend) {
       i_corn = xend - 1;
       t_x = 1.0;

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 {
@@ -80,20 +80,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;
   }