Neutrals flux limiters

include/braginskii_collisions.hxx

@@ -57,6 +57,7 @@ private:
       neutral_neutral;
 
   BoutReal ei_multiplier; // Arbitrary user-set multiplier on electron-ion collisions
+  BoutReal density_floor;
 
   /// Calculated collision rates saved for post-processing and use by other components
   /// Saved in options, the BOUT++ dictionary-like object

include/evolve_density.hxx

@@ -86,6 +86,9 @@ private:
   bool diagnose; ///< Output additional diagnostics?
   Field3D flow_xlow, flow_ylow; ///< Particle flow diagnostics
 
+  bool initialize_from_mesh; ///< Initilize the Field3D N from 2D profiles stored in the
+                             ///< mesh file.
+
   /// This sets in the state
   /// - species
   ///   - <name>

include/evolve_pressure.hxx

@@ -100,6 +100,8 @@ private:
   bool fix_momentum_boundary_flux; ///< Fix momentum flux to boundary condition?
   Field3D Sp_nvh; ///< Pressure source due to artificial viscosity
   Field3D E_PdivV, E_VgradP; ///< Diagnostic energy source terms for p*Div(V) and V*Grad(P)
+  bool initialize_from_mesh; ///< Initilize the Field3D P from 2D profiles stored in the
+                             ///< mesh file.
 
   /// Inputs
   /// - species

include/fixed_density.hxx

@@ -4,6 +4,8 @@
 
 #include "component.hxx"
 
+using bout::globals::mesh;
+
 /// Set ion density to a fixed value
 ///
 struct FixedDensity : public Component {
@@ -26,6 +28,19 @@ struct FixedDensity : public Component {
 
     // Get the density and normalise
     N = options["density"].as<Field3D>() / Nnorm;
+
+    // Initilize the Field3D N from 2D profiles stored in the mesh file.
+    initialize_from_mesh =
+        options["initialize_from_mesh"]
+            .doc("Initilize field from 2D profiles stored in the mesh file?")
+            .withDefault<bool>(false);
+
+    if (initialize_from_mesh) {
+      // Try to read the initial field from the mesh
+      mesh->get(N, std::string("N") + name + "_init"); // Units: [m^-3/s]
+      N /= Nnorm;                                      // Normalization
+    }
+
     substitutePermissions("name", {name});
     substitutePermissions("vars", {"AA", "charge", "density"});
   }
@@ -51,6 +66,9 @@ private:
 
   Field3D N; ///< Species density (normalised)
 
+  bool initialize_from_mesh; ///< Initilize the Field3D N from 2D profiles stored in the
+                             ///< mesh file.
+
   /// Sets in the state the density, mass and charge of the species
   ///
   /// - species

include/isothermal.hxx

@@ -16,6 +16,9 @@ private:
   BoutReal T; ///< The normalised temperature
   Field3D P; ///< The normalised pressure
 
+  bool initialize_from_mesh; ///< Initilize the Field3D T from 2D profiles stored in the
+                             ///< mesh file.
+
   bool diagnose; ///< Output additional diagnostics?
 
   /// Inputs

include/neutral_mixed.hxx

@@ -37,17 +37,30 @@ private:
   Field3D Tn;                     ///< Neutral temperature
   Field3D Nnlim, Pnlim, logPnlim; // Limited in regions of low density
 
+  Field3D NVn_err;    ///< Difference from momentum as input from solver
+  Field3D NVn_solver; ///< Momentum as calculated in the solver
+
   BoutReal AA; ///< Atomic mass (proton = 1)
 
   std::vector<std::string> collision_names; ///< Collisions used for collisionality
   std::string
       diffusion_collisions_mode; ///< Collision selection, either afn or multispecies
   Field3D nu;                    ///< Collisionality to use for diffusion
-  Field3D Dnn;                   ///< Diffusion coefficient
+  Field3D nu_pseudo_mfp;         ///< Pseudo-collision frequency based on mean free path
+  Field3D nu_total; ///< Total collision frequency used for diffusion, including
+                    ///< pseudo-collisions
+  Field3D Dnn, Dnn_unlimited, Dmax;    ///< Diffusion coefficient
   Field3D DnnNn, DnnPn, DnnTn, DnnNVn; ///< Used for operators
-  BoutReal flux_limit;                 ///< Diffusive flux limit
-  BoutReal diffusion_limit;            ///< Maximum diffusion coefficient
+  Field3D kappa_n_unlimited, kappa_n_max_par, kappa_n_max_perp;
+  BoutReal flux_limit_adv;       ///< Diffusive flux limit
+  BoutReal flux_limit_cond_par;  ///< Limit for parallel conductive flux
+  BoutReal flux_limit_cond_perp; ///< Limit for perpendicular conductive flux
+  BoutReal flux_limit_visc_par;  ///< Limit for parallel viscous flux
+  BoutReal flux_limit_visc_perp; ///< Limit for perpendicular viscous flux
+
+  BoutReal diffusion_limit; ///< Maximum diffusion coefficient
   BoutReal neutral_lmax;
+  BoutReal flux_limiter_sharpness; ///< Sharpness of flux limiter transition
 
   bool sheath_ydown, sheath_yup;
 
@@ -64,8 +77,18 @@ private:
   bool neutral_conduction; ///< Include heat conduction?
   bool evolve_momentum;    ///< Evolve parallel momentum?
   bool normalise_sources;  ///< Normalise input sources?
+  bool perp_ion_coupling;  ///< Include coupling to ion perpendicular velocity?
+
+  // Temporary variables
+  Field3D debug;          ///< Debug variable FIXME: remove
+  bool double_count_lmax; ///< Include neutral_lmax in Dmax and kappa_max as well as Dnn?
+  bool legacy_thermal_speed; ///< Use legacy definition of thermal speed in flux limiter?
+  bool legacy_limiter_form;  ///< Use legacy form of flux limiter rather than SOLPS-style
 
-  Field3D kappa_n, eta_n; ///< Neutral conduction and viscosity
+  Field3D kappa_n, eta_n_unlimited;      ///< Neutral conduction and viscosity
+  Field3D kappa_n_perp, eta_n_perp;      ///< Neutral conduction and viscosity
+  Field3D kappa_n_par, eta_n_par;        ///< Neutral conduction and viscosity
+  Field3D eta_n_max_par, eta_n_max_perp; ///< Viscosity reduction factor for flux-limiting
 
   bool nonorthogonal_operators;   ///< Use nonorthogonal operators for radial transport?
   bool precondition{true};        ///< Enable preconditioner?
@@ -77,6 +100,8 @@ private:
   Field3D Sn, Sp, Snv; ///< Particle, pressure and momentum source
   Field3D sound_speed; ///< Sound speed for use with Lax flux
 
+  bool zero_timederivs; ///< Set the time derivatives to zero?
+
   bool output_ddt; ///< Save time derivatives?
   bool diagnose;   ///< Save additional diagnostics?
 

include/zero_current.hxx

@@ -30,9 +30,11 @@ struct ZeroCurrent : public Component {
 private:
   std::string name; ///< Name of this species
   BoutReal charge;  ///< The charge of this species
+  BoutReal atomic_mass; // Atomic mass
 
   Field3D velocity; ///< Species velocity (for writing to output)
-
+  Field3D momentum; ///< Species momentum (for writing to output)
+
   /// Required inputs
   /// - species
   ///   - <name>

src/anomalous_diffusion.cxx

@@ -137,6 +137,8 @@ void AnomalousDiffusion::transform_impl(GuardedOptions& state) {
                                      flow_xlow, flow_ylow));
     add(species["energy_flow_xlow"], flow_xlow);
     add(species["energy_flow_ylow"], flow_ylow);
+
+    set(species["anomalous_D"], anomalous_D);
   }
 
   if (include_chi) {

src/braginskii_collisions.cxx

@@ -59,6 +59,8 @@ BraginskiiCollisions::BraginskiiCollisions(const std::string& name, Options& all
                       .doc("User-set arbitrary multiplier on electron-ion collision rate")
                       .withDefault<BoutReal>(1.0);
 
+  density_floor = options["density_floor"].doc("Minimum density floor").withDefault(1e-8);
+
   diagnose =
       options["diagnose"].doc("Output additional diagnostics?").withDefault<bool>(false);
 
@@ -144,7 +146,7 @@ void BraginskiiCollisions::collide(GuardedOptions& species1, GuardedOptions& spe
     const Field3D density2 = GET_NOBOUNDARY(Field3D, species2["density"]);
 
     const Field3D nu = filledFrom(nu_12, [&](auto& i) {
-      return nu_12[i] * (A1 / A2) * density1[i] / softFloor(density2[i], 1e-5);
+      return nu_12[i] * (A1 / A2) * density1[i] / softFloor(density2[i], density_floor);
     });
 
     add(species2["collision_frequency"], nu); // Total collision frequency

src/evolve_density.cxx

@@ -143,6 +143,18 @@ EvolveDensity::EvolveDensity(std::string name, Options& alloptions, Solver* solv
   }
   substitutePermissions("name", {name});
   substitutePermissions("outputs", outputs);
+
+  // Initilize the Field3D N from 2D profiles stored in the mesh file.
+  initialize_from_mesh =
+      n_options["initialize_from_mesh"]
+          .doc("Initilize field from 2D profiles stored in the mesh file?")
+          .withDefault<bool>(false);
+
+  if (initialize_from_mesh) {
+    // Try to read the initial field from the mesh
+    mesh->get(N, std::string("N") + name + "_init"); // Units: [m^-3/s]
+    N /= Nnorm;                                      // Normalization
+  }
 }
 
 void EvolveDensity::transform_impl(GuardedOptions& state) {

src/evolve_pressure.cxx

@@ -100,10 +100,10 @@ EvolvePressure::EvolvePressure(std::string name, Options& alloptions, Solver* so
   const BoutReal Nnorm = units["inv_meters_cubed"];
   const BoutReal Tnorm = units["eV"];
   const BoutReal Omega_ci = 1. / units["seconds"].as<BoutReal>();
+  const BoutReal Pnorm = SI::qe * Tnorm * Nnorm; // Pressure normalisation
 
   auto& p_options = alloptions[std::string("P") + name];
-  source_normalisation =
-      SI::qe * Nnorm * Tnorm * Omega_ci; // [Pa/s] or [W/m^3] if converted to energy
+  source_normalisation = Pnorm * Omega_ci; // [Pa/s] or [W/m^3] if converted to energy
   time_normalisation = 1. / Omega_ci;    // [s]
 
   // Try to read the pressure source from the mesh
@@ -164,6 +164,18 @@ EvolvePressure::EvolvePressure(std::string name, Options& alloptions, Solver* so
                            .doc("Include parallel heat conduction?")
                            .withDefault<bool>(true);
 
+  // Initilize the Field3D P from 2D profiles stored in the mesh file.
+  initialize_from_mesh =
+      p_options["initialize_from_mesh"]
+          .doc("Initilize field from 2D profiles stored in the mesh file?")
+          .withDefault<bool>(false);
+
+  if (initialize_from_mesh) {
+    // Try to read the initial field from the mesh
+    mesh->get(P, std::string("P") + name + "_init"); // Units: Pascal [Pa]
+    P /= Pnorm;                                      // Normalization
+  }
+
   if (source_time_dependent) {
     setPermissions(readOnly("time"));
   }

src/isothermal.cxx

@@ -1,7 +1,11 @@
 
+#include "../include/isothermal.hxx"
 #include <bout/constants.hxx>
+#include <bout/mesh.hxx>
 
-#include "../include/isothermal.hxx"
+using bout::globals::mesh;
+
+using bout::globals::mesh;
 
 Isothermal::Isothermal(std::string name, Options& alloptions, Solver* UNUSED(solver))
     : Component({readIfSet(fmt::format("species:{}:density", name), Regions::Interior),
@@ -15,6 +19,18 @@ Isothermal::Isothermal(std::string name, Options& alloptions, Solver* UNUSED(sol
   T = options["temperature"].doc("Constant temperature [eV]").as<BoutReal>()
       / Tnorm; // Normalise
 
+  // Initilize the Field3D T from 2D profiles stored in the mesh file.
+  initialize_from_mesh =
+      options["initialize_from_mesh"]
+          .doc("Initilize field from 2D profiles stored in the mesh file?")
+          .withDefault<bool>(false);
+
+  if (initialize_from_mesh) {
+    // Try to read the initial field from the mesh
+    mesh->get(T, std::string("T") + name + "_init"); // Units: [eV]
+    T /= Tnorm;                                      // Normalization
+  }
+
   diagnose = options["diagnose"]
     .doc("Save additional output diagnostics")
     .withDefault<bool>(false);