Implement 2D neutral equation fixes originally part of AFN flux limiter branch

include/div_ops.hxx

@@ -23,12 +23,12 @@
 
  */
 
-#ifndef __DIV_OPS_H__
-#define __DIV_OPS_H__
+#ifndef HERMES_DIV_OPS_H
+#define HERMES_DIV_OPS_H
 
 #include <bout/field3d.hxx>
-#include <bout/vector3d.hxx>
 #include <bout/fv_ops.hxx>
+#include <bout/vector3d.hxx>
 
 /*!
  * Diffusion in index space
@@ -58,7 +58,8 @@ const Field2D Laplace_FV(const Field2D& k, const Field2D& f);
 /// Perpendicular diffusion including X and Y directions
 const Field3D Div_a_Grad_perp_upwind(const Field3D& a, const Field3D& f);
 /// Version of function that returns flows
-const Field3D Div_a_Grad_perp_upwind_flows(const Field3D& a, const Field3D& f, Field3D& flux_xlow, Field3D& flux_ylow);
+const Field3D Div_a_Grad_perp_upwind_flows(const Field3D& a, const Field3D& f,
+                                           Field3D& flux_xlow, Field3D& flux_ylow);
 
 namespace FV {
 
@@ -100,8 +101,7 @@ struct Superbee {
       BoutReal sign = SIGN(gL);
       gL = fabs(gL);
       gR = fabs(gR);
-      BoutReal half_slope = sign * BOUTMAX(BOUTMIN(gL, 0.5*gR),
-                                           BOUTMIN(gR, 0.5*gL));
+      BoutReal half_slope = sign * BOUTMAX(BOUTMIN(gL, 0.5 * gR), BOUTMIN(gR, 0.5 * gL));
       n.L = n.c - half_slope;
       n.R = n.c + half_slope;
     }
@@ -210,10 +210,9 @@ const Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in,
             flux = bndryval * vpar * vpar;
           } else {
             // Add flux due to difference in boundary values
-            flux =
-                s.R * vpar * sv.R
-                + BOUTMAX(wave_speed(i, j, k), fabs(v(i, j, k)), fabs(v(i, j + 1, k)))
-                  * (s.R * sv.R - bndryval * vpar);
+            flux = s.R * vpar * sv.R
+                   + BOUTMAX(wave_speed(i, j, k), fabs(v(i, j, k)), fabs(v(i, j + 1, k)))
+                         * (s.R * sv.R - bndryval * vpar);
           }
         } else {
           // Maximum wave speed in the two cells
@@ -239,10 +238,9 @@ const Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in,
             flux = bndryval * vpar * vpar;
           } else {
             // Add flux due to difference in boundary values
-            flux =
-                s.L * vpar * sv.L
-                - BOUTMAX(wave_speed(i, j, k), fabs(v(i, j, k)), fabs(v(i, j - 1, k)))
-                  * (s.L * sv.L - bndryval * vpar);
+            flux = s.L * vpar * sv.L
+                   - BOUTMAX(wave_speed(i, j, k), fabs(v(i, j, k)), fabs(v(i, j - 1, k)))
+                         * (s.L * sv.L - bndryval * vpar);
           }
         } else {
           // Maximum wave speed in the two cells
@@ -454,6 +452,235 @@ const Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in,
   return are_unaligned ? fromFieldAligned(result, "RGN_NOBNDRY") : result;
 }
 
+/// 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) * coord->g11(i, j)
+                                 + coord->J(i + 1, j) * coord->g11(i + 1, j))
+                              / (coord->dx(i, j) + coord->dx(i + 1, j));
+
+        result(i, j, k) += fout / (coord->dx(i, j) * coord->J(i, j));
+        result(i + 1, j, k) -= fout / (coord->dx(i + 1, j) * coord->J(i + 1, 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 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++) {
+
+      BoutReal coef_u =
+          0.5
+          * (coord->g_23(i, j) / SQ(coord->J(i, j) * coord->Bxy(i, j))
+             + coord->g_23(i, j + 1) / SQ(coord->J(i, j + 1) * coord->Bxy(i, j + 1)));
+
+      BoutReal coef_d =
+          0.5
+          * (coord->g_23(i, j) / SQ(coord->J(i, j) * coord->Bxy(i, j))
+             + coord->g_23(i, j - 1) / SQ(coord->J(i, j - 1) * coord->Bxy(i, j - 1)));
+
+      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;
+
+        // 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);
+
+        // Y derivative
+        BoutReal dfdy = 2. * (fup(i, j + 1, k) - fc(i, j, k))
+                        / (coord->dy(i, j + 1) + coord->dy(i, j));
+
+        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) * coord->g23(i, j)
+                           + coord->J(i, j + 1) * coord->g23(i, j + 1))
+                        * (dfdz - coef_u * dfdy);
+
+        yzresult(i, j, k) = fout / (coord->dy(i, j) * coord->J(i, j));
+
+        // 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);
+
+        dfdy = 2. * (fc(i, j, k) - fdown(i, j - 1, k))
+               / (coord->dy(i, j) + coord->dy(i, j - 1));
+
+        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) * coord->g23(i, j)
+                  + coord->J(i, j - 1) * coord->g23(i, j - 1))
+               * (dfdz - coef_d * dfdy);
+
+        yzresult(i, j, k) -= fout / (coord->dy(i, j) * coord->J(i, j));
+      }
+    }
+  }
+
+  // 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++) {
+      // Coefficient in front of df/dy term
+      BoutReal coef = coord->g_23(i, j)
+                      / (coord->dy(i, j + 1) + 2. * coord->dy(i, j) + coord->dy(i, j - 1))
+                      / SQ(coord->J(i, j) * coord->Bxy(i, j));
+
+      for (int k = 0; k < mesh->LocalNz; k++) {
+        // 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)
+
+            // - 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);
+        yzresult(i, j, kp) -= fout / coord->dz(i, j);
+      }
+    }
+  }
+  // Check if we need to transform back
+  if (f.hasParallelSlices() && a.hasParallelSlices()) {
+    result += yzresult;
+  } else {
+    result += fromFieldAligned(yzresult);
+  }
+
+  return result;
+}
+
 } // namespace FV
 
-#endif //  __DIV_OPS_H__
+#endif //  HERMES_DIV_OPS_H

include/neutral_mixed.hxx

@@ -42,26 +42,40 @@ private:
   BoutReal AA; ///< Atomic mass (proton = 1)
 
   Field3D Dnn; ///< Diffusion coefficient
-  Field3D DnnNn, DnnPn, DnnTn, DnnNVn; ///< Used for operators
+  Field3D DnnNn, DnnPn, DnnNVn;
 
   bool sheath_ydown, sheath_yup;
 
   BoutReal nn_floor; ///< Minimum Nn used when dividing NVn by Nn to get Vn.
+  BoutReal pn_floor; ///< Minimum Pn used when dividing Pn by Nn to get Tn.
 
   BoutReal flux_limit; ///< Diffusive flux limit
   BoutReal diffusion_limit;    ///< Maximum diffusion coefficient
 
   bool neutral_viscosity; ///< include viscosity?
+  bool neutral_conduction; ///< Include heat conduction?
   bool evolve_momentum; ///< Evolve parallel momentum?
+
+  Field3D kappa_n, eta_n; ///< Neutral conduction and viscosity
 
   bool precondition {true}; ///< Enable preconditioner?
+  bool lax_flux; ///< Use Lax flux for advection terms
   std::unique_ptr<Laplacian> inv; ///< Laplacian inversion used for preconditioning
 
   Field3D density_source, pressure_source; ///< External input source
   Field3D Sn, Sp, Snv; ///< Particle, pressure and momentum source
+  Field3D sound_speed; ///< Sound speed for use with Lax flux
+  Field3D perp_nn_adv_src; ///< Source due to perpendicular advection operator
+  Field3D par_nn_adv_src; ///< Source due to parallel advection operator
 
   bool output_ddt; ///< Save time derivatives?
   bool diagnose; ///< Save additional diagnostics?
+
+  // Flow diagnostics
+  Field3D particle_flow_xlow, particle_flow_ylow;
+  Field3D momentum_flow_xlow, momentum_flow_ylow;
+  Field3D energy_flow_xlow, energy_flow_ylow;
+  Field3D conduction_flow_xlow, conduction_flow_ylow;
 };
 
 namespace {