Improvements for hermitespline interpolation

.git-blame-ignore-revs

@@ -11,3 +11,5 @@ a71cad2dd6ace5741a754e2ca7daacd4bb094e0e
 2c2402ed59c91164eaff46dee0f79386b7347e9e
 05b7c571544c3bcb153fce67d12b9ac48947fc2d
 c8f38049359170a34c915e209276238ea66b9a1e
+65880e4af16cb95438437bf057c4b9fdb427b142
+d1cfb8abd6aa5c76e6c1a4d7ab20929c65f8afc2

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,6 +27,7 @@
 #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
@@ -38,6 +39,7 @@
 #endif
 
 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 +137,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>
+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
@@ -167,15 +181,20 @@ protected:
   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() {
+  ~XZHermiteSplineBase() override {
 #if HS_USE_PETSC
     if (isInit) {
       MatDestroy(&petscWeights);
@@ -205,61 +224,8 @@ 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 XZMonotonicHermiteSpline = XZHermiteSplineBase<true>;
+using XZHermiteSpline = XZHermiteSplineBase<false>;
 
 /// XZLagrange4pt interpolation class
 ///

include/bout/msg_stack.hxx

@@ -80,7 +80,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) {