Neutral passive momentum

.gitignore

@@ -23,4 +23,7 @@ __pycache__
 spack-build*
 spack.lock
 .spack-env
-views
+views
+
+# vscode
+src/Hermes3_branches.code-workspace

include/anomalous_diffusion.hxx

@@ -52,6 +52,7 @@ private:
   Field2D anomalous_nu; ///< Anomalous momentum diffusion coefficient
 
   bool anomalous_sheath_flux; ///< Allow anomalous diffusion into sheath?
+  bool anomalous_transport_momentum;
 };
 
 

include/neutral_mixed.hxx

@@ -61,6 +61,22 @@ private:
   bool neutral_viscosity; ///< include viscosity?
   bool neutral_conduction; ///< Include heat conduction?
   bool evolve_momentum; ///< Evolve parallel momentum?
+  bool passive_momentum; ///< only evolve density, passive NVn and Tn=Ti
+  std::string temperature_from;
+
+  // --- Equilibrium NVn (used when passive_momentum = true) ---
+  // Caching sentinel: filled on first call to finally(), stays non-empty thereafter.
+  std::vector<std::string> equilibrium_momentum_collision_names;
+
+  // Collision frequency names — all stored on the neutral species' localstate
+  std::string equilibrium_nu_cx_name;   ///< Charge exchange collision freq name
+  std::string equilibrium_nu_iz_name;   ///< Ionisation collision freq name
+  std::string equilibrium_nu_rec_name;  ///< Recombination collision freq name e.g. "d+_d_rec"
+
+  // Normalisations (stored for use in finally())
+  BoutReal Nnorm, Tnorm, FreqNorm;
+  //BoutReal Nnorm_nu, Tnorm_nu, FreqNorm_nu;
+  Field3D nu_cx_out, nu_iz_out, nu_rec_out, Nn_eq_out; ///< Saved for output
 
   Field3D kappa_n, eta_n; ///< Neutral conduction and viscosity
 

src/anomalous_diffusion.cxx

@@ -17,6 +17,10 @@ AnomalousDiffusion::AnomalousDiffusion(std::string name, Options& alloptions, So
 
   Options& options = alloptions[name];
 
+  // easy way to fix anomalous + EM
+  anomalous_transport_momentum = options["anomalous_transport_momentum"]
+                    .doc("Anomalous transport of parallel momentum?")
+                    .withDefault(true);
   // Set in the mesh or options (or both)
   anomalous_D = 0.0;
   include_D = (mesh->get(anomalous_D, std::string("D_") + name) == 0)
@@ -68,8 +72,10 @@ void AnomalousDiffusion::transform(Options& state) {
                         : 0.0;
   Field2D T2D = DC(T);
 
+  //const Field3D V =
+  //    species.isSet("velocity") ? GET_NOBOUNDARY(Field3D, species["velocity"]) : 0.0;
   const Field3D V =
-      species.isSet("velocity") ? GET_NOBOUNDARY(Field3D, species["velocity"]) : 0.0;
+      (anomalous_transport_momentum && species.isSet("velocity")) ? GET_NOBOUNDARY(Field3D, species["velocity"]) : 0.0;
   Field2D V2D = DC(V);
 
   if (!anomalous_sheath_flux) {