Fluid equation of state with density floor

.git-blame-ignore-revs

@@ -3,3 +3,4 @@ f5de620399e88f139e236e40675facfd0cfd6b13
 37efa12abe0961add8715b797d3926664977d453
 8e49d5933131ddcb2f97bb2ac7095a1be80ad7b7
 8b96b9318f166a6cab2ae19d7ee456abcee51458
+225df6a0a9371363e3b0ea22ca5d3d3f7e3bb181

include/evolve_momentum.hxx

@@ -28,7 +28,6 @@ private:
   std::string name;     ///< Short name of species e.g "e"
 
   Field3D NV;           ///< Species parallel momentum (normalised, evolving)
-  Field3D NV_err;       ///< Difference from momentum as input from solver
   Field3D NV_solver;    ///< Momentum as calculated in the solver
   Field3D V;            ///< Species parallel velocity
 

include/evolve_pressure.hxx

@@ -26,7 +26,6 @@ struct EvolvePressure : public Component {
   ///   - poloidal_flows       Include poloidal ExB flows? Default is true
   ///   - precon               Enable preconditioner? Note: solver may not use it even if
   ///                          enabled.
-  ///   - p_div_v              Use p * Div(v) form? Default is v * Grad(p) form
   ///   - thermal_conduction   Include parallel heat conduction? Default is true
   ///
   /// - P<name>  e.g. "Pe", "Pd+"
@@ -59,16 +58,15 @@ struct EvolvePressure : public Component {
 private:
   std::string name; ///< Short name of the species e.g. h+
 
-  Field3D P;    ///< Pressure (normalised)
-  Field3D T, N; ///< Temperature, density
+  Field3D P;        ///< Pressure (normalised)
+  Field3D P_solver; ///< Save to restore at the end
+  Field3D T, N;     ///< Temperature, density
 
   bool bndry_flux;
   bool neumann_boundary_average_z; ///< Apply neumann boundary with Z average?
   bool poloidal_flows;
   bool thermal_conduction; ///< Include thermal conduction?
 
-  bool p_div_v; ///< Use p*Div(v) form? False -> v * Grad(p)
-
   bool evolve_log; ///< Evolve logarithm of P?
   Field3D logP;    ///< Natural logarithm of P
 

include/hermes_utils.hxx

@@ -13,6 +13,12 @@ inline BoutReal floor(BoutReal value, BoutReal min) { return std::max(value, min
 /// The intention is to keep the RHS function differentiable
 ///
 /// Note: This function cannot be used with min = 0!
+///
+/// Uses value + min * exp(-value / min)
+///
+/// Alternatives include:
+///  - sqrt(value^2 + min^2)
+///
 inline BoutReal softFloor(BoutReal value, BoutReal min) {
   value = std::max(value, 0.0);
   return value + min * exp(-value / min);

include/neutral_full_velocity.hxx

@@ -34,6 +34,7 @@ private:
 
   Field2D Nn2D;                // Neutral gas density (evolving)
   Field2D Pn2D;                // Neutral gas pressure (evolving)
+  Field2D Pn2D_solver;         // Pn2D from the solver
   Vector2D Vn2D;               // Neutral gas velocity
   Vector2D Vn2D_contravariant; ///< Neutral gas velocity v^x, v^y, v^z
   Field2D Tn2D;

include/neutral_mixed.hxx

@@ -33,6 +33,7 @@ private:
   std::string name; ///< Species name
 
   Field3D Nn, Pn, NVn;            // Density, pressure and parallel momentum
+  Field3D Pn_solver;              ///< Saved to restore in finally
   Field3D Vn;                     ///< Neutral parallel velocity
   Field3D Tn;                     ///< Neutral temperature
   Field3D Nnlim, Pnlim, logPnlim; // Limited in regions of low density

src/evolve_momentum.cxx

@@ -60,7 +60,6 @@ EvolveMomentum::EvolveMomentum(std::string name, Options& alloptions, Solver* so
 
   // Set to zero so set for output
   momentum_source = 0.0;
-  NV_err = 0.0;
 
   substitutePermissions("name", {name});
   substitutePermissions("outputs", {"velocity", "momentum"});
@@ -83,7 +82,6 @@ void EvolveMomentum::transform_impl(GuardedOptions& state) {
   NV_solver = NV; // Save the momentum as calculated by the solver
   NV = AA * N * V; // Re-calculate consistent with V and N
   // Note: Now NV and NV_solver will differ when N < density_floor
-  NV_err = NV - NV_solver; // This is used in the finally() function
   set(species["momentum"], NV);
 }
 
@@ -123,7 +121,10 @@ void EvolveMomentum::finally(const Options& state) {
 
       const Field3D phi = get<Field3D>(state["fields"]["phi"]);
 
-      ddt(NV) = -Div_n_bxGrad_f_B_XPPM(NV, phi, bndry_flux, poloidal_flows,
+      Field3D NVint = AA * Nlim * V; // Evolving momentum with Nlim rather than N
+      NVint.setBoundaryTo(NV);
+
+      ddt(NV) = -Div_n_bxGrad_f_B_XPPM(NVint, phi, bndry_flux, poloidal_flows,
                                        true); // ExB drift
 
       // Parallel electric field
@@ -171,7 +172,7 @@ void EvolveMomentum::finally(const Options& state) {
   //    otherwise energy conservation is affected
   //  - using the same operator as in density and pressure equations doesn't work
   ddt(NV) -= AA * FV::Div_par_fvv<hermes::Limiter>(Nlim, V, fastest_wave, fix_momentum_boundary_flux);
-  
+
   // Parallel pressure gradient
   if (species.isSet("pressure")) {
     Field3D P = get<Field3D>(species["pressure"]);
@@ -216,12 +217,6 @@ void EvolveMomentum::finally(const Options& state) {
     ddt(NV) += momentum_source;
   }
 
-  // If N < density_floor then NV and NV_solver may differ
-  // -> Add term to force NV_solver towards NV
-  // Note: This correction is calculated in transform()
-  //       because NV may be modified by electromagnetic terms
-  ddt(NV) += NV_err;
-
   // Scale time derivatives
   if (state.isSet("scale_timederivs")) {
     ddt(NV) *= get<Field3D>(state["scale_timederivs"]);

src/evolve_pressure.cxx

@@ -78,10 +78,6 @@ EvolvePressure::EvolvePressure(std::string name, Options& alloptions, Solver* so
   poloidal_flows =
       options["poloidal_flows"].doc("Include poloidal ExB flow").withDefault<bool>(true);
 
-  p_div_v = options["p_div_v"]
-                .doc("Use p*Div(v) form? Default, false => v * Grad(p) form")
-                .withDefault<bool>(false);
-
   hyper_z = options["hyper_z"].doc("Hyper-diffusion in Z").withDefault(-1.0);
 
   hyper_z_T = options["hyper_z_T"]
@@ -221,11 +217,20 @@ void EvolvePressure::transform_impl(GuardedOptions& state) {
   // Not using density boundary condition
   N = getNoBoundary<Field3D>(species["density"]);
 
-  Field3D Pfloor = floor(P, 0.0);
-  T = Pfloor / softFloor(N, density_floor);
-  Pfloor = N * T; // Ensure consistency
-
-  set(species["pressure"], Pfloor);
+  // Nnlim is used where division by neutral density is needed
+  // The equation of state is modified at low density:
+  //
+  // e = Cv T Nlim / N    <- Specific internal energy
+  // p = N T
+  //
+  // The internal energy evolution of (N * e) is therefore
+  // evolving P_solver = Nlim * T
+  // rather than pressure P = N * T
+  T = floor(P, 0.0) / softFloor(N, density_floor);
+  P_solver = P; // Save solver variable to restore later
+  P = N * T;    // Equation of state
+
+  set(species["pressure"], P);
   set(species["temperature"], T);
 }
 
@@ -235,21 +240,20 @@ void EvolvePressure::finally(const Options& state) {
   const auto& species = state["species"][name];
 
   // Get updated pressure and temperature with boundary conditions
-  // Note: Retain pressures which fall below zero
-  P.clearParallelSlices();
-  P.setBoundaryTo(get<Field3D>(species["pressure"]));
-  Field3D Pfloor = floor(P, 0.0); // Restricted to never go below zero
-
+  P = get<Field3D>(species["pressure"]);
   T = get<Field3D>(species["temperature"]);
   N = get<Field3D>(species["density"]);
 
+  Field3D Nlim = softFloor(N, density_floor);
+  Field3D Pint = Nlim * T; // Internal energy uses limited density
+
   if (species.isSet("charge") and (fabs(get<BoutReal>(species["charge"])) > 1e-5)
       and state.isSection("fields") and state["fields"].isSet("phi")) {
     // Electrostatic potential set and species is charged -> include ExB flow
 
     Field3D phi = get<Field3D>(state["fields"]["phi"]);
 
-    ddt(P) = -Div_n_bxGrad_f_B_XPPM(P, phi, bndry_flux, poloidal_flows, true);
+    ddt(P) = -Div_n_bxGrad_f_B_XPPM(Pint, phi, bndry_flux, poloidal_flows, true);
   } else {
     ddt(P) = 0.0;
   }
@@ -266,33 +270,22 @@ void EvolvePressure::finally(const Options& state) {
       fastest_wave = sqrt(T / AA);
     }
 
-    if (p_div_v) {
-      // Use the P * Div(V) form
-      ddt(P) -= FV::Div_par_mod<hermes::Limiter>(P, V, fastest_wave, flow_ylow_advection);
+    // Use V * Grad(P) form
+    //
+    // The internal energy term Pint = Nlim * T
+    ddt(P) -= FV::Div_par_mod<hermes::Limiter>(Pint + (2. / 3) * P, V, fastest_wave,
+                                               flow_ylow_advection);
 
-      // Work done. This balances energetically a term in the momentum equation
-      E_PdivV = -Pfloor * Div_par(V);
-      ddt(P) += (2. / 3) * E_PdivV;
+    E_VgradP = V * Grad_par(P);
+    ddt(P) += (2. / 3) * E_VgradP;
 
-    } else {
-      // Use V * Grad(P) form
-      // Note: A mixed form has been tried (on 1D neon example)
-      //       -(4/3)*FV::Div_par(P,V) + (1/3)*(V * Grad_par(P) - P * Div_par(V))
-      //       Caused heating of charged species near sheath like p_div_v
-      ddt(P) -=
-          (5. / 3)
-          * FV::Div_par_mod<hermes::Limiter>(P, V, fastest_wave, flow_ylow_advection);
-
-      E_VgradP = V * Grad_par(P);
-      ddt(P) += (2. / 3) * E_VgradP;
-    }
     flow_ylow_advection *= 5. / 2; // Energy flow
     flow_ylow = flow_ylow_advection;
 
     if (state.isSection("fields") and state["fields"].isSet("Apar_flutter")) {
       // Magnetic flutter term
       const Field3D Apar_flutter = get<Field3D>(state["fields"]["Apar_flutter"]);
-      ddt(P) -= (5. / 3) * Div_n_g_bxGrad_f_B_XZ(P, V, -Apar_flutter);
+      ddt(P) -= Div_n_g_bxGrad_f_B_XZ(Pint + (2. / 3) * P, V, -Apar_flutter);
       ddt(P) += (2. / 3) * V * bracket(P, Apar_flutter, BRACKET_ARAKAWA);
     }
 
@@ -406,6 +399,11 @@ void EvolvePressure::finally(const Options& state) {
       flow_ylow += get<Field3D>(species["energy_flow_ylow"]);
     }
   }
+
+  // Restore solver pressure.
+  // Keep boundary conditions for post-processing.
+  P_solver.setBoundaryTo(P);
+  P = P_solver;
 }
 
 void EvolvePressure::outputVars(Options& state) {
@@ -465,30 +463,16 @@ void EvolvePressure::outputVars(Options& state) {
                     {"species", name},
                     {"source", "evolve_pressure"}});
 
-    if (p_div_v) {
-      if (E_PdivV.isAllocated()) {
-        set_with_attrs(state["E" + name + "_PdivV"], E_PdivV,
-                       {{"time_dimension", "t"},
-                        {"units", "W / m^-3"},
-                        {"conversion", Pnorm * Omega_ci},
-                        {"standard_name", "energy source"},
-                        {"long_name", name + " energy source due to pressure gradient"},
-                        {"species", name},
-                        {"source", "evolve_pressure"}});
-      }
-    } else {
-      if (E_VgradP.isAllocated()) {
-        set_with_attrs(state["E" + name + "_VgradP"], E_VgradP,
-                       {{"time_dimension", "t"},
-                        {"units", "W / m^-3"},
-                        {"conversion", Pnorm * Omega_ci},
-                        {"standard_name", "energy source"},
-                        {"long_name", name + " energy source due to pressure gradient"},
-                        {"species", name},
-                        {"source", "evolve_pressure"}});
-      }
+    if (E_VgradP.isAllocated()) {
+      set_with_attrs(state["E" + name + "_VgradP"], E_VgradP,
+                     {{"time_dimension", "t"},
+                      {"units", "W / m^-3"},
+                      {"conversion", Pnorm * Omega_ci},
+                      {"standard_name", "energy source"},
+                      {"long_name", name + " energy source due to pressure gradient"},
+                      {"species", name},
+                      {"source", "evolve_pressure"}});
     }
-
     if (flow_xlow.isAllocated()) {
       set_with_attrs(
           state[fmt::format("ef{}_tot_xlow", name)], flow_xlow,

src/neutral_full_velocity.cxx

@@ -298,12 +298,6 @@ void NeutralFullVelocity::transform_impl(GuardedOptions& state) {
   Nn2D = floor(Nn2D, 0.0);
   Pn2D = floor(Pn2D, 0.0);
 
-  // Non-zero floor can use a differentiable soft floor
-  Field2D Nnlim = softFloor(Nn2D, density_floor);
-
-  Tn2D = Pn2D / Nnlim;
-  Tn2D.applyBoundary("neumann");
-
   //////////////////////////////////////////////////////
   // 2D (X-Y) full velocity model
   //
@@ -347,6 +341,23 @@ void NeutralFullVelocity::transform_impl(GuardedOptions& state) {
   // V_{||n} = b dot V_n = Vn2D.y / (JB)
   Vnpar = Vn2D.y / (coord->J * coord->Bxy);
 
+  // Non-zero floor can use a differentiable soft floor
+  Field2D Nnlim = softFloor(Nn2D, density_floor);
+
+  // Equation of state at low density
+  //
+  // e = Cv T Nlim / N
+  // p = N T
+  //
+  // The internal energy evolution of (N * e) is therefore
+  // evolving Pn2D_solver = Nlim * T
+  // rather than pressure Pn2D = N * T
+
+  Tn2D = Pn2D / Nnlim;
+  Tn2D.applyBoundary("neumann");
+  Pn2D_solver = Pn2D; // Save solver variable to restore later
+  Pn2D = Tn2D * Nn2D; // Equation of state
+
   // Set values in the state
   auto localstate = state["species"][name];
   set(localstate["density"], Nn2D);
@@ -624,7 +635,7 @@ void NeutralFullVelocity::finally(const Options& state) {
 
   //////////////////////////////////////////////////////
   // Neutral pressure
-  ddt(Pn2D) = -adiabatic_index * Div(Vn2D, Pn2D)
+  ddt(Pn2D) = -Div(Vn2D, Nn2D_floor * Tn2D + (adiabatic_index - 1.) * Pn2D)
               + (adiabatic_index - 1.) * (Vn2D_contravariant * GradPn2D);
 
   if (constant_transport_coef) {
@@ -725,6 +736,9 @@ void NeutralFullVelocity::finally(const Options& state) {
   ddt(Vn2D).toContravariant();
   Vn2D.toContravariant();
 
+  // Restore solver Pn2D so that restart files have the correcte value
+  Pn2D = Pn2D_solver;
+
   // Set time derivatives to zero
   if (zero_timederivs) {
     Field2D zero{0.0};