From d26b825f2c65fbd946623526e79009bddbd72da0 Mon Sep 17 00:00:00 2001
From: David Bold <dave@ipp.mpg.de>
Date: Tue, 18 Mar 2025 15:32:02 +0100
Subject: [PATCH 01/24] Track yoffset in SpecificInd

Useful for Field3D with parallel slices, so we can write `f[i.yp()]`
instead of `f.yup()[i.yp()]`
---
 include/bout/region.hxx | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/include/bout/region.hxx b/include/bout/region.hxx
index f441b3edd7..a249aae272 100644
--- a/include/bout/region.hxx
+++ b/include/bout/region.hxx
@@ -168,9 +168,11 @@ template <IND_TYPE N>
 struct SpecificInd {
   int ind = -1;         ///< 1D index into Field
   int ny = -1, nz = -1; ///< Sizes of y and z dimensions
+  int yoffset = 0;
 
   SpecificInd() = default;
   SpecificInd(int i, int ny, int nz) : ind(i), ny(ny), nz(nz){};
+  SpecificInd(int i, int ny, int nz, int yoffset) : ind(i), ny(ny), nz(nz), yoffset(yoffset) {};
   explicit SpecificInd(int i) : ind(i){};
 
   /// Allow explicit conversion to an int
@@ -286,7 +288,7 @@ struct SpecificInd {
     }
 #endif
     ASSERT3(std::abs(dy) < ny);
-    return {ind + (dy * nz), ny, nz};
+    return {ind + (dy * nz), ny, nz, yoffset + dy};
   }
   /// The index one point -1 in y
   inline SpecificInd ym(int dy = 1) const { return yp(-dy); }

From b27d19ddf1df1331a94e68610b61bd9393c609a5 Mon Sep 17 00:00:00 2001
From: David Bold <dave@ipp.mpg.de>
Date: Tue, 18 Mar 2025 15:33:52 +0100
Subject: [PATCH 02/24] Track yoffset in Field3D

This is needed wo we can both write `f[i.yp()]` while allowing old code
to call `f.yup()[i.yp()]` as well.
---
 include/bout/field3d.hxx                     |  3 ++
 src/field/field3d.cxx                        | 19 +++++++++-
 src/field/gen_fieldops.jinja                 |  4 ++
 src/field/generated_fieldops.cxx             | 40 ++++++++++++++++++--
 src/mesh/interpolation/hermite_spline_xz.cxx |  3 ++
 src/mesh/parallel/fci.cxx                    |  4 +-
 src/mesh/parallel/identity.cxx               |  2 +
 src/mesh/parallel/shiftedmetricinterp.cxx    |  1 +
 8 files changed, 69 insertions(+), 7 deletions(-)

diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx
index 0643efbe0a..950e3ba150 100644
--- a/include/bout/field3d.hxx
+++ b/include/bout/field3d.hxx
@@ -547,6 +547,9 @@ private:
   template <class T>
   Options* track(const T& change, std::string operation);
   Options* track(const BoutReal& change, std::string operation);
+public:
+  int yoffset{0};
+  void setParallelRegions();
 };
 
 // Non-member overloaded operators
diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx
index 74a6f0853c..fbf339e7e6 100644
--- a/src/field/field3d.cxx
+++ b/src/field/field3d.cxx
@@ -66,7 +66,7 @@ Field3D::Field3D(Mesh* localmesh, CELL_LOC location_in, DirectionTypes direction
 /// later)
 Field3D::Field3D(const Field3D& f)
     : Field(f), data(f.data), yup_fields(f.yup_fields), ydown_fields(f.ydown_fields),
-      regionID(f.regionID) {
+      regionID(f.regionID), yoffset(f.yoffset) {
 
   TRACE("Field3D(Field3D&)");
 
@@ -157,7 +157,9 @@ void Field3D::splitParallelSlices() {
     // Note the fields constructed here will be fully overwritten by the
     // ParallelTransform, so we don't need a full constructor
     yup_fields.emplace_back(fieldmesh);
+    yup_fields[i].yoffset = i + 1;
     ydown_fields.emplace_back(fieldmesh);
+    ydown_fields[i].yoffset = -1 - i;
     if (isFci()) {
       yup_fields[i].setRegion(fmt::format("RGN_YPAR_{:+d}", i + 1));
       ydown_fields[i].setRegion(fmt::format("RGN_YPAR_{:+d}", -i - 1));
@@ -224,6 +226,16 @@ bool Field3D::requiresTwistShift(bool twist_shift_enabled) {
                                                                      getDirectionY());
 }
 
+void Field3D::setParallelRegions() {
+  for (int i = 0; i < fieldmesh->ystart; ++i ) {
+    yup_fields[i].yoffset = i + 1;
+    yup_fields[i].setRegion(fmt::format("RGN_YPAR_{:+d}", i + 1));
+    ydown_fields[i].yoffset = -1 - i;
+    ydown_fields[i].setRegion(fmt::format("RGN_YPAR_{:+d}", -1 - i));
+  }
+}
+
+
 // Not in header because we need to access fieldmesh
 BoutReal& Field3D::operator()(const IndPerp& d, int jy) {
   return operator[](fieldmesh->indPerpto3D(d, jy));
@@ -277,6 +289,7 @@ Field3D& Field3D::operator=(const Field3D& rhs) {
   yup_fields = rhs.yup_fields;
   ydown_fields = rhs.ydown_fields;
   regionID = rhs.regionID;
+  yoffset = rhs.yoffset;
 
   // Copy the data and data sizes
   nx = rhs.nx;
@@ -295,12 +308,13 @@ Field3D& Field3D::operator=(Field3D&& rhs) {
   // Move parallel slices or delete existing ones.
   yup_fields = std::move(rhs.yup_fields);
   ydown_fields = std::move(rhs.ydown_fields);
+  regionID = rhs.regionID;
+  yoffset = rhs.yoffset;
 
   // Move the data and data sizes
   nx = rhs.nx;
   ny = rhs.ny;
   nz = rhs.nz;
-  regionID = rhs.regionID;
 
   data = std::move(rhs.data);
 
@@ -896,6 +910,7 @@ void swap(Field3D& first, Field3D& second) noexcept {
   swap(first.deriv, second.deriv);
   swap(first.yup_fields, second.yup_fields);
   swap(first.ydown_fields, second.ydown_fields);
+  swap(first.yoffset, second.yoffset);
 }
 
 const Region<Ind3D>&
diff --git a/src/field/gen_fieldops.jinja b/src/field/gen_fieldops.jinja
index 88e877c197..4cea18f4c8 100644
--- a/src/field/gen_fieldops.jinja
+++ b/src/field/gen_fieldops.jinja
@@ -12,6 +12,8 @@
     {% if lhs == rhs == "Field3D" %}
       {{out.name}}.setRegion({{lhs.name}}.getMesh()->getCommonRegion({{lhs.name}}.getRegionID(),
                                                                      {{rhs.name}}.getRegionID()));
+      ASSERT2({{lhs.name}}.yoffset == {{rhs.name}}.yoffset);
+      {{out.name}}.yoffset = {{lhs.name}}.yoffset;
 #if BOUT_USE_FCI_AUTOMAGIC
       if ({{lhs.name}}.isFci() and {{lhs.name}}.hasParallelSlices() and {{rhs.name}}.hasParallelSlices()) {
         {{out.name}}.splitParallelSlices();
@@ -23,6 +25,7 @@
 #endif
     {% elif lhs == "Field3D" %}
       {{out.name}}.setRegion({{lhs.name}}.getRegionID());
+      {{out.name}}.yoffset = {{lhs.name}}.yoffset;
       {% if rhs == "BoutReal" %}
 #if BOUT_USE_FCI_AUTOMAGIC
       if ({{lhs.name}}.isFci() and {{lhs.name}}.hasParallelSlices()) {
@@ -36,6 +39,7 @@
       {% endif %}
     {% elif rhs == "Field3D" %}
       {{out.name}}.setRegion({{rhs.name}}.getRegionID());
+      {{out.name}}.yoffset = {{rhs.name}}.yoffset;
       {% if lhs == "BoutReal" %}
 #if BOUT_USE_FCI_AUTOMAGIC
       if ({{rhs.name}}.isFci() and {{rhs.name}}.hasParallelSlices()) {
diff --git a/src/field/generated_fieldops.cxx b/src/field/generated_fieldops.cxx
index 75d2ede82d..4e141f7f4c 100644
--- a/src/field/generated_fieldops.cxx
+++ b/src/field/generated_fieldops.cxx
@@ -15,6 +15,8 @@ Field3D operator*(const Field3D& lhs, const Field3D& rhs) {
   checkData(rhs);
 
   result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID()));
+  ASSERT2(lhs.yoffset == rhs.yoffset);
+  result.yoffset = lhs.yoffset;
 #if BOUT_USE_FCI_AUTOMAGIC
   if (lhs.isFci() and lhs.hasParallelSlices() and rhs.hasParallelSlices()) {
     result.splitParallelSlices();
@@ -87,6 +89,8 @@ Field3D operator/(const Field3D& lhs, const Field3D& rhs) {
   checkData(rhs);
 
   result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID()));
+  ASSERT2(lhs.yoffset == rhs.yoffset);
+  result.yoffset = lhs.yoffset;
 #if BOUT_USE_FCI_AUTOMAGIC
   if (lhs.isFci() and lhs.hasParallelSlices() and rhs.hasParallelSlices()) {
     result.splitParallelSlices();
@@ -159,6 +163,8 @@ Field3D operator+(const Field3D& lhs, const Field3D& rhs) {
   checkData(rhs);
 
   result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID()));
+  ASSERT2(lhs.yoffset == rhs.yoffset);
+  result.yoffset = lhs.yoffset;
 #if BOUT_USE_FCI_AUTOMAGIC
   if (lhs.isFci() and lhs.hasParallelSlices() and rhs.hasParallelSlices()) {
     result.splitParallelSlices();
@@ -231,6 +237,8 @@ Field3D operator-(const Field3D& lhs, const Field3D& rhs) {
   checkData(rhs);
 
   result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID()));
+  ASSERT2(lhs.yoffset == rhs.yoffset);
+  result.yoffset = lhs.yoffset;
 #if BOUT_USE_FCI_AUTOMAGIC
   if (lhs.isFci() and lhs.hasParallelSlices() and rhs.hasParallelSlices()) {
     result.splitParallelSlices();
@@ -303,6 +311,7 @@ Field3D operator*(const Field3D& lhs, const Field2D& rhs) {
   checkData(rhs);
 
   result.setRegion(lhs.getRegionID());
+  result.yoffset = lhs.yoffset;
 
   Mesh* localmesh = lhs.getMesh();
 
@@ -329,7 +338,9 @@ Field3D& Field3D::operator*=(const Field2D& rhs) {
 
     // Delete existing parallel slices. We don't copy parallel slices, so any
     // that currently exist will be incorrect.
-    { clearParallelSlices(); }
+    {
+      clearParallelSlices();
+    }
 
     checkData(*this);
     checkData(rhs);
@@ -364,6 +375,7 @@ Field3D operator/(const Field3D& lhs, const Field2D& rhs) {
   checkData(rhs);
 
   result.setRegion(lhs.getRegionID());
+  result.yoffset = lhs.yoffset;
 
   Mesh* localmesh = lhs.getMesh();
 
@@ -391,7 +403,9 @@ Field3D& Field3D::operator/=(const Field2D& rhs) {
 
     // Delete existing parallel slices. We don't copy parallel slices, so any
     // that currently exist will be incorrect.
-    { clearParallelSlices(); }
+    {
+      clearParallelSlices();
+    }
 
     checkData(*this);
     checkData(rhs);
@@ -427,6 +441,7 @@ Field3D operator+(const Field3D& lhs, const Field2D& rhs) {
   checkData(rhs);
 
   result.setRegion(lhs.getRegionID());
+  result.yoffset = lhs.yoffset;
 
   Mesh* localmesh = lhs.getMesh();
 
@@ -453,7 +468,9 @@ Field3D& Field3D::operator+=(const Field2D& rhs) {
 
     // Delete existing parallel slices. We don't copy parallel slices, so any
     // that currently exist will be incorrect.
-    { clearParallelSlices(); }
+    {
+      clearParallelSlices();
+    }
 
     checkData(*this);
     checkData(rhs);
@@ -488,6 +505,7 @@ Field3D operator-(const Field3D& lhs, const Field2D& rhs) {
   checkData(rhs);
 
   result.setRegion(lhs.getRegionID());
+  result.yoffset = lhs.yoffset;
 
   Mesh* localmesh = lhs.getMesh();
 
@@ -514,7 +532,9 @@ Field3D& Field3D::operator-=(const Field2D& rhs) {
 
     // Delete existing parallel slices. We don't copy parallel slices, so any
     // that currently exist will be incorrect.
-    { clearParallelSlices(); }
+    {
+      clearParallelSlices();
+    }
 
     checkData(*this);
     checkData(rhs);
@@ -640,6 +660,7 @@ Field3D operator*(const Field3D& lhs, const BoutReal rhs) {
   checkData(rhs);
 
   result.setRegion(lhs.getRegionID());
+  result.yoffset = lhs.yoffset;
 #if BOUT_USE_FCI_AUTOMAGIC
   if (lhs.isFci() and lhs.hasParallelSlices()) {
     result.splitParallelSlices();
@@ -708,6 +729,7 @@ Field3D operator/(const Field3D& lhs, const BoutReal rhs) {
   checkData(rhs);
 
   result.setRegion(lhs.getRegionID());
+  result.yoffset = lhs.yoffset;
 #if BOUT_USE_FCI_AUTOMAGIC
   if (lhs.isFci() and lhs.hasParallelSlices()) {
     result.splitParallelSlices();
@@ -778,6 +800,7 @@ Field3D operator+(const Field3D& lhs, const BoutReal rhs) {
   checkData(rhs);
 
   result.setRegion(lhs.getRegionID());
+  result.yoffset = lhs.yoffset;
 #if BOUT_USE_FCI_AUTOMAGIC
   if (lhs.isFci() and lhs.hasParallelSlices()) {
     result.splitParallelSlices();
@@ -846,6 +869,7 @@ Field3D operator-(const Field3D& lhs, const BoutReal rhs) {
   checkData(rhs);
 
   result.setRegion(lhs.getRegionID());
+  result.yoffset = lhs.yoffset;
 #if BOUT_USE_FCI_AUTOMAGIC
   if (lhs.isFci() and lhs.hasParallelSlices()) {
     result.splitParallelSlices();
@@ -915,6 +939,7 @@ Field3D operator*(const Field2D& lhs, const Field3D& rhs) {
   checkData(rhs);
 
   result.setRegion(rhs.getRegionID());
+  result.yoffset = rhs.yoffset;
 
   Mesh* localmesh = lhs.getMesh();
 
@@ -941,6 +966,7 @@ Field3D operator/(const Field2D& lhs, const Field3D& rhs) {
   checkData(rhs);
 
   result.setRegion(rhs.getRegionID());
+  result.yoffset = rhs.yoffset;
 
   Mesh* localmesh = lhs.getMesh();
 
@@ -967,6 +993,7 @@ Field3D operator+(const Field2D& lhs, const Field3D& rhs) {
   checkData(rhs);
 
   result.setRegion(rhs.getRegionID());
+  result.yoffset = rhs.yoffset;
 
   Mesh* localmesh = lhs.getMesh();
 
@@ -993,6 +1020,7 @@ Field3D operator-(const Field2D& lhs, const Field3D& rhs) {
   checkData(rhs);
 
   result.setRegion(rhs.getRegionID());
+  result.yoffset = rhs.yoffset;
 
   Mesh* localmesh = lhs.getMesh();
 
@@ -2209,6 +2237,7 @@ Field3D operator*(const BoutReal lhs, const Field3D& rhs) {
   checkData(rhs);
 
   result.setRegion(rhs.getRegionID());
+  result.yoffset = rhs.yoffset;
 #if BOUT_USE_FCI_AUTOMAGIC
   if (rhs.isFci() and rhs.hasParallelSlices()) {
     result.splitParallelSlices();
@@ -2238,6 +2267,7 @@ Field3D operator/(const BoutReal lhs, const Field3D& rhs) {
   checkData(rhs);
 
   result.setRegion(rhs.getRegionID());
+  result.yoffset = rhs.yoffset;
 #if BOUT_USE_FCI_AUTOMAGIC
   if (rhs.isFci() and rhs.hasParallelSlices()) {
     result.splitParallelSlices();
@@ -2267,6 +2297,7 @@ Field3D operator+(const BoutReal lhs, const Field3D& rhs) {
   checkData(rhs);
 
   result.setRegion(rhs.getRegionID());
+  result.yoffset = rhs.yoffset;
 #if BOUT_USE_FCI_AUTOMAGIC
   if (rhs.isFci() and rhs.hasParallelSlices()) {
     result.splitParallelSlices();
@@ -2296,6 +2327,7 @@ Field3D operator-(const BoutReal lhs, const Field3D& rhs) {
   checkData(rhs);
 
   result.setRegion(rhs.getRegionID());
+  result.yoffset = rhs.yoffset;
 #if BOUT_USE_FCI_AUTOMAGIC
   if (rhs.isFci() and rhs.hasParallelSlices()) {
     result.splitParallelSlices();
diff --git a/src/mesh/interpolation/hermite_spline_xz.cxx b/src/mesh/interpolation/hermite_spline_xz.cxx
index 03a062df42..394002424e 100644
--- a/src/mesh/interpolation/hermite_spline_xz.cxx
+++ b/src/mesh/interpolation/hermite_spline_xz.cxx
@@ -378,6 +378,7 @@ Field3D XZHermiteSplineBase<monotonic>::interpolate(const Field3D& f,
 
   ASSERT1(f.getMesh() == localmesh);
   Field3D f_interp{emptyFrom(f)};
+  f_interp.yoffset = y_offset;
 
   const auto region2 = y_offset != 0 ? fmt::format("RGN_YPAR_{:+d}", y_offset) : region;
 
@@ -395,6 +396,7 @@ Field3D XZHermiteSplineBase<monotonic>::interpolate(const Field3D& f,
   MatMult(petscWeights, rhs, result);
   VecGetArrayRead(result, &cptr);
   BOUT_FOR(i, f.getRegion(region2)) {
+    i.yoffset = y_offset;
     f_interp[i] = cptr[int(i)];
     if constexpr (monotonic) {
       const auto iyp = i;
@@ -508,6 +510,7 @@ Field3D XZMonotonicHermiteSplineLegacy::interpolate(const Field3D& f,
   ASSERT1(f.getMesh() == localmesh);
   Field3D f_interp(f.getMesh());
   f_interp.allocate();
+  f_interp.yoffset = y_offset;
 
   // Derivatives are used for tension and need to be on dimensionless
   // coordinates
diff --git a/src/mesh/parallel/fci.cxx b/src/mesh/parallel/fci.cxx
index 580897f47a..0c3d2625ca 100644
--- a/src/mesh/parallel/fci.cxx
+++ b/src/mesh/parallel/fci.cxx
@@ -99,6 +99,7 @@ bool load_parallel_metric_component(std::string name, Field3D& component, int of
   pcom.setRegion(fmt::format("RGN_YPAR_{:+d}", offset));
   pcom.name = name;
   BOUT_FOR(i, component.getRegion("RGN_NOBNDRY")) { pcom[i.yp(offset)] = tmp[i]; }
+  ASSERT2(pcom.yoffset == offset);
   return isValid;
 }
 #endif
@@ -430,8 +431,8 @@ void FCITransform::calcParallelSlices(Field3D& f) {
   // Interpolate f onto yup and ydown fields
   for (const auto& map : field_line_maps) {
     f.ynext(map.offset) = map.interpolate(f);
-    f.ynext(map.offset).setRegion(fmt::format("RGN_YPAR_{:+d}", map.offset));
   }
+  f.setParallelRegions();
 }
 
 void FCITransform::integrateParallelSlices(Field3D& f) {
@@ -449,6 +450,7 @@ void FCITransform::integrateParallelSlices(Field3D& f) {
   for (const auto& map : field_line_maps) {
     f.ynext(map.offset) = map.integrate(f);
   }
+  f.setParallelRegions();
 }
 
 void FCITransform::loadParallelMetrics(Coordinates* coords) {
diff --git a/src/mesh/parallel/identity.cxx b/src/mesh/parallel/identity.cxx
index 90fcaef45e..b9869a22eb 100644
--- a/src/mesh/parallel/identity.cxx
+++ b/src/mesh/parallel/identity.cxx
@@ -23,7 +23,9 @@ void ParallelTransformIdentity::calcParallelSlices(Field3D& f) {
   f.splitParallelSlices();
 
   for (int i = 0; i < f.getMesh()->ystart; ++i) {
+    f_copy.yoffset = i + 1;
     f.yup(i) = f_copy;
+    f_copy.yoffset = -1 - i;
     f.ydown(i) = f_copy;
   }
 }
diff --git a/src/mesh/parallel/shiftedmetricinterp.cxx b/src/mesh/parallel/shiftedmetricinterp.cxx
index dfb397c626..bfe53b38e6 100644
--- a/src/mesh/parallel/shiftedmetricinterp.cxx
+++ b/src/mesh/parallel/shiftedmetricinterp.cxx
@@ -220,6 +220,7 @@ void ShiftedMetricInterp::calcParallelSlices(Field3D& f) {
   for (const auto& interp : parallel_slice_interpolators) {
     f.ynext(interp->y_offset) = interp->interpolate(f);
   }
+  f.setParallelRegions();
 }
 
 /*!

From aaf102d7a96cc9888b96fdd0370b1ed8abb5a02f Mon Sep 17 00:00:00 2001
From: David Bold <dave@ipp.mpg.de>
Date: Tue, 18 Mar 2025 15:34:19 +0100
Subject: [PATCH 03/24] Add some checks to ensure the parallel offset is
 correct

---
 include/bout/field3d.hxx | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx
index 950e3ba150..8864b52630 100644
--- a/include/bout/field3d.hxx
+++ b/include/bout/field3d.hxx
@@ -271,23 +271,27 @@ public:
   /// Return reference to yup field
   Field3D& yup(std::vector<Field3D>::size_type index = 0) {
     ASSERT2(index < yup_fields.size());
+    ASSERT2(yup_fields[index].yoffset == static_cast<int>(index) + 1);
     return yup_fields[index];
   }
   /// Return const reference to yup field
   const Field3D& yup(std::vector<Field3D>::size_type index = 0) const {
     ASSERT2(index < yup_fields.size());
+    ASSERT2(yup_fields[index].yoffset == static_cast<int>(index) + 1);
     return yup_fields[index];
   }
 
   /// Return reference to ydown field
   Field3D& ydown(std::vector<Field3D>::size_type index = 0) {
     ASSERT2(index < ydown_fields.size());
+    ASSERT2(ydown_fields[index].yoffset == -1 - static_cast<int>(index));
     return ydown_fields[index];
   }
 
   /// Return const reference to ydown field
   const Field3D& ydown(std::vector<Field3D>::size_type index = 0) const {
     ASSERT2(index < ydown_fields.size());
+    ASSERT2(ydown_fields[index].yoffset == -1 - static_cast<int>(index));
     return ydown_fields[index];
   }
 

From feec777a13439de2819e0a35612dc778ac0e9dfa Mon Sep 17 00:00:00 2001
From: David Bold <dave@ipp.mpg.de>
Date: Tue, 18 Mar 2025 16:29:01 +0100
Subject: [PATCH 04/24] Automatically select the appropriate parallel slice

---
 include/bout/field3d.hxx | 19 +++++++++++++++++--
 1 file changed, 17 insertions(+), 2 deletions(-)

diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx
index 8864b52630..5bcd70cfeb 100644
--- a/include/bout/field3d.hxx
+++ b/include/bout/field3d.hxx
@@ -355,8 +355,23 @@ public:
     return std::end(getRegion("RGN_ALL"));
   };
 
-  BoutReal& operator[](const Ind3D& d) { return data[d.ind]; }
-  const BoutReal& operator[](const Ind3D& d) const { return data[d.ind]; }
+  BoutReal& operator[](const Ind3D& d) {
+    if (d.yoffset) {
+      if (yoffset == 0) {
+	if (hasParallelSlices()) {
+	  return ynext(d.yoffset)[d];
+	}
+#if CHECK >= 2
+	else if (isFci()) { // We probably should assert here that this is field aligned
+	  throw BoutException("Tried to access parallel slices, but they are not calculated!");
+	}
+#endif
+      } else {
+	ASSERT2(d.yoffset == yoffset);
+      }
+    }
+    return data[d.ind]; }
+  const BoutReal& operator[](const Ind3D& d) const { return (*const_cast<Field3D*>(this))[d]; }
 
   BoutReal& operator()(const IndPerp& d, int jy);
   const BoutReal& operator()(const IndPerp& d, int jy) const;

From 49bae2127128e03bad38d9a170a57ceef2554cac Mon Sep 17 00:00:00 2001
From: dschwoerer <5637662+dschwoerer@users.noreply.github.com>
Date: Tue, 18 Mar 2025 15:31:50 +0000
Subject: [PATCH 05/24] Apply clang-format changes

---
 include/bout/field3d.hxx         | 23 ++++++++++++++---------
 include/bout/region.hxx          |  3 ++-
 src/field/field3d.cxx            |  3 +--
 src/field/generated_fieldops.cxx | 16 ++++------------
 4 files changed, 21 insertions(+), 24 deletions(-)

diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx
index 5bcd70cfeb..39c4f5780a 100644
--- a/include/bout/field3d.hxx
+++ b/include/bout/field3d.hxx
@@ -358,20 +358,24 @@ public:
   BoutReal& operator[](const Ind3D& d) {
     if (d.yoffset) {
       if (yoffset == 0) {
-	if (hasParallelSlices()) {
-	  return ynext(d.yoffset)[d];
-	}
+        if (hasParallelSlices()) {
+          return ynext(d.yoffset)[d];
+        }
 #if CHECK >= 2
-	else if (isFci()) { // We probably should assert here that this is field aligned
-	  throw BoutException("Tried to access parallel slices, but they are not calculated!");
-	}
+        else if (isFci()) { // We probably should assert here that this is field aligned
+          throw BoutException(
+              "Tried to access parallel slices, but they are not calculated!");
+        }
 #endif
       } else {
-	ASSERT2(d.yoffset == yoffset);
+        ASSERT2(d.yoffset == yoffset);
       }
     }
-    return data[d.ind]; }
-  const BoutReal& operator[](const Ind3D& d) const { return (*const_cast<Field3D*>(this))[d]; }
+    return data[d.ind];
+  }
+  const BoutReal& operator[](const Ind3D& d) const {
+    return (*const_cast<Field3D*>(this))[d];
+  }
 
   BoutReal& operator()(const IndPerp& d, int jy);
   const BoutReal& operator()(const IndPerp& d, int jy) const;
@@ -566,6 +570,7 @@ private:
   template <class T>
   Options* track(const T& change, std::string operation);
   Options* track(const BoutReal& change, std::string operation);
+
 public:
   int yoffset{0};
   void setParallelRegions();
diff --git a/include/bout/region.hxx b/include/bout/region.hxx
index a249aae272..7edd63f912 100644
--- a/include/bout/region.hxx
+++ b/include/bout/region.hxx
@@ -172,7 +172,8 @@ struct SpecificInd {
 
   SpecificInd() = default;
   SpecificInd(int i, int ny, int nz) : ind(i), ny(ny), nz(nz){};
-  SpecificInd(int i, int ny, int nz, int yoffset) : ind(i), ny(ny), nz(nz), yoffset(yoffset) {};
+  SpecificInd(int i, int ny, int nz, int yoffset)
+      : ind(i), ny(ny), nz(nz), yoffset(yoffset){};
   explicit SpecificInd(int i) : ind(i){};
 
   /// Allow explicit conversion to an int
diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx
index fbf339e7e6..fb76244016 100644
--- a/src/field/field3d.cxx
+++ b/src/field/field3d.cxx
@@ -227,7 +227,7 @@ bool Field3D::requiresTwistShift(bool twist_shift_enabled) {
 }
 
 void Field3D::setParallelRegions() {
-  for (int i = 0; i < fieldmesh->ystart; ++i ) {
+  for (int i = 0; i < fieldmesh->ystart; ++i) {
     yup_fields[i].yoffset = i + 1;
     yup_fields[i].setRegion(fmt::format("RGN_YPAR_{:+d}", i + 1));
     ydown_fields[i].yoffset = -1 - i;
@@ -235,7 +235,6 @@ void Field3D::setParallelRegions() {
   }
 }
 
-
 // Not in header because we need to access fieldmesh
 BoutReal& Field3D::operator()(const IndPerp& d, int jy) {
   return operator[](fieldmesh->indPerpto3D(d, jy));
diff --git a/src/field/generated_fieldops.cxx b/src/field/generated_fieldops.cxx
index 4e141f7f4c..c1c450b2a9 100644
--- a/src/field/generated_fieldops.cxx
+++ b/src/field/generated_fieldops.cxx
@@ -338,9 +338,7 @@ Field3D& Field3D::operator*=(const Field2D& rhs) {
 
     // Delete existing parallel slices. We don't copy parallel slices, so any
     // that currently exist will be incorrect.
-    {
-      clearParallelSlices();
-    }
+    { clearParallelSlices(); }
 
     checkData(*this);
     checkData(rhs);
@@ -403,9 +401,7 @@ Field3D& Field3D::operator/=(const Field2D& rhs) {
 
     // Delete existing parallel slices. We don't copy parallel slices, so any
     // that currently exist will be incorrect.
-    {
-      clearParallelSlices();
-    }
+    { clearParallelSlices(); }
 
     checkData(*this);
     checkData(rhs);
@@ -468,9 +464,7 @@ Field3D& Field3D::operator+=(const Field2D& rhs) {
 
     // Delete existing parallel slices. We don't copy parallel slices, so any
     // that currently exist will be incorrect.
-    {
-      clearParallelSlices();
-    }
+    { clearParallelSlices(); }
 
     checkData(*this);
     checkData(rhs);
@@ -532,9 +526,7 @@ Field3D& Field3D::operator-=(const Field2D& rhs) {
 
     // Delete existing parallel slices. We don't copy parallel slices, so any
     // that currently exist will be incorrect.
-    {
-      clearParallelSlices();
-    }
+    { clearParallelSlices(); }
 
     checkData(*this);
     checkData(rhs);

From 247c2fb16c979081f0762efaa8d5a35e407530fc Mon Sep 17 00:00:00 2001
From: David Bold <dave@ipp.mpg.de>
Date: Wed, 19 Mar 2025 10:44:49 +0100
Subject: [PATCH 06/24] Fix: preserve yoffset with other operations

---
 include/bout/region.hxx | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/include/bout/region.hxx b/include/bout/region.hxx
index 7edd63f912..6151bc6f1a 100644
--- a/include/bout/region.hxx
+++ b/include/bout/region.hxx
@@ -278,7 +278,9 @@ struct SpecificInd {
     }
   }
 
-  inline SpecificInd xp(int dx = 1) const { return {ind + (dx * ny * nz), ny, nz}; }
+  inline SpecificInd xp(int dx = 1) const {
+    return {ind + (dx * ny * nz), ny, nz, yoffset};
+  }
   /// The index one point -1 in x
   inline SpecificInd xm(int dx = 1) const { return xp(-dx); }
   /// The index one point +1 in y
@@ -300,7 +302,7 @@ struct SpecificInd {
   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};
+    return {(ind + dz) % nz < dz ? ind - nz + dz : ind + dz, ny, nz, yoffset};
   }
   /// The index one point -1 in z. Wraps around zstart to zend
   /// An alternative, non-branching calculation is :
@@ -309,7 +311,7 @@ struct SpecificInd {
   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};
+    return {(ind) % nz < dz ? ind + nz - dz : ind - dz, ny, nz, yoffset};
   }
   /// Automatically select zm or zp depending on sign
   inline SpecificInd zpm(int dz) const { return dz > 0 ? zp(dz) : zm(-dz); }

From 6e834db4880ff3185222c00fc3c8ee200a7bb0e4 Mon Sep 17 00:00:00 2001
From: David Bold <dave@ipp.mpg.de>
Date: Wed, 19 Mar 2025 10:45:43 +0100
Subject: [PATCH 07/24] Add a version of `.i.yp()` that does not shift the
 parallel field

Interpolation methods might need that, but in general it is unlike to be
of use.
---
 include/bout/region.hxx | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/include/bout/region.hxx b/include/bout/region.hxx
index 6151bc6f1a..e7750747ed 100644
--- a/include/bout/region.hxx
+++ b/include/bout/region.hxx
@@ -293,6 +293,16 @@ struct SpecificInd {
     ASSERT3(std::abs(dy) < ny);
     return {ind + (dy * nz), ny, nz, yoffset + dy};
   }
+  inline SpecificInd yp_no_parallel_shift(int dy = 1) const {
+#if CHECK >= 4
+    if (y() + dy < 0 or y() + dy >= ny) {
+      throw BoutException("Offset in y ({:d}) would go out of bounds at {:d}", dy, ind);
+    }
+#endif
+    ASSERT3(std::abs(dy) < ny);
+    return {ind + (dy * nz), ny, nz, yoffset};
+  }
+
   /// 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

From 242b9b2be958e87dccc1c38ca9149a666bed5ae1 Mon Sep 17 00:00:00 2001
From: David Bold <dave@ipp.mpg.de>
Date: Wed, 19 Mar 2025 10:51:32 +0100
Subject: [PATCH 08/24] Add y-shifted regions by default

---
 src/mesh/mesh.cxx                             |  9 +++++++
 src/mesh/parallel/fci.cxx                     |  9 -------
 .../unit/mesh/parallel/test_shiftedmetric.cxx | 27 +++----------------
 3 files changed, 13 insertions(+), 32 deletions(-)

diff --git a/src/mesh/mesh.cxx b/src/mesh/mesh.cxx
index cb6cb8410c..8cc204ff41 100644
--- a/src/mesh/mesh.cxx
+++ b/src/mesh/mesh.cxx
@@ -685,6 +685,15 @@ void Mesh::createDefaultRegions() {
                                 + getRegion3D("RGN_YGUARDS") + getRegion3D("RGN_ZGUARDS"))
                                    .unique());
 
+  // Add region for parallel slices
+  for (int i = 1; i <= ystart; ++i) {
+    for (const int offset : {-i, i}) {
+      const auto region = fmt::format("RGN_YPAR_{:+d}", offset);
+      addRegion3D(region, Region<Ind3D>(xstart, xend, ystart + offset, yend + offset, 0,
+                                        LocalNz - 1, LocalNy, LocalNz));
+    }
+  }
+
   //2D regions
   addRegion2D("RGN_ALL", Region<Ind2D>(0, LocalNx - 1, 0, LocalNy - 1, 0, 0, LocalNy, 1,
                                        maxregionblocksize));
diff --git a/src/mesh/parallel/fci.cxx b/src/mesh/parallel/fci.cxx
index 0c3d2625ca..9d6997f660 100644
--- a/src/mesh/parallel/fci.cxx
+++ b/src/mesh/parallel/fci.cxx
@@ -340,15 +340,6 @@ FCIMap::FCIMap(Mesh& mesh, const Coordinates::FieldMetric& UNUSED(dy), Options&
   region_no_boundary = region_no_boundary.mask(to_remove);
 
   interp->setRegion(region_no_boundary);
-
-  const auto region = fmt::format("RGN_YPAR_{:+d}", offset);
-  if (not map_mesh.hasRegion3D(region)) {
-    // The valid region for this slice
-    map_mesh.addRegion3D(
-        region, Region<Ind3D>(map_mesh.xstart, map_mesh.xend, map_mesh.ystart + offset,
-                              map_mesh.yend + offset, 0, map_mesh.LocalNz - 1,
-                              map_mesh.LocalNy, map_mesh.LocalNz));
-  }
 }
 
 Field3D FCIMap::integrate(Field3D& f) const {
diff --git a/tests/unit/mesh/parallel/test_shiftedmetric.cxx b/tests/unit/mesh/parallel/test_shiftedmetric.cxx
index 7e111e9c4e..e1ad088404 100644
--- a/tests/unit/mesh/parallel/test_shiftedmetric.cxx
+++ b/tests/unit/mesh/parallel/test_shiftedmetric.cxx
@@ -284,25 +284,6 @@ TEST_F(ShiftedMetricTest, ToFromFieldAlignedFieldPerp) {
 }
 
 TEST_F(ShiftedMetricTest, CalcParallelSlices) {
-  // We don't shift in the guard cells, and the parallel slices are
-  // stored offset in y, therefore we need to make new regions that we
-  // can compare the expected and actual outputs over
-  output_info.disable();
-  mesh->addRegion3D("RGN_YUP",
-                    Region<Ind3D>(0, mesh->LocalNx - 1, mesh->ystart + 1, mesh->yend + 1,
-                                  0, mesh->LocalNz - 1, mesh->LocalNy, mesh->LocalNz));
-  mesh->addRegion3D("RGN_YUP2",
-                    Region<Ind3D>(0, mesh->LocalNx - 1, mesh->ystart + 2, mesh->yend + 2,
-                                  0, mesh->LocalNz - 1, mesh->LocalNy, mesh->LocalNz));
-
-  mesh->addRegion3D("RGN_YDOWN",
-                    Region<Ind3D>(0, mesh->LocalNx - 1, mesh->ystart - 1, mesh->yend - 1,
-                                  0, mesh->LocalNz - 1, mesh->LocalNy, mesh->LocalNz));
-  mesh->addRegion3D("RGN_YDOWN2",
-                    Region<Ind3D>(0, mesh->LocalNx - 1, mesh->ystart - 2, mesh->yend - 2,
-                                  0, mesh->LocalNz - 1, mesh->LocalNy, mesh->LocalNz));
-  output_info.enable();
-
   // Actual interesting bit here!
   input.getCoordinates()->getParallelTransform().calcParallelSlices(input);
   // Expected output values
@@ -412,9 +393,9 @@ TEST_F(ShiftedMetricTest, CalcParallelSlices) {
                               {0., 0., 0., 0., 0.},
                               {0., 0., 0., 0., 0.}}});
 
-  EXPECT_TRUE(IsFieldEqual(input.ynext(1), expected_up_1, "RGN_YUP", FFTTolerance));
-  EXPECT_TRUE(IsFieldEqual(input.ynext(2), expected_up_2, "RGN_YUP2", FFTTolerance));
-  EXPECT_TRUE(IsFieldEqual(input.ynext(-1), expected_down_1, "RGN_YDOWN", FFTTolerance));
-  EXPECT_TRUE(IsFieldEqual(input.ynext(-2), expected_down2, "RGN_YDOWN2", FFTTolerance));
+  EXPECT_TRUE(IsFieldEqual(input.ynext(1), expected_up_1, "RGN_YPAR_+1", FFTTolerance));
+  EXPECT_TRUE(IsFieldEqual(input.ynext(2), expected_up_2, "RGN_YPAR_+2", FFTTolerance));
+  EXPECT_TRUE(IsFieldEqual(input.ynext(-1), expected_down_1, "RGN_YPAR_-1", FFTTolerance));
+  EXPECT_TRUE(IsFieldEqual(input.ynext(-2), expected_down2, "RGN_YPAR_-2", FFTTolerance));
 }
 #endif

From 53c9e31d6831f7b5ee06e286f39ff0fd4db78c08 Mon Sep 17 00:00:00 2001
From: David Bold <dave@ipp.mpg.de>
Date: Wed, 19 Mar 2025 10:51:56 +0100
Subject: [PATCH 09/24] Fix hermite_spline_z for auto-shift

---
 src/mesh/interpolation/hermite_spline_z.cxx | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/mesh/interpolation/hermite_spline_z.cxx b/src/mesh/interpolation/hermite_spline_z.cxx
index c4c44cb7b4..02d8d45a61 100644
--- a/src/mesh/interpolation/hermite_spline_z.cxx
+++ b/src/mesh/interpolation/hermite_spline_z.cxx
@@ -151,6 +151,7 @@ Field3D ZHermiteSpline::interpolate(const Field3D& f,
 
   ASSERT1(f.getMesh() == localmesh);
   Field3D f_interp{emptyFrom(f)};
+  f_interp.yoffset = y_offset;
 
   std::string local_fz_region;
   if (region_str == "DEFAULT") {
@@ -168,7 +169,7 @@ Field3D ZHermiteSpline::interpolate(const Field3D& f,
   Field3D fz = bout::derivatives::index::DDZ(f, CELL_DEFAULT, "DEFAULT", local_fz_region);
 
   BOUT_FOR(i, local_region) {
-    const auto corner = k_corner[i.ind].yp(y_offset);
+    const auto corner = k_corner[i.ind].yp_no_parallel_shift(y_offset);
     const auto corner_zp1 = corner.zp();
 
     // Interpolate in Z

From fe888aeba67974fe07402b34c96a988d834a15e2 Mon Sep 17 00:00:00 2001
From: David Bold <dave@ipp.mpg.de>
Date: Wed, 19 Mar 2025 10:52:09 +0100
Subject: [PATCH 10/24] Avoid warning with 3D metrics

---
 src/physics/smoothing.cxx | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/physics/smoothing.cxx b/src/physics/smoothing.cxx
index 0a1391907f..036c475cd7 100644
--- a/src/physics/smoothing.cxx
+++ b/src/physics/smoothing.cxx
@@ -350,7 +350,7 @@ BoutReal Average_XY(const Field2D& var) {
   return Vol_Glb;
 }
 
-BoutReal Vol_Integral(const Field2D& var) {
+BoutReal Vol_Integral([[maybe_unused]] const Field2D& var) {
 #if BOUT_USE_METRIC_3D
   AUTO_TRACE();
   throw BoutException("Vol_Intregral currently incompatible with 3D metrics");

From 76dc241a041b9d5cf75636af47c97d6aa99b34c6 Mon Sep 17 00:00:00 2001
From: dschwoerer <5637662+dschwoerer@users.noreply.github.com>
Date: Wed, 19 Mar 2025 09:53:00 +0000
Subject: [PATCH 11/24] Apply clang-format changes

---
 tests/unit/mesh/parallel/test_shiftedmetric.cxx | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/tests/unit/mesh/parallel/test_shiftedmetric.cxx b/tests/unit/mesh/parallel/test_shiftedmetric.cxx
index e1ad088404..3d7a52607f 100644
--- a/tests/unit/mesh/parallel/test_shiftedmetric.cxx
+++ b/tests/unit/mesh/parallel/test_shiftedmetric.cxx
@@ -395,7 +395,8 @@ TEST_F(ShiftedMetricTest, CalcParallelSlices) {
 
   EXPECT_TRUE(IsFieldEqual(input.ynext(1), expected_up_1, "RGN_YPAR_+1", FFTTolerance));
   EXPECT_TRUE(IsFieldEqual(input.ynext(2), expected_up_2, "RGN_YPAR_+2", FFTTolerance));
-  EXPECT_TRUE(IsFieldEqual(input.ynext(-1), expected_down_1, "RGN_YPAR_-1", FFTTolerance));
+  EXPECT_TRUE(
+      IsFieldEqual(input.ynext(-1), expected_down_1, "RGN_YPAR_-1", FFTTolerance));
   EXPECT_TRUE(IsFieldEqual(input.ynext(-2), expected_down2, "RGN_YPAR_-2", FFTTolerance));
 }
 #endif

From eeb7c14eafa7cdc29171ea8248c5159a71001f29 Mon Sep 17 00:00:00 2001
From: David Bold <dave@ipp.mpg.de>
Date: Mon, 2 Jun 2025 10:39:40 +0200
Subject: [PATCH 12/24] Ajust tests for auto parallel slices

The tests are probably not meant to work with parallel slices.
---
 .../unit/include/bout/test_petsc_indexer.cxx  | 26 ++++++++++---------
 tests/unit/include/bout/test_petsc_matrix.cxx |  6 ++---
 2 files changed, 17 insertions(+), 15 deletions(-)

diff --git a/tests/unit/include/bout/test_petsc_indexer.cxx b/tests/unit/include/bout/test_petsc_indexer.cxx
index 0e9f597eae..70745e7203 100644
--- a/tests/unit/include/bout/test_petsc_indexer.cxx
+++ b/tests/unit/include/bout/test_petsc_indexer.cxx
@@ -474,9 +474,9 @@ TEST_P(ParallelIndexerTest, TestConvertIndex3D) {
         const int global = index->getGlobal(i.xm()),
                   otherXind = (this->pe_xind == 0) ? this->pe_xind : this->pe_xind - 1,
                   otherYind = this->pe_yind - 1;
-        const Ind3D otherInd =
-            i.offset((this->pe_xind == 0) ? -1 : this->xend - this->xstart,
-                     this->yend - this->ystart + 1, 0);
+        Ind3D otherInd = i.offset((this->pe_xind == 0) ? -1 : this->xend - this->xstart,
+                                  this->yend - this->ystart + 1, 0);
+        otherInd.yoffset = 0;
         const int otherGlobal = this->getIndexer<Field3D>(indexers, otherXind, otherYind)
                                     ->getGlobal(otherInd);
         EXPECT_NE(global, -1);
@@ -486,7 +486,7 @@ TEST_P(ParallelIndexerTest, TestConvertIndex3D) {
       int global = index->getGlobal(i);
       int otherGlobal =
           this->getIndexer<Field3D>(indexers, this->pe_xind, this->pe_yind - 1)
-              ->getGlobal(i.yp(this->yend - this->ystart + 1));
+              ->getGlobal(i.yp_no_parallel_shift(this->yend - this->ystart + 1));
       EXPECT_EQ(global, otherGlobal);
 
       if (i.x() == this->xend) {
@@ -494,9 +494,10 @@ TEST_P(ParallelIndexerTest, TestConvertIndex3D) {
                   otherXind = (this->pe_xind == this->nxpe - 1) ? this->pe_xind
                                                                 : this->pe_xind + 1,
                   otherYind = this->pe_yind - 1;
-        const Ind3D otherInd =
+        Ind3D otherInd =
             i.offset((this->pe_xind == this->nxpe - 1) ? 1 : this->xstart - this->xend,
                      this->yend - this->ystart + 1, 0);
+        otherInd.yoffset = 0;
         const int otherGlobal = this->getIndexer<Field3D>(indexers, otherXind, otherYind)
                                     ->getGlobal(otherInd);
         EXPECT_NE(global, -1);
@@ -509,9 +510,9 @@ TEST_P(ParallelIndexerTest, TestConvertIndex3D) {
         const int global = index->getGlobal(i.xm()),
                   otherXind = (this->pe_xind == 0) ? this->pe_xind : this->pe_xind - 1,
                   otherYind = this->pe_yind + 1;
-        const Ind3D otherInd =
-            i.offset((this->pe_xind == 0) ? -1 : this->xend - this->xstart,
-                     this->ystart - this->yend - 1, 0);
+        Ind3D otherInd = i.offset((this->pe_xind == 0) ? -1 : this->xend - this->xstart,
+                                  this->ystart - this->yend - 1, 0);
+        otherInd.yoffset = 0;
         const int otherGlobal = this->getIndexer<Field3D>(indexers, otherXind, otherYind)
                                     ->getGlobal(otherInd);
         EXPECT_NE(global, -1);
@@ -521,7 +522,7 @@ TEST_P(ParallelIndexerTest, TestConvertIndex3D) {
       int global = index->getGlobal(i);
       int otherGlobal =
           this->getIndexer<Field3D>(indexers, this->pe_xind, this->pe_yind + 1)
-              ->getGlobal(i.ym(this->yend - this->ystart + 1));
+              ->getGlobal(i.yp_no_parallel_shift(-(this->yend - this->ystart + 1)));
       EXPECT_EQ(global, otherGlobal);
 
       if (i.x() == this->xend) {
@@ -529,9 +530,10 @@ TEST_P(ParallelIndexerTest, TestConvertIndex3D) {
                   otherXind = (this->pe_xind == this->nxpe - 1) ? this->pe_xind
                                                                 : this->pe_xind + 1,
                   otherYind = this->pe_yind + 1;
-        const Ind3D otherInd =
+        Ind3D otherInd =
             i.offset((this->pe_xind == this->nxpe - 1) ? 1 : this->xstart - this->xend,
                      this->ystart - this->yend - 1, 0);
+        otherInd.yoffset = 0;
         const int otherGlobal = this->getIndexer<Field3D>(indexers, otherXind, otherYind)
                                     ->getGlobal(otherInd);
         EXPECT_NE(global, -1);
@@ -577,7 +579,7 @@ TEST_P(ParallelIndexerTest, TestConvertIndex2D) {
       int global = index->getGlobal(i);
       int otherGlobal =
           this->getIndexer<Field2D>(indexers, this->pe_xind, this->pe_yind - 1)
-              ->getGlobal(i.yp(this->yend - this->ystart + 1));
+              ->getGlobal(i.yp_no_parallel_shift(this->yend - this->ystart + 1));
       EXPECT_EQ(global, otherGlobal);
 
       if (i.x() == this->xend) {
@@ -612,7 +614,7 @@ TEST_P(ParallelIndexerTest, TestConvertIndex2D) {
       int global = index->getGlobal(i);
       int otherGlobal =
           this->getIndexer<Field2D>(indexers, this->pe_xind, this->pe_yind + 1)
-              ->getGlobal(i.ym(this->yend - this->ystart + 1));
+              ->getGlobal(i.yp_no_parallel_shift(-(this->yend - this->ystart + 1)));
       EXPECT_EQ(global, otherGlobal);
 
       if (i.x() == this->xend) {
diff --git a/tests/unit/include/bout/test_petsc_matrix.cxx b/tests/unit/include/bout/test_petsc_matrix.cxx
index 4f8a4291d6..7200504aa8 100644
--- a/tests/unit/include/bout/test_petsc_matrix.cxx
+++ b/tests/unit/include/bout/test_petsc_matrix.cxx
@@ -55,7 +55,7 @@ class PetscMatrixTest : public FakeMeshFixture {
     if constexpr (std::is_same_v<F, FieldPerp>) {
       indexB = indexA.zp();
     } else {
-      indexB = indexA.yp();
+      indexB = indexA.yp_no_parallel_shift();
     }
     iWU0 = indexB.xm();
     iWU1 = indexB;
@@ -65,9 +65,9 @@ class PetscMatrixTest : public FakeMeshFixture {
       iWD1 = indexB;
       iWD2 = indexB.zp();
     } else {
-      iWD0 = indexB.ym();
+      iWD0 = indexB.yp_no_parallel_shift(-1);
       iWD1 = indexB;
-      iWD2 = indexB.yp();
+      iWD2 = indexB.yp_no_parallel_shift();
     }
     std::unique_ptr<MockTransform> transform =
         bout::utils::make_unique<MockTransform>(*bout::globals::mesh);

From fa1bcb4a8e359468f153378bb0f61620958a3b78 Mon Sep 17 00:00:00 2001
From: David Bold <dschwoerer@users.noreply.github.com>
Date: Mon, 21 Jul 2025 10:14:13 +0200
Subject: [PATCH 13/24] Apply suggestions from code review

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>
---
 include/bout/field3d.hxx | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx
index 39c4f5780a..0dea9b7a7b 100644
--- a/include/bout/field3d.hxx
+++ b/include/bout/field3d.hxx
@@ -356,7 +356,7 @@ public:
   };
 
   BoutReal& operator[](const Ind3D& d) {
-    if (d.yoffset) {
+    if (d.yoffset != 0) {
       if (yoffset == 0) {
         if (hasParallelSlices()) {
           return ynext(d.yoffset)[d];
@@ -364,7 +364,9 @@ public:
 #if CHECK >= 2
         else if (isFci()) { // We probably should assert here that this is field aligned
           throw BoutException(
-              "Tried to access parallel slices, but they are not calculated!");
+	if (isFci()) { // We probably should assert here that this is field aligned
+	  throw BoutException("Tried to access parallel slices, but they are not calculated!");
+	}
         }
 #endif
       } else {

From 79e2579944884c64ef3bb55b2c0d87824b9c0d47 Mon Sep 17 00:00:00 2001
From: dschwoerer <5637662+dschwoerer@users.noreply.github.com>
Date: Mon, 21 Jul 2025 08:15:31 +0000
Subject: [PATCH 14/24] Apply clang-format changes

---
 include/bout/field3d.hxx | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx
index 0dea9b7a7b..05505df689 100644
--- a/include/bout/field3d.hxx
+++ b/include/bout/field3d.hxx
@@ -365,7 +365,8 @@ public:
         else if (isFci()) { // We probably should assert here that this is field aligned
           throw BoutException(
 	if (isFci()) { // We probably should assert here that this is field aligned
-	  throw BoutException("Tried to access parallel slices, but they are not calculated!");
+            throw BoutException(
+                "Tried to access parallel slices, but they are not calculated!");
 	}
         }
 #endif

From b714f7a615d8140fce0671b03637a904b8619192 Mon Sep 17 00:00:00 2001
From: David Bold <dave@ipp.mpg.de>
Date: Mon, 21 Jul 2025 10:17:05 +0200
Subject: [PATCH 15/24] Add header explicitly for clang-tidy

---
 include/bout/field3d.hxx         | 1 +
 src/field/gen_fieldops.py        | 1 +
 src/field/generated_fieldops.cxx | 1 +
 3 files changed, 3 insertions(+)

diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx
index 05505df689..2c93d5ba72 100644
--- a/include/bout/field3d.hxx
+++ b/include/bout/field3d.hxx
@@ -29,6 +29,7 @@ class Field3D;
 #include "bout/array.hxx"
 #include "bout/assert.hxx"
 #include "bout/bout_types.hxx"
+#include "bout/boutexception.hxx"
 #include "bout/field.hxx"
 #include "bout/field2d.hxx"
 #include "bout/fieldperp.hxx"
diff --git a/src/field/gen_fieldops.py b/src/field/gen_fieldops.py
index 29631ff7aa..bf71a37671 100755
--- a/src/field/gen_fieldops.py
+++ b/src/field/gen_fieldops.py
@@ -65,6 +65,7 @@ def smart_open(filename, mode="r"):
 )
 
 header = """// This file is autogenerated - see gen_fieldops.py
+#include <bout/assert.hxx>
 #include <bout/field2d.hxx>
 #include <bout/field3d.hxx>
 #include <bout/globals.hxx>
diff --git a/src/field/generated_fieldops.cxx b/src/field/generated_fieldops.cxx
index c1c450b2a9..eeb32225b4 100644
--- a/src/field/generated_fieldops.cxx
+++ b/src/field/generated_fieldops.cxx
@@ -1,4 +1,5 @@
 // This file is autogenerated - see gen_fieldops.py
+#include <bout/assert.hxx>
 #include <bout/field2d.hxx>
 #include <bout/field3d.hxx>
 #include <bout/globals.hxx>

From 742e6d051e205988702a1f99a30fd1a5a923f68c Mon Sep 17 00:00:00 2001
From: David Bold <dave@ipp.mpg.de>
Date: Mon, 21 Jul 2025 10:24:21 +0200
Subject: [PATCH 16/24] Fixup `suggestions from code-review`

---
 include/bout/field3d.hxx | 8 +++-----
 1 file changed, 3 insertions(+), 5 deletions(-)

diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx
index 2c93d5ba72..e53094e762 100644
--- a/include/bout/field3d.hxx
+++ b/include/bout/field3d.hxx
@@ -363,13 +363,11 @@ public:
           return ynext(d.yoffset)[d];
         }
 #if CHECK >= 2
-        else if (isFci()) { // We probably should assert here that this is field aligned
+        if (isFci()) { // We probably should assert here that this is field aligned
           throw BoutException(
-	if (isFci()) { // We probably should assert here that this is field aligned
-            throw BoutException(
-                "Tried to access parallel slices, but they are not calculated!");
-	}
+              "Tried to access parallel slices, but they are not calculated!");
         }
+      }
 #endif
       } else {
         ASSERT2(d.yoffset == yoffset);

From 58caf56070f1aac88b1eb175abf4278b2ce3b5fa Mon Sep 17 00:00:00 2001
From: dschwoerer <5637662+dschwoerer@users.noreply.github.com>
Date: Mon, 21 Jul 2025 08:26:23 +0000
Subject: [PATCH 17/24] Apply clang-format changes

---
 include/bout/field3d.hxx | 13 ++++++-------
 1 file changed, 6 insertions(+), 7 deletions(-)

diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx
index e53094e762..5c5538b3ac 100644
--- a/include/bout/field3d.hxx
+++ b/include/bout/field3d.hxx
@@ -369,15 +369,14 @@ public:
         }
       }
 #endif
-      } else {
-        ASSERT2(d.yoffset == yoffset);
-      }
+    } else {
+      ASSERT2(d.yoffset == yoffset);
     }
-    return data[d.ind];
-  }
-  const BoutReal& operator[](const Ind3D& d) const {
-    return (*const_cast<Field3D*>(this))[d];
   }
+  return data[d.ind];
+} const BoutReal& operator[](const Ind3D& d) const {
+  return (*const_cast<Field3D*>(this))[d];
+}
 
   BoutReal& operator()(const IndPerp& d, int jy);
   const BoutReal& operator()(const IndPerp& d, int jy) const;

From 6f8a4c9a30db6ea7d965ea7bf428ec93f9b23998 Mon Sep 17 00:00:00 2001
From: David Bold <dave@ipp.mpg.de>
Date: Mon, 21 Jul 2025 10:51:05 +0200
Subject: [PATCH 18/24] Remove extra } from code review

---
 include/bout/field3d.hxx | 1 -
 1 file changed, 1 deletion(-)

diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx
index 5c5538b3ac..8bec0e3646 100644
--- a/include/bout/field3d.hxx
+++ b/include/bout/field3d.hxx
@@ -367,7 +367,6 @@ public:
           throw BoutException(
               "Tried to access parallel slices, but they are not calculated!");
         }
-      }
 #endif
     } else {
       ASSERT2(d.yoffset == yoffset);

From 159c26e31db0866fab34c6f69927ff1bb1225022 Mon Sep 17 00:00:00 2001
From: dschwoerer <5637662+dschwoerer@users.noreply.github.com>
Date: Mon, 21 Jul 2025 08:56:26 +0000
Subject: [PATCH 19/24] Apply clang-format changes

---
 include/bout/field3d.hxx | 13 +++++++------
 1 file changed, 7 insertions(+), 6 deletions(-)

diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx
index 8bec0e3646..a9896af49e 100644
--- a/include/bout/field3d.hxx
+++ b/include/bout/field3d.hxx
@@ -368,14 +368,15 @@ public:
               "Tried to access parallel slices, but they are not calculated!");
         }
 #endif
-    } else {
-      ASSERT2(d.yoffset == yoffset);
+      } else {
+        ASSERT2(d.yoffset == yoffset);
+      }
     }
+    return data[d.ind];
+  }
+  const BoutReal& operator[](const Ind3D& d) const {
+    return (*const_cast<Field3D*>(this))[d];
   }
-  return data[d.ind];
-} const BoutReal& operator[](const Ind3D& d) const {
-  return (*const_cast<Field3D*>(this))[d];
-}
 
   BoutReal& operator()(const IndPerp& d, int jy);
   const BoutReal& operator()(const IndPerp& d, int jy) const;

From 9496211bfc72c24bc3e3aea3eb8cae356889a86e Mon Sep 17 00:00:00 2001
From: David Bold <dave@ipp.mpg.de>
Date: Mon, 21 Jul 2025 12:52:36 +0200
Subject: [PATCH 20/24] Add yoffset to emptyFrom

---
 include/bout/field.hxx     | 3 ++-
 include/bout/field2d.hxx   | 2 +-
 include/bout/field3d.hxx   | 4 +++-
 include/bout/fieldperp.hxx | 2 +-
 src/field/field2d.cxx      | 2 +-
 src/field/field3d.cxx      | 4 ++--
 src/field/fieldperp.cxx    | 2 +-
 7 files changed, 11 insertions(+), 8 deletions(-)

diff --git a/include/bout/field.hxx b/include/bout/field.hxx
index 27835ecbd7..77a24ea96a 100644
--- a/include/bout/field.hxx
+++ b/include/bout/field.hxx
@@ -133,6 +133,7 @@ public:
   virtual void setRegion(const std::string& UNUSED(region_name)) {}
   virtual void resetRegion() {}
   virtual std::optional<size_t> getRegionID() const { return {}; }
+  virtual int getYoffset() const { return 0; }
 
 private:
   /// Labels for the type of coordinate system this field is defined over
@@ -186,7 +187,7 @@ template <typename T>
 inline T emptyFrom(const T& f) {
   static_assert(bout::utils::is_Field_v<T>, "emptyFrom only works on Fields");
   return T(f.getMesh(), f.getLocation(), {f.getDirectionY(), f.getDirectionZ()},
-           f.getRegionID())
+           f.getRegionID(), f.getYoffset())
       .allocate();
 }
 
diff --git a/include/bout/field2d.hxx b/include/bout/field2d.hxx
index 5eab330e8e..f51424f5bd 100644
--- a/include/bout/field2d.hxx
+++ b/include/bout/field2d.hxx
@@ -67,7 +67,7 @@ public:
   Field2D(Mesh* localmesh = nullptr, CELL_LOC location_in = CELL_CENTRE,
           DirectionTypes directions_in = {YDirectionType::Standard,
                                           ZDirectionType::Average},
-          std::optional<size_t> region = {});
+          std::optional<size_t> region = {}, int yoffset = 0);
 
   /*!
    * Copy constructor. After this both fields
diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx
index a9896af49e..1a70fd32f1 100644
--- a/include/bout/field3d.hxx
+++ b/include/bout/field3d.hxx
@@ -168,7 +168,7 @@ public:
   Field3D(Mesh* localmesh = nullptr, CELL_LOC location_in = CELL_CENTRE,
           DirectionTypes directions_in = {YDirectionType::Standard,
                                           ZDirectionType::Standard},
-          std::optional<size_t> regionID = {});
+          std::optional<size_t> regionID = {}, int yoffset = 0);
 
   /*!
    * Copy constructor
@@ -344,6 +344,8 @@ public:
   void setRegion(std::optional<size_t> id) override;
   std::optional<size_t> getRegionID() const override { return regionID; };
 
+  int getYoffset() const override { return yoffset; };
+
   /// Return a Region<Ind2D> reference to use to iterate over the x- and
   /// y-indices of this field
   const Region<Ind2D>& getRegion2D(REGION region) const;
diff --git a/include/bout/fieldperp.hxx b/include/bout/fieldperp.hxx
index b50eef1991..8cb649a29d 100644
--- a/include/bout/fieldperp.hxx
+++ b/include/bout/fieldperp.hxx
@@ -59,7 +59,7 @@ public:
             int yindex_in = -1,
             DirectionTypes directions_in = {YDirectionType::Standard,
                                             ZDirectionType::Standard},
-            std::optional<size_t> regionID = {});
+            std::optional<size_t> regionID = {}, int yoffset = 0);
 
   /*!
    * Copy constructor. After this the data
diff --git a/src/field/field2d.cxx b/src/field/field2d.cxx
index 4591d228f7..e6674eb7e9 100644
--- a/src/field/field2d.cxx
+++ b/src/field/field2d.cxx
@@ -49,7 +49,7 @@
 #include <bout/assert.hxx>
 
 Field2D::Field2D(Mesh* localmesh, CELL_LOC location_in, DirectionTypes directions_in,
-                 std::optional<size_t> UNUSED(regionID))
+                 std::optional<size_t> UNUSED(regionID), int UNUSED(yoffset))
     : Field(localmesh, location_in, directions_in) {
 
   if (fieldmesh) {
diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx
index fb76244016..658f3fd8ce 100644
--- a/src/field/field3d.cxx
+++ b/src/field/field3d.cxx
@@ -49,8 +49,8 @@
 
 /// Constructor
 Field3D::Field3D(Mesh* localmesh, CELL_LOC location_in, DirectionTypes directions_in,
-                 std::optional<size_t> regionID)
-    : Field(localmesh, location_in, directions_in), regionID{regionID} {
+                 std::optional<size_t> regionID, const int yoffset)
+  : Field(localmesh, location_in, directions_in), regionID{regionID}, yoffset{yoffset} {
 #if BOUT_USE_TRACK
   name = "<F3D>";
 #endif
diff --git a/src/field/fieldperp.cxx b/src/field/fieldperp.cxx
index c7a196a57a..0ab3343d7b 100644
--- a/src/field/fieldperp.cxx
+++ b/src/field/fieldperp.cxx
@@ -35,7 +35,7 @@
 #include <bout/utils.hxx>
 
 FieldPerp::FieldPerp(Mesh* localmesh, CELL_LOC location_in, int yindex_in,
-                     DirectionTypes directions, std::optional<size_t> UNUSED(regionID))
+                     DirectionTypes directions, std::optional<size_t> UNUSED(regionID), int UNUSED(yoffset))
     : Field(localmesh, location_in, directions), yindex(yindex_in) {
   if (fieldmesh) {
     nx = fieldmesh->LocalNx;

From 9b052eb5415d3047a68fc2035fcdc7ca486c1108 Mon Sep 17 00:00:00 2001
From: David Bold <dave@ipp.mpg.de>
Date: Mon, 21 Jul 2025 12:53:23 +0200
Subject: [PATCH 21/24] Use correct yoffset for accessing parallel fields

---
 src/mesh/coordinates_accessor.cxx | 11 ++++++++---
 1 file changed, 8 insertions(+), 3 deletions(-)

diff --git a/src/mesh/coordinates_accessor.cxx b/src/mesh/coordinates_accessor.cxx
index aff546c2b0..d77445fb95 100644
--- a/src/mesh/coordinates_accessor.cxx
+++ b/src/mesh/coordinates_accessor.cxx
@@ -55,9 +55,14 @@ CoordinatesAccessor::CoordinatesAccessor(const Coordinates* coords) {
     COPY_STRIPE(J);
 
     data[stripe_size * ind.ind + static_cast<int>(Offset::B)] = coords->Bxy[ind];
-    data[stripe_size * ind.ind + static_cast<int>(Offset::Byup)] = coords->Bxy.yup()[ind];
-    data[stripe_size * ind.ind + static_cast<int>(Offset::Bydown)] =
-        coords->Bxy.ydown()[ind];
+    {
+      auto indc = ind;
+      indc.yoffset = 1;
+      data[stripe_size * ind.ind + static_cast<int>(Offset::Byup)] = coords->Bxy.yup()[indc];
+      indc.yoffset = -1;
+      data[stripe_size * ind.ind + static_cast<int>(Offset::Bydown)] =
+        coords->Bxy.ydown()[indc];
+    }
 
     COPY_STRIPE(G1, G3);
     COPY_STRIPE(g11, g12, g13, g22, g23, g33);

From 085948c2e8f5b0156b5bfb13544ebd386b7c93f7 Mon Sep 17 00:00:00 2001
From: David Bold <dave@ipp.mpg.de>
Date: Mon, 21 Jul 2025 12:53:40 +0200
Subject: [PATCH 22/24] Fix mock fields to have realistic yoffset

---
 tests/unit/fake_mesh_fixture.hxx              |  7 ++++++-
 tests/unit/mesh/test_coordinates_accessor.cxx | 21 ++++++++++++++++---
 2 files changed, 24 insertions(+), 4 deletions(-)

diff --git a/tests/unit/fake_mesh_fixture.hxx b/tests/unit/fake_mesh_fixture.hxx
index 2758dbe416..06f865b66c 100644
--- a/tests/unit/fake_mesh_fixture.hxx
+++ b/tests/unit/fake_mesh_fixture.hxx
@@ -58,8 +58,13 @@ public:
     test_coords->d1_dx = test_coords->d1_dy = 0.2;
     test_coords->d1_dz = 0.0;
 #if BOUT_USE_METRIC_3D
+    auto Bup = test_coords->Bxy;
+    Bup.yoffset = 1;
+    auto Bdown = test_coords->Bxy;
+    Bdown.yoffset = -1;
     test_coords->Bxy.splitParallelSlices();
-    test_coords->Bxy.yup() = test_coords->Bxy.ydown() = test_coords->Bxy;
+    test_coords->Bxy.yup() = Bup;
+    test_coords->Bxy.ydown() = Bdown;
 #endif
 
     // No call to Coordinates::geometry() needed here
diff --git a/tests/unit/mesh/test_coordinates_accessor.cxx b/tests/unit/mesh/test_coordinates_accessor.cxx
index 16985981bb..1b96bb731e 100644
--- a/tests/unit/mesh/test_coordinates_accessor.cxx
+++ b/tests/unit/mesh/test_coordinates_accessor.cxx
@@ -81,8 +81,13 @@ TEST_F(CoordinatesAccessorTest, ClearBoth) {
   coords.non_uniform = true;
   coords.d1_dx = coords.d1_dy = coords.d1_dz = 0.1;
 #if BOUT_USE_METRIC_3D
+  auto Bup = coords.Bxy;
+  Bup.yoffset = 1;
+  auto Bdown = coords.Bxy;
+  Bdown.yoffset = -1;
   coords.Bxy.splitParallelSlices();
-  coords.Bxy.yup() = coords.Bxy.ydown() = coords.Bxy;
+  coords.Bxy.yup() = Bup;
+  coords.Bxy.ydown() = Bdown;
 #endif
 
   CoordinatesAccessor acc(mesh->getCoordinates());
@@ -121,8 +126,13 @@ TEST_F(CoordinatesAccessorTest, ClearOneTwo) {
   coords.non_uniform = true;
   coords.d1_dx = coords.d1_dy = coords.d1_dz = 0.1;
 #if BOUT_USE_METRIC_3D
+  auto Bup = coords.Bxy;
+  Bup.yoffset = 1;
+  auto Bdown = coords.Bxy;
+  Bdown.yoffset = -1;
   coords.Bxy.splitParallelSlices();
-  coords.Bxy.yup() = coords.Bxy.ydown() = coords.Bxy;
+  coords.Bxy.yup() = Bup;
+  coords.Bxy.ydown() = Bdown;
 #endif
 
   CoordinatesAccessor acc(mesh->getCoordinates());
@@ -163,8 +173,13 @@ TEST_F(CoordinatesAccessorTest, ClearTwoOneNone) {
   coords.non_uniform = true;
   coords.d1_dx = coords.d1_dy = coords.d1_dz = 0.1;
 #if BOUT_USE_METRIC_3D
+  auto Bup = coords.Bxy;
+  Bup.yoffset = 1;
+  auto Bdown = coords.Bxy;
+  Bdown.yoffset = -1;
   coords.Bxy.splitParallelSlices();
-  coords.Bxy.yup() = coords.Bxy.ydown() = coords.Bxy;
+  coords.Bxy.yup() = Bup;
+  coords.Bxy.ydown() = Bdown;
 #endif
 
   CoordinatesAccessor acc(mesh->getCoordinates());

From 288eb775ea16151d8530a0ed621b70f5f070043c Mon Sep 17 00:00:00 2001
From: dschwoerer <5637662+dschwoerer@users.noreply.github.com>
Date: Mon, 21 Jul 2025 10:55:34 +0000
Subject: [PATCH 23/24] Apply clang-format changes

---
 src/field/field3d.cxx             | 2 +-
 src/field/fieldperp.cxx           | 3 ++-
 src/mesh/coordinates_accessor.cxx | 5 +++--
 3 files changed, 6 insertions(+), 4 deletions(-)

diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx
index 658f3fd8ce..b66aa9c58d 100644
--- a/src/field/field3d.cxx
+++ b/src/field/field3d.cxx
@@ -50,7 +50,7 @@
 /// Constructor
 Field3D::Field3D(Mesh* localmesh, CELL_LOC location_in, DirectionTypes directions_in,
                  std::optional<size_t> regionID, const int yoffset)
-  : Field(localmesh, location_in, directions_in), regionID{regionID}, yoffset{yoffset} {
+    : Field(localmesh, location_in, directions_in), regionID{regionID}, yoffset{yoffset} {
 #if BOUT_USE_TRACK
   name = "<F3D>";
 #endif
diff --git a/src/field/fieldperp.cxx b/src/field/fieldperp.cxx
index 0ab3343d7b..a3458543c9 100644
--- a/src/field/fieldperp.cxx
+++ b/src/field/fieldperp.cxx
@@ -35,7 +35,8 @@
 #include <bout/utils.hxx>
 
 FieldPerp::FieldPerp(Mesh* localmesh, CELL_LOC location_in, int yindex_in,
-                     DirectionTypes directions, std::optional<size_t> UNUSED(regionID), int UNUSED(yoffset))
+                     DirectionTypes directions, std::optional<size_t> UNUSED(regionID),
+                     int UNUSED(yoffset))
     : Field(localmesh, location_in, directions), yindex(yindex_in) {
   if (fieldmesh) {
     nx = fieldmesh->LocalNx;
diff --git a/src/mesh/coordinates_accessor.cxx b/src/mesh/coordinates_accessor.cxx
index d77445fb95..aca0c319e8 100644
--- a/src/mesh/coordinates_accessor.cxx
+++ b/src/mesh/coordinates_accessor.cxx
@@ -58,10 +58,11 @@ CoordinatesAccessor::CoordinatesAccessor(const Coordinates* coords) {
     {
       auto indc = ind;
       indc.yoffset = 1;
-      data[stripe_size * ind.ind + static_cast<int>(Offset::Byup)] = coords->Bxy.yup()[indc];
+      data[stripe_size * ind.ind + static_cast<int>(Offset::Byup)] =
+          coords->Bxy.yup()[indc];
       indc.yoffset = -1;
       data[stripe_size * ind.ind + static_cast<int>(Offset::Bydown)] =
-        coords->Bxy.ydown()[indc];
+          coords->Bxy.ydown()[indc];
     }
 
     COPY_STRIPE(G1, G3);

From f3fd2420930446a44402052a66fa389662f0966a Mon Sep 17 00:00:00 2001
From: David Bold <dave@ipp.mpg.de>
Date: Mon, 21 Jul 2025 10:54:45 +0200
Subject: [PATCH 24/24] Disable hypre3d test

---
 tests/integrated/test-laplace-hypre3d/CMakeLists.txt | 1 +
 1 file changed, 1 insertion(+)

diff --git a/tests/integrated/test-laplace-hypre3d/CMakeLists.txt b/tests/integrated/test-laplace-hypre3d/CMakeLists.txt
index 2645c18c67..b09b416feb 100644
--- a/tests/integrated/test-laplace-hypre3d/CMakeLists.txt
+++ b/tests/integrated/test-laplace-hypre3d/CMakeLists.txt
@@ -7,4 +7,5 @@ bout_add_integrated_test(test-laplace-hypre3d
     data_slab_sol/BOUT.inp
   USE_RUNTEST
   REQUIRES BOUT_HAS_HYPRE
+  REQUIRES BOUT_RUN_ALL_TESTS
   )