Revert to Div_a_Grad_perp, add flows

Upwinding and limiting with MC/MinMod are causing checkerboarding and/or crashes in the neutral_mixed example and provide seemingly no benefit in production runs. Reverting to original operator but with added flow diagnostics. Left warnings on the other operators, maybe we'll get back to them at some point.

Squashed with commit to update operator choice.

include/div_ops.hxx

@@ -61,8 +61,14 @@ const Field3D D4DZ4_Index(const Field3D& f);
 const Field2D Laplace_FV(const Field2D& k, const Field2D& f);
 
 /// Perpendicular diffusion including X and Y directions
+/// Takes Div_a_Grad_perp from BOUT++ and adds flows
+const Field3D Div_a_Grad_perp_flows(const Field3D& a, const Field3D& f,
+                                           Field3D& flux_xlow, Field3D& flux_ylow);
+/// Same but with upwinding
+/// WARNING: Causes checkerboarding in neutral_mixed integrated test
 const Field3D Div_a_Grad_perp_upwind(const Field3D& a, const Field3D& f);
-/// Version of function that returns flows
+/// Same but with upwinding and flows
+/// WARNING: Causes checkerboarding in neutral_mixed integrated test
 const Field3D Div_a_Grad_perp_upwind_flows(const Field3D& a, const Field3D& f,
                                            Field3D& flux_xlow, Field3D& flux_ylow);
 

src/div_ops.cxx

@@ -792,7 +792,237 @@ const Field2D Laplace_FV(const Field2D &k, const Field2D &f) {
   return result;
 }
 
+const Field3D Div_a_Grad_perp_flows(const Field3D& a, const Field3D& f,
+                              Field3D &flow_xlow, 
+                              Field3D &flow_ylow) {
+  ASSERT2(a.getLocation() == f.getLocation());
+
+  Mesh* mesh = a.getMesh();
+
+  Field3D result{zeroFrom(f)};
+
+  Coordinates* coord = f.getCoordinates();
+
+  // Zero all flows
+  flow_xlow = 0.0;
+  flow_ylow = 0.0;
+
+  // Flux in x
+
+  int xs = mesh->xstart - 1;
+  int xe = mesh->xend;
+
+  /*
+    if(mesh->firstX())
+    xs += 1;
+  */
+  /*
+    if(mesh->lastX())
+    xe -= 1;
+  */
+
+  for (int i = xs; i <= xe; 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
+
+        BoutReal fout = 0.5 * (a(i, j, k) + a(i + 1, j, k))
+                        * (coord->J(i, j, k) * coord->g11(i, j, k)
+                           + coord->J(i + 1, j, k) * coord->g11(i + 1, j, k))
+                        * (f(i + 1, j, k) - f(i, 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));
+
+        // Flow will be positive in the positive coordinate direction
+        flow_xlow(i + 1, j, k) = -1.0 * fout * coord->dy(i, j) * coord->dz(i, j);
+      }
+    }
+  }
+
+  // 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 g23up(mesh), g23down(mesh);
+  Field3D g_23up(mesh), g_23down(mesh);
+  Field3D Jup(mesh), Jdown(mesh);
+  Field3D dyup(mesh), dydown(mesh);
+  Field3D dzup(mesh), dzdown(mesh);
+  Field3D Bxyup(mesh), Bxydown(mesh);
+
+  // Values on this y slice (centre).
+  // This is needed because toFieldAligned may modify the field
+  Field3D fc = f;
+  Field3D ac = a;
+
+  Field3D g23c = coord->g23;
+  Field3D g_23c = coord->g_23;
+  Field3D Jc = coord->J;
+  Field3D dyc = coord->dy;
+  Field3D dzc = coord->dz;
+  Field3D Bxyc = coord->Bxy;
+
+  // Result of the Y and Z fluxes
+  Field3D yzresult(mesh);
+  yzresult.allocate();
+
+  if (f.hasParallelSlices() && a.hasParallelSlices()) {
+    // Both inputs have yup and ydown
+
+    fup = f.yup();
+    fdown = f.ydown();
+
+    aup = a.yup();
+    adown = a.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);
+
+    yzresult.setDirectionY(YDirectionType::Aligned);
+  }
+
+  if (bout::build::use_metric_3d) {
+    // 3D Metric, need yup/ydown fields.
+    // Requires previous communication of metrics
+    // -- should insert communication here?
+    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?");
+    }
+
+    g23up = coord->g23.yup();
+    g23down = coord->g23.ydown();
+
+    g_23up = coord->g_23.yup();
+    g_23down = coord->g_23.ydown();
+
+    Jup = coord->J.yup();
+    Jdown = coord->J.ydown();
+
+    dyup = coord->dy.yup();
+    dydown = coord->dy.ydown();
+
+    dzup = coord->dz.yup();
+    dzdown = coord->dz.ydown();
+
+    Bxyup = coord->Bxy.yup();
+    Bxydown = coord->Bxy.ydown();
+
+  } else {
+    // No 3D metrics
+    // Need to shift to/from field aligned coordinates
+    g23up = g23down = g23c = toFieldAligned(coord->g23);
+    g_23up = g_23down = g_23c = toFieldAligned(coord->g_23);
+    Jup = Jdown = Jc = toFieldAligned(coord->J);
+    dyup = dydown = dyc = toFieldAligned(coord->dy);
+    dzup = dzdown = dzc = toFieldAligned(coord->dz);
+    Bxyup = Bxydown = Bxyc = toFieldAligned(coord->Bxy);
+  }
+
+  // Y flux
+
+  for (int i = mesh->xstart; i <= mesh->xend; i++) {
+    for (int j = mesh->ystart; j <= mesh->yend; j++) {
+      for (int k = 0; k < mesh->LocalNz; k++) {
+        // Calculate flux between j and j+1
+        int kp = (k + 1) % mesh->LocalNz;
+        int km = (k - 1 + mesh->LocalNz) % mesh->LocalNz;
+
+        BoutReal coef =
+            0.5
+            * (g_23c(i, j, k) / SQ(Jc(i, j, k) * Bxyc(i, j, k))
+               + g_23up(i, j + 1, k) / SQ(Jup(i, j + 1, k) * Bxyup(i, j + 1, k)));
+
+        // Calculate Z derivative at y boundary
+        BoutReal dfdz =
+            0.5 * (fc(i, j, kp) - fc(i, j, km) + fup(i, j + 1, kp) - fup(i, j + 1, km))
+            / (dzc(i, j, k) + dzup(i, j + 1, k));
+
+        // Y derivative
+        BoutReal dfdy =
+            2. * (fup(i, j + 1, k) - fc(i, j, k)) / (dyup(i, j + 1, k) + dyc(i, j, k));
+
+        BoutReal fout =
+            0.25 * (ac(i, j, k) + aup(i, j + 1, k))
+            * (Jc(i, j, k) * g23c(i, j, k) + Jup(i, j + 1, k) * g23up(i, j + 1, k))
+            * (dfdz - coef * dfdy);
+
+        yzresult(i, j, k) = fout / (dyc(i, j, k) * Jc(i, j, k));
+
+        // Calculate flux between j and j-1
+        coef =
+            0.5
+            * (g_23c(i, j, k) / SQ(Jc(i, j, k) * Bxyc(i, j, k))
+               + g_23down(i, j - 1, k) / SQ(Jdown(i, j - 1, k) * Bxydown(i, j - 1, k)));
+
+        dfdz = 0.5
+               * (fc(i, j, kp) - fc(i, j, km) + fdown(i, j - 1, kp) - fdown(i, j - 1, km))
+               / (dzc(i, j, k) + dzdown(i, j - 1, k));
+
+        dfdy = 2. * (fc(i, j, k) - fdown(i, j - 1, k))
+               / (dyc(i, j, k) + dydown(i, j - 1, k));
+
+        fout = 0.25 * (ac(i, j, k) + adown(i, j - 1, k))
+               * (Jc(i, j, k) * g23c(i, j, k) + Jdown(i, j - 1, k) * g23down(i, j - 1, k))
+               * (dfdz - coef * dfdy);
+
+        yzresult(i, j, k) -= fout / (dyc(i, j, k) * Jc(i, j, k));
+
+        // Flow will be positive in the positive coordinate direction
+        flow_ylow(i, j, k) = -1.0 * fout * coord->dx(i, j) * coord->dz(i, j);
+      }
+    }
+  }
+
+  // Z flux
+
+  for (int i = mesh->xstart; i <= mesh->xend; i++) {
+    for (int j = mesh->ystart; j <= mesh->yend; j++) {
+      for (int k = 0; k < mesh->LocalNz; k++) {
+        // Calculate flux between k and k+1
+        int kp = (k + 1) % mesh->LocalNz;
+
+        // Coefficient in front of df/dy term
+        BoutReal coef = g_23c(i, j, k)
+                        / (dyup(i, j + 1, k) + 2. * dyc(i, j, k) + dydown(i, j - 1, k))
+                        / SQ(Jc(i, j, k) * Bxyc(i, j, k));
+
+        BoutReal fout =
+            0.25 * (ac(i, j, k) + ac(i, j, kp))
+            * (Jc(i, j, k) * coord->g33(i, j, k) + Jc(i, j, kp) * coord->g33(i, j, kp))
+            * ( // df/dz
+                (fc(i, j, kp) - fc(i, j, k)) / dzc(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)));
+
+        yzresult(i, j, k) += fout / (Jc(i, j, k) * dzc(i, j, k));
+        yzresult(i, j, kp) -= fout / (Jc(i, j, kp) * dzc(i, j, kp));
+      }
+    }
+  }
+  // Check if we need to transform back
+  if (f.hasParallelSlices() && a.hasParallelSlices()) {
+    result += yzresult;
+  } else {
+    result += fromFieldAligned(yzresult);
+    flow_ylow = fromFieldAligned(flow_ylow);
+  }
+
+  return result;
+}
+
 // Div ( a Grad_perp(f) )  -- diffusion
+/// WARNING: Causes checkerboarding in neutral_mixed integrated test
 const Field3D Div_a_Grad_perp_upwind(const Field3D& a, const Field3D& f) {
   ASSERT2(a.getLocation() == f.getLocation());
 
@@ -969,6 +1199,7 @@ const Field3D Div_a_Grad_perp_upwind(const Field3D& a, const Field3D& f) {
 ///           ylow(i,j)
 ///
 ///
+/// WARNING: Causes checkerboarding in neutral_mixed integrated test
 const Field3D Div_a_Grad_perp_upwind_flows(const Field3D& a, const Field3D& f,
                                            Field3D &flow_xlow,
                                            Field3D &flow_ylow) {

src/neutral_mixed.cxx

@@ -350,7 +350,7 @@ void NeutralMixed::finally(const Options& state) {
   // Neutral density
   TRACE("Neutral density");
 
-  perp_nn_adv_src = Div_a_Grad_perp_upwind_flows(DnnNn, logPnlim,
+  perp_nn_adv_src = Div_a_Grad_perp_flows(DnnNn, logPnlim,
                                    particle_flow_xlow,
                                    particle_flow_ylow); // Perpendicular advection
 
@@ -376,7 +376,7 @@ void NeutralMixed::finally(const Options& state) {
     ddt(NVn) =
         -AA * FV::Div_par_fvv<ParLimiter>(Nnlim, Vn, sound_speed) // Momentum flow
         - Grad_par(Pn) // Pressure gradient
-      + Div_a_Grad_perp_upwind_flows(DnnNVn, logPnlim,
+      + Div_a_Grad_perp_flows(DnnNVn, logPnlim,
                                      momentum_flow_xlow,
                                      momentum_flow_ylow) // Perpendicular advection
       ;
@@ -414,17 +414,16 @@ void NeutralMixed::finally(const Options& state) {
 
   ddt(Pn) = - FV::Div_par_mod<ParLimiter>(Pn, Vn, sound_speed) // Parallel advection
             - (2. / 3) * Pn * Div_par(Vn)                      // Compression
-    + Div_a_Grad_perp_upwind_flows(DnnPn, logPnlim,
+    + Div_a_Grad_perp_flows(DnnPn, logPnlim,
                                    energy_flow_xlow, energy_flow_ylow) // Perpendicular advection
      ;
 
-  // The factor here is likely 5/2 as we're advecting internal energy and pressure.
-  // Doing this still leaves a heat imbalance factor of 0.11 in the cells, but better than 0.33 with 3/2.
+  // The factor here is 5/2 as we're advecting internal energy and pressure.
   energy_flow_xlow *= 5/2; 
   energy_flow_ylow *= 5/2;
 
   if (neutral_conduction) {
-    ddt(Pn) += Div_a_Grad_perp_upwind_flows(DnnNn, Tn,
+    ddt(Pn) += Div_a_Grad_perp_flows(DnnNn, Tn,
                         conduction_flow_xlow, conduction_flow_ylow)    // Perpendicular conduction
       + FV::Div_par_K_Grad_par(DnnNn, Tn)        // Parallel conduction
       ;