Merge pull request #3298 from boutproject/hermitespline

Improvements for hermitespline interpolation

Author: bendudson

.git-blame-ignore-revs

@@ -11,6 +11,8 @@ a71cad2dd6ace5741a754e2ca7daacd4bb094e0e
 2c2402ed59c91164eaff46dee0f79386b7347e9e
 05b7c571544c3bcb153fce67d12b9ac48947fc2d
 c8f38049359170a34c915e209276238ea66b9a1e
+65880e4af16cb95438437bf057c4b9fdb427b142
+d1cfb8abd6aa5c76e6c1a4d7ab20929c65f8afc2
 8d5cb39e03c2644715a50684f8cd0318b4e82676
 ec69e8838be2dde140a915e50506f8e1ce3cb534
 f2bc0488a298f136294c523bc5ab4086d090059b

CMakeLists.txt

@@ -296,11 +296,12 @@ set(BOUT_SOURCES
     ./src/mesh/interpolation/hermite_spline_z.cxx
     ./src/mesh/interpolation/interpolation_z.cxx
     ./src/mesh/interpolation/lagrange_4pt_xz.cxx
-    ./src/mesh/interpolation/monotonic_hermite_spline_xz.cxx
     ./src/mesh/invert3x3.hxx
     ./src/mesh/mesh.cxx
     ./src/mesh/parallel/fci.cxx
     ./src/mesh/parallel/fci.hxx
+    ./src/mesh/parallel/fci_comm.cxx
+    ./src/mesh/parallel/fci_comm.hxx
     ./src/mesh/parallel/identity.cxx
     ./src/mesh/parallel/shiftedmetric.cxx
     ./src/mesh/parallel/shiftedmetricinterp.cxx

include/bout/interpolation_xz.hxx

@@ -27,17 +27,18 @@
 #include <bout/bout_types.hxx>
 #include <bout/generic_factory.hxx>
 #include <bout/mask.hxx>
+#include <array>
 
-#define USE_NEW_WEIGHTS 1
 #if BOUT_HAS_PETSC
-#define HS_USE_PETSC 1
-#endif
-
-#ifdef HS_USE_PETSC
 #include "bout/petsclib.hxx"
 #endif
 
+namespace {
+enum class implementation_type { new_weights, petsc, legacy };
+}
+
 class Options;
+class GlobalField3DAccess;
 
 /// Interpolate a field onto a perturbed set of points
 const Field3D interpolate(const Field3D& f, const Field3D& delta_x,
@@ -135,14 +136,26 @@ public:
   }
 };
 
-class XZHermiteSpline : public XZInterpolation {
+/// Monotonic Hermite spline interpolator
+///
+/// Similar to XZHermiteSpline, so uses most of the same code.
+/// Forces the interpolated result to be in the range of the
+/// neighbouring cell values. This prevents unphysical overshoots,
+/// but also degrades accuracy near maxima and minima.
+/// Perhaps should only impose near boundaries, since that is where
+/// problems most obviously occur.
+template <bool monotonic, implementation_type imp_type>
+class XZHermiteSplineBase : public XZInterpolation {
 protected:
   /// This is protected rather than private so that it can be
   /// extended and used by HermiteSplineMonotonic
 
   Tensor<SpecificInd<IND_TYPE::IND_3D>> i_corner; // index of bottom-left grid point
   Tensor<int> k_corner;                           // z-index of bottom-left grid point
 
+  std::unique_ptr<GlobalField3DAccess> gf3daccess;
+  Tensor<std::array<int, 4>> g3dinds;
+
   // 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
@@ -160,23 +173,28 @@ protected:
 
   std::vector<Field3D> newWeights;
 
-#if HS_USE_PETSC
+#if BOUT_HAS_PETSC
   PetscLib* petsclib;
   bool isInit{false};
-  Mat petscWeights;
-  Vec rhs, result;
+  Mat petscWeights{};
+  Vec rhs{}, result{};
 #endif
 
+  /// Factors to allow for some wiggleroom
+  BoutReal abs_fac_monotonic{1e-2};
+  BoutReal rel_fac_monotonic{1e-1};
+
 public:
-  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) {
+  XZHermiteSplineBase(Mesh* mesh = nullptr, [[maybe_unused]] Options* options = nullptr)
+      : XZHermiteSplineBase(0, mesh, options) {}
+  XZHermiteSplineBase(int y_offset = 0, Mesh* mesh = nullptr, Options* options = nullptr);
+  XZHermiteSplineBase(const BoutMask& mask, int y_offset = 0, Mesh* mesh = nullptr,
+                      Options* options = nullptr)
+      : XZHermiteSplineBase(y_offset, mesh, options) {
     setRegion(regionFromMask(mask, localmesh));
   }
-  ~XZHermiteSpline() {
-#if HS_USE_PETSC
+  ~XZHermiteSplineBase() override {
+#if BOUT_HAS_PETSC
     if (isInit) {
       MatDestroy(&petscWeights);
       VecDestroy(&rhs);
@@ -205,61 +223,21 @@ public:
   getWeightsForYApproximation(int i, int j, int k, int yoffset) override;
 };
 
-/// Monotonic Hermite spline interpolator
-///
-/// Similar to XZHermiteSpline, so uses most of the same code.
-/// Forces the interpolated result to be in the range of the
-/// neighbouring cell values. This prevents unphysical overshoots,
-/// 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, 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");
-    }
-  }
-  XZMonotonicHermiteSpline(int y_offset = 0, Mesh* mesh = nullptr)
-      : XZHermiteSpline(y_offset, mesh) {
-    if (localmesh->getNXPE() > 1) {
-      throw BoutException("Do not support MPI splitting in X");
-    }
-  }
-  XZMonotonicHermiteSpline(const BoutMask& mask, int y_offset = 0, Mesh* mesh = nullptr)
-      : XZHermiteSpline(mask, y_offset, mesh) {
-    if (localmesh->getNXPE() > 1) {
-      throw BoutException("Do not support MPI splitting in X");
-    }
-  }
-
-  using XZHermiteSpline::interpolate;
-  /// Interpolate using precalculated weights.
-  /// This function is called by the other interpolate functions
-  /// in the base class XZHermiteSpline.
-  Field3D interpolate(const Field3D& f,
-                      const std::string& region = "RGN_NOBNDRY") const override;
-};
+using XZMonotonicHermiteSplineSerial =
+    XZHermiteSplineBase<true, implementation_type::new_weights>;
+using XZHermiteSplineSerial =
+    XZHermiteSplineBase<false, implementation_type::new_weights>;
+using XZMonotonicHermiteSplineLegacy =
+    XZHermiteSplineBase<true, implementation_type::legacy>;
+using XZHermiteSplineLegacy = XZHermiteSplineBase<false, implementation_type::legacy>;
+#if BOUT_HAS_PETSC
+using XZMonotonicHermiteSpline = XZHermiteSplineBase<true, implementation_type::petsc>;
+using XZHermiteSpline = XZHermiteSplineBase<false, implementation_type::petsc>;
+#else
+using XZMonotonicHermiteSpline =
+    XZHermiteSplineBase<true, implementation_type::new_weights>;
+using XZHermiteSpline = XZHermiteSplineBase<false, implementation_type::new_weights>;
+#endif
 
 /// XZLagrange4pt interpolation class
 ///

include/bout/mesh.hxx

@@ -360,6 +360,8 @@ public:
   virtual int getXProcIndex() const = 0; ///< This processor's index in X direction
   virtual int getYProcIndex() const = 0; ///< This processor's index in Y direction
   virtual int getZProcIndex() const = 0; ///< This processor's index in Z direction
+  /// The index of a processor with given X, Y, and Z index
+  virtual int getProcIndex(int X, int Y, int Z) const = 0;
 
   // X communications
   virtual bool firstX()

include/bout/msg_stack.hxx

@@ -82,7 +82,7 @@ public:
   }
 
   void pop() {}
-  void pop(int [[maybe_unused]] id) {}
+  void pop([[maybe_unused]] int id) {}
   void clear() {}
 
   void dump() {}

include/bout/region.hxx

@@ -139,7 +139,7 @@ class BoutMask;
   BOUT_FOR_OMP(index, (region), for schedule(BOUT_OPENMP_SCHEDULE) nowait)
 // NOLINTEND(cppcoreguidelines-macro-usage,bugprone-macro-parentheses)
 
-enum class IND_TYPE { IND_3D = 0, IND_2D = 1, IND_PERP = 2 };
+enum class IND_TYPE { IND_3D = 0, IND_2D = 1, IND_PERP = 2, IND_GLOBAL_3D = 3 };
 
 /// Indices base class for Fields -- Regions are dereferenced into these
 ///
@@ -386,6 +386,7 @@ inline SpecificInd<N> operator-(SpecificInd<N> lhs, const SpecificInd<N>& rhs) {
 using Ind3D = SpecificInd<IND_TYPE::IND_3D>;
 using Ind2D = SpecificInd<IND_TYPE::IND_2D>;
 using IndPerp = SpecificInd<IND_TYPE::IND_PERP>;
+using IndG3D = SpecificInd<IND_TYPE::IND_GLOBAL_3D>;
 
 /// Get string representation of Ind3D
 inline std::string toString(const Ind3D& i) {

src/mesh/impls/bout/boutmesh.cxx

@@ -555,6 +555,8 @@ int BoutMesh::load() {
   PE_XIND = MYPE % NXPE;
   PE_ZIND = 0;
 
+  ASSERT2(MYPE == getProcIndex(PE_XIND, PE_YIND, PE_ZIND));
+
   // Set the other grid sizes from nx, ny, nz
   setDerivedGridSizes();
 
@@ -1036,6 +1038,8 @@ void BoutMesh::createXBoundaries() {
   }
 }
 
+int BoutMesh::getProcIndex(int X, int Y, int Z) const { return Y * NXPE + X; }
+
 void BoutMesh::createYBoundaries() {
   if (MYG <= 0) {
     return;

src/mesh/impls/bout/boutmesh.hxx

@@ -63,6 +63,7 @@ public:
   int getXProcIndex() const override; ///< This processor's index in X direction
   int getYProcIndex() const override; ///< This processor's index in Y direction
   int getZProcIndex() const override; ///< This processor's index in Z direction
+  int getProcIndex(int X, int Y, int Z) const override;
 
   /////////////////////////////////////////////
   // X communications