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/globalindexer.hxx

@@ -17,21 +17,25 @@ using IndexerPtr = std::shared_ptr<GlobalIndexer<T>>;
 
 using InterpolationWeights = std::vector<ParallelTransform::PositionsAndWeights>;
 
-/*!
- * An object which accepts index objects produced by iterating over
- * fields and returns a global index. This index can be used when
- * constructing PETSc arrays. Guard regions used for communication
- * between processes will have the indices of the part of the interior
- * region they are mirroring. Boundaries required by the stencil will
- * have unique indices. If no stencil is provided then only the
- * interior region will be assigned global indices. By default, the
- * indexer is fully initialised so that guard cells are communicated
- * (ensuring they hold the appropriate global indices). However, by
- * passing `autoInitialise = false` this behaviour can be turned off
- * and the user can then manually call the `initialise()` method
- * later. This can be useful for mocking/faking the class when
- * testing.
- */
+/// Convert local indices on current process to consistent, unique
+/// global flat indices across all processes.
+///
+/// Note that this can only convert indices both local to *and owned by* the
+/// current process!
+///
+/// An object which accepts index objects produced by iterating over
+/// fields and returns a global index. This index can be used when
+/// constructing PETSc arrays. Guard regions used for communication
+/// between processes will have the indices of the part of the interior
+/// region they are mirroring. Boundaries required by the stencil will
+/// have unique indices. If no stencil is provided then only the
+/// interior region will be assigned global indices. By default, the
+/// indexer is fully initialised so that guard cells are communicated
+/// (ensuring they hold the appropriate global indices). However, by
+/// passing `autoInitialise = false` this behaviour can be turned off
+/// and the user can then manually call the `initialise()` method
+/// later. This can be useful for mocking/faking the class when
+/// testing.
 template <class T>
 class GlobalIndexer {
 public:
@@ -98,7 +102,7 @@ public:
     }
   }
 
-  virtual ~GlobalIndexer() {}
+  virtual ~GlobalIndexer() = default;
 
   /// Call this immediately after construction when running unit tests.
   void initialiseTest() {}
@@ -123,7 +127,7 @@ public:
 
   /// Convert the local index object to a global index which can be
   /// used in PETSc vectors and matrices.
-  int getGlobal(const ind_type& ind) const {
+  virtual int getGlobal(const ind_type& ind) const {
     return static_cast<int>(std::round(indices[ind]));
   }
 
@@ -224,7 +228,8 @@ private:
   Array<int> int_indices;
   /// The first and last global index on this processor (inclusive in
   /// both cases)
-  int globalStart, globalEnd;
+  int globalStart = 0;
+  int globalEnd = 0;
 
   /// Stencil for which this indexer has been configured
   OperatorStencil<ind_type> stencils;
@@ -237,4 +242,49 @@ private:
   mutable std::vector<int> numDiagonal, numOffDiagonal;
 };
 
+namespace bout {
+/// A simplified local-to-global index converter: includes all boundaries, but
+/// not guard cells, and assumes a simple rectangular division of processes.
+///
+/// The conversion indices are not precalculated or cached, so this may be less
+/// efficient than `GlobalIndexer`. However, it can convert "local" indices
+/// outside of the range owned by the current process, for example a local index
+/// with `x = -1` which would be owned by the process to the left.
+class SimpleGlobalIndexer : public GlobalIndexer<Field3D> {
+  // Stencil to ensure x-boundaries (to a depth of 2) are included in the base
+  // indexer
+  //
+  // Actually, I'm not 100% sure what this does, but this appears to the minimum
+  // necessary to get things to work
+  static auto XZstencil(const Mesh& mesh) {
+    using ind_type = Field3D::ind_type;
+    const IndexOffset<ind_type> zero;
+    OperatorStencil<ind_type> stencil;
+
+    stencil.add(
+        [&mesh](ind_type ind) -> bool {
+          return (mesh.xstart <= ind.x() and ind.x() <= mesh.xend)
+                 and (mesh.ystart <= ind.y() and ind.y() <= mesh.yend)
+                 and (mesh.zstart <= ind.z() and ind.z() <= mesh.zend);
+        },
+        {zero.xm(2), zero.xp(2), zero.zm(2), zero.zp(2)});
+
+    return stencil;
+  }
+
+public:
+  SimpleGlobalIndexer() = default;
+  explicit SimpleGlobalIndexer(Mesh* localmesh, bool autoInitialise = true)
+      : GlobalIndexer<Field3D>(localmesh, XZstencil(*localmesh), autoInitialise) {}
+
+  int getGlobal(const ind_type& ind) const override {
+    const auto& mesh = *getMesh();
+    return ((mesh.getGlobalXIndex(ind.x()) * mesh.GlobalNy)
+            + mesh.getGlobalYIndex(ind.y()))
+               * mesh.GlobalNz
+           + mesh.getGlobalZIndex(ind.z());
+  }
+};
+} // namespace bout
+
 #endif // BOUT_GLOBALINDEXER_H

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,12 @@
 #ifndef BOUT_INTERP_XZ_H
 #define BOUT_INTERP_XZ_H
 
+#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 +99,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 +176,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 +224,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 +248,68 @@ 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;
+  IndexerPtr<Field3D> indexer;
+  PetscMatrix<Field3D> weights;
+
+  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 = default;
+
+  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 +325,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);
-}