From a333188125e25094e00f9f11dc8377d862c56340 Mon Sep 17 00:00:00 2001
From: David Bold <dave@ipp.mpg.de>
Date: Thu, 2 Jul 2026 12:55:06 +0200
Subject: [PATCH 01/11] Port Div_par_fvv_heating from hermes-3

---
 include/bout/fv_ops_impl.hxx | 261 +++++++++++++++++++++++++++++++++++
 1 file changed, 261 insertions(+)

diff --git a/include/bout/fv_ops_impl.hxx b/include/bout/fv_ops_impl.hxx
index 3755467027..f99ec961e3 100644
--- a/include/bout/fv_ops_impl.hxx
+++ b/include/bout/fv_ops_impl.hxx
@@ -981,5 +981,266 @@ Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in,
   }
   return fromFieldAligned(result, "RGN_NOBNDRY");
 }
+
+// Calculates viscous heating due to numerical momentum fluxes
+// and flow of kinetic energy (in flow_ylow)
+template <typename CellEdges = MC>
+Field3D Div_par_fvv_heating(const Field3D& f_in, const Field3D& v_in,
+                            const Field3D& wave_speed_in, Field3D& flow_ylow,
+                            bool fixflux = true) {
+
+  ASSERT1(areFieldsCompatible(f_in, v_in));
+  ASSERT1(areFieldsCompatible(f_in, wave_speed_in));
+
+  Mesh* mesh = f_in.getMesh();
+  Coordinates* coord = f_in.getCoordinates();
+  CellEdges cellboundary;
+
+  if (f_in.isFci()) {
+    // FCI version, using yup/down fields
+    ASSERT1(f_in.hasParallelSlices());
+    ASSERT1(v_in.hasParallelSlices());
+
+    const auto B = coord->Bxy;
+    const auto B_up = coord->Bxy.yup();
+    const auto B_down = coord->Bxy.ydown();
+
+    const auto& f_up = f_in.yup();
+    const auto& f_down = f_in.ydown();
+
+    const auto& v_up = v_in.yup();
+    const auto& v_down = v_in.ydown();
+
+    const auto g_22 = coord->g_22;
+    const auto dy = coord->dy;
+
+    Field3D result{emptyFrom(f_in)};
+    flow_ylow = zeroFrom(f_in);
+
+    BOUT_FOR(i, f_in.getRegion("RGN_NOBNDRY")) {
+      const auto iyp = i.yp();
+      const auto iym = i.ym();
+
+      //Maximum local wave speed
+      const BoutReal amax =
+          BOUTMAX(wave_speed_in[i], fabs(v_in[i]), fabs(v_up[iyp]), fabs(v_down[iym]));
+
+      result[i] =
+          B[i]
+          * ((f_up[iyp] * v_up[iyp] * v_up[iyp] / B_up[iyp])
+             - (f_down[iym] * v_down[iym] * v_down[iym] / B_down[iym])
+             // Penalty terms. This implementation is very dissipative.
+             // Note: This version adds a viscosity that damps gradients of velocity
+             + amax * (f_in[i] + f_up[iyp]) * (v_in[i] - v_up[iyp]) / (B[i] + B_up[iyp])
+             + amax * (f_in[i] + f_down[iym]) * (v_in[i] - v_down[iym])
+                   / (B[i] + B_down[iym]))
+          / (2 * dy[i] * sqrt(g_22[i]));
+
+#if CHECK > 0
+      if (!std::isfinite(result[i])) {
+        throw BoutException("Non-finite value in Div_par_fvv at {}\n"
+                            "fup {} vup {} fdown {} vdown {} amax {}\n",
+                            "B {} Bup {} Bdown {} dy {} sqrt(g_22} {}", i, f_up[i],
+                            v_up[i], f_down[i], v_down[i], amax, B[i], B_up[i], B_down[i],
+                            dy[i], sqrt(g_22[i]));
+      }
+#endif
+    }
+    return result;
+  }
+
+  /// Ensure that f, v and wave_speed are field aligned
+  Field3D f = toFieldAligned(f_in, "RGN_NOX");
+  Field3D v = toFieldAligned(v_in, "RGN_NOX");
+  Field3D wave_speed = toFieldAligned(wave_speed_in, "RGN_NOX");
+
+  // result and flow_ylow are field-aligned.
+  // Will be converted to non-aligned before return.
+  Field3D result{zeroFrom(f)};
+  flow_ylow = zeroFrom(f);
+
+  // Only need one guard cell, so no need to communicate fluxes
+  // Instead calculate in guard cells to preserve fluxes
+  int ys = mesh->ystart - 1;
+  int ye = mesh->yend + 1;
+
+  for (int i = mesh->xstart; i <= mesh->xend; i++) {
+
+    if (!mesh->firstY(i) || mesh->periodicY(i)) {
+      // Calculate in guard cell to get fluxes consistent between processors
+      ys = mesh->ystart - 1;
+    } else {
+      // Don't include the boundary cell. Note that this implies special
+      // handling of boundaries later
+      ys = mesh->ystart;
+    }
+
+    if (!mesh->lastY(i) || mesh->periodicY(i)) {
+      // Calculate in guard cells
+      ye = mesh->yend + 1;
+    } else {
+      // Not in boundary cells
+      ye = mesh->yend;
+    }
+
+    for (int j = ys; j <= ye; j++) {
+      for (int k = 0; k < mesh->LocalNz; k++) {
+        // Pre-calculate factors which multiply fluxes
+        // Note: In 3D metric geometries these quantities can depend on (i,j,k),
+        // so calculate inside the k loop.
+
+        // For right cell boundaries
+        BoutReal common_factor =
+            (coord->J(i, j, k) + coord->J(i, j + 1, k))
+            / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j + 1, k)));
+
+        const BoutReal flux_factor_rc =
+            common_factor / (coord->dy(i, j, k) * coord->J(i, j, k));
+        const BoutReal area_rp =
+            common_factor * coord->dx(i, j + 1, k) * coord->dz(i, j + 1, k);
+
+        // For left cell boundaries
+        common_factor = (coord->J(i, j, k) + coord->J(i, j - 1, k))
+                        / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j - 1, k)));
+
+        const BoutReal flux_factor_lc =
+            common_factor / (coord->dy(i, j, k) * coord->J(i, j, k));
+        const BoutReal area_lc = common_factor * coord->dx(i, j, k) * coord->dz(i, j, k);
+
+        ////////////////////////////////////////////
+        // Reconstruct f at the cell faces
+        // This calculates s.R and s.L for the Right and Left
+        // face values on this cell
+
+        // Reconstruct f at the cell faces
+        Stencil1D s;
+        s.c = f(i, j, k);
+        s.m = f(i, j - 1, k);
+        s.p = f(i, j + 1, k);
+
+        cellboundary(s); // Calculate s.R and s.L
+
+        // Reconstruct v at the cell faces
+        Stencil1D sv;
+        sv.c = v(i, j, k);
+        sv.m = v(i, j - 1, k);
+        sv.p = v(i, j + 1, k);
+
+        cellboundary(sv);
+
+        ////////////////////////////////////////////
+        // Right boundary
+
+        // Calculate velocity at right boundary (y+1/2)
+        BoutReal v_mid = 0.5 * (sv.c + sv.p);
+        // And mid-point density at right boundary
+        BoutReal n_mid = 0.5 * (s.c + s.p);
+
+        if (mesh->lastY(i) && (j == mesh->yend) && !mesh->periodicY(i)) {
+          // Last point in domain
+
+          // Expected loss of kinetic energy into boundary
+          // This is used in the sheath boundary condition to calculate
+          // energy losses.
+          const BoutReal expected_ke = 0.5 * n_mid * v_mid * v_mid * v_mid;
+
+          BoutReal flux_mom;
+          if (fixflux) {
+            // Mid-point consistent with boundary conditions
+            // but kinetic energy loss will not match expected
+            // -> Adjust energy balance in pressure equation
+            flux_mom = n_mid * v_mid * v_mid;
+          } else {
+            flux_mom = s.R * sv.R * sv.R
+                       + BOUTMAX(wave_speed(i, j, k), fabs(sv.c), fabs(sv.p))
+                             * (s.R * sv.R - n_mid * v_mid);
+          }
+
+          // Assume that particle flux is fixed to boundary value
+          const BoutReal flux_part = n_mid * v_mid;
+
+          // d/dt(1/2 m n v^2) = v * d/dt(mnv) - 1/2 m v^2 * dn/dt
+          const BoutReal actual_ke = sv.c * flux_mom - 0.5 * sv.c * sv.c * flux_part;
+
+          // Note: If the actual loss was higher than expected, then
+          //       plasma heating is needed to compensate
+          result(i, j, k) += (actual_ke - expected_ke) * flux_factor_rc;
+
+          // Final flow through boundary is the expected value
+          flow_ylow(i, j + 1, k) += expected_ke * area_rp; //expected_ke * area_rp;
+
+        } else {
+          // Maximum wave speed in the two cells
+          const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j + 1, k),
+                                        fabs(sv.c), fabs(sv.p));
+
+          // Viscous heating due to relaxation of velocity towards midpoint
+          result(i, j, k) +=
+              (amax + 0.5 * sv.R) * s.R * (sv.c - sv.p) * (sv.R - v_mid) * flux_factor_rc;
+
+          // Kinetic energy flow into next cell.
+          // Note: Different from flow out of this cell; the difference
+          //       is in the viscous heating.
+          const BoutReal flux_part = s.R * 0.5 * (sv.R + amax);
+          const BoutReal flux_mom = flux_part * sv.R;
+
+          flow_ylow(i, j + 1, k) +=
+              (sv.p * flux_mom - 0.5 * SQ(sv.p) * flux_part) * area_rp;
+        }
+
+        ////////////////////////////////////////////
+        // Calculate at left boundary
+
+        v_mid = 0.5 * (sv.c + sv.m);
+        n_mid = 0.5 * (s.c + s.m);
+
+        // Expected KE loss. Note minus sign because negative v into boundary
+        const BoutReal expected_ke = -0.5 * n_mid * v_mid * v_mid * v_mid;
+
+        if (mesh->firstY(i) && (j == mesh->ystart) && !mesh->periodicY(i)) {
+          // First point in domain
+          BoutReal flux_mom;
+          if (fixflux) {
+            // Use mid-point to be consistent with boundary conditions
+            flux_mom = n_mid * v_mid * v_mid;
+          } else {
+            // Add flux due to difference in boundary values
+            flux_mom = s.L * sv.L * sv.L
+                       - BOUTMAX(wave_speed(i, j, k), fabs(sv.c), fabs(sv.m))
+                             * (s.L * sv.L - n_mid * v_mid);
+          }
+
+          // Assume that density flux is fixed to boundary value
+          const BoutReal flux_part = n_mid * v_mid;
+
+          // d/dt(1/2 m n v^2) = v * d/dt(mnv) - 1/2 m v^2 * dn/dt
+          const BoutReal actual_ke = -sv.c * flux_mom + 0.5 * sv.c * sv.c * flux_part;
+
+          result(i, j, k) += (actual_ke - expected_ke) * flux_factor_lc;
+
+          flow_ylow(i, j, k) -= expected_ke * area_lc;
+        } else {
+          // Maximum wave speed in the two cells
+          const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j - 1, k),
+                                        fabs(sv.c), fabs(sv.m));
+
+          // Viscous heating due to relaxation
+          result(i, j, k) +=
+              (amax - 0.5 * sv.L) * s.L * (sv.c - sv.m) * (sv.L - v_mid) * flux_factor_lc;
+
+          // Kinetic energy flow into this cell.
+          // Note: Different from flow out of left cell; the difference
+          //       is in the viscous heating.
+          const BoutReal flux_part = s.L * 0.5 * (sv.L - amax);
+          const BoutReal flux_mom = flux_part * sv.L;
+
+          flow_ylow(i, j, k) += (sv.c * flux_mom - 0.5 * SQ(sv.c) * flux_part) * area_lc;
+        }
+      }
+    }
+  }
+  flow_ylow = fromFieldAligned(flow_ylow, "RGN_NOBNDRY");
+  return fromFieldAligned(result, "RGN_NOBNDRY");
+}
 } // namespace FV
 #endif // BOUT_FV_OPS_H

From 4b8f73fb0e871a2e897cc59ba384511ee8ac307c Mon Sep 17 00:00:00 2001
From: David Bold <dave@ipp.mpg.de>
Date: Thu, 2 Jul 2026 12:56:02 +0200
Subject: [PATCH 02/11] Port Div_a_Grad_perp_limit from hermes-3

---
 include/bout/fv_ops_impl.hxx | 252 +++++++++++++++++++++++++++++++++++
 1 file changed, 252 insertions(+)

diff --git a/include/bout/fv_ops_impl.hxx b/include/bout/fv_ops_impl.hxx
index f99ec961e3..ac9041a19b 100644
--- a/include/bout/fv_ops_impl.hxx
+++ b/include/bout/fv_ops_impl.hxx
@@ -1242,5 +1242,257 @@ Field3D Div_par_fvv_heating(const Field3D& f_in, const Field3D& v_in,
   flow_ylow = fromFieldAligned(flow_ylow, "RGN_NOBNDRY");
   return fromFieldAligned(result, "RGN_NOBNDRY");
 }
+
+/// Div ( a g Grad_perp(f) )  -- Perpendicular gradient-driven advection
+///
+/// This version uses a slope limiter to calculate cell edge values of g in X,
+/// the advects the upwind cell edge.
+///
+/// 1st order upwinding is used in Y.
+template <typename CellEdges = MC>
+const Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g,
+                                    const Field3D& f) {
+  ASSERT2(a.getLocation() == f.getLocation());
+
+  Mesh* mesh = a.getMesh();
+
+  // Requires at least 2 communication guard cells in X, 1 in Y
+  ASSERT1(mesh->xstart >= 2);
+  ASSERT1(mesh->ystart >= 1);
+
+  CellEdges cellboundary;
+
+  Field3D result{zeroFrom(f)};
+
+  Coordinates* coord = f.getCoordinates();
+
+  // Flux in x
+
+  for (int i = mesh->xstart - 1; i <= mesh->xend; i++) {
+    for (int j = mesh->ystart; j <= mesh->yend; j++) {
+      for (int k = 0; k < mesh->LocalNz; k++) {
+        // Calculate flux from i to i+1
+
+        const BoutReal gradient = f(i + 1, j, k) - f(i, j, k);
+
+        // Mid-point average boundary value of 'a'
+        const BoutReal aedge = 0.5 * (a(i + 1, j, k) + a(i, j, k));
+        BoutReal gedge;
+        if (((i == mesh->xstart - 1) and mesh->firstX())
+            or ((i == mesh->xend) and mesh->lastX())) {
+          // Mid-point average boundary value of 'g'
+          gedge = 0.5 * (g(i + 1, j, k) + g(i, j, k));
+        } else if (gradient > 0) {
+          // Flux is from (i+1) to (i)
+          // Reconstruct `g` at left of (i+1, j, k)
+
+          Stencil1D sg;
+          sg.m = g(i, j, k);
+          sg.c = g(i + 1, j, k);
+          sg.p = g(i + 2, j, k);
+          cellboundary(sg); // Calculate sg.R and sg.L
+
+          gedge = sg.L;
+        } else {
+          // Flux is from (i) to (i+1)
+          // Reconstruct `g` at right of (i, j, k)
+
+          Stencil1D sg;
+          sg.m = g(i - 1, j, k);
+          sg.c = g(i, j, k);
+          sg.p = g(i + 1, j, k);
+          cellboundary(sg); // Calculate sg.R and sg.L
+
+          gedge = sg.R;
+        }
+
+        // Flux across cell edge
+        const BoutReal fout = gradient * aedge * gedge
+                              * (coord->J(i, j, k) * coord->g11(i, j, k)
+                                 + coord->J(i + 1, j, k) * coord->g11(i + 1, j, k))
+                              / (coord->dx(i, j, k) + coord->dx(i + 1, j, k));
+
+        result(i, j, k) += fout / (coord->dx(i, j, k) * coord->J(i, j, k));
+        result(i + 1, j, k) -= fout / (coord->dx(i + 1, j, k) * coord->J(i + 1, j, k));
+      }
+    }
+  }
+
+  // Y and Z fluxes require Y derivatives
+
+  // Fields containing values along the magnetic field
+  Field3D fup(mesh), fdown(mesh);
+  Field3D aup(mesh), adown(mesh);
+  Field3D gup(mesh), gdown(mesh);
+
+  // Values on this y slice (centre).
+  // This is needed because toFieldAligned may modify the field
+  Field3D ac = a;
+  Field3D gc = g;
+  Field3D fc = f;
+
+  // Result of the Y and Z fluxes
+  Field3D yzresult(mesh);
+  yzresult.allocate();
+
+  if (f.hasParallelSlices() && a.hasParallelSlices() && g.hasParallelSlices()) {
+    // All inputs have yup and ydown
+
+    fup = f.yup();
+    fdown = f.ydown();
+
+    aup = a.yup();
+    adown = a.ydown();
+
+    gup = g.yup();
+    gdown = g.ydown();
+  } else {
+    // At least one input doesn't have yup/ydown fields.
+    // Need to shift to/from field aligned coordinates
+
+    fup = fdown = fc = toFieldAligned(f);
+    aup = adown = ac = toFieldAligned(a);
+    gup = gdown = gc = toFieldAligned(g);
+    yzresult.setDirectionY(YDirectionType::Aligned);
+  }
+
+  // Y flux
+
+  for (int i = mesh->xstart; i <= mesh->xend; i++) {
+    for (int j = mesh->ystart; j <= mesh->yend; j++) {
+#if BOUT_USE_METRIC_3D
+      for (int k = 0; k < mesh->LocalNz; k++)
+#else
+      int k = 0;
+#endif
+      {
+        BoutReal coef_u =
+            0.5
+            * (coord->g_23(i, j, k) / SQ(coord->J(i, j, k) * coord->Bxy(i, j, k))
+               + coord->g_23(i, j + 1, k)
+                     / SQ(coord->J(i, j + 1, k) * coord->Bxy(i, j + 1, k)));
+
+        BoutReal coef_d =
+            0.5
+            * (coord->g_23(i, j, k) / SQ(coord->J(i, j, k) * coord->Bxy(i, j, k))
+               + coord->g_23(i, j - 1, k)
+                     / SQ(coord->J(i, j - 1, k) * coord->Bxy(i, j - 1, k)));
+
+#if not BOUT_USE_METRIC_3D
+        for (int k = 0; k < mesh->LocalNz; k++)
+#endif
+        {
+          // Calculate flux between j and j+1
+          int kp = (k + 1) % mesh->LocalNz;
+          int km = (k - 1 + mesh->LocalNz) % mesh->LocalNz;
+
+          // Calculate Z derivative at y boundary
+          BoutReal dfdz =
+              0.25 * (fc(i, j, kp) - fc(i, j, km) + fup(i, j + 1, kp) - fup(i, j + 1, km))
+              / coord->dz(i, j, k);
+
+          // Y derivative
+          BoutReal dfdy = 2. * (fup(i, j + 1, k) - fc(i, j, k))
+                          / (coord->dy(i, j + 1, k) + coord->dy(i, j, k));
+
+          BoutReal aedge = 0.5 * (ac(i, j, k) + aup(i, j + 1, k));
+          BoutReal gedge;
+          if ((j == mesh->yend) and mesh->lastY(i)) {
+            // Midpoint boundary value
+            gedge = 0.5 * (gc(i, j, k) + gup(i, j + 1, k));
+          } else if (dfdy > 0) {
+            // Flux from (j+1) to (j)
+            gedge = gup(i, j + 1, k);
+          } else {
+            // Flux from (j) to (j+1)
+            gedge = gc(i, j, k);
+          }
+
+          BoutReal fout = aedge * gedge * 0.5
+                          * (coord->J(i, j, k) * coord->g23(i, j, k)
+                             + coord->J(i, j + 1, k) * coord->g23(i, j + 1, k))
+                          * (dfdz - coef_u * dfdy);
+
+          yzresult(i, j, k) = fout / (coord->dy(i, j, k) * coord->J(i, j, k));
+
+          // Calculate flux between j and j-1
+          dfdz =
+              0.25
+              * (fc(i, j, kp) - fc(i, j, km) + fdown(i, j - 1, kp) - fdown(i, j - 1, km))
+              / coord->dz(i, j, k);
+
+          dfdy = 2. * (fc(i, j, k) - fdown(i, j - 1, k))
+                 / (coord->dy(i, j, k) + coord->dy(i, j - 1, k));
+
+          aedge = 0.5 * (ac(i, j, k) + adown(i, j - 1, k));
+          if ((j == mesh->ystart) and mesh->firstY(i)) {
+            gedge = 0.5 * (gc(i, j, k) + gdown(i, j - 1, k));
+          } else if (dfdy > 0) {
+            gedge = gc(i, j, k);
+          } else {
+            gedge = gdown(i, j - 1, k);
+          }
+
+          fout = aedge * gedge * 0.5
+                 * (coord->J(i, j, k) * coord->g23(i, j, k)
+                    + coord->J(i, j - 1, k) * coord->g23(i, j - 1, k))
+                 * (dfdz - coef_d * dfdy);
+
+          yzresult(i, j, k) -= fout / (coord->dy(i, j, k) * coord->J(i, j, k));
+        }
+      }
+    }
+  }
+
+  // Z flux
+  // Easier since all metrics constant in Z
+
+  for (int i = mesh->xstart; i <= mesh->xend; i++) {
+    for (int j = mesh->ystart; j <= mesh->yend; j++) {
+#if BOUT_USE_METRIC_3D
+      for (int k = 0; k < mesh->LocalNz; k++)
+#else
+      int k = 0;
+#endif
+      {
+        // Coefficient in front of df/dy term
+        BoutReal coef =
+            coord->g_23(i, j, k)
+            / (coord->dy(i, j + 1, k) + 2. * coord->dy(i, j, k) + coord->dy(i, j - 1, k))
+            / SQ(coord->J(i, j, k) * coord->Bxy(i, j, k));
+#if not BOUT_USE_METRIC_3D
+        for (int k = 0; k < mesh->LocalNz; k++)
+#endif
+        {
+          // Calculate flux between k and k+1
+          int kp = (k + 1) % mesh->LocalNz;
+
+          BoutReal gradient =
+              // df/dz
+              (fc(i, j, kp) - fc(i, j, k)) / coord->dz(i, j, k)
+
+              // - g_yz * df/dy / SQ(J*B)
+              - coef
+                    * (fup(i, j + 1, k) + fup(i, j + 1, kp) - fdown(i, j - 1, k)
+                       - fdown(i, j - 1, kp));
+
+          BoutReal fout = gradient * 0.5 * (ac(i, j, kp) + ac(i, j, k))
+                          * ((gradient > 0) ? gc(i, j, kp) : gc(i, j, k));
+
+          yzresult(i, j, k) += fout / coord->dz(i, j, k);
+          yzresult(i, j, kp) -= fout / coord->dz(i, j, kp);
+        }
+      }
+    }
+  }
+  // Check if we need to transform back
+  if (f.hasParallelSlices() && a.hasParallelSlices()) {
+    result += yzresult;
+  } else {
+    result += fromFieldAligned(yzresult);
+  }
+  return result;
+}
+
 } // namespace FV
 #endif // BOUT_FV_OPS_H

From d42c5f5f9bc9c2f06f4b8f1874bb7b4259e52349 Mon Sep 17 00:00:00 2001
From: David Bold <dave@ipp.mpg.de>
Date: Thu, 2 Jul 2026 12:58:49 +0200
Subject: [PATCH 03/11] Declare communicateFluxes only once

---
 include/bout/fv_ops.hxx      | 17 +++++++++++++++++
 include/bout/fv_ops_impl.hxx | 12 +++---------
 2 files changed, 20 insertions(+), 9 deletions(-)

diff --git a/include/bout/fv_ops.hxx b/include/bout/fv_ops.hxx
index b2b0e49120..c7a6f6774c 100644
--- a/include/bout/fv_ops.hxx
+++ b/include/bout/fv_ops.hxx
@@ -165,5 +165,22 @@ Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in,
 template <typename CellEdges = MC>
 Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in,
                     const Field3D& wave_speed_in, bool fixflux = true);
+
+/// Calculates viscous heating due to numerical momentum fluxes
+/// and flow of kinetic energy (in flow_ylow)
+template <typename CellEdges = MC>
+Field3D Div_par_fvv_heating(const Field3D& f_in, const Field3D& v_in,
+                            const Field3D& wave_speed_in, Field3D& flow_ylow,
+                            bool fixflux = true);
+
+/// Div ( a g Grad_perp(f) )  -- Perpendicular gradient-driven advection
+///
+/// This version uses a slope limiter to calculate cell edge values of g in X,
+/// the advects the upwind cell edge.
+///
+/// 1st order upwinding is used in Y.
+template <typename CellEdges = MC>
+const Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, const Field3D& f);
+
 } // namespace FV
 #endif // BOUT_FV_OPS_H
diff --git a/include/bout/fv_ops_impl.hxx b/include/bout/fv_ops_impl.hxx
index ac9041a19b..ecbfe4add4 100644
--- a/include/bout/fv_ops_impl.hxx
+++ b/include/bout/fv_ops_impl.hxx
@@ -270,12 +270,6 @@ struct WENO3 {
   }
 };
 
-/*!
-   * Communicate fluxes between processors
-   * Takes values in guard cells, and adds them to cells
-   */
-void communicateFluxes(Field3D& f);
-
 /// Finite volume parallel divergence
 ///
 /// Preserves the sum of f*J*dx*dy*dz over the domain
@@ -984,10 +978,10 @@ Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in,
 
 // Calculates viscous heating due to numerical momentum fluxes
 // and flow of kinetic energy (in flow_ylow)
-template <typename CellEdges = MC>
+template <typename CellEdges>
 Field3D Div_par_fvv_heating(const Field3D& f_in, const Field3D& v_in,
                             const Field3D& wave_speed_in, Field3D& flow_ylow,
-                            bool fixflux = true) {
+                            bool fixflux) {
 
   ASSERT1(areFieldsCompatible(f_in, v_in));
   ASSERT1(areFieldsCompatible(f_in, wave_speed_in));
@@ -1249,7 +1243,7 @@ Field3D Div_par_fvv_heating(const Field3D& f_in, const Field3D& v_in,
 /// the advects the upwind cell edge.
 ///
 /// 1st order upwinding is used in Y.
-template <typename CellEdges = MC>
+template <typename CellEdges>
 const Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g,
                                     const Field3D& f) {
   ASSERT2(a.getLocation() == f.getLocation());

From d9277946570e15bc6dff8125924fa3daa40a2acf Mon Sep 17 00:00:00 2001
From: David Bold <dave@ipp.mpg.de>
Date: Thu, 2 Jul 2026 13:09:41 +0200
Subject: [PATCH 04/11] Add explicit instantiations for moved functions

---
 src/mesh/fv_ops.cxx | 36 +++++++++++++++++++++++++++++++++++-
 1 file changed, 35 insertions(+), 1 deletion(-)

diff --git a/src/mesh/fv_ops.cxx b/src/mesh/fv_ops.cxx
index 1bd00d34a9..01c7ca8352 100644
--- a/src/mesh/fv_ops.cxx
+++ b/src/mesh/fv_ops.cxx
@@ -573,6 +573,11 @@ template Field3D Div_f_v<Upwind>(const Field3D& n_in, const Vector3D& v, bool bn
 template Field3D Div_par_mod<Upwind>(const Field3D& f_in, const Field3D& v_in,
                                      const Field3D& wave_speed_in, Field3D& flow_ylow,
                                      bool fixflux = true);
+template Field3D Div_par_fvv_heating<Upwind>(const Field3D& f_in, const Field3D& v_in,
+                                             const Field3D& wave_speed_in,
+                                             Field3D& flow_ylow, bool fixflux = true);
+template Field3D Div_a_Grad_perp_limit<Upwind>(const Field3D& a, const Field3D& g,
+                                               const Field3D& f);
 
 template Field3D Div_par_fvv<Fromm>(const Field3D& f_in, const Field3D& v_in,
                                     const Field3D& wave_speed_in, bool fixflux = true);
@@ -582,6 +587,11 @@ template Field3D Div_f_v<Fromm>(const Field3D& n_in, const Vector3D& v, bool bnd
 template Field3D Div_par_mod<Fromm>(const Field3D& f_in, const Field3D& v_in,
                                     const Field3D& wave_speed_in, Field3D& flow_ylow,
                                     bool fixflux = true);
+template Field3D Div_par_fvv_heating<Fromm>(const Field3D& f_in, const Field3D& v_in,
+                                            const Field3D& wave_speed_in,
+                                            Field3D& flow_ylow, bool fixflux = true);
+template Field3D Div_a_Grad_perp_limit<Fromm>(const Field3D& a, const Field3D& g,
+                                              const Field3D& f);
 
 template Field3D Div_par_fvv<MinMod>(const Field3D& f_in, const Field3D& v_in,
                                      const Field3D& wave_speed_in, bool fixflux = true);
@@ -591,6 +601,11 @@ template Field3D Div_f_v<MinMod>(const Field3D& n_in, const Vector3D& v, bool bn
 template Field3D Div_par_mod<MinMod>(const Field3D& f_in, const Field3D& v_in,
                                      const Field3D& wave_speed_in, Field3D& flow_ylow,
                                      bool fixflux = true);
+template Field3D Div_par_fvv_heating<MinMod>(const Field3D& f_in, const Field3D& v_in,
+                                             const Field3D& wave_speed_in,
+                                             Field3D& flow_ylow, bool fixflux = true);
+template Field3D Div_a_Grad_perp_limit<MinMod>(const Field3D& a, const Field3D& g,
+                                               const Field3D& f);
 
 template Field3D Div_par_fvv<MC>(const Field3D& f_in, const Field3D& v_in,
                                  const Field3D& wave_speed_in, bool fixflux = true);
@@ -600,6 +615,11 @@ template Field3D Div_f_v<MC>(const Field3D& n_in, const Vector3D& v, bool bndry_
 template Field3D Div_par_mod<MC>(const Field3D& f_in, const Field3D& v_in,
                                  const Field3D& wave_speed_in, Field3D& flow_ylow,
                                  bool fixflux = true);
+template Field3D Div_par_fvv_heating<MC>(const Field3D& f_in, const Field3D& v_in,
+                                         const Field3D& wave_speed_in, Field3D& flow_ylow,
+                                         bool fixflux = true);
+template Field3D Div_a_Grad_perp_limit<MC>(const Field3D& a, const Field3D& g,
+                                           const Field3D& f);
 
 template Field3D Div_par_fvv<Superbee>(const Field3D& f_in, const Field3D& v_in,
                                        const Field3D& wave_speed_in, bool fixflux = true);
@@ -610,6 +630,11 @@ template Field3D Div_f_v<Superbee>(const Field3D& n_in, const Vector3D& v,
 template Field3D Div_par_mod<Superbee>(const Field3D& f_in, const Field3D& v_in,
                                        const Field3D& wave_speed_in, Field3D& flow_ylow,
                                        bool fixflux = true);
+template Field3D Div_par_fvv_heating<Superbee>(const Field3D& f_in, const Field3D& v_in,
+                                               const Field3D& wave_speed_in,
+                                               Field3D& flow_ylow, bool fixflux = true);
+template Field3D Div_a_Grad_perp_limit<Superbee>(const Field3D& a, const Field3D& g,
+                                                 const Field3D& f);
 
 template Field3D Div_par_fvv<VanAlbada>(const Field3D& f_in, const Field3D& v_in,
                                         const Field3D& wave_speed_in,
@@ -621,6 +646,11 @@ template Field3D Div_f_v<VanAlbada>(const Field3D& n_in, const Vector3D& v,
 template Field3D Div_par_mod<VanAlbada>(const Field3D& f_in, const Field3D& v_in,
                                         const Field3D& wave_speed_in, Field3D& flow_ylow,
                                         bool fixflux = true);
+template Field3D Div_par_fvv_heating<VanAlbada>(const Field3D& f_in, const Field3D& v_in,
+                                                const Field3D& wave_speed_in,
+                                                Field3D& flow_ylow, bool fixflux = true);
+template Field3D Div_a_Grad_perp_limit<VanAlbada>(const Field3D& a, const Field3D& g,
+                                                  const Field3D& f);
 
 template Field3D Div_par_fvv<WENO3>(const Field3D& f_in, const Field3D& v_in,
                                     const Field3D& wave_speed_in, bool fixflux = true);
@@ -630,5 +660,9 @@ template Field3D Div_f_v<WENO3>(const Field3D& n_in, const Vector3D& v, bool bnd
 template Field3D Div_par_mod<WENO3>(const Field3D& f_in, const Field3D& v_in,
                                     const Field3D& wave_speed_in, Field3D& flow_ylow,
                                     bool fixflux = true);
-
+template Field3D Div_par_fvv_heating<WENO3>(const Field3D& f_in, const Field3D& v_in,
+                                            const Field3D& wave_speed_in,
+                                            Field3D& flow_ylow, bool fixflux = true);
+template Field3D Div_a_Grad_perp_limit<WENO3>(const Field3D& a, const Field3D& g,
+                                              const Field3D& f);
 } // Namespace FV

From 3c4383d12ec98ede8b6154c427ae0a4c00604e77 Mon Sep 17 00:00:00 2001
From: David Bold <dave@ipp.mpg.de>
Date: Thu, 2 Jul 2026 13:12:13 +0200
Subject: [PATCH 05/11] Return should be non-const

---
 include/bout/fv_ops.hxx      | 2 +-
 include/bout/fv_ops_impl.hxx | 3 +--
 2 files changed, 2 insertions(+), 3 deletions(-)

diff --git a/include/bout/fv_ops.hxx b/include/bout/fv_ops.hxx
index c7a6f6774c..98dc689bfd 100644
--- a/include/bout/fv_ops.hxx
+++ b/include/bout/fv_ops.hxx
@@ -180,7 +180,7 @@ Field3D Div_par_fvv_heating(const Field3D& f_in, const Field3D& v_in,
 ///
 /// 1st order upwinding is used in Y.
 template <typename CellEdges = MC>
-const Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, const Field3D& f);
+Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, const Field3D& f);
 
 } // namespace FV
 #endif // BOUT_FV_OPS_H
diff --git a/include/bout/fv_ops_impl.hxx b/include/bout/fv_ops_impl.hxx
index ecbfe4add4..9e436ad95d 100644
--- a/include/bout/fv_ops_impl.hxx
+++ b/include/bout/fv_ops_impl.hxx
@@ -1244,8 +1244,7 @@ Field3D Div_par_fvv_heating(const Field3D& f_in, const Field3D& v_in,
 ///
 /// 1st order upwinding is used in Y.
 template <typename CellEdges>
-const Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g,
-                                    const Field3D& f) {
+Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, const Field3D& f) {
   ASSERT2(a.getLocation() == f.getLocation());
 
   Mesh* mesh = a.getMesh();

From 665798b525bf49d3f728dee14c9db7408edd515e Mon Sep 17 00:00:00 2001
From: David Bold <dave@ipp.mpg.de>
Date: Mon, 6 Jul 2026 15:25:17 +0200
Subject: [PATCH 06/11] Apply clang-tidy fixes

---
 include/bout/fv_ops_impl.hxx | 65 +++++++++++++++++++-----------------
 1 file changed, 34 insertions(+), 31 deletions(-)

diff --git a/include/bout/fv_ops_impl.hxx b/include/bout/fv_ops_impl.hxx
index 9e436ad95d..9400f5deb1 100644
--- a/include/bout/fv_ops_impl.hxx
+++ b/include/bout/fv_ops_impl.hxx
@@ -986,7 +986,7 @@ Field3D Div_par_fvv_heating(const Field3D& f_in, const Field3D& v_in,
   ASSERT1(areFieldsCompatible(f_in, v_in));
   ASSERT1(areFieldsCompatible(f_in, wave_speed_in));
 
-  Mesh* mesh = f_in.getMesh();
+  const Mesh* mesh = f_in.getMesh();
   Coordinates* coord = f_in.getCoordinates();
   CellEdges cellboundary;
 
@@ -1138,23 +1138,23 @@ Field3D Div_par_fvv_heating(const Field3D& f_in, const Field3D& v_in,
           // energy losses.
           const BoutReal expected_ke = 0.5 * n_mid * v_mid * v_mid * v_mid;
 
-          BoutReal flux_mom;
+          BoutReal flux_mom = NAN;
           if (fixflux) {
             // Mid-point consistent with boundary conditions
             // but kinetic energy loss will not match expected
             // -> Adjust energy balance in pressure equation
             flux_mom = n_mid * v_mid * v_mid;
           } else {
-            flux_mom = s.R * sv.R * sv.R
-                       + BOUTMAX(wave_speed(i, j, k), fabs(sv.c), fabs(sv.p))
-                             * (s.R * sv.R - n_mid * v_mid);
+            flux_mom = (s.R * sv.R * sv.R)
+                       + (BOUTMAX(wave_speed(i, j, k), fabs(sv.c), fabs(sv.p))
+                          * (s.R * sv.R - n_mid * v_mid));
           }
 
           // Assume that particle flux is fixed to boundary value
           const BoutReal flux_part = n_mid * v_mid;
 
           // d/dt(1/2 m n v^2) = v * d/dt(mnv) - 1/2 m v^2 * dn/dt
-          const BoutReal actual_ke = sv.c * flux_mom - 0.5 * sv.c * sv.c * flux_part;
+          const BoutReal actual_ke = (sv.c * flux_mom) - (0.5 * sv.c * sv.c * flux_part);
 
           // Note: If the actual loss was higher than expected, then
           //       plasma heating is needed to compensate
@@ -1193,22 +1193,22 @@ Field3D Div_par_fvv_heating(const Field3D& f_in, const Field3D& v_in,
 
         if (mesh->firstY(i) && (j == mesh->ystart) && !mesh->periodicY(i)) {
           // First point in domain
-          BoutReal flux_mom;
+          BoutReal flux_mom = NAN;
           if (fixflux) {
             // Use mid-point to be consistent with boundary conditions
             flux_mom = n_mid * v_mid * v_mid;
           } else {
             // Add flux due to difference in boundary values
-            flux_mom = s.L * sv.L * sv.L
-                       - BOUTMAX(wave_speed(i, j, k), fabs(sv.c), fabs(sv.m))
-                             * (s.L * sv.L - n_mid * v_mid);
+            flux_mom = (s.L * sv.L * sv.L)
+                       - (BOUTMAX(wave_speed(i, j, k), fabs(sv.c), fabs(sv.m))
+                          * (s.L * sv.L - n_mid * v_mid));
           }
 
           // Assume that density flux is fixed to boundary value
           const BoutReal flux_part = n_mid * v_mid;
 
           // d/dt(1/2 m n v^2) = v * d/dt(mnv) - 1/2 m v^2 * dn/dt
-          const BoutReal actual_ke = -sv.c * flux_mom + 0.5 * sv.c * sv.c * flux_part;
+          const BoutReal actual_ke = (-sv.c * flux_mom) + (0.5 * sv.c * sv.c * flux_part);
 
           result(i, j, k) += (actual_ke - expected_ke) * flux_factor_lc;
 
@@ -1270,7 +1270,7 @@ Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, const Field3D&
 
         // Mid-point average boundary value of 'a'
         const BoutReal aedge = 0.5 * (a(i + 1, j, k) + a(i, j, k));
-        BoutReal gedge;
+        BoutReal gedge = NAN;
         if (((i == mesh->xstart - 1) and mesh->firstX())
             or ((i == mesh->xend) and mesh->lastX())) {
           // Mid-point average boundary value of 'g'
@@ -1314,9 +1314,12 @@ Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, const Field3D&
   // Y and Z fluxes require Y derivatives
 
   // Fields containing values along the magnetic field
-  Field3D fup(mesh), fdown(mesh);
-  Field3D aup(mesh), adown(mesh);
-  Field3D gup(mesh), gdown(mesh);
+  Field3D fup(mesh);
+  Field3D fdown(mesh);
+  Field3D aup(mesh);
+  Field3D adown(mesh);
+  Field3D gup(mesh);
+  Field3D gdown(mesh);
 
   // Values on this y slice (centre).
   // This is needed because toFieldAligned may modify the field
@@ -1356,16 +1359,16 @@ Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, const Field3D&
 #if BOUT_USE_METRIC_3D
       for (int k = 0; k < mesh->LocalNz; k++)
 #else
-      int k = 0;
+      const int k = 0;
 #endif
       {
-        BoutReal coef_u =
+        const BoutReal coef_u =
             0.5
             * (coord->g_23(i, j, k) / SQ(coord->J(i, j, k) * coord->Bxy(i, j, k))
                + coord->g_23(i, j + 1, k)
                      / SQ(coord->J(i, j + 1, k) * coord->Bxy(i, j + 1, k)));
 
-        BoutReal coef_d =
+        const BoutReal coef_d =
             0.5
             * (coord->g_23(i, j, k) / SQ(coord->J(i, j, k) * coord->Bxy(i, j, k))
                + coord->g_23(i, j - 1, k)
@@ -1376,8 +1379,8 @@ Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, const Field3D&
 #endif
         {
           // Calculate flux between j and j+1
-          int kp = (k + 1) % mesh->LocalNz;
-          int km = (k - 1 + mesh->LocalNz) % mesh->LocalNz;
+          const int kp = (k + 1) % mesh->LocalNz;
+          const int km = (k - 1 + mesh->LocalNz) % mesh->LocalNz;
 
           // Calculate Z derivative at y boundary
           BoutReal dfdz =
@@ -1389,7 +1392,7 @@ Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, const Field3D&
                           / (coord->dy(i, j + 1, k) + coord->dy(i, j, k));
 
           BoutReal aedge = 0.5 * (ac(i, j, k) + aup(i, j + 1, k));
-          BoutReal gedge;
+          BoutReal gedge = NAN;
           if ((j == mesh->yend) and mesh->lastY(i)) {
             // Midpoint boundary value
             gedge = 0.5 * (gc(i, j, k) + gup(i, j + 1, k));
@@ -1445,11 +1448,11 @@ Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, const Field3D&
 #if BOUT_USE_METRIC_3D
       for (int k = 0; k < mesh->LocalNz; k++)
 #else
-      int k = 0;
+      const int k = 0;
 #endif
       {
         // Coefficient in front of df/dy term
-        BoutReal coef =
+        const BoutReal coef =
             coord->g_23(i, j, k)
             / (coord->dy(i, j + 1, k) + 2. * coord->dy(i, j, k) + coord->dy(i, j - 1, k))
             / SQ(coord->J(i, j, k) * coord->Bxy(i, j, k));
@@ -1458,19 +1461,19 @@ Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, const Field3D&
 #endif
         {
           // Calculate flux between k and k+1
-          int kp = (k + 1) % mesh->LocalNz;
+          const int kp = (k + 1) % mesh->LocalNz;
 
-          BoutReal gradient =
+          const BoutReal gradient =
               // df/dz
-              (fc(i, j, kp) - fc(i, j, k)) / coord->dz(i, j, k)
+              ((fc(i, j, kp) - fc(i, j, k)) / coord->dz(i, j, k))
 
               // - g_yz * df/dy / SQ(J*B)
-              - coef
-                    * (fup(i, j + 1, k) + fup(i, j + 1, kp) - fdown(i, j - 1, k)
-                       - fdown(i, j - 1, kp));
+              - (coef
+                 * (fup(i, j + 1, k) + fup(i, j + 1, kp) - fdown(i, j - 1, k)
+                    - fdown(i, j - 1, kp)));
 
-          BoutReal fout = gradient * 0.5 * (ac(i, j, kp) + ac(i, j, k))
-                          * ((gradient > 0) ? gc(i, j, kp) : gc(i, j, k));
+          const BoutReal fout = gradient * 0.5 * (ac(i, j, kp) + ac(i, j, k))
+                                * ((gradient > 0) ? gc(i, j, kp) : gc(i, j, k));
 
           yzresult(i, j, k) += fout / coord->dz(i, j, k);
           yzresult(i, j, kp) -= fout / coord->dz(i, j, kp);

From 6eeaacb1d4404b2f3cdf0c7dfe72ba68c5c34915 Mon Sep 17 00:00:00 2001
From: David Bold <dave@ipp.mpg.de>
Date: Mon, 6 Jul 2026 15:27:31 +0200
Subject: [PATCH 07/11] Prefer std::abs

This is a template that always uses the correct datatype, so we do not
implicitly cast as double.
---
 include/bout/fv_ops_impl.hxx | 45 ++++++++++++++++++------------------
 1 file changed, 23 insertions(+), 22 deletions(-)

diff --git a/include/bout/fv_ops_impl.hxx b/include/bout/fv_ops_impl.hxx
index 9400f5deb1..09b4a2f87f 100644
--- a/include/bout/fv_ops_impl.hxx
+++ b/include/bout/fv_ops_impl.hxx
@@ -91,7 +91,7 @@ private:
       return 0.0;
     }
 
-    if (fabs(a) < fabs(b)) {
+    if (std::abs(a) < std::abs(b)) {
       return a;
     }
     return b;
@@ -126,7 +126,7 @@ private:
     }
 
     // Return the minimum absolute value
-    return SIGN(a) * BOUTMIN(fabs(a), fabs(b), fabs(c));
+    return SIGN(a) * BOUTMIN(std::abs(a), std::abs(b), std::abs(c));
   }
 };
 
@@ -167,8 +167,8 @@ struct Superbee {
       n.L = n.R = n.c;
     } else {
       const BoutReal sign = SIGN(gL);
-      const BoutReal abs_gL = fabs(gL);
-      const BoutReal abs_gR = fabs(gR);
+      const BoutReal abs_gL = std::abs(gL);
+      const BoutReal abs_gR = std::abs(gR);
       const BoutReal half_slope =
           sign * BOUTMAX(BOUTMIN(abs_gL, 0.5 * abs_gR), BOUTMIN(abs_gR, 0.5 * abs_gL));
       n.L = n.c - half_slope;
@@ -741,7 +741,7 @@ Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in,
         } else {
           // Maximum wave speed in the two cells
           const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j + 1, k),
-                                        fabs(v(i, j, k)), fabs(v(i, j + 1, k)));
+                                        std::abs(v(i, j, k)), std::abs(v(i, j + 1, k)));
 
           flux = s.R * 0.5 * (sv.R + amax);
         }
@@ -769,7 +769,7 @@ Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in,
 
           // Maximum wave speed in the two cells
           const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j - 1, k),
-                                        fabs(v(i, j, k)), fabs(v(i, j - 1, k)));
+                                        std::abs(v(i, j, k)), std::abs(v(i, j - 1, k)));
 
           flux = s.L * 0.5 * (sv.L - amax);
         }
@@ -825,8 +825,8 @@ Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in,
       const auto iym = i.ym();
 
       // Maximum local wave speed
-      const BoutReal amax =
-          BOUTMAX(wave_speed_in[i], fabs(v_in[i]), fabs(v_up[iyp]), fabs(v_down[iym]));
+      const BoutReal amax = BOUTMAX(wave_speed_in[i], std::abs(v_in[i]),
+                                    std::abs(v_up[iyp]), std::abs(v_down[iym]));
 
       const BoutReal term = (f_up[iyp] * v_up[iyp] * v_up[iyp] / B_up[iyp])
                             - (f_down[iym] * v_down[iym] * v_down[iym] / B_down[iym]);
@@ -837,7 +837,7 @@ Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in,
           + (amax * (f_in[i] * v_in[i] - f_down[iym] * v_down[iym])
              / (B[i] + B_down[iym]));
 
-      if (fabs(penalty) > fabs(term) and penalty * v_in[i] > 0) {
+      if (std::abs(penalty) > std::abs(term) and penalty * v_in[i] > 0) {
         if (term * penalty > 0) {
           penalty = term;
         } else {
@@ -928,14 +928,15 @@ Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in,
             flux = n_mid_r * v_mid_r * v_mid_r;
           } else {
             // Add flux due to difference in boundary values
-            flux = (s.R * sv.R * sv.R) // Use right cell edge values
-                   + (BOUTMAX(wave_speed(i, j, k), fabs(sv.c), fabs(sv.p)) * n_mid_r
-                      * (sv.R - v_mid_r)); // Damp differences in velocity, not flux
+            flux =
+                (s.R * sv.R * sv.R) // Use right cell edge values
+                + (BOUTMAX(wave_speed(i, j, k), std::abs(sv.c), std::abs(sv.p)) * n_mid_r
+                   * (sv.R - v_mid_r)); // Damp differences in velocity, not flux
           }
         } else {
           // Maximum wave speed in the two cells
           const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j + 1, k),
-                                        fabs(sv.c), fabs(sv.p));
+                                        std::abs(sv.c), std::abs(sv.p));
 
           flux = s.R * 0.5 * (sv.R + amax) * sv.R;
         }
@@ -957,13 +958,13 @@ Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in,
           } else {
             // Add flux due to difference in boundary values
             flux = (s.L * sv.L * sv.L)
-                   - (BOUTMAX(wave_speed(i, j, k), fabs(sv.c), fabs(sv.m)) * n_mid_l
-                      * (sv.L - v_mid_l));
+                   - (BOUTMAX(wave_speed(i, j, k), std::abs(sv.c), std::abs(sv.m))
+                      * n_mid_l * (sv.L - v_mid_l));
           }
         } else {
           // Maximum wave speed in the two cells
           const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j - 1, k),
-                                        fabs(sv.c), fabs(sv.m));
+                                        std::abs(sv.c), std::abs(sv.m));
 
           flux = s.L * 0.5 * (sv.L - amax) * sv.L;
         }
@@ -1016,8 +1017,8 @@ Field3D Div_par_fvv_heating(const Field3D& f_in, const Field3D& v_in,
       const auto iym = i.ym();
 
       //Maximum local wave speed
-      const BoutReal amax =
-          BOUTMAX(wave_speed_in[i], fabs(v_in[i]), fabs(v_up[iyp]), fabs(v_down[iym]));
+      const BoutReal amax = BOUTMAX(wave_speed_in[i], std::abs(v_in[i]),
+                                    std::abs(v_up[iyp]), std::abs(v_down[iym]));
 
       result[i] =
           B[i]
@@ -1146,7 +1147,7 @@ Field3D Div_par_fvv_heating(const Field3D& f_in, const Field3D& v_in,
             flux_mom = n_mid * v_mid * v_mid;
           } else {
             flux_mom = (s.R * sv.R * sv.R)
-                       + (BOUTMAX(wave_speed(i, j, k), fabs(sv.c), fabs(sv.p))
+                       + (BOUTMAX(wave_speed(i, j, k), std::abs(sv.c), std::abs(sv.p))
                           * (s.R * sv.R - n_mid * v_mid));
           }
 
@@ -1166,7 +1167,7 @@ Field3D Div_par_fvv_heating(const Field3D& f_in, const Field3D& v_in,
         } else {
           // Maximum wave speed in the two cells
           const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j + 1, k),
-                                        fabs(sv.c), fabs(sv.p));
+                                        std::abs(sv.c), std::abs(sv.p));
 
           // Viscous heating due to relaxation of velocity towards midpoint
           result(i, j, k) +=
@@ -1200,7 +1201,7 @@ Field3D Div_par_fvv_heating(const Field3D& f_in, const Field3D& v_in,
           } else {
             // Add flux due to difference in boundary values
             flux_mom = (s.L * sv.L * sv.L)
-                       - (BOUTMAX(wave_speed(i, j, k), fabs(sv.c), fabs(sv.m))
+                       - (BOUTMAX(wave_speed(i, j, k), std::abs(sv.c), std::abs(sv.m))
                           * (s.L * sv.L - n_mid * v_mid));
           }
 
@@ -1216,7 +1217,7 @@ Field3D Div_par_fvv_heating(const Field3D& f_in, const Field3D& v_in,
         } else {
           // Maximum wave speed in the two cells
           const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j - 1, k),
-                                        fabs(sv.c), fabs(sv.m));
+                                        std::abs(sv.c), std::abs(sv.m));
 
           // Viscous heating due to relaxation
           result(i, j, k) +=

From fcd84b3b1625330d3391a7a2536f3bb59a287e58 Mon Sep 17 00:00:00 2001
From: Ben Dudson <dudson2@llnl.gov>
Date: Tue, 7 Jul 2026 14:54:23 -0700
Subject: [PATCH 08/11] Div_a_Grad_perp_limit: Use slices like Div_a_Grad_perp

Functions both now use `makeslices` to handle field-aligned
or FCI metrics.
---
 include/bout/fv_ops_impl.hxx | 169 +++++++++++++++++++----------------
 1 file changed, 93 insertions(+), 76 deletions(-)

diff --git a/include/bout/fv_ops_impl.hxx b/include/bout/fv_ops_impl.hxx
index 9d169223b6..7ab4c2a042 100644
--- a/include/bout/fv_ops_impl.hxx
+++ b/include/bout/fv_ops_impl.hxx
@@ -270,6 +270,22 @@ struct WENO3 {
   }
 };
 
+template <class T>
+struct Slices {
+  T c;
+  T up;
+  T down;
+
+  Slices(bool use_slices, const T& field)
+      : c(use_slices ? field : toFieldAligned(field)), up(use_slices ? field.yup() : c),
+        down(use_slices ? field.ydown() : c) {}
+};
+
+template <class T>
+Slices<T> makeslices(bool use_slices, const T& field) {
+  return Slices<T>(use_slices, field);
+}
+
 /// Finite volume parallel divergence
 ///
 /// Preserves the sum of f*J*dx*dy*dz over the domain
@@ -1313,44 +1329,44 @@ Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, const Field3D&
     }
   }
 
-  // Y and Z fluxes require Y derivatives
+  const bool fci =
+      f.hasParallelSlices() && a.hasParallelSlices() && g.hasParallelSlices();
 
-  // Fields containing values along the magnetic field
-  Field3D fup(mesh);
-  Field3D fdown(mesh);
-  Field3D aup(mesh);
-  Field3D adown(mesh);
-  Field3D gup(mesh);
-  Field3D gdown(mesh);
+#if BOUT_USE_METRIC_3D
+  if (fci) {
+    // 3D Metric, need yup/ydown fields.
+    // Requires previous communication of metrics.
+    if (!coord->g23.hasParallelSlices() || !coord->g_23.hasParallelSlices()
+        || !coord->dy.hasParallelSlices() || !coord->dz.hasParallelSlices()
+        || !coord->Bxy.hasParallelSlices() || !coord->J.hasParallelSlices()) {
+      throw BoutException("metrics have no yup/down: Maybe communicate in init?");
+    }
+  }
+#endif
+
+  // Y and Z fluxes require Y derivatives
 
   // Values on this y slice (centre).
   // This is needed because toFieldAligned may modify the field
-  Field3D ac = a;
-  Field3D gc = g;
-  Field3D fc = f;
+  const auto f_slice = makeslices(fci, f);
+  const auto a_slice = makeslices(fci, a);
+  const auto g_slice = makeslices(fci, g);
 
-  // Result of the Y and Z fluxes
-  Field3D yzresult(mesh);
-  yzresult.allocate();
-
-  if (f.hasParallelSlices() && a.hasParallelSlices() && g.hasParallelSlices()) {
-    // All inputs have yup and ydown
-
-    fup = f.yup();
-    fdown = f.ydown();
-
-    aup = a.yup();
-    adown = a.ydown();
-
-    gup = g.yup();
-    gdown = g.ydown();
-  } else {
-    // At least one input doesn't have yup/ydown fields.
-    // Need to shift to/from field aligned coordinates
+#if BOUT_USE_METRIC_3D
+  const bool metric_fci = fci;
+#else
+  constexpr bool metric_fci = false;
+#endif
+  const auto g23 = makeslices(metric_fci, coord->g23);
+  const auto g_23 = makeslices(metric_fci, coord->g_23);
+  const auto J = makeslices(metric_fci, coord->J);
+  const auto dy = makeslices(metric_fci, coord->dy);
+  const auto dz = makeslices(metric_fci, coord->dz);
+  const auto Bxy = makeslices(metric_fci, coord->Bxy);
 
-    fup = fdown = fc = toFieldAligned(f);
-    aup = adown = ac = toFieldAligned(a);
-    gup = gdown = gc = toFieldAligned(g);
+  // Result of the Y and Z fluxes
+  Field3D yzresult(0.0, mesh);
+  if (!fci) {
     yzresult.setDirectionY(YDirectionType::Aligned);
   }
 
@@ -1366,15 +1382,14 @@ Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, const Field3D&
       {
         const BoutReal coef_u =
             0.5
-            * (coord->g_23(i, j, k) / SQ(coord->J(i, j, k) * coord->Bxy(i, j, k))
-               + coord->g_23(i, j + 1, k)
-                     / SQ(coord->J(i, j + 1, k) * coord->Bxy(i, j + 1, k)));
+            * (g_23.c(i, j, k) / SQ(J.c(i, j, k) * Bxy.c(i, j, k))
+               + g_23.up(i, j + 1, k) / SQ(J.up(i, j + 1, k) * Bxy.up(i, j + 1, k)));
 
         const BoutReal coef_d =
             0.5
-            * (coord->g_23(i, j, k) / SQ(coord->J(i, j, k) * coord->Bxy(i, j, k))
-               + coord->g_23(i, j - 1, k)
-                     / SQ(coord->J(i, j - 1, k) * coord->Bxy(i, j - 1, k)));
+            * (g_23.c(i, j, k) / SQ(J.c(i, j, k) * Bxy.c(i, j, k))
+               + g_23.down(i, j - 1, k)
+                     / SQ(J.down(i, j - 1, k) * Bxy.down(i, j - 1, k)));
 
 #if not BOUT_USE_METRIC_3D
         for (int k = 0; k < mesh->LocalNz; k++)
@@ -1385,58 +1400,59 @@ Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, const Field3D&
           const int km = (k - 1 + mesh->LocalNz) % mesh->LocalNz;
 
           // Calculate Z derivative at y boundary
-          BoutReal dfdz =
-              0.25 * (fc(i, j, kp) - fc(i, j, km) + fup(i, j + 1, kp) - fup(i, j + 1, km))
-              / coord->dz(i, j, k);
+          BoutReal dfdz = 0.5
+                          * (f_slice.c(i, j, kp) - f_slice.c(i, j, km)
+                             + f_slice.up(i, j + 1, kp) - f_slice.up(i, j + 1, km))
+                          / (dz.c(i, j, k) + dz.up(i, j + 1, k));
 
           // Y derivative
-          BoutReal dfdy = 2. * (fup(i, j + 1, k) - fc(i, j, k))
-                          / (coord->dy(i, j + 1, k) + coord->dy(i, j, k));
+          BoutReal dfdy = 2. * (f_slice.up(i, j + 1, k) - f_slice.c(i, j, k))
+                          / (dy.up(i, j + 1, k) + dy.c(i, j, k));
 
-          BoutReal aedge = 0.5 * (ac(i, j, k) + aup(i, j + 1, k));
+          BoutReal aedge = 0.5 * (a_slice.c(i, j, k) + a_slice.up(i, j + 1, k));
           BoutReal gedge = NAN;
           if ((j == mesh->yend) and mesh->lastY(i)) {
             // Midpoint boundary value
-            gedge = 0.5 * (gc(i, j, k) + gup(i, j + 1, k));
+            gedge = 0.5 * (g_slice.c(i, j, k) + g_slice.up(i, j + 1, k));
           } else if (dfdy > 0) {
             // Flux from (j+1) to (j)
-            gedge = gup(i, j + 1, k);
+            gedge = g_slice.up(i, j + 1, k);
           } else {
             // Flux from (j) to (j+1)
-            gedge = gc(i, j, k);
+            gedge = g_slice.c(i, j, k);
           }
 
-          BoutReal fout = aedge * gedge * 0.5
-                          * (coord->J(i, j, k) * coord->g23(i, j, k)
-                             + coord->J(i, j + 1, k) * coord->g23(i, j + 1, k))
-                          * (dfdz - coef_u * dfdy);
+          BoutReal fout =
+              aedge * gedge * 0.5
+              * (J.c(i, j, k) * g23.c(i, j, k) + J.up(i, j + 1, k) * g23.up(i, j + 1, k))
+              * (dfdz - coef_u * dfdy);
 
-          yzresult(i, j, k) = fout / (coord->dy(i, j, k) * coord->J(i, j, k));
+          yzresult(i, j, k) += fout / (dy.c(i, j, k) * J.c(i, j, k));
 
           // Calculate flux between j and j-1
-          dfdz =
-              0.25
-              * (fc(i, j, kp) - fc(i, j, km) + fdown(i, j - 1, kp) - fdown(i, j - 1, km))
-              / coord->dz(i, j, k);
+          dfdz = 0.5
+                 * (f_slice.c(i, j, kp) - f_slice.c(i, j, km) + f_slice.down(i, j - 1, kp)
+                    - f_slice.down(i, j - 1, km))
+                 / (dz.c(i, j, k) + dz.down(i, j - 1, k));
 
-          dfdy = 2. * (fc(i, j, k) - fdown(i, j - 1, k))
-                 / (coord->dy(i, j, k) + coord->dy(i, j - 1, k));
+          dfdy = 2. * (f_slice.c(i, j, k) - f_slice.down(i, j - 1, k))
+                 / (dy.c(i, j, k) + dy.down(i, j - 1, k));
 
-          aedge = 0.5 * (ac(i, j, k) + adown(i, j - 1, k));
+          aedge = 0.5 * (a_slice.c(i, j, k) + a_slice.down(i, j - 1, k));
           if ((j == mesh->ystart) and mesh->firstY(i)) {
-            gedge = 0.5 * (gc(i, j, k) + gdown(i, j - 1, k));
+            gedge = 0.5 * (g_slice.c(i, j, k) + g_slice.down(i, j - 1, k));
           } else if (dfdy > 0) {
-            gedge = gc(i, j, k);
+            gedge = g_slice.c(i, j, k);
           } else {
-            gedge = gdown(i, j - 1, k);
+            gedge = g_slice.down(i, j - 1, k);
           }
 
           fout = aedge * gedge * 0.5
-                 * (coord->J(i, j, k) * coord->g23(i, j, k)
-                    + coord->J(i, j - 1, k) * coord->g23(i, j - 1, k))
+                 * (J.c(i, j, k) * g23.c(i, j, k)
+                    + J.down(i, j - 1, k) * g23.down(i, j - 1, k))
                  * (dfdz - coef_d * dfdy);
 
-          yzresult(i, j, k) -= fout / (coord->dy(i, j, k) * coord->J(i, j, k));
+          yzresult(i, j, k) -= fout / (dy.c(i, j, k) * J.c(i, j, k));
         }
       }
     }
@@ -1455,9 +1471,9 @@ Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, const Field3D&
       {
         // Coefficient in front of df/dy term
         const BoutReal coef =
-            coord->g_23(i, j, k)
-            / (coord->dy(i, j + 1, k) + 2. * coord->dy(i, j, k) + coord->dy(i, j - 1, k))
-            / SQ(coord->J(i, j, k) * coord->Bxy(i, j, k));
+            g_23.c(i, j, k)
+            / (dy.up(i, j + 1, k) + 2. * dy.c(i, j, k) + dy.down(i, j - 1, k))
+            / SQ(J.c(i, j, k) * Bxy.c(i, j, k));
 #if not BOUT_USE_METRIC_3D
         for (int k = 0; k < mesh->LocalNz; k++)
 #endif
@@ -1467,24 +1483,25 @@ Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, const Field3D&
 
           const BoutReal gradient =
               // df/dz
-              ((fc(i, j, kp) - fc(i, j, k)) / coord->dz(i, j, k))
+              ((f_slice.c(i, j, kp) - f_slice.c(i, j, k)) / dz.c(i, j, k))
 
               // - g_yz * df/dy / SQ(J*B)
               - (coef
-                 * (fup(i, j + 1, k) + fup(i, j + 1, kp) - fdown(i, j - 1, k)
-                    - fdown(i, j - 1, kp)));
+                 * (f_slice.up(i, j + 1, k) + f_slice.up(i, j + 1, kp)
+                    - f_slice.down(i, j - 1, k) - f_slice.down(i, j - 1, kp)));
 
-          const BoutReal fout = gradient * 0.5 * (ac(i, j, kp) + ac(i, j, k))
-                                * ((gradient > 0) ? gc(i, j, kp) : gc(i, j, k));
+          const BoutReal fout =
+              gradient * 0.5 * (a_slice.c(i, j, kp) + a_slice.c(i, j, k))
+              * ((gradient > 0) ? g_slice.c(i, j, kp) : g_slice.c(i, j, k));
 
-          yzresult(i, j, k) += fout / coord->dz(i, j, k);
-          yzresult(i, j, kp) -= fout / coord->dz(i, j, kp);
+          yzresult(i, j, k) += fout / dz.c(i, j, k);
+          yzresult(i, j, kp) -= fout / dz.c(i, j, kp);
         }
       }
     }
   }
   // Check if we need to transform back
-  if (f.hasParallelSlices() && a.hasParallelSlices()) {
+  if (fci) {
     result += yzresult;
   } else {
     result += fromFieldAligned(yzresult);

From 5bc750a13e47225bcba91eeb74d9f8fdee6b9bf6 Mon Sep 17 00:00:00 2001
From: Ben Dudson <dudson2@llnl.gov>
Date: Tue, 7 Jul 2026 21:59:32 -0700
Subject: [PATCH 09/11] Document new FV ops

---
 manual/sphinx/user_docs/coordinates.rst       | 35 +++++------
 .../user_docs/differential_operators.rst      | 59 +++++++++++++++++++
 2 files changed, 77 insertions(+), 17 deletions(-)

diff --git a/manual/sphinx/user_docs/coordinates.rst b/manual/sphinx/user_docs/coordinates.rst
index 2ac086e05e..d270dda468 100644
--- a/manual/sphinx/user_docs/coordinates.rst
+++ b/manual/sphinx/user_docs/coordinates.rst
@@ -77,7 +77,7 @@ The cross products are:
 .. math::
    :label: eq:psithetazetacrossproducts
 
-   \boldsymbol{e}_\psi\times\boldsymbol{e}_\theta = J_{\psi\theta\zeta} \nabla\zeta \qquad 
+   \boldsymbol{e}_\psi\times\boldsymbol{e}_\theta = J_{\psi\theta\zeta} \nabla\zeta \qquad
    \boldsymbol{e}_\psi\times\boldsymbol{e}_\zeta = -J_{\psi\theta\zeta} \nabla\theta \qquad
    \boldsymbol{e}_\theta\times\boldsymbol{e}_\zeta = J_{\psi\theta\zeta} \nabla\psi
 
@@ -122,7 +122,7 @@ To align to the magnetic field we define a local field line pitch `\nu`:
    \frac{{B_{\text{tor}}}{h_\theta}}{{B_{\text{pol}}}R}
    \end{aligned}
 
-The sign of the poloidal field `{B_{\text{pol}}}` and toroidal field 
+The sign of the poloidal field `{B_{\text{pol}}}` and toroidal field
 `{B_{\text{tor}}}` can be either + or -.
 
 The field-aligned coordinates `\left(x,y,z\right)` are defined by:
@@ -158,7 +158,7 @@ The reciprocal basis vectors are
    \nabla z = \nabla\zeta - \sigma_{B\text{pol}}\left[\int_{\theta_0}^\theta{\frac{\partial \nu\left(\psi,\theta\right)}{\partial \psi}} d\theta\right] \nabla\psi
    - \sigma_{B\text{pol}}\nu\left(\psi, \theta\right)\nabla\theta
    \end{aligned}
-  
+
 The term in square brackets is the integrated local shear:
 
 .. math::
@@ -167,18 +167,18 @@ The term in square brackets is the integrated local shear:
    \begin{aligned}
    I = \int_{y_0}^y\frac{\partial\nu\left(x, y\right)}{\partial\psi}dy\end{aligned}
 
-  
+
 The basis vectors are:
 
 .. math::
    :label: eq:basisvectors
-   
+
    \begin{aligned}
    \boldsymbol{e}_x =& J_{xyz}\left(\nabla y \times \nabla z\right) = {\sigma_{B\text{pol}}} {\boldsymbol{e}}_\psi + I{\boldsymbol{e}}_\zeta \\
    \boldsymbol{e}_y =& J_{xyz}\left(\nabla z \times \nabla x\right) = {\boldsymbol{e}}_\theta + \nu{\boldsymbol{e}}_\phi \\
    \boldsymbol{e}_z =& J_{xyz}\left(\nabla x \times \nabla y\right) = {\boldsymbol{e}}_\zeta
    \end{aligned}
- 
+
 where `{\boldsymbol{e}}_\phi =
 {\sigma_{B\text{pol}}}{\boldsymbol{e}}_\zeta` is always anticlockwise when
 seen from above the tokamak looking down. The direction of
@@ -458,7 +458,7 @@ The reciprocal basis vectors are
 
 .. math::
    :label: eq:reciprocalbasisvectorsrighthanded
-   
+
    \begin{aligned}
    \nabla x =& {\sigma_{B\text{pol}}} \nabla \psi \\
    \nabla \eta =& {\sigma_{B\text{pol}}} \nabla \theta \\
@@ -469,7 +469,7 @@ and basis vectors
 
 .. math::
    :label: eq:basisvectorsrighthanded
-   
+
    \begin{aligned}
    \boldsymbol{e}_x =& J_{x\eta z}\left(\nabla y \times \nabla z\right) = {\sigma_{B\text{pol}}} {\boldsymbol{e}}_\psi + I{\boldsymbol{e}}_\zeta \\
    \boldsymbol{e}_\eta =& J_{x\eta z}\left(\nabla z \times \nabla x\right) = {\sigma_{B\text{pol}}} {\boldsymbol{e}}_\theta + \nu{\boldsymbol{e}}_\zeta \\
@@ -853,7 +853,7 @@ Components of `\nabla\times{\boldsymbol{b}}` are [#curvature]_:
 
    \begin{aligned}
    \left(\nabla\times{\boldsymbol{b}}\right)^x =& {\sigma_y}\frac{{B_{\text{pol}}}}{{h_\theta}}{\frac{\partial }{\partial y}}\left(\frac{{B_{\text{tor}}}
-   R}{B}\right) \\ 
+   R}{B}\right) \\
    \left(\nabla\times{\boldsymbol{b}}\right)^y =& -{\sigma_y}\frac{{B_{\text{pol}}}}{{h_\theta}}{\frac{\partial }{\partial x}}\left(\frac{{B_{\text{tor}}}R}{B}\right) \\
    \left(\nabla\times{\boldsymbol{b}}\right)^z =& \frac{{B_{\text{pol}}}}{{h_\theta}}{\frac{\partial }{\partial x}}\left(\frac{B{h_\theta}}{{B_{\text{pol}}}}\right) - \frac{{B_{\text{pol}}}{B_{\text{tor}}} R}{{h_\theta}B}{\frac{\partial \nu}{\partial x}} - {\sigma_y}I\frac{{|B_{\text{pol}}|}}{{h_\theta}}{\frac{\partial }{\partial y}}\left(\frac{{B_{\text{tor}}} R}{B}\right) \\
    \end{aligned}
@@ -919,7 +919,7 @@ we can re-write the above components as:
                 flip, the `x`-coordinate stays the same and the
                 `z`-coordinate flips sign. The `y`-coordinate stays
                 the same, or the `\eta`-coordinate flips sign.
-                Therefore the x-component of `\nabla\times\boldsymbol{b}` 
+                Therefore the x-component of `\nabla\times\boldsymbol{b}`
                 should flip sign, the `z`-component should not flip
                 sign (product of two sign flips), and the
                 `y`-component should flip sign if '`y`' is `y` and not
@@ -1383,7 +1383,7 @@ and
            \frac{{B_{\text{pol}}}}{B{h_\theta}}{\frac{\partial }{\partial y}}\left(\frac{{B_{\text{tor}}}R}{B}\right)\end{aligned}
 
 The second and third terms partly cancel, and using
-`\sigma_y\sigma_{B\text{pol}}{\frac{\partial I}{\partial y}} = 
+`\sigma_y\sigma_{B\text{pol}}{\frac{\partial I}{\partial y}} =
 {\frac{\partial \nu}{\partial x}}`
 
 .. math::
@@ -1750,10 +1750,10 @@ Differential geometry
 =====================
 
 .. warning:: The following are notes from [haeseleer]_. If you are new to this
-	     topic it is strongly suggested to read [haeseleer]_ chapter 2
-	     instead, as here not all terms are defined, and the discussion of
-	     co- and contra variant components is incomplete. Similiarly, the
-	     notation is based on [haeseleer]_ and not explained in detail.
+             topic it is strongly suggested to read [haeseleer]_ chapter 2
+             instead, as here not all terms are defined, and the discussion of
+             co- and contra variant components is incomplete. Similiarly, the
+             notation is based on [haeseleer]_ and not explained in detail.
 
 Sets of vectors `\left\{\mathbf{A, B, C}\right\}` and
 `\left\{\mathbf{a, b, c}\right\}` are reciprocal if
@@ -2229,8 +2229,9 @@ The perpendicular Laplacian can therefore be written in divergence form as:
    +& \frac{1}{J}\frac{\partial}{\partial z}\left[Jg^{xz}\left(\frac{\partial}{\partial x} - \frac{g_{xy}}{g_{yy}}\frac{\partial}{\partial y}\right) + Jg^{zz}\left(\frac{\partial}{\partial z} - \frac{g_{yz}}{g_{yy}}\frac{\partial}{\partial y}\right)\right]
    \end{aligned}
 
-This form is currently implemented in ``FV::Div_a_Grad_perp``
-(``bout/fv_ops.hxx``) but that operator currently assumes that the
+This form is currently implemented in ``FV::Div_a_Grad_perp`` and the
+limited variant ``FV::Div_a_Grad_perp_limit`` (both in
+``bout/fv_ops.hxx``), but these operators currently assume that the
 off-diagonal terms `g^{xz}` and `g^{xy}` are zero, which is the case
 for orthogonal grids with shifted metrics but not in general.
 
diff --git a/manual/sphinx/user_docs/differential_operators.rst b/manual/sphinx/user_docs/differential_operators.rst
index bcc04dcc29..e3a2d03d79 100644
--- a/manual/sphinx/user_docs/differential_operators.rst
+++ b/manual/sphinx/user_docs/differential_operators.rst
@@ -613,6 +613,65 @@ reconstructed values of :math:`f` and :math:`v` remain consistent with the
 other finite-volume advection operators.
 
 
+Parallel momentum-flux heating ``Div_par_fvv_heating``
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+This operator is a companion to ``FV::Div_par_fvv``. It calculates the heating
+associated with the numerical momentum-flux discretisation, and returns a
+diagnostic for the kinetic-energy flow through the lower :math:`y` cell face:
+
+::
+
+   template<typename CellEdges = MC>
+   Field3D Div_par_fvv_heating(const Field3D &f_in, const Field3D &v_in,
+                               const Field3D &a, Field3D &flow_ylow,
+                               bool fixflux=true);
+
+
+The returned field is intended for use in parallel momentum and energy systems
+where ``FV::Div_par_fvv`` is used in the momentum equation. The extra output
+argument ``flow_ylow`` stores the lower-face kinetic-energy flow, including the
+area factor. For FCI fields this diagnostic is currently returned as zero.
+
+As with the other finite-volume parallel operators, ``fixflux=true`` adjusts
+boundary fluxes to remain consistent with the imposed boundary values.
+
+
+Perpendicular finite-volume operators
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The finite-volume operators in ``bout/fv_ops.hxx`` also include conservative
+perpendicular operators. The baseline operator
+``FV::Div_a_Grad_perp(a, f)`` calculates
+:math:`\nabla \cdot \left( a \nabla_\perp f \right)` in flux-conservative form:
+
+::
+
+   Field3D Div_a_Grad_perp(const Field3D &a, const Field3D &f);
+
+
+There is also a limited variant:
+
+::
+
+   template<typename CellEdges = MC>
+   Field3D Div_a_Grad_perp_limit(const Field3D &a, const Field3D &g,
+                                 const Field3D &f);
+
+
+This calculates :math:`\nabla \cdot \left( a g \nabla_\perp f \right)` using a
+limited reconstruction of :math:`g` on cell faces. The arguments are:
+
+* ``a`` - the transport coefficient
+* ``g`` - the factor reconstructed and upwinded at cell faces
+* ``f`` - the field whose perpendicular gradient is used
+
+In the :math:`x` direction, the chosen ``CellEdges`` limiter is used to
+reconstruct cell-edge values of :math:`g`. In the :math:`y` direction,
+first-order upwinding is used. As with ``FV::Div_par``, the template parameter
+selects the limiter family described in :ref:`sec-slope-limiters`.
+
+
 Example and convergence test
 ++++++++++++++++++++++++++++
 

From 0f34609247361f6857fda40fdd68ea5dd309792d Mon Sep 17 00:00:00 2001
From: David Bold <dave@ipp.mpg.de>
Date: Mon, 6 Jul 2026 15:29:33 +0200
Subject: [PATCH 10/11] Prefer BoutNaN

---
 include/bout/fv_ops_impl.hxx | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/include/bout/fv_ops_impl.hxx b/include/bout/fv_ops_impl.hxx
index 7ab4c2a042..1b7a36f874 100644
--- a/include/bout/fv_ops_impl.hxx
+++ b/include/bout/fv_ops_impl.hxx
@@ -383,7 +383,7 @@ Field3D Div_par(const Field3D& f_in, const Field3D& v_in, const Field3D& wave_sp
 
         // Calculate velocity at right boundary (y+1/2)
         BoutReal vpar = 0.5 * (v(i, j, k) + v(i, j + 1, k));
-        BoutReal flux = NAN;
+        BoutReal flux = BoutNaN;
 
         if (is_last_y && (j == mesh->yend) && !is_periodic_y) {
           // Last point in domain
@@ -509,7 +509,7 @@ Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux) {
     if ((i.x() == mesh->xend) && (mesh->lastX())) {
       // At right boundary in X
       if (bndry_flux) {
-        BoutReal flux = NAN;
+        BoutReal flux = BoutNaN;
         if (vR > 0.0) {
           // Flux to boundary
           flux = vR * s.R;
@@ -536,7 +536,7 @@ Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux) {
       // At left boundary in X
 
       if (bndry_flux) {
-        BoutReal flux = NAN;
+        BoutReal flux = BoutNaN;
         if (vL < 0.0) {
           // Flux to boundary
           flux = vL * s.L;
@@ -935,7 +935,7 @@ Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in,
         const BoutReal v_mid_r = 0.5 * (sv.c + sv.p);
         // And mid-point density at right boundary
         const BoutReal n_mid_r = 0.5 * (s.c + s.p);
-        BoutReal flux = NAN;
+        BoutReal flux = BoutNaN;
 
         if (mesh->lastY(i) && (j == mesh->yend) && !mesh->periodicY(i)) {
           // Last point in domain
@@ -1156,7 +1156,7 @@ Field3D Div_par_fvv_heating(const Field3D& f_in, const Field3D& v_in,
           // energy losses.
           const BoutReal expected_ke = 0.5 * n_mid * v_mid * v_mid * v_mid;
 
-          BoutReal flux_mom = NAN;
+          BoutReal flux_mom = BoutNaN;
           if (fixflux) {
             // Mid-point consistent with boundary conditions
             // but kinetic energy loss will not match expected
@@ -1211,7 +1211,7 @@ Field3D Div_par_fvv_heating(const Field3D& f_in, const Field3D& v_in,
 
         if (mesh->firstY(i) && (j == mesh->ystart) && !mesh->periodicY(i)) {
           // First point in domain
-          BoutReal flux_mom = NAN;
+          BoutReal flux_mom = BoutNaN;
           if (fixflux) {
             // Use mid-point to be consistent with boundary conditions
             flux_mom = n_mid * v_mid * v_mid;
@@ -1288,7 +1288,7 @@ Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, const Field3D&
 
         // Mid-point average boundary value of 'a'
         const BoutReal aedge = 0.5 * (a(i + 1, j, k) + a(i, j, k));
-        BoutReal gedge = NAN;
+        BoutReal gedge = BoutNaN;
         if (((i == mesh->xstart - 1) and mesh->firstX())
             or ((i == mesh->xend) and mesh->lastX())) {
           // Mid-point average boundary value of 'g'
@@ -1410,7 +1410,7 @@ Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, const Field3D&
                           / (dy.up(i, j + 1, k) + dy.c(i, j, k));
 
           BoutReal aedge = 0.5 * (a_slice.c(i, j, k) + a_slice.up(i, j + 1, k));
-          BoutReal gedge = NAN;
+          BoutReal gedge = BoutNaN;
           if ((j == mesh->yend) and mesh->lastY(i)) {
             // Midpoint boundary value
             gedge = 0.5 * (g_slice.c(i, j, k) + g_slice.up(i, j + 1, k));

From 7f8adf5b3f461aba2903590669349a283b8246e7 Mon Sep 17 00:00:00 2001
From: David Bold <dave@ipp.mpg.de>
Date: Wed, 8 Jul 2026 10:56:05 +0200
Subject: [PATCH 11/11] Remove unused headers

---
 include/bout/fv_ops.hxx | 6 ------
 1 file changed, 6 deletions(-)

diff --git a/include/bout/fv_ops.hxx b/include/bout/fv_ops.hxx
index 17b4748ae4..2ae0d3d547 100644
--- a/include/bout/fv_ops.hxx
+++ b/include/bout/fv_ops.hxx
@@ -5,17 +5,11 @@
 #ifndef BOUT_FV_OPS_H
 #define BOUT_FV_OPS_H
 
-#include "bout/assert.hxx"
 #include "bout/bout_types.hxx"
-#include "bout/build_defines.hxx"
-#include "bout/coordinates.hxx"
 #include "bout/field.hxx"
 #include "bout/field3d.hxx"
-#include "bout/globals.hxx"
 #include "bout/mesh.hxx"
 #include "bout/output_bout_types.hxx" // NOLINT(unused-includes, misc-include-cleaner)
-#include "bout/region.hxx"
-#include "bout/utils.hxx"
 #include "bout/vector2d.hxx"
 
 namespace FV {