WIP: Add communication in Z

The actual significant changes are quite small, but these have big knock-on
effects in the tests. Main changes:

- Z is still periodic, but this is now implemented using guard cells instead of
  modulus in `SpecificInd::zp/zm`
- Adds a new `Mesh` method `sendZ`
- `NZPE` is now a user-settable option
- Domain decomposition is done in Z before the usual X-Y decomposition
- The X/Y _boundary_ regions do not include the Z guards
  - These corners should be filled in by `Mesh::communicate`

Changes to the tests:
- `FakeMesh` now has 1 Z guard cell, which changes the size of the interior
  domain, and requires many of the `Field*` tests to be updated
- `IndexOffsetTest` needs to skip the Z guards, just like X, Y
- Derivative tests need extra point for interior + 2 * guards
- The Petsc and Hypre Laplace unit tests require the Z guards to be communicated
  in order to fill in the global indices. This requires `FakeMesh` to actually
  implement `sendZ`. This just copies the relevant points into the Z guards
  - Also, we should not check the result in the Z guards
  - I've also pulled out a bunch of common machinery between these two tests
- `BoutMeshTest` now needs to understand Z guards and Z processors

Author: ZedThree

include/bout/mesh.hxx

@@ -349,6 +349,9 @@ public:
   /// Send only the y-guard cells
   virtual comm_handle sendY(FieldGroup& g, comm_handle handle = nullptr) = 0;
 
+  /// Send only the z-guard cells
+  virtual comm_handle sendZ(FieldGroup& g, comm_handle handle = nullptr) = 0;
+
   /// Wait for the handle, return error code
   virtual int wait(comm_handle handle) = 0; ///< Wait for the handle, return error code
 
@@ -812,9 +815,9 @@ protected:
   /// Read a 1D array of integers
   const std::vector<int> readInts(const std::string& name, int n);
 
-  /// Calculates the size of a message for a given x and y range
-  int msg_len(const std::vector<FieldData*>& var_list, int xge, int xlt, int yge,
-              int ylt);
+  /// Calculates the size of a message for a given (x, y, z) range
+  int msg_len(const std::vector<FieldData*>& var_list, int xge, int xlt, int yge, int ylt,
+              int zge, int zlt);
 
   /// Initialise derivatives
   void derivs_init(Options* options);

include/bout/region.hxx

@@ -170,8 +170,8 @@ struct SpecificInd {
   int ny = -1, nz = -1; ///< Sizes of y and z dimensions
 
   SpecificInd() = default;
-  SpecificInd(int i, int ny, int nz) : ind(i), ny(ny), nz(nz){};
-  explicit SpecificInd(int i) : ind(i){};
+  SpecificInd(int i, int ny, int nz) : ind(i), ny(ny), nz(nz) {};
+  explicit SpecificInd(int i) : ind(i) {};
 
   /// Allow explicit conversion to an int
   explicit operator int() const { return ind; }
@@ -290,26 +290,15 @@ struct SpecificInd {
   }
   /// The index one point -1 in y
   inline SpecificInd ym(int dy = 1) const { return yp(-dy); }
-  /// The index one point +1 in z. Wraps around zend to zstart
-  /// An alternative, non-branching calculation is :
-  /// ind + dz - nz * ((ind + dz) / nz  - ind / nz)
-  /// but this appears no faster (and perhaps slower).
+  /// The index one point +1 in z
   inline SpecificInd zp(int dz = 1) const {
-    ASSERT3(dz >= 0);
-    dz = dz <= nz ? dz : dz % nz; //Fix in case dz > nz, if not force it to be in range
-    return {(ind + dz) % nz < dz ? ind - nz + dz : ind + dz, ny, nz};
-  }
-  /// The index one point -1 in z. Wraps around zstart to zend
-  /// An alternative, non-branching calculation is :
-  /// ind - dz + nz * ( (nz + ind) / nz - (nz + ind - dz) / nz)
-  /// but this appears no faster (and perhaps slower).
-  inline SpecificInd zm(int dz = 1) const {
-    dz = dz <= nz ? dz : dz % nz; //Fix in case dz > nz, if not force it to be in range
-    ASSERT3(dz >= 0);
-    return {(ind) % nz < dz ? ind + nz - dz : ind - dz, ny, nz};
+    if constexpr (N == IND_TYPE::IND_2D) {
+      return *this;
+    }
+    return {ind + dz, ny, nz};
   }
-  /// Automatically select zm or zp depending on sign
-  inline SpecificInd zpm(int dz) const { return dz > 0 ? zp(dz) : zm(-dz); }
+  /// The index one point -1 in z
+  inline SpecificInd zm(int dz = 1) const { return zp(-dz); }
 
   // and for 2 cells
   inline SpecificInd xpp() const { return xp(2); }
@@ -320,9 +309,7 @@ struct SpecificInd {
   inline SpecificInd zmm() const { return zm(2); }
 
   /// Generic offset of \p index in multiple directions simultaneously
-  inline SpecificInd offset(int dx, int dy, int dz) const {
-    return zpm(dz).yp(dy).xp(dx);
-  }
+  inline SpecificInd offset(int dx, int dy, int dz) const { return zp(dz).yp(dy).xp(dx); }
 };
 
 /// Relational operators
@@ -490,10 +477,9 @@ template <typename T = Ind3D>
 class Region {
   // Following prevents a Region being created with anything other
   // than Ind2D, Ind3D or IndPerp as template type
-  static_assert(
-      std::is_base_of_v<
-          Ind2D, T> || std::is_base_of_v<Ind3D, T> || std::is_base_of_v<IndPerp, T>,
-      "Region must be templated with one of IndPerp, Ind2D or Ind3D");
+  static_assert(std::is_base_of_v<Ind2D, T> || std::is_base_of_v<Ind3D, T>
+                    || std::is_base_of_v<IndPerp, T>,
+                "Region must be templated with one of IndPerp, Ind2D or Ind3D");
 
 public:
   using data_type = T;
@@ -569,7 +555,7 @@ public:
   };
 
   Region(RegionIndices& indices, int maxregionblocksize = MAXREGIONBLOCKSIZE)
-      : indices(indices), blocks(getContiguousBlocks(maxregionblocksize)){};
+      : indices(indices), blocks(getContiguousBlocks(maxregionblocksize)) {};
 
   // We need to first set the blocks, and only after that call getRegionIndices.
   // Do not put in the member initialisation

src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx

@@ -436,8 +436,7 @@ OperatorStencil<Ind3D> LaplaceHypre3d::getStencil(Mesh* localmesh,
   std::transform(pw.begin(), pw.end(), std::back_inserter(interpPattern),
                  [localmesh](ParallelTransform::PositionsAndWeights p) -> OffsetInd3D {
                    return {localmesh->xstart - p.i, localmesh->ystart - p.j,
-                           localmesh->LocalNz - p.k < p.k ? p.k - localmesh->LocalNz
-                                                          : p.k};
+                           localmesh->zstart - p.k};
                  });
 
   OffsetInd3D zero;

src/invert/laplace/impls/petsc3damg/petsc3damg.cxx

@@ -490,9 +490,7 @@ OperatorStencil<Ind3D> LaplacePetsc3dAmg::getStencil(Mesh* localmesh,
       std::back_inserter(interpPattern),
       [localmesh](ParallelTransform::PositionsAndWeights position) -> OffsetInd3D {
         return {localmesh->xstart - position.i, localmesh->ystart - position.j,
-                ((localmesh->LocalNz - position.k) < position.k)
-                    ? position.k - localmesh->LocalNz
-                    : position.k};
+                localmesh->zstart - position.k};
       });
 
   const OffsetInd3D zero;