Merge pull request #3206 from boutproject/add-assert-non-parallel

Add some asserts for non parallelised XZ interpolation

Author: ZedThree

include/bout/interpolation_xz.hxx

@@ -238,6 +238,9 @@ public:
                       const std::string& region = "RGN_NOBNDRY") const override;
 };
 
+/// XZLagrange4pt interpolation class
+///
+/// Does not support MPI splitting in X
 class XZLagrange4pt : public XZInterpolation {
   Tensor<int> i_corner; // x-index of bottom-left grid point
   Tensor<int> k_corner; // z-index of bottom-left grid point
@@ -271,6 +274,9 @@ public:
   BoutReal lagrange_4pt(const BoutReal v[], BoutReal offset) const;
 };
 
+/// XZBilinear interpolation calss
+///
+/// Does not support MPI splitting in X.
 class XZBilinear : public XZInterpolation {
   Tensor<int> i_corner; // x-index of bottom-left grid point
   Tensor<int> k_corner; // z-index of bottom-left grid point

src/mesh/interpolation/bilinear_xz.cxx

@@ -31,6 +31,10 @@ XZBilinear::XZBilinear(int y_offset, Mesh* mesh)
     : XZInterpolation(y_offset, mesh), w0(localmesh), w1(localmesh), w2(localmesh),
       w3(localmesh) {
 
+  if (localmesh->getNXPE() > 1) {
+    throw BoutException("XZBilinear interpolation does not support MPI splitting in X");
+  }
+
   // Index arrays contain guard cells in order to get subscripts right
   i_corner.reallocate(localmesh->LocalNx, localmesh->LocalNy, localmesh->LocalNz);
   k_corner.reallocate(localmesh->LocalNx, localmesh->LocalNy, localmesh->LocalNz);

src/mesh/interpolation/hermite_spline_xz.cxx

@@ -101,8 +101,8 @@ class IndConverter {
   }
 };
 
-XZHermiteSpline::XZHermiteSpline(int y_offset, Mesh* mesh)
-    : XZInterpolation(y_offset, mesh), h00_x(localmesh), h01_x(localmesh),
+XZHermiteSpline::XZHermiteSpline(int y_offset, Mesh* meshin)
+    : XZInterpolation(y_offset, meshin), h00_x(localmesh), h01_x(localmesh),
       h10_x(localmesh), h11_x(localmesh), h00_z(localmesh), h01_z(localmesh),
       h10_z(localmesh), h11_z(localmesh) {
 
@@ -346,6 +346,9 @@ Field3D XZHermiteSpline::interpolate(const Field3D& f, const std::string& region
   ASSERT1(f.getMesh() == localmesh);
   Field3D f_interp{emptyFrom(f)};
 
+  const auto region2 =
+      y_offset == 0 ? "RGN_NOY" : fmt::format("RGN_YPAR_{:+d}", y_offset);
+
 #if USE_NEW_WEIGHTS
 #ifdef HS_USE_PETSC
   BoutReal* ptr;
@@ -355,7 +358,6 @@ Field3D XZHermiteSpline::interpolate(const Field3D& f, const std::string& region
   VecRestoreArray(rhs, &ptr);
   MatMult(petscWeights, rhs, result);
   VecGetArrayRead(result, &cptr);
-  const auto region2 = y_offset == 0 ? region : fmt::format("RGN_YPAR_{:+d}", y_offset);
   BOUT_FOR(i, f.getRegion(region2)) {
     f_interp[i] = cptr[int(i)];
     ASSERT2(std::isfinite(cptr[int(i)]));
@@ -375,11 +377,10 @@ Field3D XZHermiteSpline::interpolate(const Field3D& f, const std::string& region
     }
   }
 #endif
-  return f_interp;
 #else
   // Derivatives are used for tension and need to be on dimensionless
   // coordinates
-  const auto region2 = fmt::format("RGN_YPAR_{:+d}", y_offset);
+
   // f has been communcated, and thus we can assume that the x-boundaries are
   // also valid in the y-boundary.  Thus the differentiated field needs no
   // extra comms.
@@ -418,8 +419,10 @@ Field3D XZHermiteSpline::interpolate(const Field3D& f, const std::string& region
     ASSERT2(std::isfinite(f_interp[iyp]) || i.x() < localmesh->xstart
             || i.x() > localmesh->xend);
   }
-  return f_interp;
 #endif
+  f_interp.setRegion(region2);
+  ASSERT2(f_interp.getRegionID());
+  return f_interp;
 }
 
 Field3D XZHermiteSpline::interpolate(const Field3D& f, const Field3D& delta_x,

src/mesh/interpolation/lagrange_4pt_xz.cxx

@@ -29,6 +29,11 @@
 XZLagrange4pt::XZLagrange4pt(int y_offset, Mesh* mesh)
     : XZInterpolation(y_offset, mesh), t_x(localmesh), t_z(localmesh) {
 
+  if (localmesh->getNXPE() > 1) {
+    throw BoutException(
+        "XZLagrange4pt interpolation does not support MPI splitting in X");
+  }
+
   // Index arrays contain guard cells in order to get subscripts right
   i_corner.reallocate(localmesh->LocalNx, localmesh->LocalNy, localmesh->LocalNz);
   k_corner.reallocate(localmesh->LocalNx, localmesh->LocalNy, localmesh->LocalNz);