Merge pull request #192 from bendudson/diffusive-2d-neutrals

Allow solving in 2D/3D without neutral momentum equation

Author: bendudson

include/neutral_mixed.hxx

@@ -52,6 +52,7 @@ private:
   BoutReal diffusion_limit;    ///< Maximum diffusion coefficient
 
   bool neutral_viscosity; ///< include viscosity?
+  bool evolve_momentum; ///< Evolve parallel momentum?
 
   bool precondition {true}; ///< Enable preconditioner?
   std::unique_ptr<Laplacian> inv; ///< Laplacian inversion used for preconditioning

src/neutral_mixed.cxx

@@ -31,7 +31,19 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver*
   // Evolving variables e.g name is "h" or "h+"
   solver->add(Nn, std::string("N") + name);
   solver->add(Pn, std::string("P") + name);
-  solver->add(NVn, std::string("NV") + name);
+  
+
+  evolve_momentum = options["evolve_momentum"]
+  .doc("Evolve parallel neutral momentum?")
+  .withDefault<bool>(true);
+
+  if (evolve_momentum) {
+    solver->add(NVn, std::string("NV") + name);
+  } else {
+    output_warn.write("WARNING: Not evolving neutral parallel momentum. NVn and Vn set to zero\n");
+    NVn = 0.0;
+    Vn = 0.0;
+  }
 
   sheath_ydown = options["sheath_ydown"]
                      .doc("Enable wall boundary conditions at ydown")
@@ -328,35 +340,42 @@ void NeutralMixed::finally(const Options& state) {
   }
   ddt(Nn) += Sn; // Always add density_source
 
-  /////////////////////////////////////////////////////
-  // Neutral momentum
-  TRACE("Neutral momentum");
+  if (evolve_momentum) {
 
-  ddt(NVn) =
-      -AA * FV::Div_par_fvv<hermes::Limiter>(Nnlim, Vn, sound_speed) // Momentum flow
-      - Grad_par(Pn)                                                 // Pressure gradient
-      + FV::Div_a_Grad_perp(DnnNVn, logPnlim) // Perpendicular diffusion
-      ;
+    /////////////////////////////////////////////////////
+    // Neutral momentum
+    TRACE("Neutral momentum");
 
-  if (neutral_viscosity) {
-    // NOTE: The following viscosity terms are are not (yet) balanced
-    //       by a viscous heating term
+    ddt(NVn) =
+        -AA * FV::Div_par_fvv<hermes::Limiter>(Nnlim, Vn, sound_speed) // Momentum flow
+        - Grad_par(Pn)                                                 // Pressure gradient
+        + FV::Div_a_Grad_perp(DnnNVn, logPnlim) // Perpendicular diffusion
+        ;
 
-    // Relationship between heat conduction and viscosity for neutral
-    // gas Chapman, Cowling "The Mathematical Theory of Non-Uniform
-    // Gases", CUP 1952 Ferziger, Kaper "Mathematical Theory of
-    // Transport Processes in Gases", 1972
-    // eta_n = (2. / 5) * kappa_n;
-    //
+    if (neutral_viscosity) {
+      // NOTE: The following viscosity terms are are not (yet) balanced
+      //       by a viscous heating term
 
-    ddt(NVn) += AA * FV::Div_a_Grad_perp((2. / 5) * DnnNn, Vn)    // Perpendicular viscosity
-              + AA * FV::Div_par_K_Grad_par((2. / 5) * DnnNn, Vn) // Parallel viscosity
-      ;
-  }
+      // Relationship between heat conduction and viscosity for neutral
+      // gas Chapman, Cowling "The Mathematical Theory of Non-Uniform
+      // Gases", CUP 1952 Ferziger, Kaper "Mathematical Theory of
+      // Transport Processes in Gases", 1972
+      // eta_n = (2. / 5) * kappa_n;
+      //
+
+      ddt(NVn) += AA * FV::Div_a_Grad_perp((2. / 5) * DnnNn, Vn)    // Perpendicular viscosity
+                + AA * FV::Div_par_K_Grad_par((2. / 5) * DnnNn, Vn) // Parallel viscosity
+        ;
+    }
 
-  if (localstate.isSet("momentum_source")) {
-    Snv = get<Field3D>(localstate["momentum_source"]);
-    ddt(NVn) += Snv;
+    if (localstate.isSet("momentum_source")) {
+      Snv = get<Field3D>(localstate["momentum_source"]);
+      ddt(NVn) += Snv;
+    }
+
+  } else {
+    ddt(NVn) = 0;
+    Snv = 0;
   }
 
   /////////////////////////////////////////////////////
@@ -533,6 +552,8 @@ void NeutralMixed::precon(const Options& state, BoutReal gamma) {
   inv->setCoefD(coef);
 
   ddt(Nn) = inv->solve(ddt(Nn));
-  ddt(NVn) = inv->solve(ddt(NVn));
+  if (evolve_momentum) {
+    ddt(NVn) = inv->solve(ddt(NVn));
+  }
   ddt(Pn) = inv->solve(ddt(Pn));
 }