Add Mesh::communicate_no_slices

include/bout/mesh.hxx

@@ -262,6 +262,16 @@ public:
     communicate(g);
   }
 
+  /// Communicate Fields without calculating the parallel slices
+  template <typename... Ts>
+  void communicate_no_slices(Ts&... ts) {
+    FieldGroup g(ts...);
+    const bool old = calcParallelSlices_on_communicate;
+    calcParallelSlices_on_communicate = false;
+    communicate(g);
+    calcParallelSlices_on_communicate = old;
+  }
+
   template <typename... Ts>
   void communicateXZ(Ts&... ts) {
     FieldGroup g(ts...);
@@ -795,9 +805,11 @@ protected:
   /// Mesh options section
   Options* options{nullptr};
 
+private:
   /// Set whether to call calcParallelSlices on all communicated fields (true) or not (false)
   bool calcParallelSlices_on_communicate{true};
 
+protected:
   /// Read a 1D array of integers
   const std::vector<int> readInts(const std::string& name, int n);
 

src/mesh/coordinates.cxx

@@ -23,19 +23,7 @@
 #include "parallel/fci.hxx"
 #include "parallel/shiftedmetricinterp.hxx"
 
-// use anonymous namespace so this utility function is not available outside this file
 namespace {
-template <typename T, typename... Ts>
-// Use sendY()/sendX() and wait() instead of Mesh::communicate() to ensure we
-// don't try to calculate parallel slices as Coordinates are not constructed yet
-void communicate(T& t, Ts&... ts) {
-  FieldGroup g(t, ts...);
-  auto h = t.getMesh()->sendY(g);
-  t.getMesh()->wait(h);
-  h = t.getMesh()->sendX(g);
-  t.getMesh()->wait(h);
-}
-
 /// Interpolate a Field2D to a new CELL_LOC with interp_to.
 /// Communicates to set internal guard cells.
 /// Boundary guard cells are set by extrapolating from the grid, like
@@ -55,7 +43,7 @@ Field2D interpolateAndExtrapolate(const Field2D& f, CELL_LOC location, bool extr
   // communicate f. We will sort out result's boundary guard cells below, but
   // not f's so we don't want to change f.
   result.allocate();
-  communicate(result);
+  localmesh->communicate_no_slices(result);
 
   // Extrapolate into boundaries (if requested) so that differential geometry
   // terms can be interpolated if necessary
@@ -204,7 +192,7 @@ Field3D interpolateAndExtrapolate(const Field3D& f_, CELL_LOC location,
   // communicate f. We will sort out result's boundary guard cells below, but
   // not f's so we don't want to change f.
   result.allocate();
-  communicate(result);
+  localmesh->communicate_no_slices(result);
 
   // Extrapolate into boundaries (if requested) so that differential geometry
   // terms can be interpolated if necessary
@@ -444,7 +432,7 @@ Coordinates::Coordinates(Mesh* mesh, Options* options)
                                  transform.get());
 
   if (mesh->periodicX) {
-    communicate(dx);
+    mesh->communicate_no_slices(dx);
   }
 
   dy = interpolateAndExtrapolate(dy, location, extrapolate_x, extrapolate_y, false,
@@ -553,7 +541,7 @@ Coordinates::Coordinates(Mesh* mesh, Options* options)
     // Compare calculated and loaded values
     output_warn.write("\tMaximum difference in J is {:e}\n", max(abs(J - Jcalc)));
 
-    communicate(J);
+    mesh->communicate_no_slices(J);
 
     // Re-evaluate Bxy using new J
     Bxy = sqrt(g_22) / J;
@@ -660,7 +648,7 @@ Coordinates::Coordinates(Mesh* mesh, Options* options, const CELL_LOC loc,
                           extrapolate_y, false, transform.get());
 
     if (mesh->periodicX) {
-      communicate(dx);
+      mesh->communicate_no_slices(dx);
     }
 
     getAtLocAndFillGuards(mesh, dy, "dy", suffix, location, 1.0, extrapolate_x,
@@ -944,8 +932,8 @@ const Field2D& Coordinates::zlength() const {
 int Coordinates::geometry(bool recalculate_staggered,
                           bool force_interpolate_from_centre) {
   TRACE("Coordinates::geometry");
-  communicate(dx, dy, dz, g11, g22, g33, g12, g13, g23, g_11, g_22, g_33, g_12, g_13,
-              g_23, J, Bxy);
+  localmesh->communicate_no_slices(dx, dy, dz, g11, g22, g33, g12, g13, g23, g_11, g_22,
+                                   g_33, g_12, g_13, g_23, J, Bxy);
 
   output_progress.write("Calculating differential geometry terms\n");
 
@@ -1028,20 +1016,21 @@ int Coordinates::geometry(bool recalculate_staggered,
           + 0.5 * g33 * DDY(g_33);
 
   auto tmp = J * g12;
-  communicate(tmp);
+  localmesh->communicate_no_slices(tmp);
   G1 = (DDX(J * g11) + DDY(tmp) + DDZ(J * g13)) / J;
   tmp = J * g22;
-  communicate(tmp);
+  localmesh->communicate_no_slices(tmp);
   G2 = (DDX(J * g12) + DDY(tmp) + DDZ(J * g23)) / J;
   tmp = J * g23;
-  communicate(tmp);
+  localmesh->communicate_no_slices(tmp);
   G3 = (DDX(J * g13) + DDY(tmp) + DDZ(J * g33)) / J;
 
   // Communicate christoffel symbol terms
   output_progress.write("\tCommunicating connection terms\n");
 
-  communicate(G1_11, G1_22, G1_33, G1_12, G1_13, G1_23, G2_11, G2_22, G2_33, G2_12, G2_13,
-              G2_23, G3_11, G3_22, G3_33, G3_12, G3_13, G3_23, G1, G2, G3);
+  localmesh->communicate_no_slices(G1_11, G1_22, G1_33, G1_12, G1_13, G1_23, G2_11, G2_22,
+                                   G2_33, G2_12, G2_13, G2_23, G3_11, G3_22, G3_33, G3_12,
+                                   G3_13, G3_23, G1, G2, G3);
 
   // Set boundary guard cells of Christoffel symbol terms
   // Ideally, when location is staggered, we would set the upper/outer boundary point
@@ -1098,7 +1087,7 @@ int Coordinates::geometry(bool recalculate_staggered,
           "\tWARNING: differencing quantity 'd2x' not found. Calculating from dx\n");
       d1_dx = bout::derivatives::index::DDX(1. / dx); // d/di(1/dx)
 
-      communicate(d1_dx);
+      localmesh->communicate_no_slices(d1_dx);
       d1_dx =
           interpolateAndExtrapolate(d1_dx, location, true, true, true, transform.get());
     } else {
@@ -1115,7 +1104,7 @@ int Coordinates::geometry(bool recalculate_staggered,
           "\tWARNING: differencing quantity 'd2y' not found. Calculating from dy\n");
       d1_dy = DDY(1. / dy); // d/di(1/dy)
 
-      communicate(d1_dy);
+      localmesh->communicate_no_slices(d1_dy);
       d1_dy =
           interpolateAndExtrapolate(d1_dy, location, true, true, true, transform.get());
     } else {
@@ -1132,7 +1121,7 @@ int Coordinates::geometry(bool recalculate_staggered,
       output_warn.write(
           "\tWARNING: differencing quantity 'd2z' not found. Calculating from dz\n");
       d1_dz = bout::derivatives::index::DDZ(1. / dz);
-      communicate(d1_dz);
+      localmesh->communicate_no_slices(d1_dz);
       d1_dz =
           interpolateAndExtrapolate(d1_dz, location, true, true, true, transform.get());
     } else {
@@ -1152,7 +1141,7 @@ int Coordinates::geometry(bool recalculate_staggered,
           "\tWARNING: differencing quantity 'd2x' not found. Calculating from dx\n");
       d1_dx = bout::derivatives::index::DDX(1. / dx); // d/di(1/dx)
 
-      communicate(d1_dx);
+      localmesh->communicate_no_slices(d1_dx);
       d1_dx =
           interpolateAndExtrapolate(d1_dx, location, true, true, true, transform.get());
     } else {
@@ -1167,7 +1156,7 @@ int Coordinates::geometry(bool recalculate_staggered,
           "\tWARNING: differencing quantity 'd2y' not found. Calculating from dy\n");
       d1_dy = DDY(1. / dy); // d/di(1/dy)
 
-      communicate(d1_dy);
+      localmesh->communicate_no_slices(d1_dy);
       d1_dy =
           interpolateAndExtrapolate(d1_dy, location, true, true, true, transform.get());
     } else {
@@ -1183,7 +1172,7 @@ int Coordinates::geometry(bool recalculate_staggered,
           "\tWARNING: differencing quantity 'd2z' not found. Calculating from dz\n");
       d1_dz = bout::derivatives::index::DDZ(1. / dz);
 
-      communicate(d1_dz);
+      localmesh->communicate_no_slices(d1_dz);
       d1_dz =
           interpolateAndExtrapolate(d1_dz, location, true, true, true, transform.get());
     } else {
@@ -1196,7 +1185,7 @@ int Coordinates::geometry(bool recalculate_staggered,
     d1_dz = 0;
 #endif
   }
-  communicate(d1_dx, d1_dy, d1_dz);
+  localmesh->communicate_no_slices(d1_dx, d1_dy, d1_dz);
 
   if (location == CELL_CENTRE && recalculate_staggered) {
     // Re-calculate interpolated Coordinates at staggered locations
@@ -1366,7 +1355,7 @@ void fixZShiftGuards(Field2D& zShift) {
                                      not localmesh->sourceHasYBoundaryGuards(), false);
 
   // make sure zShift has been communicated
-  communicate(zShift);
+  localmesh->communicate_no_slices(zShift);
 
   // Correct guard cells for discontinuity of zShift at poloidal branch cut
   for (int x = 0; x < localmesh->LocalNx; x++) {

src/mesh/interpolation/monotonic_hermite_spline_xz.cxx

@@ -36,26 +36,10 @@ Field3D XZMonotonicHermiteSpline::interpolate(const Field3D& f,
   // Derivatives are used for tension and need to be on dimensionless
   // coordinates
   Field3D fx = bout::derivatives::index::DDX(f, CELL_DEFAULT, "DEFAULT");
-  localmesh->communicateXZ(fx);
-  // communicate in y, but do not calculate parallel slices
-  {
-    auto h = localmesh->sendY(fx);
-    localmesh->wait(h);
-  }
   Field3D fz = bout::derivatives::index::DDZ(f, CELL_DEFAULT, "DEFAULT", "RGN_ALL");
-  localmesh->communicateXZ(fz);
-  // communicate in y, but do not calculate parallel slices
-  {
-    auto h = localmesh->sendY(fz);
-    localmesh->wait(h);
-  }
+  localmesh->communicate_no_slices(fx, fz);
   Field3D fxz = bout::derivatives::index::DDX(fz, CELL_DEFAULT, "DEFAULT");
-  localmesh->communicateXZ(fxz);
-  // communicate in y, but do not calculate parallel slices
-  {
-    auto h = localmesh->sendY(fxz);
-    localmesh->wait(h);
-  }
+  localmesh->communicate_no_slices(fxz);
 
   const auto curregion{getRegion(region)};
   BOUT_FOR(i, curregion) {

tests/unit/fake_parallel_mesh.hxx

@@ -53,7 +53,6 @@ public:
     StaggerGrids = false;
     periodicX = false;
     IncIntShear = false;
-    calcParallelSlices_on_communicate = true;
     options = Options::getRoot();
     mpi = mpiSmart.get();
   }