Fixes for 1D examples

examples/tokamak-1D/extra/1D-hydrogen/BOUT.inp

@@ -18,25 +18,31 @@ timestep = 10000
 MXG = 0  # No guard cells in X
 
 [mesh]
+# 1D simulation, use "y" as the dimension along the fieldline
 nx = 1
 ny = 400   # Resolution along field-line
 nz = 1
 
 length = 30           # Length of the domain in meters
-length_xpt = 10   # Length from midplane to X-point [m]
+length_xpt = 10   # Length from midplane to X-point [m] (i.e. this is where the source ends)
 
 dymin = 0.1  # Minimum grid spacing near target, as fraction of average. Must be > 0 and < 1
 
-# Parallel grid spacing
+# Parallel grid spacing — grid refinement near the divertor target
 dy = (length / ny) * (1 + (1-dymin)*(1-y/pi))
 
-# Calculate where the source ends in grid index
+# Calculate where the source ends in grid index (i.e. at the X-point)
 source = length_xpt / length
 y_xpt = pi * ( 2 - dymin - sqrt( (2-dymin)^2 - 4*(1-dymin)*source ) ) / (1 - dymin)
 
+# Separatrix index set to -1 for 1D simulations
 ixseps1 = -1
 ixseps2 = -1
 
+# In 1D, the Jacobian is inversely proportional to volume
+# Total flux expansion can be set up with a non-constant J profile
+J = 1
+
 [hermes]
 # Notes:
 #  - electrons after other species, so the density can be set by quasineutrality
@@ -57,22 +63,34 @@ Tnorm = 100
 [solver]
 type = beuler  # Backward Euler steady-state solver
 
-snes_type = newtonls
+timestep = 1e-2  # Initial timestep
+
+#### FIRST ORDER IN TIME, FASTER SOLVER BEULER
+type = beuler                    # Backward Euler steady-state solver
+snes_type = newtonls             # Newton iterations
+ksp_type = gmres                 # GMRES method
+pc_type = hypre                  # Preconditioner type
+pc_hypre_type = euclid           # Hypre preconditioner type
+max_nonlinear_iterations = 10    # Max Newton iterations per timestep
+lag_jacobian = 500               # How long to wait before Jacobian recalculation per Newton iteration
+####
+
+# Solver tolerances
+atol = 1e-7                      # Absolute tolerance (controls small numbers)
+rtol = 1e-5                      # Relative tolerance (primary convergence control)
 
-max_nonlinear_iterations = 10
+matrix_free_operator = true
 
 diagnose = false # Print diagnostic information?
 
-atol = 1e-7
-rtol = 1e-5
 
 [sheath_boundary]
 
 lower_y = false
 upper_y = true
 
 [neutral_parallel_diffusion]
-
+diffusion_collisions_mode = multispecies  # "afn" or "multispecies"
 dneut = 10   # (B / Bpol)^2 in neutral diffusion terms
 
 ####################################
@@ -108,15 +126,18 @@ function = 1
 source_shape = H(mesh:y_xpt - y)*1e20 # For upstream_density_feedback
 
 [Pd+]
+
+# Initial condition for ion pressure (in terms of hermes:Nnorm * hermes:Tnorm)
 function = 1
 
 powerflux = 2.5e7  # Input power flux in W/m^2
 
+# Divide by length of source region to get heat source in W/m3
 source = (powerflux*2/3 / (mesh:length_xpt))*H(mesh:y_xpt - y)  # Input power as function of y
 
 [NVd+]
 
-function = 0
+function = 0  # Initialise with Vi=0
 
 ####################################
 
@@ -130,38 +151,35 @@ AA = 2
 thermal_conduction = true
 
 [Nd]
-
 function = 0.001
 
 [Pd]
-
 function = 0.0001
 
 ####################################
 
 [e] # Electrons
-type = quasineutral, zero_current, evolve_pressure, noflow_boundary
+type = quasineutral, evolve_pressure, zero_current, noflow_boundary
 
 noflow_upper_y = false
 
 charge = -1
 AA = 1/1836
-
 thermal_conduction = true  # in evolve_pressure
+diagnose = true
 
 [Pe]
 
 function = `Pd+:function`  # Same as ion pressure initially
-
-source = `Pd+:source`  # Same as ion pressure source
+source = `Pd+:source`      # Same as ion pressure source
 
 ####################################
 
 [recycling]
-
 species = d+
 
 [reactions]
+diagnose = true
 type = (
         d + e -> d+ + 2e,     # Deuterium ionisation
         d+ + e -> d,          # Deuterium recombination

include/component.hxx

@@ -354,7 +354,7 @@ bool hermesDataInvalid([[maybe_unused]] const T& value) {
 /// Doesn't check boundary cells
 template<>
 inline bool hermesDataInvalid(const Field3D& value) {
-  for (auto& i : value.getRegion("RGN_NOBNDRY")) {
+  for (const auto& i : value.getRegion("RGN_NOBNDRY")) {
     if (!std::isfinite(value[i])) {
       return true;
     }

src/braginskii_collisions.cxx

@@ -18,7 +18,9 @@
 
 #include "../include/braginskii_collisions.hxx"
 #include "../include/component.hxx"
+#include "../include/guarded_options.hxx"
 #include "../include/hermes_utils.hxx"
+#include "../include/permissions.hxx"
 
 BraginskiiCollisions::BraginskiiCollisions(const std::string& name, Options& alloptions,
                                            Solver*)
@@ -370,7 +372,7 @@ void BraginskiiCollisions::transform_impl(GuardedOptions& state) {
             BoutReal const coulomb_log =
                 29.91
                 - log((Z1 * Z2 * (AA1 + AA2)) / (AA1 * Tlim2 + AA2 * Tlim1)
-                      * sqrt(Nlim1 * SQ(Z1) / Tlim1 + Nlim2 * SQ(Z2) / Tlim2));
+                      * sqrt((Nlim1 * SQ(Z1) / Tlim1) + (Nlim2 * SQ(Z2) / Tlim2)));
 
             // Calculate v_a^2, v_b^2
             const BoutReal v1sq = 2 * Tlim1 * SI::qe / mass1;
@@ -390,6 +392,10 @@ void BraginskiiCollisions::transform_impl(GuardedOptions& state) {
         } else {
           // species1 charged, species2 neutral
 
+          if (!ion_neutral) {
+            continue;
+          }
+
           // Scattering of charged species 1
           // Neutral density
           Field3D const Nn = GET_NOBOUNDARY(Field3D, species2["density"]);

src/neutral_parallel_diffusion.cxx

@@ -1,8 +1,11 @@
+#include "../include/component.hxx"
 #include "../include/neutral_parallel_diffusion.hxx"
 #include "../include/hermes_utils.hxx"
+#include "../include/guarded_options.hxx"
 
 #include <bout/constants.hxx>
 #include <bout/fv_ops.hxx>
+#include <bout/output_bout_types.hxx>
 
 using bout::globals::mesh;
 
@@ -89,6 +92,8 @@ void NeutralParallelDiffusion::transform_impl(GuardedOptions& state) {
     for (const auto& collision_name : collision_names[species_name]) {
       nu += GET_VALUE(Field3D, species["collision_frequencies"][collision_name]);
     }
+    // If nu is zero then Dn will be NaN
+    nu = floor(nu, 1e-5);
 
     // Diffusion coefficient
     BoutReal advection_factor = 0;