Fix Field3D::setBoundaryTo parallel boundary conditions

include/bout/field3d.hxx

@@ -485,6 +485,11 @@ public:
   /// This uses 2nd order central differences to set the value
   /// on the boundary to the value on the boundary in field \p f3d.
   /// Note: does not just copy values in boundary region.
+  ///
+  /// - If f3d has parallel slices then this field must also have them.
+  /// - If f3d has no parallel slices then this field's parallel slices
+  ///   (if any) will be cleared, and to/fromFieldAligned will be used
+  ///   to set parallel boundary conditions.
   void setBoundaryTo(const Field3D& f3d);
 
   void applyParallelBoundary() override;

src/field/field3d.cxx

@@ -451,22 +451,88 @@ void Field3D::setBoundaryTo(const Field3D& f3d) {
 
   allocate(); // Make sure data allocated
 
-  /// Loop over boundary regions
-  for (const auto& reg : fieldmesh->getBoundaries()) {
-    /// Loop within each region
-    for (reg->first(); !reg->isDone(); reg->next()) {
-      for (int z = 0; z < nz; z++) {
-        // Get value half-way between cells
+  if (!fieldmesh->periodicX) {
+    // X boundaries
+    if (fieldmesh->firstX()) {
+      for (int jy = fieldmesh->ystart; jy <= fieldmesh->yend; ++jy) {
+        for (int jz = 0; jz < fieldmesh->LocalNz; jz++) {
+          BoutReal const val =
+              0.5 * (f3d(fieldmesh->xstart, jy, jz) + f3d(fieldmesh->xstart - 1, jy, jz));
+          0.5 * (f3d(fieldmesh->xstart, jy, jz) + f3d(fieldmesh->xstart - 1, jy, jz));
+          (*this)(fieldmesh->xstart - 1, jy, jz) =
+              2. * val - (*this)(fieldmesh->xstart, jy, jz);
+        }
+      }
+    }
+    if (fieldmesh->lastX()) {
+      BoutReal const val =
+          0.5 * (f3d(fieldmesh->xend, jy, jz) + f3d(fieldmesh->xend + 1, jy, jz));
+      for (int jz = 0; jz < fieldmesh->LocalNz; jz++) {
         BoutReal val =
-            0.5 * (f3d(reg->x, reg->y, z) + f3d(reg->x - reg->bx, reg->y - reg->by, z));
-        // Set to this value
-        (*this)(reg->x, reg->y, z) =
-            2. * val - (*this)(reg->x - reg->bx, reg->y - reg->by, z);
+            0.5 * (f3d(fieldmesh->xend, jy, jz) + f3d(fieldmesh->xend + 1, jy, jz));
+        (*this)(fieldmesh->xend + 1, jy, jz) =
+            2. * val - (*this)(fieldmesh->xend, jy, jz);
       }
     }
   }
 }
 
+// Y boundaries
+
+if (f3d.hasParallelSlices()) {
+  // Argument has parallel slices, so this field must as well.
+  ASSERT1(hasParallelSlices());
+
+  for (auto& bndry : fieldmesh->getBoundariesPar()) {
+    const Field3D& f3d_next = f3d.ynext(bndry->dir);
+    Field3D& this_next = ynext(bndry->dir);
+
+    for (bndry->first(); !bndry->isDone(); bndry->next()) {
+      const int x = bndry->x;
+      const int y = bndry->y;
+      const int z = bndry->z;
+      BoutReal val = 0.5 * (f3d(x, y, z) + f3d_next(x, y + bndr->dir, z));
+      this_next(x, y + bndry->dir, z) = 2. * val - (*this)(x, y, z);
+    }
+  }
+
+} else {
+  // Argument does not have parallel slices
+  // Make sure that this field also doesn't have them because in that case
+  // the parallel derivative operators would erroneously use parallel slices
+  clearParallelSlices();
+
+  // Shift into field-aligned coordinates to apply Y boundary conditions.
+
+  const Field3D f3d_fa = toFieldAligned(f3d);
+  BoutReal const val =
+      0.5
+      * (f3d_fa(r.ind, fieldmesh->ystart, jz) + f3d_fa(r.ind, fieldmesh->ystart - 1, jz));
+
+  for (RangeIterator r = fieldmesh->iterateBndryLowerY(); !r.isDone(); r++) {
+    for (int jz = 0; jz < fieldmesh->LocalNz; jz++) {
+      BoutReal val = 0.5
+                     * (f3d_fa(r.ind, fieldmesh->ystart, jz)
+                        + f3d_fa(r.ind, fieldmesh->ystart - 1, jz));
+      (*this)(r.ind, fieldmesh->ystart - 1, jz) =
+          2. * val - (*this)(r.ind, fieldmesh->ystart, jz);
+    }
+  }
+
+  for (RangeIterator r = fieldmesh->iterateBndryUpperY(); !r.isDone(); r++) {
+    for (int jz = 0; jz < fieldmesh->LocalNz; jz++) {
+      BoutReal val =
+          0.5
+          * (f3d_fa(r.ind, fieldmesh->yend, jz) + f3d_fa(r.ind, fieldmesh->yend + 1, jz));
+      (*this)(r.ind, fieldmesh->yend + 1, jz) =
+          2. * val - (*this)(r.ind, fieldmesh->yend, jz);
+    }
+  }
+
+  (*this) = fromFieldAligned(*this);
+}
+}
+
 void Field3D::applyParallelBoundary() {
 
   TRACE("Field3D::applyParallelBoundary()");