Merge branch 'master' into E-AFN-neutral-advection

docs/sphinx/code_structure.rst

@@ -54,7 +54,7 @@ for new variables which are added to the state:
     * `charge`  Charge, in units of proton charge (i.e. electron=-1)
 
     * `density`
-    * `momentum`
+    * `momentum` Parallel momentum
     * `pressure`
     * `velocity` Parallel velocity
     * `temperature`
@@ -74,12 +74,14 @@ for new variables which are added to the state:
 * `fields`
 
   * `vorticity`
-  * `phi`       Electrostatic potential
-  * `DivJdia`   Divergence of diamagnetic current
-  * `DivJcol`   Divergence of collisional current
-  * `DivJextra` Divergence of current, including 2D parallel current
-    closures.  Not including diamagnetic, parallel current due to
-    flows, or polarisation currents
+  * `phi`           Electrostatic potential
+  * `Apar`          Electromagnetic potential b dot A in induction terms
+  * `Apar_flutter`  The electromagnetic potential (b dot A) in flutter terms
+  * `DivJdia`       Divergence of diamagnetic current
+  * `DivJcol`       Divergence of collisional current
+  * `DivJextra`     Divergence of current, including 2D parallel current
+                    closures.  Not including diamagnetic, parallel current due to
+                    flows, or polarisation currents
 
 For example to get the electron density::
 

docs/sphinx/components.rst

@@ -2044,11 +2044,19 @@ electromagnetic
 ~~~~~~~~~~~~~~~
 
 **Notes**: When using this module,
+
 1. Set ``sound_speed:alfven_wave=true`` so that the shear Alfven wave
    speed is included in the calculation of the fastest parallel wave
    speed for numerical dissipation.
-2. For tokamak simulations use zero-Laplacian boundary conditions
-   by setting ``electromagnetic:apar_boundary_neumann=false``.
+2. For tokamak simulations use Neumann boundary condition on the core
+   and Dirichlet on SOL and PF boundaries by setting
+   ``electromagnetic:apar_core_neumann=true`` (this is the default).
+3. Set the potential core boundary to be constant in Y by setting
+   ``vorticity:phi_core_averagey = true``
+4. Magnetic flutter terms must be enabled to be active
+   (``electromagnetic:magnetic_flutter=true``).  They use an
+   ``Apar_flutter`` field, not the ``Apar`` field that is calculated
+   from the induction terms.
 
 This component modifies the definition of momentum of all species, to
 include the contribution from the electromagnetic potential
@@ -2062,8 +2070,17 @@ Assumes that "momentum" :math:`p_s` calculated for all species
    p_s = m_s n_s v_{||s} + Z_s e n_s A_{||}
 
 which arises once the electromagnetic contribution to the force on
-each species is included in the momentum equation. This is normalised
-so that in dimensionless quantities
+each species is included in the momentum equation. This requires
+an additional term in the momentum equation:
+
+.. math::
+
+   \frac{\partial p_s}{\partial t} = \cdots + Z_s e A_{||} \frac{\partial n_s}{\partial t}
+
+This is implemented so that the density time-derivative is calculated using the lowest order
+terms (parallel flow, ExB drift and a low density numerical diffusion term).
+
+The above equations are normalised so that in dimensionless quantities:
 
 .. math::
 
@@ -2098,5 +2115,22 @@ To convert the species momenta into a current, we take the sum of
 
    - \frac{1}{\beta_{em}} \nabla_\perp^2 A_{||} + \sum_s \frac{Z^2 n_s}{A}A_{||} = \sum_s \frac{Z}{A} p_s
 
+The toroidal variation of density :math:`n_s` must be kept in this
+equation.  By default the iterative "Naulin" solver is used to do
+this: A fast FFT-based method is used in a fixed point iteration,
+correcting for the density variation.
+
+Magnetic flutter terms are disabled by default, and can be enabled by setting
+
+.. code-block:: ini
+
+   [electromagnetic]
+   magnetic_flutter = true
+
+This writes an ``Apar_flutter`` field to the state, which then enables perturbed
+parallel derivative terms in the ``evolve_density``, ``evolve_pressure``, ``evolve_energy`` and
+``evolve_momentum`` components. Parallel flow terms are modified, and parallel heat
+conduction.
+
 .. doxygenstruct:: Electromagnetic
    :members:

include/div_ops.hxx

@@ -23,8 +23,8 @@
 
  */
 
-#ifndef HERMES_DIV_OPS_H
-#define HERMES_DIV_OPS_H
+#ifndef DIV_OPS_H
+#define DIV_OPS_H
 
 #include <bout/field3d.hxx>
 #include <bout/fv_ops.hxx>
@@ -44,6 +44,11 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D& n, const Field3D& f,
                                     bool bndry_flux = true, bool poloidal = false,
                                     bool positive = false);
 
+/// This version has an extra coefficient 'g' that is linearly interpolated
+/// onto cell faces
+const Field3D Div_n_g_bxGrad_f_B_XZ(const Field3D &n, const Field3D &g, const Field3D &f, 
+                                    bool bndry_flux = true, bool positive = false);
+
 const Field3D Div_Perp_Lap_FV_Index(const Field3D& a, const Field3D& f, bool xflux);
 
 const Field3D Div_Z_FV_Index(const Field3D& a, const Field3D& f);
@@ -216,9 +221,10 @@ 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 * sv.R * 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
@@ -244,9 +250,10 @@ 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 * sv.L * 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
@@ -690,4 +697,4 @@ const Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, const Fi
 
 } // namespace FV
 
-#endif //  HERMES_DIV_OPS_H
+#endif //  DIV_OPS_H

include/electromagnetic.hxx

@@ -64,6 +64,13 @@ private:
 
   std::unique_ptr<Laplacian> aparSolver; // Laplacian solver in X-Z
 
+  bool const_gradient; // Set neumann boundaries by extrapolation
+  BoutReal apar_boundary_timescale; // Relaxation timescale
+  BoutReal last_time;  // The last time the boundaries were updated
+
+  bool magnetic_flutter; ///< Set the magnetic flutter term?
+  Field3D Apar_flutter;
+
   bool diagnose; ///< Output additional diagnostics?
 };
 

include/vorticity.hxx

@@ -37,6 +37,8 @@ struct Vorticity : public Component {
   ///     Relax radial phi boundaries towards zero-gradient?
   ///   - phi_boundary_timescale: float, 1e-4
   ///     Timescale for phi boundary relaxation [seconds]
+  ///   - phi_core_averagey: bool, default false
+  ///     Average phi core boundary in Y? (if phi_boundary_relax)
   ///   - phi_dissipation: bool, default true
   ///     Parallel dissipation of potential (Recommended)
   ///   - poloidal_flows: bool, default true
@@ -128,7 +130,8 @@ private:
   bool phi_boundary_relax; ///< Relax boundary to zero-gradient
   BoutReal phi_boundary_timescale; ///< Relaxation timescale [normalised]
   BoutReal phi_boundary_last_update; ///< Time when last updated
-
+  bool phi_core_averagey; ///< Average phi core boundary in Y?
+
   bool split_n0; // Split phi into n=0 and n!=0 components
   LaplaceXY* laplacexy; // Laplacian solver in X-Y (n=0)
 

src/div_ops.cxx

@@ -1372,3 +1372,160 @@ const Field3D Div_a_Grad_perp_flows(const Field3D& a, const Field3D& f,
   return result;
 }
 
+const Field3D Div_n_g_bxGrad_f_B_XZ(const Field3D &n, const Field3D &g, const Field3D &f, 
+                                    bool bndry_flux, bool positive) {
+  Field3D result{0.0};
+
+  Coordinates *coord = mesh->getCoordinates();
+
+  //////////////////////////////////////////
+  // X-Z advection.
+  //
+  //             Z
+  //             |
+  //
+  //    fmp --- vU --- fpp
+  //     |      nU      |
+  //     |               |
+  //    vL nL        nR vR    -> X
+  //     |               |
+  //     |      nD       |
+  //    fmm --- vD --- fpm
+  //
+
+  int nz = mesh->LocalNz;
+  for (int i = mesh->xstart; i <= mesh->xend; i++) {
+    for (int j = mesh->ystart; j <= mesh->yend; j++) {
+      for (int k = 0; k < nz; k++) {
+        int kp = (k + 1) % nz;
+        int kpp = (kp + 1) % nz;
+        int km = (k - 1 + nz) % nz;
+        int kmm = (km - 1 + nz) % nz;
+
+        // 1) Interpolate stream function f onto corners fmp, fpp, fpm
+
+        BoutReal fmm = 0.25 * (f(i, j, k) + f(i - 1, j, k) + f(i, j, km) +
+                               f(i - 1, j, km));
+        BoutReal fmp = 0.25 * (f(i, j, k) + f(i, j, kp) + f(i - 1, j, k) +
+                               f(i - 1, j, kp)); // 2nd order accurate
+        BoutReal fpp = 0.25 * (f(i, j, k) + f(i, j, kp) + f(i + 1, j, k) +
+                               f(i + 1, j, kp));
+        BoutReal fpm = 0.25 * (f(i, j, k) + f(i + 1, j, k) + f(i, j, km) +
+                               f(i + 1, j, km));
+
+        // 2) Calculate velocities on cell faces
+
+        BoutReal vU = coord->J(i, j) * (fmp - fpp) / coord->dx(i, j); // -J*df/dx
+        BoutReal vD = coord->J(i, j) * (fmm - fpm) / coord->dx(i, j); // -J*df/dx
+
+        BoutReal vR = 0.5 * (coord->J(i, j) + coord->J(i + 1, j)) * (fpp - fpm) /
+                      coord->dz(i, j); // J*df/dz
+        BoutReal vL = 0.5 * (coord->J(i, j) + coord->J(i - 1, j)) * (fmp - fmm) /
+                      coord->dz(i, j); // J*df/dz
+
+        // 3) Calculate g on cell faces
+
+        BoutReal gU = 0.5 * (g(i, j, kp) + g(i, j, k));
+        BoutReal gD = 0.5 * (g(i, j, km) + g(i, j, k));
+        BoutReal gR = 0.5 * (g(i + 1, j, k) + g(i, j, k));
+        BoutReal gL = 0.5 * (g(i - 1, j, k) + g(i, j, k));
+
+        // 4) Calculate n on the cell faces. The sign of the
+        //    velocity determines which side is used.
+
+        // X direction
+        Stencil1D s;
+        s.c = n(i, j, k);
+        s.m = n(i - 1, j, k);
+        s.mm = n(i - 2, j, k);
+        s.p = n(i + 1, j, k);
+        s.pp = n(i + 2, j, k);
+
+        MC(s, coord->dx(i, j));
+
+        // Right side
+        if ((i == mesh->xend) && (mesh->lastX())) {
+          // At right boundary in X
+
+          if (bndry_flux) {
+            BoutReal flux;
+            if (vR > 0.0) {
+              // Flux to boundary
+              flux = vR * s.R * gR;
+            } else {
+              // Flux in from boundary
+              flux = vR * 0.5 * (n(i + 1, j, k) + n(i, j, k)) * gR;
+            }
+
+            result(i, j, k) += flux / (coord->dx(i, j) * coord->J(i, j));
+            result(i + 1, j, k) -=
+                flux / (coord->dx(i + 1, j) * coord->J(i + 1, j));
+          }
+        } else {
+          // Not at a boundary
+          if (vR > 0.0) {
+            // Flux out into next cell
+            BoutReal flux = vR * s.R * gR;
+            result(i, j, k) += flux / (coord->dx(i, j) * coord->J(i, j));
+            result(i + 1, j, k) -=
+                flux / (coord->dx(i + 1, j) * coord->J(i + 1, j));
+          }
+        }
+
+        // Left side
+
+        if ((i == mesh->xstart) && (mesh->firstX())) {
+          // At left boundary in X
+
+          if (bndry_flux) {
+            BoutReal flux;
+
+            if (vL < 0.0) {
+              // Flux to boundary
+              flux = vL * s.L * gL;
+
+            } else {
+              // Flux in from boundary
+              flux = vL * 0.5 * (n(i - 1, j, k) + n(i, j, k)) * gL;
+            }
+            result(i, j, k) -= flux / (coord->dx(i, j) * coord->J(i, j));
+            result(i - 1, j, k) +=
+                flux / (coord->dx(i - 1, j) * coord->J(i - 1, j));
+          }
+        } else {
+          // Not at a boundary
+
+          if (vL < 0.0) {
+            BoutReal flux = vL * s.L * gL;
+            result(i, j, k) -= flux / (coord->dx(i, j) * coord->J(i, j));
+            result(i - 1, j, k) +=
+                flux / (coord->dx(i - 1, j) * coord->J(i - 1, j));
+          }
+        }
+
+        /// NOTE: Need to communicate fluxes
+
+        // Z direction
+        s.m = n(i, j, km);
+        s.mm = n(i, j, kmm);
+        s.p = n(i, j, kp);
+        s.pp = n(i, j, kpp);
+
+        MC(s, coord->dz(i, j));
+
+        if (vU > 0.0) {
+          BoutReal flux = vU * s.R * gU/ (coord->J(i, j) * coord->dz(i, j));
+          result(i, j, k) += flux;
+          result(i, j, kp) -= flux;
+        }
+        if (vD < 0.0) {
+          BoutReal flux = vD * s.L * gD / (coord->J(i, j) * coord->dz(i, j));
+          result(i, j, k) -= flux;
+          result(i, j, km) += flux;
+        }
+      }
+    }
+  }
+  FV::communicateFluxes(result);
+  return result;
+}