diff --git a/.git-blame-ignore-revs b/.git-blame-ignore-revs index 5bc6bab8..ffb9768f 100644 --- a/.git-blame-ignore-revs +++ b/.git-blame-ignore-revs @@ -3,6 +3,7 @@ f5de620399e88f139e236e40675facfd0cfd6b13 37efa12abe0961add8715b797d3926664977d453 8e49d5933131ddcb2f97bb2ac7095a1be80ad7b7 8b96b9318f166a6cab2ae19d7ee456abcee51458 +225df6a0a9371363e3b0ea22ca5d3d3f7e3bb181 1069be8bda209261843b3bfa8583e3ee8545e152 171ae56b72fb7c8aff5cd02e66e9eac7315dfd4e 793edfad44d1b59be4cb141bd30a2b086e989369 diff --git a/include/evolve_momentum.hxx b/include/evolve_momentum.hxx index 02711488..37c09292 100644 --- a/include/evolve_momentum.hxx +++ b/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 diff --git a/include/evolve_pressure.hxx b/include/evolve_pressure.hxx index 32b2ef9b..5746476e 100644 --- a/include/evolve_pressure.hxx +++ b/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 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 diff --git a/include/hermes_utils.hxx b/include/hermes_utils.hxx index 8b17c2be..3f48cec5 100644 --- a/include/hermes_utils.hxx +++ b/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); diff --git a/include/neutral_full_velocity.hxx b/include/neutral_full_velocity.hxx index ccca3853..775df3fe 100644 --- a/include/neutral_full_velocity.hxx +++ b/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; diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index ca4f8b7c..f4c9753d 100644 --- a/include/neutral_mixed.hxx +++ b/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 diff --git a/src/evolve_momentum.cxx b/src/evolve_momentum.cxx index 68863212..4e216498 100644 --- a/src/evolve_momentum.cxx +++ b/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(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(Nlim, V, fastest_wave, fix_momentum_boundary_flux); - + // Parallel pressure gradient if (species.isSet("pressure")) { Field3D P = get(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(state["scale_timederivs"]); diff --git a/src/evolve_pressure.cxx b/src/evolve_pressure.cxx index ebfb606f..314a8662 100644 --- a/src/evolve_pressure.cxx +++ b/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(true); - p_div_v = options["p_div_v"] - .doc("Use p*Div(v) form? Default, false => v * Grad(p) form") - .withDefault(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(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(species["pressure"])); - Field3D Pfloor = floor(P, 0.0); // Restricted to never go below zero - + P = get(species["pressure"]); T = get(species["temperature"]); N = get(species["density"]); + Field3D Nlim = softFloor(N, density_floor); + Field3D Pint = Nlim * T; // Internal energy uses limited density + if (species.isSet("charge") and (fabs(get(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(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(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(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(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(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(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, diff --git a/src/neutral_full_velocity.cxx b/src/neutral_full_velocity.cxx index ff690f9a..85948f34 100644 --- a/src/neutral_full_velocity.cxx +++ b/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}; diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 6f632c8f..acc9aabf 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -222,10 +222,20 @@ void NeutralMixed::transform_impl(GuardedOptions& state) { Nn = floor(Nn, 0.0); Pn = floor(Pn, 0.0); - // Nnlim Used where division by neutral density is needed + // 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 Pn_solver = Nlim * Tn + // rather than pressure Pn = Nn * Tn Nnlim = softFloor(Nn, density_floor); - Tn = Pn / Nnlim; + Tn = Pn / Nnlim; // Internal energy Tn.applyBoundary(); + Pn_solver = Pn; // Save solver variable to restore later + Pn = Tn * Nn; // Equation of state, so Pn is now pressure Vn = NVn / (AA * Nnlim); Vn.applyBoundary("neumann"); @@ -329,6 +339,7 @@ void NeutralMixed::finally(const Options& state) { Pnlim = softFloor(Pn, pressure_floor); logPnlim = log(Pnlim); logPnlim.applyBoundary(); + /////////////////////////////////////////////////////// // Calculate cross-field diffusion from collision frequency // @@ -339,6 +350,7 @@ void NeutralMixed::finally(const Options& state) { if (collisionality_override > 0.0) { // user has set an override for collision frequency Dnn = (Tn / AA) / collisionality_override; + } else { if (localstate.isSet("collision_frequency")) { // Collisionality @@ -508,20 +520,21 @@ void NeutralMixed::finally(const Options& state) { // Neutral pressure TRACE("Neutral pressure"); - ddt(Pn) = -(5. / 3) - * FV::Div_par_mod( // Parallel advection - Pn, Vn, sound_speed, ef_adv_par_ylow) + // The equation of state is modified by the density floor, + // so the advection of internal energy and work done are combined as: + Field3D e_plus_p = Nnlim * Tn + (2. / 3) * Pn; + + ddt(Pn) = -FV::Div_par_mod( // Parallel advection + e_plus_p, Vn, sound_speed, ef_adv_par_ylow) + (2. / 3) * Vn * Grad_par(Pn); // Work done // Perpendicular advection of pressure if (nonorthogonal_operators) { - ddt(Pn) += - (5. / 3) - * Div_a_Grad_perp_nonorthog(DnnPn, logPnlim, ef_adv_perp_xlow, ef_adv_perp_ylow); + ddt(Pn) += Div_a_Grad_perp_nonorthog(Dnn * e_plus_p, logPnlim, ef_adv_perp_xlow, + ef_adv_perp_ylow); } else { - ddt(Pn) += - (5. / 3) - * Div_a_Grad_perp_flows(DnnPn, logPnlim, ef_adv_perp_xlow, ef_adv_perp_ylow); + ddt(Pn) += Div_a_Grad_perp_flows(Dnn * e_plus_p, logPnlim, ef_adv_perp_xlow, + ef_adv_perp_ylow); } // The factor here is 5/2 as we're advecting internal energy and pressure. @@ -672,6 +685,9 @@ void NeutralMixed::finally(const Options& state) { } } #endif + + // Restore solver Pn + Pn = Pn_solver; } void NeutralMixed::outputVars(Options& state) {