Added floors to several components.

.git-blame-ignore-revs

@@ -2,3 +2,4 @@ f5de620399e88f139e236e40675facfd0cfd6b13
 07616298519738fad0c1bb259e4d021cd6f51d58
 37efa12abe0961add8715b797d3926664977d453
 8e49d5933131ddcb2f97bb2ac7095a1be80ad7b7
+eaf52851d149ca10f69371e3c78e413f2e3e60cf

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/braginskii_electron_viscosity.hxx

@@ -37,6 +37,9 @@ private:
   bool diagnose;            ///< Output viscosity diagnostic?
   Field3D viscosity;        ///< The viscosity momentum source
 
+  BoutReal density_floor;
+  BoutReal pressure_floor; ///< When non-zero pressure is needed
+
   /// Inputs
   /// - species
   ///   - e

include/braginskii_ion_viscosity.hxx

@@ -61,6 +61,7 @@ private:
                                      ///< frequency change
   BoutReal bounce_frequency_R;       ///< Input major radius
   BoutReal density_floor;            ///< Minimum density used in calculating Pi_ciperp
+  BoutReal pressure_floor;           ///< When non-zero pressure is needed
   bool diagnose;                     ///< Output additional diagnostics?
 
   /// Per-species diagnostics

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,8 @@ 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) * softFloor(density1[i], density_floor)
+             / softFloor(density2[i], density_floor);
     });
 
     add(species2["collision_frequency"], nu); // Total collision frequency

src/braginskii_electron_viscosity.cxx

@@ -29,6 +29,14 @@ BraginskiiElectronViscosity::BraginskiiElectronViscosity(const std::string& name
                         .doc("Viscosity flux limiter coefficient. <0 = turned off")
                         .withDefault(-1.0);
 
+  density_floor = options["density_floor"].doc("Minimum density floor").withDefault(1e-8);
+
+  BoutReal temperature_floor =
+      options["temperature_floor"].doc("Low temperature scale").withDefault<BoutReal>(0.1)
+      / get<BoutReal>(alloptions["units"]["eV"]);
+
+  pressure_floor = density_floor * temperature_floor;
+
   diagnose = options["diagnose"].doc("Output diagnostics?").withDefault<bool>(false);
 }
 
@@ -44,8 +52,8 @@ void BraginskiiElectronViscosity::transform_impl(GuardedOptions& state) {
     throw BoutException("No electron velocity => Can't calculate electron viscosity");
   }
 
-  const Field3D tau = 1. / get<Field3D>(species["collision_frequency"]);
-  const Field3D P = get<Field3D>(species["pressure"]);
+  const Field3D tau = 1. / softFloor(get<Field3D>(species["collision_frequency"]), 1e-15);
+  const Field3D P = floor(get<Field3D>(species["pressure"]), pressure_floor);
   const Field3D V = get<Field3D>(species["velocity"]);
 
   Coordinates* coord = P.getCoordinates();
@@ -59,10 +67,10 @@ void BraginskiiElectronViscosity::transform_impl(GuardedOptions& state) {
     // SOLPS-style flux limiter
     // Values of alpha ~ 0.5 typically
 
-    const Field3D q_cl = eta * Grad_par(V);   // Collisional value
+    const Field3D q_cl = eta * abs(Grad_par(V)); // Collisional value
     const Field3D q_fl = eta_limit_alpha * P; // Flux limit
 
-    eta = eta / (1. + abs(q_cl / q_fl));
+    eta = eta / (1. + softFloor(q_cl, 1e-15) / softFloor(q_fl, 1e-15));
 
     eta.getMesh()->communicate(eta);
     eta.applyBoundary("neumann");

src/braginskii_friction.cxx

@@ -105,7 +105,6 @@ void BraginskiiFriction::transform_impl(GuardedOptions& state) {
       }
 
       const Field3D nu = GET_VALUE(Field3D, species1["collision_frequencies"][coll_name]);
-      const Field3D density2 = GET_VALUE(Field3D, species2["density"]);
       const Field3D velocity2 = species2.isSet("velocity")
                                     ? GET_NOBOUNDARY(Field3D, species2["velocity"])
                                     : 0.0;

src/braginskii_ion_viscosity.cxx

@@ -48,6 +48,14 @@ BraginskiiIonViscosity::BraginskiiIonViscosity(const std::string& name,
                         .doc("Viscosity flux limiter coefficient. <0 = turned off")
                         .withDefault(-1.0);
 
+  density_floor = options["density_floor"].doc("Minimum density floor").withDefault(1e-8);
+
+  BoutReal temperature_floor =
+      options["temperature_floor"].doc("Low temperature scale").withDefault<BoutReal>(0.1)
+      / get<BoutReal>(alloptions["units"]["eV"]);
+
+  pressure_floor = density_floor * temperature_floor;
+
   diagnose =
       options["diagnose"].doc("Output additional diagnostics?").withDefault<bool>(false);
 
@@ -84,8 +92,6 @@ BraginskiiIonViscosity::BraginskiiIonViscosity(const std::string& name,
                                 "frequency modification to viscosity")
                            .withDefault(2.0);
 
-  density_floor = options["density_floor"].doc("Minimum density floor").withDefault(1e-7);
-
   if (perpendicular) {
     // Read curvature vector
     try {
@@ -230,7 +236,7 @@ void BraginskiiIonViscosity::transform_impl(GuardedOptions& state) {
     }
 
     const Field3D tau = 1. / nu;
-    const Field3D P = get<Field3D>(species["pressure"]);
+    const Field3D P = floor(get<Field3D>(species["pressure"]), pressure_floor);
     const Field3D T = get<Field3D>(species["temperature"]);
     const Field3D N = get<Field3D>(species["density"]);
 
@@ -272,10 +278,10 @@ void BraginskiiIonViscosity::transform_impl(GuardedOptions& state) {
         // SOLPS-style flux limiter
         // Values of alpha ~ 0.5 typically
 
-        const Field3D q_cl = eta * Grad_par(V);   // Collisional value
+        const Field3D q_cl = eta * abs(Grad_par(V)); // Collisional value
         const Field3D q_fl = eta_limit_alpha * P; // Flux limit
 
-        eta = eta / (1. + abs(q_cl / q_fl));
+        eta = eta / (1. + (q_cl / q_fl));
 
         eta.getMesh()->communicate(eta);
         eta.applyBoundary("neumann");

src/braginskii_thermal_force.cxx

@@ -38,7 +38,7 @@ void BraginskiiThermalForce::transform_impl(GuardedOptions& state) {
 
       const BoutReal Z = get<BoutReal>(species["charge"]);
       // Don't need density boundary
-      const Field3D nz = GET_NOBOUNDARY(Field3D, species["density"]);
+      const Field3D nz = floor(GET_NOBOUNDARY(Field3D, species["density"]), 0.0);
 
       const Field3D ion_force = nz * (0.71 * SQ(Z)) * Grad_Te;
 
@@ -123,7 +123,8 @@ void BraginskiiThermalForce::transform_impl(GuardedOptions& state) {
 
         const BoutReal mz = get<BoutReal>((*heavy)["AA"]);
         const BoutReal Z = get<BoutReal>((*heavy)["charge"]);
-        const Field3D nz = GET_NOBOUNDARY(Field3D, (*heavy)["density"]);
+        const Field3D nz = floor(GET_NOBOUNDARY(Field3D, (*heavy)["density"]), 0.0);
+        ;
 
         if (Z == 0.0) {
           continue; // Check that the charge is not zero

src/hydrogen_charge_exchange.cxx

@@ -20,7 +20,7 @@ void HydrogenChargeExchange::calculate_rates(GuardedOptions&& atom1, GuardedOpti
 
   // Calculate effective temperature in eV
   Field3D Teff = (Tatom / Aatom + Tion / Aion) * Tnorm;
-  for (auto& i : Teff.getRegion("RGN_NOBNDRY")) {
+  for (auto& i : Teff.getRegion("RGN_ALL")) {
     if (Teff[i] < 0.01) {
       Teff[i] = 0.01;
     } else if (Teff[i] > 10000) {
@@ -44,8 +44,8 @@ void HydrogenChargeExchange::calculate_rates(GuardedOptions&& atom1, GuardedOpti
   // Optionally multiply by arbitrary multiplier
   const Field3D sigmav = exp(ln_sigmav) * (1e-6 * Nnorm / FreqNorm) * rate_multiplier;
 
-  const Field3D Natom = get<Field3D>(atom1["density"]);
-  const Field3D Nion = get<Field3D>(ion1["density"]);
+  const Field3D Natom = softFloor(get<Field3D>(atom1["density"]), 0.0);
+  const Field3D Nion = softFloor(get<Field3D>(ion1["density"]), 0.0);
 
   R = Natom * Nion * sigmav; // Rate coefficient in [m^-3 s^-1]
 

src/relax_potential.cxx

@@ -7,6 +7,7 @@
 #include <bout/derivs.hxx>
 #include <bout/difops.hxx>
 #include <bout/fv_ops.hxx>
+#include <bout/output_bout_types.hxx>
 
 using bout::globals::mesh;
 
@@ -843,6 +844,20 @@ void RelaxPotential::finally(const Options& state) {
 
     ddt(phi1) = lambda_1 * (phi_vort - Vort);
   }
+
+#if CHECKLEVEL >= 1
+  for (auto& i : Vort.getRegion("RGN_NOBNDRY")) {
+    if (!std::isfinite(ddt(Vort)[i])) {
+      throw BoutException("ddt(Vort) non-finite at {}\n", i);
+    }
+  }
+
+  for (auto& i : phi1.getRegion("RGN_NOBNDRY")) {
+    if (!std::isfinite(ddt(phi1)[i])) {
+      throw BoutException("ddt(phi1) non-finite at {}\n", i);
+    }
+  }
+#endif
 }
 
 void RelaxPotential::outputVars(Options& state) {

src/sheath_boundary_simple.cxx

@@ -43,9 +43,6 @@ BoutReal limitFree(BoutReal fm, BoutReal fc, BoutReal mode) {
     throw BoutException("Unknown boundary mode");
   }
 
-  return fp;  // Extrapolation
-
-
 #if CHECKLEVEL >= 2
   if (!std::isfinite(fp)) {
     throw BoutException("SheathBoundary limitFree: {}, {} -> {}", fm, fc, fp);

src/vorticity.cxx

@@ -4,11 +4,12 @@
 #include "../include/hermes_utils.hxx"
 
 #include <bout/constants.hxx>
-#include <bout/fv_ops.hxx>
-#include <bout/invert/laplacexy.hxx>
 #include <bout/derivs.hxx>
 #include <bout/difops.hxx>
+#include <bout/fv_ops.hxx>
+#include <bout/invert/laplacexy.hxx>
 #include <bout/invert_laplace.hxx>
+#include <bout/output_bout_types.hxx>
 
 using bout::globals::mesh;
 
@@ -901,6 +902,14 @@ void Vorticity::finally(const Options& state) {
       }
     }
   }
+
+#if CHECKLEVEL >= 1
+  for (auto& i : Vort.getRegion("RGN_NOBNDRY")) {
+    if (!std::isfinite(ddt(Vort)[i])) {
+      throw BoutException("ddt(Vort) non-finite at {}\n", i);
+    }
+  }
+#endif
 }
 
 void Vorticity::outputVars(Options& state) {