Merge branch 'turbulence-cases' into paper-inputs

Author: bendudson

include/diamagnetic_drift.hxx

@@ -22,6 +22,7 @@ struct DiamagneticDrift : public Component {
 private:
   Vector2D Curlb_B;
   bool bndry_flux;
+  Field2D diamag_form;
 };
 
 namespace {

include/evolve_density.hxx

@@ -27,6 +27,8 @@ struct EvolveDensity : public Component {
   /// - N<name>   e.g. "Ne", "Nd+"
   ///   - source    Source of particles [/m^3/s]
   ///               NOTE: This overrides mesh input N<name>_src
+  ///   - source_only_in_core         Zero the source outside the closed field-line region?
+  ///   - neumann_boundary_average_z  Apply Neumann boundaries with Z average?
   ///
   EvolveDensity(std::string name, Options &options, Solver *solver);
 
@@ -65,6 +67,7 @@ private:
 
   bool bndry_flux;      ///< Allow flows through boundaries?
   bool poloidal_flows;  ///< Include ExB flow in Y direction?
+  bool neumann_boundary_average_z; ///< Apply neumann boundary with Z average?
 
   BoutReal density_floor;
   bool low_n_diffuse;   ///< Parallel diffusion at low density

include/evolve_momentum.hxx

@@ -48,6 +48,7 @@ private:
   BoutReal hyper_z;  ///< Hyper-diffusion
 
   bool diagnose; ///< Output additional diagnostics?
+  bool fix_momentum_boundary_flux; ///< Fix momentum flux to boundary condition?
 };
 
 namespace {

include/evolve_pressure.hxx

@@ -33,6 +33,8 @@ struct EvolvePressure : public Component {
   /// - P<name>  e.g. "Pe", "Pd+"
   ///   - source     Source of pressure [Pa / s].
   ///                NOTE: This overrides mesh input P<name>_src
+  ///   - source_only_in_core         Zero the source outside the closed field-line region?
+  ///   - neumann_boundary_average_z  Apply Neumann boundaries with Z average?
   ///
   EvolvePressure(std::string name, Options& options, Solver* solver);
 
@@ -74,6 +76,7 @@ private:
   Field3D T, N; ///< Temperature, density
 
   bool bndry_flux;
+  bool neumann_boundary_average_z; ///< Apply neumann boundary with Z average?
   bool poloidal_flows;
   bool thermal_conduction;    ///< Include thermal conduction?
   BoutReal kappa_coefficient; ///< Leading numerical coefficient in parallel heat flux calculation

include/sheath_boundary.hxx

@@ -25,6 +25,7 @@ struct SheathBoundary : public Component {
   ///   - lower_y                  Boundary on lower y?
   ///   - upper_y                  Boundary on upper y?
   ///   - wall_potential           Voltage of the wall [Volts]
+  ///   - floor_potential          Apply floor to sheath potential?
   ///   - secondary_electron_coef  Effective secondary electron emission coefficient
   ///   - sin_alpha                Sine of the angle between magnetic field line and wall surface (0 to 1)
   ///   - always_set_phi           Always set phi field? Default is to only modify if already set
@@ -84,6 +85,8 @@ private:
   bool always_set_phi; ///< Set phi field?
 
   Field3D wall_potential; ///< Voltage at the wall. Normalised units.
+
+  bool floor_potential; ///< Apply floor to sheath potential?
 };
 
 namespace {

include/sound_speed.hxx

@@ -4,6 +4,8 @@
 
 #include "component.hxx"
 
+#include <bout/constants.hxx>
+
 /// Calculate the system sound speed
 ///
 /// This uses the sum of all species pressures and mass densities
@@ -14,6 +16,25 @@ struct SoundSpeed : public Component {
     electron_dynamics = options["electron_dynamics"]
       .doc("Include electron sound speed?")
       .withDefault<bool>(true);
+
+    alfven_wave = options["alfven_wave"]
+      .doc("Include Alfven wave speed?")
+      .withDefault<bool>(false);
+    if (alfven_wave) {
+      // Calculate normalisation factor
+      const auto& units = alloptions["units"];
+      const auto Bnorm = get<BoutReal>(units["Tesla"]);
+      const auto Nnorm = get<BoutReal>(units["inv_meters_cubed"]);
+      const auto Cs0 = get<BoutReal>(units["meters"]) / get<BoutReal>(units["seconds"]);
+      beta_norm = Bnorm / sqrt(SI::mu0 * Nnorm * SI::Mp) / Cs0;
+    }
+
+    temperature_floor = options["temperature_floor"]
+      .doc("Minimum temperature when calculating sound speeds [eV]")
+      .withDefault(0.0);
+    if (temperature_floor > 0.0) {
+      temperature_floor /= get<BoutReal>(alloptions["units"]["eV"]);
+    }
   }
 
   /// This sets in the state
@@ -31,6 +52,9 @@ struct SoundSpeed : public Component {
 
 private:
   bool electron_dynamics; ///< Include electron sound speed?
+  bool alfven_wave; ///< Include Alfven wave speed?
+  BoutReal beta_norm{0.0}; ///< Normalisation factor for Alfven speed
+  BoutReal temperature_floor; ///< Minimum temperature when calculating speed
 };
 
 namespace {

include/vorticity.hxx

@@ -91,6 +91,7 @@ private:
 
   bool vort_dissipation; ///< Parallel dissipation of vorticity
   bool phi_dissipation;  ///< Parallel dissipation of potential
+  bool phi_sheath_dissipation; ///< Dissipation at the sheath if phi < 0
 
   bool phi_boundary_relax; ///< Relax boundary to zero-gradient
   BoutReal phi_boundary_timescale; ///< Relaxation timescale [normalised]
@@ -102,6 +103,7 @@ private:
   Field2D Bsq; // SQ(coord->Bxy)
   Vector2D Curlb_B; // Curvature vector Curl(b/B)
   BoutReal hyper_z; ///< Hyper-viscosity in Z
+  Field2D viscosity; ///< Kinematic viscosity
 
   // Diagnostic outputs
   Field3D DivJdia, DivJcol; // Divergence of diamagnetic and collisional current

include/zero_current.hxx

@@ -37,6 +37,12 @@ struct ZeroCurrent : public Component {
   /// 
   void transform(Options &state) override;
 
+  void finally(const Options &state) override {
+    // Get the velocity with boundary condition applied.
+    // This is for output only
+    velocity = get<Field3D>(state["species"][name]["velocity"]);
+  }
+
   void outputVars(Options &state) override;
 private:
   std::string name; ///< Name of this species

src/diamagnetic_drift.cxx

@@ -1,4 +1,5 @@
 #include <bout/fv_ops.hxx>
+#include <vecops.hxx>
 
 #include "../include/diamagnetic_drift.hxx"
 
@@ -13,6 +14,10 @@ DiamagneticDrift::DiamagneticDrift(std::string name, Options& alloptions,
   bndry_flux =
       options["bndry_flux"].doc("Allow fluxes through boundary?").withDefault<bool>(true);
 
+  diamag_form = options["diamag_form"]
+    .doc("Form of diamagnetic drift: 0 = gradient; 1 = divergence")
+    .withDefault(Field2D(1.0));
+
   // Read curvature vector
   Curlb_B.covariant = false; // Contravariant
   if (mesh->get(Curlb_B, "bxcv")) {
@@ -41,6 +46,15 @@ DiamagneticDrift::DiamagneticDrift(std::string name, Options& alloptions,
   Curlb_B.z *= SQ(Lnorm);
 
   Curlb_B *= 2. / mesh->getCoordinates()->Bxy;
+
+  // Set drift to zero through sheath boundaries.
+  // Flux through those cell faces should be set by sheath.
+  for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) {
+    Curlb_B.y(r.ind, mesh->ystart - 1) = -Curlb_B.y(r.ind, mesh->ystart);
+  }
+  for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) {
+    Curlb_B.y(r.ind, mesh->yend + 1) = -Curlb_B.y(r.ind, mesh->yend);
+  }
 }
 
 void DiamagneticDrift::transform(Options& state) {
@@ -62,17 +76,28 @@ void DiamagneticDrift::transform(Options& state) {
 
     if (IS_SET(species["density"])) {
       auto N = GET_VALUE(Field3D, species["density"]);
-      subtract(species["density_source"], FV::Div_f_v(N, vD, bndry_flux));
+
+      // Divergence form: Div(n v_D)
+      Field3D div_form = FV::Div_f_v(N, vD, bndry_flux);
+      // Gradient form: Curlb_B dot Grad(N T / q)
+      Field3D grad_form = Curlb_B * Grad(N * T / q);
+
+      subtract(species["density_source"], diamag_form * div_form + (1. - diamag_form) * grad_form);
     }
 
     if (IS_SET(species["pressure"])) {
       auto P = get<Field3D>(species["pressure"]);
-      subtract(species["energy_source"], (5. / 2) * FV::Div_f_v(P, vD, bndry_flux));
+
+      Field3D div_form = FV::Div_f_v(P, vD, bndry_flux);
+      Field3D grad_form = Curlb_B * Grad(P * T / q);
+      subtract(species["energy_source"], (5. / 2) * (diamag_form * div_form + (1. - diamag_form) * grad_form));
     }
 
     if (IS_SET(species["momentum"])) {
       auto NV = get<Field3D>(species["momentum"]);
-      subtract(species["momentum_source"], FV::Div_f_v(NV, vD, bndry_flux));
+      Field3D div_form = FV::Div_f_v(NV, vD, bndry_flux);
+      Field3D grad_form = Curlb_B * Grad(NV * T / q);
+      subtract(species["momentum_source"], diamag_form * div_form + (1. - diamag_form) * grad_form);
     }
   }
 }

src/evolve_density.cxx

@@ -81,6 +81,25 @@ EvolveDensity::EvolveDensity(std::string name, Options& alloptions, Solver* solv
                .doc("Source term in ddt(N" + name + std::string("). Units [m^-3/s]"))
                .withDefault(source)
            / (Nnorm * Omega_ci);
+
+  if (alloptions[std::string("N") + name]["source_only_in_core"]
+      .doc("Zero the source outside the closed field-line region?")
+      .withDefault<bool>(false)) {
+    for (int x = mesh->xstart; x <= mesh->xend; x++) {
+      if (!mesh->periodicY(x)) {
+        // Not periodic, so not in core
+        for (int y = mesh->ystart; y <= mesh->yend; y++) {
+          for (int z = mesh->zstart; z <= mesh->zend; z++) {
+            source(x, y, z) = 0.0;
+          }
+        }
+      }
+    }
+  }
+
+  neumann_boundary_average_z = alloptions[std::string("N") + name]["neumann_boundary_average_z"]
+    .doc("Apply neumann boundary with Z average?")
+    .withDefault<bool>(false);
 }
 
 void EvolveDensity::transform(Options& state) {
@@ -93,6 +112,40 @@ void EvolveDensity::transform(Options& state) {
 
   mesh->communicate(N);
 
+  if (neumann_boundary_average_z) {
+    // Take Z (usually toroidal) average and apply as X (radial) boundary condition
+    if (mesh->firstX()) {
+      for (int j = mesh->ystart; j <= mesh->yend; j++) {
+        BoutReal Navg = 0.0; // Average N in Z
+        for (int k = 0; k < mesh->LocalNz; k++) {
+          Navg += N(mesh->xstart, j, k);
+        }
+        Navg /= mesh->LocalNz;
+
+        // Apply boundary condition
+        for (int k = 0; k < mesh->LocalNz; k++) {
+          N(mesh->xstart - 1, j, k) = 2. * Navg - N(mesh->xstart, j, k);
+          N(mesh->xstart - 2, j, k) = N(mesh->xstart - 1, j, k);
+        }
+      }
+    }
+
+    if (mesh->lastX()) {
+      for (int j = mesh->ystart; j <= mesh->yend; j++) {
+        BoutReal Navg = 0.0; // Average N in Z
+        for (int k = 0; k < mesh->LocalNz; k++) {
+          Navg += N(mesh->xend, j, k);
+        }
+        Navg /= mesh->LocalNz;
+
+        for (int k = 0; k < mesh->LocalNz; k++) {
+          N(mesh->xend + 1, j, k) = 2. * Navg - N(mesh->xend, j, k);
+          N(mesh->xend + 2, j, k) = N(mesh->xend + 1, j, k);
+        }
+      }
+    }
+  }
+
   auto& species = state["species"][name];
   set(species["density"], floor(N, 0.0)); // Density in state always >= 0
   set(species["AA"], AA);                 // Atomic mass

src/evolve_momentum.cxx

@@ -43,6 +43,10 @@ EvolveMomentum::EvolveMomentum(std::string name, Options &alloptions, Solver *so
   diagnose = options["diagnose"]
     .doc("Output additional diagnostics?")
     .withDefault<bool>(false);
+
+  fix_momentum_boundary_flux = options["fix_momentum_boundary_flux"]
+    .doc("Fix Y boundary momentum flux to boundary midpoint value?")
+    .withDefault<bool>(false);
 }
 
 void EvolveMomentum::transform(Options &state) {
@@ -115,7 +119,7 @@ void EvolveMomentum::finally(const Options &state) {
 
   // Note: Density floor should be consistent with calculation of V
   //       otherwise energy conservation is affected
-  ddt(NV) -= AA * FV::Div_par_fvv(Nlim, V, fastest_wave);
+  ddt(NV) -= AA * FV::Div_par_fvv(Nlim, V, fastest_wave, fix_momentum_boundary_flux);
 
   // Parallel pressure gradient
   if (species.isSet("pressure")) {

src/evolve_pressure.cxx

@@ -90,6 +90,25 @@ EvolvePressure::EvolvePressure(std::string name, Options& alloptions, Solver* so
                     + std::string("). Units [N/m^2/s]"))
                .withDefault(source)
            / (SI::qe * Nnorm * Tnorm * Omega_ci);
+
+  if (alloptions[std::string("P") + name]["source_only_in_core"]
+      .doc("Zero the source outside the closed field-line region?")
+      .withDefault<bool>(false)) {
+    for (int x = mesh->xstart; x <= mesh->xend; x++) {
+      if (!mesh->periodicY(x)) {
+        // Not periodic, so not in core
+        for (int y = mesh->ystart; y <= mesh->yend; y++) {
+          for (int z = mesh->zstart; z <= mesh->zend; z++) {
+            source(x, y, z) = 0.0;
+          }
+        }
+      }
+    }
+  }
+
+  neumann_boundary_average_z = alloptions[std::string("P") + name]["neumann_boundary_average_z"]
+    .doc("Apply neumann boundary with Z average?")
+    .withDefault<bool>(false);
 }
 
 void EvolvePressure::transform(Options& state) {
@@ -102,6 +121,39 @@ void EvolvePressure::transform(Options& state) {
 
   mesh->communicate(P);
 
+  if (neumann_boundary_average_z) {
+    // Take Z (usually toroidal) average and apply as X (radial) boundary condition
+    if (mesh->firstX()) {
+      for (int j = mesh->ystart; j <= mesh->yend; j++) {
+        BoutReal Pavg = 0.0; // Average P in Z
+        for (int k = 0; k < mesh->LocalNz; k++) {
+          Pavg += P(mesh->xstart, j, k);
+        }
+        Pavg /= mesh->LocalNz;
+
+        // Apply boundary condition
+        for (int k = 0; k < mesh->LocalNz; k++) {
+          P(mesh->xstart - 1, j, k) = 2. * Pavg - P(mesh->xstart, j, k);
+          P(mesh->xstart - 2, j, k) = P(mesh->xstart - 1, j, k);
+        }
+      }
+    }
+
+    if (mesh->lastX()) {
+      for (int j = mesh->ystart; j <= mesh->yend; j++) {
+        BoutReal Pavg = 0.0; // Average P in Z
+        for (int k = 0; k < mesh->LocalNz; k++) {
+          Pavg += P(mesh->xend, j, k);
+        }
+        Pavg /= mesh->LocalNz;
+
+        for (int k = 0; k < mesh->LocalNz; k++) {
+          P(mesh->xend + 1, j, k) = 2. * Pavg - P(mesh->xend, j, k);
+          P(mesh->xend + 2, j, k) = P(mesh->xend + 1, j, k);
+        }
+      }
+    }
+  }
 
   auto& species = state["species"][name];
 

src/sheath_boundary.cxx

@@ -100,6 +100,10 @@ SheathBoundary::SheathBoundary(std::string name, Options &alloptions, Solver *)
   // not the value at the boundary half-way between cells. This is due
   // to how twist-shift boundary conditions and non-aligned inputs are
   // treated; using the cell boundary gives incorrect results.
+
+  floor_potential = options["floor_potential"]
+                        .doc("Apply a floor to wall potential when calculating Ve?")
+                        .withDefault<bool>(true);
 }
 
 void SheathBoundary::transform(Options &state) {
@@ -323,11 +327,12 @@ void SheathBoundary::transform(Options &state) {
         const BoutReal tesheath = 0.5 * (Te[im] + Te[i]);  // electron temperature
         const BoutReal phi_wall = wall_potential[i];
 
-        const BoutReal phisheath = floor(
-            0.5 * (phi[im] + phi[i]), phi_wall); // Electron saturation at phi = phi_wall
+        const BoutReal phisheath = floor_potential ? floor(
+            0.5 * (phi[im] + phi[i]), phi_wall) // Electron saturation at phi = phi_wall
+	    : 0.5 * (phi[im] + phi[i]);
 
         // Electron sheath heat transmission
-        const BoutReal gamma_e = 2 / (1. - Ge) + (phisheath - phi_wall) / floor(tesheath, 1e-5);
+        const BoutReal gamma_e = floor(2 / (1. - Ge) + (phisheath - phi_wall) / floor(tesheath, 1e-5), 0.0);
 
         // Electron velocity into sheath (< 0)
         const BoutReal vesheath = (tesheath < 1e-10) ?
@@ -385,10 +390,11 @@ void SheathBoundary::transform(Options &state) {
         const BoutReal nesheath = 0.5 * (Ne[ip] + Ne[i]);
         const BoutReal tesheath = 0.5 * (Te[ip] + Te[i]);  // electron temperature
         const BoutReal phi_wall = wall_potential[i];
-        const BoutReal phisheath = floor(0.5 * (phi[ip] + phi[i]), phi_wall); // Electron saturation at phi = phi_wall
+        const BoutReal phisheath = floor_potential ? floor(0.5 * (phi[ip] + phi[i]), phi_wall) // Electron saturation at phi = phi_wall
+                                                   : 0.5 * (phi[ip] + phi[i]);
 
         // Electron sheath heat transmission
-        const BoutReal gamma_e = 2 / (1. - Ge) + (phisheath - phi_wall) / floor(tesheath, 1e-5);
+        const BoutReal gamma_e = floor(2 / (1. - Ge) + (phisheath - phi_wall) / floor(tesheath, 1e-5), 0.0);
 
         // Electron velocity into sheath (> 0)
         const BoutReal vesheath = (tesheath < 1e-10) ?