PETSc XZHermiteSpline

CMakeLists.txt

@@ -270,6 +270,7 @@ set(BOUT_SOURCES
   ./src/mesh/interpolation/interpolation_z.cxx
   ./src/mesh/interpolation/lagrange_4pt_xz.cxx
   ./src/mesh/interpolation/monotonic_hermite_spline_xz.cxx
+  ./src/mesh/interpolation/petsc_hermite_spline_xz.cxx
   ./src/mesh/mesh.cxx
   ./src/mesh/parallel/fci.cxx
   ./src/mesh/parallel/fci.hxx

include/bout/interpolation_xz.hxx

@@ -1,8 +1,7 @@
 /**************************************************************************
- * Copyright 2010-2020 B.D.Dudson, S.Farley, P. Hill, J.T. Omotani, J.T. Parker,
- * M.V.Umansky, X.Q.Xu
+ * Copyright 2010-2024 BOUT++ contributors
  *
- * Contact: Ben Dudson, [email protected]
+ * Contact: Ben Dudson, [email protected]
  *
  * This file is part of BOUT++.
  *
@@ -24,7 +23,13 @@
 #ifndef BOUT_INTERP_XZ_H
 #define BOUT_INTERP_XZ_H
 
+#include "bout/bout_types.hxx"
+#include "bout/build_defines.hxx"
+#include "bout/field3d.hxx"
+#include "bout/globalindexer.hxx"
 #include "bout/mask.hxx"
+#include "bout/petsc_interface.hxx"
+#include "bout/petsclib.hxx"
 
 class Options;
 
@@ -95,20 +100,30 @@ public:
     ASSERT1(has_region);
     return getRegion();
   }
+  /// Calculate weights using given offsets in X and Z
   virtual void calcWeights(const Field3D& delta_x, const Field3D& delta_z,
                            const std::string& region = "RGN_NOBNDRY") = 0;
-  virtual void calcWeights(const Field3D& delta_x, const Field3D& delta_z,
-                           const BoutMask& mask,
-                           const std::string& region = "RGN_NOBNDRY") = 0;
+  void calcWeights(const Field3D& delta_x, const Field3D& delta_z, const BoutMask& mask,
+                   const std::string& region = "RGN_NOBNDRY") {
+    setMask(mask);
+    calcWeights(delta_x, delta_z, region);
+  }
 
+  /// Use pre-calculated weights
   virtual Field3D interpolate(const Field3D& f,
                               const std::string& region = "RGN_NOBNDRY") const = 0;
-  virtual Field3D interpolate(const Field3D& f, const Field3D& delta_x,
-                              const Field3D& delta_z,
-                              const std::string& region = "RGN_NOBNDRY") = 0;
-  virtual Field3D interpolate(const Field3D& f, const Field3D& delta_x,
-                              const Field3D& delta_z, const BoutMask& mask,
-                              const std::string& region = "RGN_NOBNDRY") = 0;
+
+  /// Calculate weights then interpolate
+  Field3D interpolate(const Field3D& f, const Field3D& delta_x, const Field3D& delta_z,
+                      const std::string& region = "RGN_NOBNDRY") {
+    calcWeights(delta_x, delta_z, region);
+    return interpolate(f, region);
+  }
+  Field3D interpolate(const Field3D& f, const Field3D& delta_x, const Field3D& delta_z,
+                      const BoutMask& mask, const std::string& region = "RGN_NOBNDRY") {
+    calcWeights(delta_x, delta_z, mask, region);
+    return interpolate(f, region);
+  }
 
   // Interpolate using the field at (x,y+y_offset,z), rather than (x,y,z)
   void setYOffset(int offset) { y_offset = offset; }
@@ -162,18 +177,10 @@ public:
 
   void calcWeights(const Field3D& delta_x, const Field3D& delta_z,
                    const std::string& region = "RGN_NOBNDRY") override;
-  void calcWeights(const Field3D& delta_x, const Field3D& delta_z, const BoutMask& mask,
-                   const std::string& region = "RGN_NOBNDRY") override;
 
-  // Use precalculated weights
   Field3D interpolate(const Field3D& f,
                       const std::string& region = "RGN_NOBNDRY") const override;
-  // Calculate weights and interpolate
-  Field3D interpolate(const Field3D& f, const Field3D& delta_x, const Field3D& delta_z,
-                      const std::string& region = "RGN_NOBNDRY") override;
-  Field3D interpolate(const Field3D& f, const Field3D& delta_x, const Field3D& delta_z,
-                      const BoutMask& mask,
-                      const std::string& region = "RGN_NOBNDRY") override;
+
   std::vector<ParallelTransform::PositionsAndWeights>
   getWeightsForYApproximation(int i, int j, int k, int yoffset) override;
 };
@@ -218,21 +225,12 @@ public:
 
   void calcWeights(const Field3D& delta_x, const Field3D& delta_z,
                    const std::string& region = "RGN_NOBNDRY") override;
-  void calcWeights(const Field3D& delta_x, const Field3D& delta_z, const BoutMask& mask,
-                   const std::string& region = "RGN_NOBNDRY") override;
+
+  using XZInterpolation::interpolate;
 
   // Use precalculated weights
   Field3D interpolate(const Field3D& f,
                       const std::string& region = "RGN_NOBNDRY") const override;
-  // Calculate weights and interpolate
-  Field3D interpolate(const Field3D& f, const Field3D& delta_x, const Field3D& delta_z,
-                      const std::string& region = "RGN_NOBNDRY") override;
-  Field3D interpolate(const Field3D& f, const Field3D& delta_x, const Field3D& delta_z,
-                      const BoutMask& mask,
-                      const std::string& region = "RGN_NOBNDRY") override;
-  BoutReal lagrange_4pt(BoutReal v2m, BoutReal vm, BoutReal vp, BoutReal v2p,
-                        BoutReal offset) const;
-  BoutReal lagrange_4pt(const BoutReal v[], BoutReal offset) const;
 };
 
 class XZBilinear : public XZInterpolation {
@@ -251,19 +249,73 @@ public:
 
   void calcWeights(const Field3D& delta_x, const Field3D& delta_z,
                    const std::string& region = "RGN_NOBNDRY") override;
-  void calcWeights(const Field3D& delta_x, const Field3D& delta_z, const BoutMask& mask,
+
+  using XZInterpolation::interpolate;
+
+  // Use precalculated weights
+  Field3D interpolate(const Field3D& f,
+                      const std::string& region = "RGN_NOBNDRY") const override;
+};
+
+#if BOUT_HAS_PETSC
+class PetscXZHermiteSpline : public XZInterpolation {
+  PetscLib petsclib;
+  bool isInit{false};
+  IndexerPtr<Field3D> indexer;
+  PetscMatrix<Field3D> weights;
+  Mat petscWeights{nullptr};
+  Vec rhs{nullptr};
+  Vec result{nullptr};
+  std::vector<Field3D> newWeights;
+
+  Tensor<Field3D::ind_type> i_corner; // x-index of bottom-left grid point
+  Tensor<int> k_corner;               // z-index of bottom-left grid point
+
+  // Basis functions for cubic Hermite spline interpolation
+  //    see http://en.wikipedia.org/wiki/Cubic_Hermite_spline
+  // The h00 and h01 basis functions are applied to the function itself
+  // and the h10 and h11 basis functions are applied to its derivative
+  // along the interpolation direction.
+  Field3D h00_x;
+  Field3D h01_x;
+  Field3D h10_x;
+  Field3D h11_x;
+  Field3D h00_z;
+  Field3D h01_z;
+  Field3D h10_z;
+  Field3D h11_z;
+
+public:
+  PetscXZHermiteSpline(int y_offset = 0, Mesh* mesh = nullptr);
+  PetscXZHermiteSpline(Mesh* mesh = nullptr) : PetscXZHermiteSpline(0, mesh) {}
+  PetscXZHermiteSpline(const BoutMask& mask, int y_offset = 0, Mesh* mesh = nullptr)
+      : PetscXZHermiteSpline(y_offset, mesh) {
+    region = regionFromMask(mask, localmesh);
+  }
+  PetscXZHermiteSpline(const PetscXZHermiteSpline& other) = delete;
+  PetscXZHermiteSpline(PetscXZHermiteSpline&& other) = default;
+  PetscXZHermiteSpline& operator=(const PetscXZHermiteSpline& other) = delete;
+  PetscXZHermiteSpline& operator=(PetscXZHermiteSpline&& other) = delete;
+  ~PetscXZHermiteSpline() override;
+
+  void calcWeights(const Field3D& delta_x, const Field3D& delta_z,
                    const std::string& region = "RGN_NOBNDRY") override;
+  void calcWeights(const Field3D& delta_x, const Field3D& delta_z, const BoutMask& mask,
+                   const std::string& region = "RGN_NOBNDRY");
 
   // Use precalculated weights
   Field3D interpolate(const Field3D& f,
                       const std::string& region = "RGN_NOBNDRY") const override;
   // Calculate weights and interpolate
   Field3D interpolate(const Field3D& f, const Field3D& delta_x, const Field3D& delta_z,
-                      const std::string& region = "RGN_NOBNDRY") override;
+                      const std::string& region = "RGN_NOBNDRY");
   Field3D interpolate(const Field3D& f, const Field3D& delta_x, const Field3D& delta_z,
-                      const BoutMask& mask,
-                      const std::string& region = "RGN_NOBNDRY") override;
+                      const BoutMask& mask, const std::string& region = "RGN_NOBNDRY");
+
+  std::vector<ParallelTransform::PositionsAndWeights>
+  getWeightsForYApproximation(int i, int j, int k, int yoffset) override;
 };
+#endif // BOUT_HAS_PETSC
 
 class XZInterpolationFactory
     : public Factory<XZInterpolation, XZInterpolationFactory, Mesh*> {
@@ -279,11 +331,30 @@ public:
   ReturnType create(const std::string& type, [[maybe_unused]] Options* options) const {
     return Factory::create(type, nullptr);
   }
-
-  static void ensureRegistered();
 };
 
 template <class DerivedType>
 using RegisterXZInterpolation = XZInterpolationFactory::RegisterInFactory<DerivedType>;
 
+using RegisterUnavailableXZInterpolation =
+    XZInterpolationFactory::RegisterUnavailableInFactory;
+
+namespace {
+const inline RegisterXZInterpolation<XZHermiteSpline> registerinterphermitespline{
+    "hermitespline"};
+const inline RegisterXZInterpolation<XZMonotonicHermiteSpline>
+    registerinterpmonotonichermitespline{"monotonichermitespline"};
+const inline RegisterXZInterpolation<XZLagrange4pt> registerinterplagrange4pt{
+    "lagrange4pt"};
+const inline RegisterXZInterpolation<XZBilinear> registerinterpbilinear{"bilinear"};
+
+#if BOUT_HAS_PETSC
+const inline RegisterXZInterpolation<PetscXZHermiteSpline> registerpetschermitespline{
+    "petschermitespline"};
+#else
+const inline RegisterUnavailableXZInterpolation register_unavailable_petsc_hermite_spline{
+    "petschermitespline", "BOUT++ was not configured with PETSc"};
+#endif
+} // namespace
+
 #endif // BOUT_INTERP_XZ_H

src/mesh/interpolation/bilinear_xz.cxx

@@ -83,12 +83,6 @@ void XZBilinear::calcWeights(const Field3D& delta_x, const Field3D& delta_z,
   }
 }
 
-void XZBilinear::calcWeights(const Field3D& delta_x, const Field3D& delta_z,
-                             const BoutMask& mask, const std::string& region) {
-  setMask(mask);
-  calcWeights(delta_x, delta_z, region);
-}
-
 Field3D XZBilinear::interpolate(const Field3D& f, const std::string& region) const {
   ASSERT1(f.getMesh() == localmesh);
   Field3D f_interp{emptyFrom(f)};
@@ -113,16 +107,3 @@ Field3D XZBilinear::interpolate(const Field3D& f, const std::string& region) con
   }
   return f_interp;
 }
-
-Field3D XZBilinear::interpolate(const Field3D& f, const Field3D& delta_x,
-                                const Field3D& delta_z, const std::string& region) {
-  calcWeights(delta_x, delta_z, region);
-  return interpolate(f, region);
-}
-
-Field3D XZBilinear::interpolate(const Field3D& f, const Field3D& delta_x,
-                                const Field3D& delta_z, const BoutMask& mask,
-                                const std::string& region) {
-  calcWeights(delta_x, delta_z, mask, region);
-  return interpolate(f, region);
-}

src/mesh/interpolation/hermite_spline_xz.cxx

@@ -111,12 +111,6 @@ void XZHermiteSpline::calcWeights(const Field3D& delta_x, const Field3D& delta_z
   }
 }
 
-void XZHermiteSpline::calcWeights(const Field3D& delta_x, const Field3D& delta_z,
-                                  const BoutMask& mask, const std::string& region) {
-  setMask(mask);
-  calcWeights(delta_x, delta_z, region);
-}
-
 /*!
  * Return position and weight of points needed to approximate the function value at the
  * point that the field line through (i,j,k) meets the (j+1)-plane. For the case where
@@ -223,16 +217,3 @@ Field3D XZHermiteSpline::interpolate(const Field3D& f, const std::string& region
   }
   return f_interp;
 }
-
-Field3D XZHermiteSpline::interpolate(const Field3D& f, const Field3D& delta_x,
-                                     const Field3D& delta_z, const std::string& region) {
-  calcWeights(delta_x, delta_z, region);
-  return interpolate(f, region);
-}
-
-Field3D XZHermiteSpline::interpolate(const Field3D& f, const Field3D& delta_x,
-                                     const Field3D& delta_z, const BoutMask& mask,
-                                     const std::string& region) {
-  calcWeights(delta_x, delta_z, mask, region);
-  return interpolate(f, region);
-}

src/mesh/interpolation/lagrange_4pt_xz.cxx

@@ -20,11 +20,29 @@
  *
  **************************************************************************/
 
+#include "bout/bout_types.hxx"
 #include "bout/globals.hxx"
 #include "bout/interpolation_xz.hxx"
 #include "bout/mesh.hxx"
 
-#include <vector>
+#include <array>
+
+namespace {
+// 4-point Lagrangian interpolation
+// offset must be between 0 and 1
+BoutReal lagrange_4pt(const BoutReal v2m, const BoutReal v1m, const BoutReal v1p,
+                      const BoutReal v2p, const BoutReal offset) {
+  return -offset * (offset - 1.0) * (offset - 2.0) * v2m / 6.0
+         + 0.5 * (offset * offset - 1.0) * (offset - 2.0) * v1m
+         - 0.5 * offset * (offset + 1.0) * (offset - 2.0) * v1p
+         + offset * (offset * offset - 1.0) * v2p / 6.0;
+}
+// Convenience helper
+BoutReal lagrange_4pt(const std::array<BoutReal, 4>& x, const BoutReal offset) {
+  return lagrange_4pt(x[0], x[1], x[2], x[3], offset);
+}
+
+} // namespace
 
 XZLagrange4pt::XZLagrange4pt(int y_offset, Mesh* mesh)
     : XZInterpolation(y_offset, mesh), t_x(localmesh), t_z(localmesh) {
@@ -75,12 +93,6 @@ void XZLagrange4pt::calcWeights(const Field3D& delta_x, const Field3D& delta_z,
   }
 }
 
-void XZLagrange4pt::calcWeights(const Field3D& delta_x, const Field3D& delta_z,
-                                const BoutMask& mask, const std::string& region) {
-  setMask(mask);
-  calcWeights(delta_x, delta_z, region);
-}
-
 Field3D XZLagrange4pt::interpolate(const Field3D& f, const std::string& region) const {
 
   ASSERT1(f.getMesh() == localmesh);
@@ -107,55 +119,25 @@ Field3D XZLagrange4pt::interpolate(const Field3D& f, const std::string& region)
     const int jz2mnew = (jz - 1 + ncz) % ncz;
 
     // Interpolate in Z first
-    BoutReal xvals[4];
-
     const int y_next = y + y_offset;
 
-    xvals[0] = lagrange_4pt(f(jx2mnew, y_next, jz2mnew), f(jx2mnew, y_next, jz),
-                            f(jx2mnew, y_next, jzpnew), f(jx2mnew, y_next, jz2pnew),
-                            t_z(x, y, z));
+    const std::array<BoutReal, 4> xvals{
+        lagrange_4pt(f(jx2mnew, y_next, jz2mnew), f(jx2mnew, y_next, jz),
+                     f(jx2mnew, y_next, jzpnew), f(jx2mnew, y_next, jz2pnew),
+                     t_z(x, y, z)),
 
-    xvals[1] = lagrange_4pt(f(jx, y_next, jz2mnew), f(jx, y_next, jz),
-                            f(jx, y_next, jzpnew), f(jx, y_next, jz2pnew), t_z(x, y, z));
+        lagrange_4pt(f(jx, y_next, jz2mnew), f(jx, y_next, jz), f(jx, y_next, jzpnew),
+                     f(jx, y_next, jz2pnew), t_z(x, y, z)),
 
-    xvals[2] =
         lagrange_4pt(f(jxpnew, y_next, jz2mnew), f(jxpnew, y_next, jz),
-                     f(jxpnew, y_next, jzpnew), f(jxpnew, y_next, jz2pnew), t_z(x, y, z));
-
-    xvals[3] = lagrange_4pt(f(jx2pnew, y_next, jz2mnew), f(jx2pnew, y_next, jz),
-                            f(jx2pnew, y_next, jzpnew), f(jx2pnew, y_next, jz2pnew),
-                            t_z(x, y, z));
+                     f(jxpnew, y_next, jzpnew), f(jxpnew, y_next, jz2pnew), t_z(x, y, z)),
 
+        lagrange_4pt(f(jx2pnew, y_next, jz2mnew), f(jx2pnew, y_next, jz),
+                     f(jx2pnew, y_next, jzpnew), f(jx2pnew, y_next, jz2pnew),
+                     t_z(x, y, z)),
+    };
     // Then in X
     f_interp(x, y_next, z) = lagrange_4pt(xvals, t_x(x, y, z));
   }
   return f_interp;
 }
-
-Field3D XZLagrange4pt::interpolate(const Field3D& f, const Field3D& delta_x,
-                                   const Field3D& delta_z, const std::string& region) {
-  calcWeights(delta_x, delta_z, region);
-  return interpolate(f, region);
-}
-
-Field3D XZLagrange4pt::interpolate(const Field3D& f, const Field3D& delta_x,
-                                   const Field3D& delta_z, const BoutMask& mask,
-                                   const std::string& region) {
-  calcWeights(delta_x, delta_z, mask, region);
-  return interpolate(f, region);
-}
-
-// 4-point Lagrangian interpolation
-// offset must be between 0 and 1
-BoutReal XZLagrange4pt::lagrange_4pt(const BoutReal v2m, const BoutReal vm,
-                                     const BoutReal vp, const BoutReal v2p,
-                                     const BoutReal offset) const {
-  return -offset * (offset - 1.0) * (offset - 2.0) * v2m / 6.0
-         + 0.5 * (offset * offset - 1.0) * (offset - 2.0) * vm
-         - 0.5 * offset * (offset + 1.0) * (offset - 2.0) * vp
-         + offset * (offset * offset - 1.0) * v2p / 6.0;
-}
-
-BoutReal XZLagrange4pt::lagrange_4pt(const BoutReal v[], const BoutReal offset) const {
-  return lagrange_4pt(v[0], v[1], v[2], v[3], offset);
-}