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,19 +451,81 @@ 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
-        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);
+  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 val =
+              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()) {
+      for (int jy = fieldmesh->ystart; jy <= fieldmesh->yend; ++jy) {
+        for (int jz = 0; jz < fieldmesh->LocalNz; jz++) {
+          BoutReal val =
+              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, z));
+        this_next(x, y, 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);
+    (*this) = toFieldAligned(*this);
+
+    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);
   }
 }