Fixes to FV operators

include/bout/fv_ops_impl.hxx

@@ -1170,7 +1170,7 @@ Field3D Div_par_fvv_heating(const Field3D& f_in, const Field3D& v_in,
           } else {
             flux_mom = (s.R * sv.R * sv.R)
                        + (BOUTMAX(wave_speed(i, j, k), std::abs(sv.c), std::abs(sv.p))
-                          * (s.R * sv.R - n_mid * v_mid));
+                          * n_mid * (sv.R - v_mid));
           }
 
           // Assume that particle flux is fixed to boundary value
@@ -1191,14 +1191,16 @@ Field3D Div_par_fvv_heating(const Field3D& f_in, const Field3D& v_in,
           const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j + 1, k),
                                         std::abs(sv.c), std::abs(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;
+          // Use the same split particle flux as Div_par_mod and Div_par_fvv.
+          const BoutReal flux_part = s.R * 0.5 * (sv.R + amax);
+
+          // Difference between kinetic energy leaving this cell and entering
+          // the next cell, using d(n*v^2/2)/dt = v*d(n*v)/dt - v^2*dn/dt/2.
+          result(i, j, k) += flux_part * (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) +=
@@ -1224,7 +1226,7 @@ Field3D Div_par_fvv_heating(const Field3D& f_in, const Field3D& v_in,
             // Add flux due to difference in boundary values
             flux_mom = (s.L * sv.L * sv.L)
                        - (BOUTMAX(wave_speed(i, j, k), std::abs(sv.c), std::abs(sv.m))
-                          * (s.L * sv.L - n_mid * v_mid));
+                          * n_mid * (sv.L - v_mid));
           }
 
           // Assume that density flux is fixed to boundary value
@@ -1241,17 +1243,17 @@ Field3D Div_par_fvv_heating(const Field3D& f_in, const Field3D& v_in,
           const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j - 1, k),
                                         std::abs(sv.c), std::abs(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;
+          const BoutReal flux_part = s.L * 0.5 * (sv.L - amax);
 
-          // Kinetic energy flow into this cell.
-          // Note: Different from flow out of left cell; the difference
+          // The left-going split flux has the opposite orientation.
+          result(i, j, k) -= flux_part * (sv.c - sv.m) * (sv.L - v_mid) * flux_factor_lc;
+
+          // Kinetic energy flow into the previous cell.
+          // Note: Different from flow out of this 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(i, j, k) += (sv.m * flux_mom - 0.5 * SQ(sv.m) * flux_part) * area_lc;
         }
       }
     }

src/mesh/fv_ops.cxx

@@ -84,7 +84,11 @@ Field3D Div_a_Grad_perp(const Field3D& a, const Field3D& f) {
     }
   }
 
-  // Y and Z fluxes require Y derivatives
+  // A singleton global Y mesh has no Y derivatives, even if guards exist.
+  const bool y_derivatives = mesh->GlobalNy > 1;
+  ASSERT1(!y_derivatives || mesh->ystart > 0);
+
+  // Y fluxes and the mixed part of the Z flux require Y derivatives.
 
   // Fields containing values along the magnetic field
 
@@ -109,7 +113,7 @@ Field3D Div_a_Grad_perp(const Field3D& a, const Field3D& f) {
   }
 
   // Y flux
-  if (mesh->ystart > 0) {
+  if (y_derivatives) {
     BOUT_FOR(i, yzresult.getRegion("RGN_NOBNDRY")) {
       auto ikp = i.zp();
       auto ikm = i.zm();
@@ -155,29 +159,27 @@ Field3D Div_a_Grad_perp(const Field3D& a, const Field3D& f) {
 
         yzresult[i] -= fout / (dy.c[i] * J.c[i]);
       }
-      // Z flux
-      {
+    }
+  }
 
-        // Coefficient in front of df/dy term
-        const BoutReal coef = g_23.c[i] / (dy.up[i.yp()] + 2. * dy.c[i] + dy.down[i.ym()])
-                              / SQ(J.c[i] * Bxy.c[i]);
-
-        const BoutReal fout =
-            0.25 * (a_slice.c[i] + a_slice.c[ikp])
-            * (J.c[i] * coord->g33()[i] + J.c[ikp] * coord->g33()[ikp])
-            * ( // df/dz
-                (f_slice.c[ikp] - f_slice.c[i]) / dz.c[i]
-                // - g_yz * df/dy / SQ(J*B)
-                - coef
-                      * (f_slice.up[i.yp()] + f_slice.up[ikp.yp()] - f_slice.down[i.ym()]
-                         - f_slice.down[ikp.ym()]));
-
-        yzresult[i] += fout / (J.c[i] * dz.c[i]);
-        yzresult[ikp] -= fout / (J.c[ikp] * dz.c[ikp]);
-      }
+  // Z flux exists independently of Y guards or Y derivatives.
+  BOUT_FOR(i, yzresult.getRegion("RGN_NOBNDRY")) {
+    const auto ikp = i.zp();
+    BoutReal y_correction = 0.0;
+    if (y_derivatives) {
+      const BoutReal coef = g_23.c[i] / (dy.up[i.yp()] + 2. * dy.c[i] + dy.down[i.ym()])
+                            / SQ(J.c[i] * Bxy.c[i]);
+      y_correction = coef
+                     * (f_slice.up[i.yp()] + f_slice.up[ikp.yp()] - f_slice.down[i.ym()]
+                        - f_slice.down[ikp.ym()]);
     }
-  } else {
-    ASSERT1(mesh->yend == 0);
+
+    const BoutReal fout = 0.25 * (a_slice.c[i] + a_slice.c[ikp])
+                          * (J.c[i] * coord->g33()[i] + J.c[ikp] * coord->g33()[ikp])
+                          * ((f_slice.c[ikp] - f_slice.c[i]) / dz.c[i] - y_correction);
+
+    yzresult[i] += fout / (J.c[i] * dz.c[i]);
+    yzresult[ikp] -= fout / (J.c[ikp] * dz.c[ikp]);
   }
 
   // Check if we need to transform back

tests/unit/fake_mesh.hxx

@@ -35,20 +35,20 @@ class Options;
 /// - This is a mesh for a single process, so the global and local
 /// indices are the same.
 ///
-/// - There is a single guard cell at each of the start/end x/y grids.
+/// - There is one X guard cell and, by default, one Y guard cell at each end.
 ///
 /// - Only the **grid** information is assumed to be used -- anything
 ///   else will likely **not** work!
 class FakeMesh : public Mesh {
 public:
-  FakeMesh(int nx, int ny, int nz, MpiWrapper& mpi_in) {
+  FakeMesh(int nx, int ny, int nz, MpiWrapper& mpi_in, int myg = 1) {
     // Mesh only on one process, so global and local indices are the
     // same
     GlobalNx = nx;
     GlobalNy = ny;
     GlobalNz = nz;
     GlobalNxNoBoundaries = nx - 2;
-    GlobalNyNoBoundaries = ny - 2;
+    GlobalNyNoBoundaries = ny - (2 * myg);
     GlobalNzNoBoundaries = nz;
     LocalNx = nx;
     LocalNy = ny;
@@ -59,20 +59,20 @@ public:
 
     // These bits only for ADIOS2, also boring due to single process
     MapCountX = nx - 2;
-    MapCountY = ny - 2;
+    MapCountY = ny - (2 * myg);
     MapCountZ = nz;
     MapGlobalX = nx;
     MapGlobalY = ny;
     MapGlobalZ = nz;
     MapLocalX = nx - 2;
-    MapLocalY = ny - 2;
+    MapLocalY = ny - (2 * myg);
     MapLocalZ = nz;
 
     // Small "inner" region
     xstart = 1;
     xend = nx - 2;
-    ystart = 1;
-    yend = ny - 2;
+    ystart = myg;
+    yend = ny - myg - 1;
     zstart = 0; // no guards
     zend = nz - 1;
 

tests/unit/fake_mesh_fixture.hxx

@@ -34,11 +34,15 @@
 /// Use this template class directly to use different sized grid:
 ///
 ///     using MyTest = FakeMeshFixture_tmpl<7, 9, 11>;
-template <int NX, int NY, int NZ, bool FCI = false>
+///
+/// Use the final template argument to retain two Y guards for reconstruction:
+///
+///     using ParallelTest = FakeMeshFixture_tmpl<5, 20, 3, false, 2>;
+template <int NX, int NY, int NZ, bool FCI = false, int MYG = 1>
 class FakeMeshFixture_tmpl : public ::testing::Test {
 public:
   FakeMeshFixture_tmpl()
-      : mesh_m(NX, NY, NZ, mpi), mesh_staggered_m(NX, NY, NZ, mpi),
+      : mesh_m(NX, NY, NZ, mpi, MYG), mesh_staggered_m(NX, NY, NZ, mpi, MYG),
         mesh_staggered(&mesh_staggered_m) {
 
     bout::globals::mpi = &mpi;
@@ -58,7 +62,8 @@ public:
     test_coords->setD1_dy(0.2);
     test_coords->setD1_dz(0.0);
 
-    if (bout::build::use_metric_3d) {
+    // Parallel slices are only allocated when Y guard cells exist.
+    if (bout::build::use_metric_3d && MYG > 0) {
       bout::FieldMetric mutable_Bxy = test_coords->Bxy();
       mutable_Bxy.splitParallelSlices();
       mutable_Bxy.yup() = test_coords->Bxy();
@@ -100,7 +105,7 @@ public:
     test_coords_staggered->setD1_dy(0.2);
     test_coords_staggered->setD1_dz(0.0);
 
-    if (bout::build::use_metric_3d) {
+    if (bout::build::use_metric_3d && MYG > 0) {
       bout::FieldMetric mutable_Bxy = test_coords_staggered->Bxy();
       mutable_Bxy.splitParallelSlices();
       test_coords_staggered->setBxy(mutable_Bxy);