diff --git a/.gitignore b/.gitignore index b95dfdea7..865391aee 100644 --- a/.gitignore +++ b/.gitignore @@ -23,4 +23,7 @@ __pycache__ spack-build* spack.lock .spack-env -views \ No newline at end of file +views + +# vscode +src/Hermes3_branches.code-workspace diff --git a/include/anomalous_diffusion.hxx b/include/anomalous_diffusion.hxx index e8c833592..03780b96f 100644 --- a/include/anomalous_diffusion.hxx +++ b/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; }; diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index 0a80aaad6..4759a30dd 100644 --- a/include/neutral_mixed.hxx +++ b/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 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 diff --git a/src/anomalous_diffusion.cxx b/src/anomalous_diffusion.cxx index cc585efaa..705059d3f 100644 --- a/src/anomalous_diffusion.cxx +++ b/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) { diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 5115c0bb4..75856f9bf 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -25,9 +25,10 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* const Options& units = alloptions["units"]; const BoutReal meters = units["meters"]; const BoutReal seconds = units["seconds"]; - const BoutReal Nnorm = units["inv_meters_cubed"]; - const BoutReal Tnorm = units["eV"]; - const BoutReal Omega_ci = 1. / units["seconds"].as(); + Nnorm = units["inv_meters_cubed"]; + Tnorm = units["eV"]; + FreqNorm = 1. / units["seconds"].as(); + const BoutReal Omega_ci = FreqNorm; // Need to take derivatives in X for cross-field diffusion terms ASSERT0(mesh->xstart > 0); @@ -41,9 +42,21 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* evolve_momentum = options["evolve_momentum"] .doc("Evolve parallel neutral momentum?") .withDefault(true); + passive_momentum = options["passive_momentum"] + .doc("Evolve only neutral density, passive momentum+pressure") + .withDefault(false); + temperature_from = options["temperature_from"] + .doc("Name of species to take temperature from. If not set, uses own species temperature.") + .withDefault(""); if (evolve_momentum) { solver->add(NVn, std::string("NV") + name); + } else if (passive_momentum) { + output_warn.write( + "WARNING: Not evolving neutral parallel momentum. " + "NVn set from diffusion + ion flow.\n"); + NVn = 0.0; + Vn = 0.0; } else { output_warn.write( "WARNING: Not evolving neutral parallel momentum. NVn and Vn set to zero\n"); @@ -181,12 +194,19 @@ void NeutralMixed::transform(Options& state) { NVn.clearParallelSlices(); Nn = floor(Nn, 0.0); - Pn = floor(Pn, 0.0); // Nnlim Used where division by neutral density is needed Nnlim = softFloor(Nn, density_floor); - Tn = Pn / Nnlim; - Tn.applyBoundary(); + // Optionally take temperature from another species (e.g. Tn = Ti) + if (!temperature_from.empty()) { + Tn = GET_NOBOUNDARY(Field3D, state["species"][temperature_from]["temperature"]); + Pn = Nn * Tn; + Pn = floor(Pn, 0.0); + } else { + Pn = floor(Pn, 0.0); + Tn = Pn / Nnlim; + Tn.applyBoundary(); + } Vn = NVn / (AA * Nnlim); Vn.applyBoundary("neumann"); @@ -275,10 +295,18 @@ void NeutralMixed::finally(const Options& state) { // Field3D logNn = log(Nn); // Field3D logTn = log(Tn); + // In passive_momentum mode, set_temperature has run its transform() by now, + // localstate["temperature"] holds the correct current Ti. Re-read Tn and + // recompute Pn -> update (Dnn, kappa_n, logPnlim) + if (passive_momentum) { + Tn = GET_VALUE(Field3D, localstate["temperature"]); + Pn = Nn * Tn; + Pnlim = softFloor(Pn, pressure_floor); + Pnlim.applyBoundary(); + } + logPnlim = log(Pnlim); logPnlim.applyBoundary(); - - /////////////////////////////////////////////////////// // Calculate cross-field diffusion from collision frequency // // @@ -454,41 +482,43 @@ 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) - + (2. / 3) * Vn * Grad_par(Pn) // Work done - + (5. / 3) * Div_a_Grad_perp_flows( // Perpendicular advection - DnnPn, logPnlim, - ef_adv_perp_xlow, ef_adv_perp_ylow) - ; - - // The factor here is 5/2 as we're advecting internal energy and pressure. - ef_adv_par_ylow *= 5./2; - ef_adv_perp_xlow *= 5./2; - ef_adv_perp_ylow *= 5./2; - - if (neutral_conduction) { - ddt(Pn) += (2. / 3) * Div_a_Grad_perp_flows( - kappa_n, Tn, // Perpendicular conduction - ef_cond_perp_xlow, ef_cond_perp_ylow) - - + (2. / 3) * Div_par_K_Grad_par_mod(kappa_n, Tn, // Parallel conduction - ef_cond_par_ylow, - false) // No conduction through target boundary - ; - // The factor here is likely 3/2 as this is pure energy flow, but needs checking. - ef_cond_perp_xlow *= 3./2; - ef_cond_perp_ylow *= 3./2; - ef_cond_par_ylow *= 3./2; - } + if (!passive_momentum) { + TRACE("Neutral pressure"); + + ddt(Pn) = - (5. / 3) * FV::Div_par_mod( // Parallel advection + Pn, Vn, sound_speed, ef_adv_par_ylow) + + (2. / 3) * Vn * Grad_par(Pn) // Work done + + (5. / 3) * Div_a_Grad_perp_flows( // Perpendicular advection + DnnPn, logPnlim, + ef_adv_perp_xlow, ef_adv_perp_ylow) + ; + + // The factor here is 5/2 as we're advecting internal energy and pressure. + ef_adv_par_ylow *= 5./2; + ef_adv_perp_xlow *= 5./2; + ef_adv_perp_ylow *= 5./2; + + if (neutral_conduction) { + ddt(Pn) += (2. / 3) * Div_a_Grad_perp_flows( + kappa_n, Tn, // Perpendicular conduction + ef_cond_perp_xlow, ef_cond_perp_ylow) + + + (2. / 3) * Div_par_K_Grad_par_mod(kappa_n, Tn, // Parallel conduction + ef_cond_par_ylow, + false) // No conduction through target boundary + ; + // The factor here is likely 3/2 as this is pure energy flow, but needs checking. + ef_cond_perp_xlow *= 3./2; + ef_cond_perp_ylow *= 3./2; + ef_cond_par_ylow *= 3./2; + } - Sp = pressure_source; - if (localstate.isSet("energy_source")) { - Sp += (2. / 3) * get(localstate["energy_source"]); + Sp = pressure_source; + if (localstate.isSet("energy_source")) { + Sp += (2. / 3) * get(localstate["energy_source"]); + } + ddt(Pn) += Sp; } - ddt(Pn) += Sp; if (evolve_momentum) { @@ -539,6 +569,77 @@ void NeutralMixed::finally(const Options& state) { Snv = 0; } + } else if (passive_momentum) { + // NVn = Nn_eq * Vi [equilibrium flow with ions] + // Ion parallel velocity + Field3D U = 0.0; + const std::string ion_name = std::string(name) + "+"; + if (state["species"].isSet(ion_name) && + state["species"][ion_name].isSet("velocity")) { + U = GET_VALUE(Field3D, state["species"][ion_name]["velocity"]); + } + Field3D Ne = GET_VALUE(Field3D, state["species"]["e"]["density"]); + //Field3D Ni = GET_VALUE(Field3D, state["species"][ion_name]["density"]); // for now assuming Ni=Ne + Field3D Te = GET_VALUE(Field3D, state["species"]["e"]["temperature"]); + + // Fetch collision frequencies; fall back to zero if not found + Field3D nu_cx = 0.0; + Field3D nu_iz = 0.0; + Field3D nu_rec = 0.0; + + // similar to Collisionality part + if (equilibrium_momentum_collision_names.empty()) { + if (localstate.isSet("collision_frequency")) { + for (const auto& collision : localstate["collision_frequencies"].getChildren()) { + const std::string coll_name = collision.second.name(); + + if (equilibrium_nu_cx_name.empty() && + collisionSpeciesMatch(coll_name, name, "+", "cx", "partial")) { + equilibrium_nu_cx_name = coll_name; + } + if (equilibrium_nu_iz_name.empty() && + collisionSpeciesMatch(coll_name, name, "+", "iz", "partial")) { + equilibrium_nu_iz_name = coll_name; + } + if (equilibrium_nu_rec_name.empty() && + collisionSpeciesMatch(coll_name, std::string(name) + "+", name, "rec", "partial")) { + equilibrium_nu_rec_name = coll_name; + } + //output_info.write( + // "\t{:s} passive_momentum: cx='{:s}' iz='{:s}' rec='{:s}'\n", + // name, + // equilibrium_nu_cx_name, equilibrium_nu_iz_name, equilibrium_nu_rec_name); + } + } + + // Cache a dummy name so the empty() check above is not re-triggered + equilibrium_momentum_collision_names.push_back("initialised"); + } + + if (!equilibrium_nu_cx_name.empty()) { + nu_cx = GET_VALUE(Field3D, localstate["collision_frequencies"][equilibrium_nu_cx_name]); + } + if (!equilibrium_nu_iz_name.empty()) { + nu_iz = GET_VALUE(Field3D, localstate["collision_frequencies"][equilibrium_nu_iz_name]); + } + if (!equilibrium_nu_rec_name.empty()) { + nu_rec = GET_VALUE(Field3D, localstate["collision_frequencies"][equilibrium_nu_rec_name]); + } + + // Equilibrium neutral density: + // Nn_eq = Ni * (Nn * nu_cx + Ne * nu_rec) / (Ni * nu_cx + Ne * nu_iz) + // Quasineutrality: Ni = Ne + Field3D Nn_eq = (Nnlim * nu_cx + Ne * nu_rec) + / softFloor(nu_cx + nu_iz, density_floor); + NVn = Nn_eq * U; + + // Save for output diagnostics + if (diagnose) { + nu_cx_out = nu_cx; + nu_iz_out = nu_iz; + nu_rec_out = nu_rec; + Nn_eq_out = Nn_eq; + } } else { ddt(NVn) = 0; Snv = 0; @@ -548,8 +649,11 @@ void NeutralMixed::finally(const Options& state) { if (state.isSet("scale_timederivs")) { Field3D scale_timederivs = get(state["scale_timederivs"]); ddt(Nn) *= scale_timederivs; - ddt(Pn) *= scale_timederivs; - ddt(NVn) *= scale_timederivs; + + if (!passive_momentum) { + ddt(Pn) *= scale_timederivs; + ddt(NVn) *= scale_timederivs; + } } if (freeze_low_density) { @@ -580,8 +684,10 @@ void NeutralMixed::finally(const Options& state) { const BoutReal meanNn = (1./6) * (2 * Nn[i] + Nn[i.xp()] + Nn[i.xm()] + Nn[i.yp()] + Nn[i.ym()]); const BoutReal factor = exp(- density_floor / meanNn); ddt(Nn)[i] = factor * ddt(Nn)[i] + (1. - factor) * Nn_s[i]; - ddt(Pn)[i] = factor * ddt(Pn)[i] + (1. - factor) * Pn_s[i]; - ddt(NVn)[i] = factor * ddt(NVn)[i] + (1. - factor) * NVn_s[i]; + if (!passive_momentum) { + ddt(Pn)[i] = factor * ddt(Pn)[i] + (1. - factor) * Pn_s[i]; + ddt(NVn)[i] = factor * ddt(NVn)[i] + (1. - factor) * NVn_s[i]; + } } } @@ -590,11 +696,13 @@ void NeutralMixed::finally(const Options& state) { if (!std::isfinite(ddt(Nn)[i])) { throw BoutException("ddt(N{}) non-finite at {}\n", name, i); } - if (!std::isfinite(ddt(Pn)[i])) { - throw BoutException("ddt(P{}) non-finite at {}\n", name, i); - } - if (!std::isfinite(ddt(NVn)[i])) { - throw BoutException("ddt(NV{}) non-finite at {}\n", name, i); + if (!passive_momentum) { + if (!std::isfinite(ddt(Pn)[i])) { + throw BoutException("ddt(P{}) non-finite at {}\n", name, i); + } + if (!std::isfinite(ddt(NVn)[i])) { + throw BoutException("ddt(NV{}) non-finite at {}\n", name, i); + } } } #endif @@ -650,16 +758,18 @@ void NeutralMixed::outputVars(Options& state) { {"conversion", Nnorm * Omega_ci}, {"long_name", std::string("Rate of change of ") + name + " number density"}, {"source", "neutral_mixed"}}); - set_with_attrs(state[std::string("ddt(P") + name + std::string(")")], ddt(Pn), - {{"time_dimension", "t"}, - {"units", "Pa s^-1"}, - {"conversion", Pnorm * Omega_ci}, - {"source", "neutral_mixed"}}); - set_with_attrs(state[std::string("ddt(NV") + name + std::string(")")], ddt(NVn), + if (!passive_momentum) { + set_with_attrs(state[std::string("ddt(P") + name + std::string(")")], ddt(Pn), + {{"time_dimension", "t"}, + {"units", "Pa s^-1"}, + {"conversion", Pnorm * Omega_ci}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("ddt(NV") + name + std::string(")")], ddt(NVn), {{"time_dimension", "t"}, {"units", "kg m^-2 s^-2"}, {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci}, {"source", "neutral_mixed"}}); + } // end if (!passive_momentum) } if (diagnose) { set_with_attrs(state[std::string("T") + name], Tn, @@ -683,20 +793,22 @@ void NeutralMixed::outputVars(Options& state) { {"standard_name", "density source"}, {"long_name", name + " number density source"}, {"source", "neutral_mixed"}}); - set_with_attrs(state[std::string("SP") + name], Sp, - {{"time_dimension", "t"}, - {"units", "Pa s^-1"}, - {"conversion", SI::qe * Tnorm * Nnorm * Omega_ci}, - {"standard_name", "pressure source"}, - {"long_name", name + " pressure source"}, - {"source", "neutral_mixed"}}); - set_with_attrs(state[std::string("SNV") + name], Snv, - {{"time_dimension", "t"}, - {"units", "kg m^-2 s^-2"}, - {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci}, - {"standard_name", "momentum source"}, - {"long_name", name + " momentum source"}, - {"source", "neutral_mixed"}}); + if (!passive_momentum) { + set_with_attrs(state[std::string("SP") + name], Sp, + {{"time_dimension", "t"}, + {"units", "Pa s^-1"}, + {"conversion", SI::qe * Tnorm * Nnorm * Omega_ci}, + {"standard_name", "pressure source"}, + {"long_name", name + " pressure source"}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("SNV") + name], Snv, + {{"time_dimension", "t"}, + {"units", "kg m^-2 s^-2"}, + {"conversion", SI::Mp * Nnorm * Cs0 * Omega_ci}, + {"standard_name", "momentum source"}, + {"long_name", name + " momentum source"}, + {"source", "neutral_mixed"}}); + } set_with_attrs(state[std::string("S") + name + std::string("_src")], density_source, {{"time_dimension", "t"}, {"units", "m^-3 s^-1"}, @@ -880,6 +992,43 @@ void NeutralMixed::outputVars(Options& state) { {"species", name}, {"source", "evolve_pressure"}}); } + if (passive_momentum) { + set_with_attrs(state[std::string("nu_cx_") + name], nu_cx_out, + {{"time_dimension", "t"}, + {"units", "s^-1"}, + {"conversion", Omega_ci}, + {"standard_name", "charge exchange frequency"}, + {"long_name", name + " charge exchange collision frequency"}, + {"species", name}, + {"source", "neutral_mixed"}}); + + set_with_attrs(state[std::string("nu_iz_") + name], nu_iz_out, + {{"time_dimension", "t"}, + {"units", "s^-1"}, + {"conversion", Omega_ci}, + {"standard_name", "ionisation frequency"}, + {"long_name", name + " ionisation collision frequency"}, + {"species", name}, + {"source", "neutral_mixed"}}); + + set_with_attrs(state[std::string("nu_rec_") + name], nu_rec_out, + {{"time_dimension", "t"}, + {"units", "s^-1"}, + {"conversion", Omega_ci}, + {"standard_name", "recombination frequency"}, + {"long_name", name + " recombination collision frequency"}, + {"species", name}, + {"source", "neutral_mixed"}}); + + set_with_attrs(state[std::string("Nn_eq_") + name], Nn_eq_out, + {{"time_dimension", "t"}, + {"units", "m^-3"}, + {"conversion", Nnorm}, + {"standard_name", "equilibrium neutral density"}, + {"long_name", name + " equilibrium neutral density"}, + {"species", name}, + {"source", "neutral_mixed"}}); + } } } @@ -887,6 +1036,11 @@ void NeutralMixed::precon([[maybe_unused]] const Options& state, BoutReal gamma) if (!precondition) { return; } + // Preconditioner only applies when evolving Pn and NVn. + // will get ddt(Nn) later + if (passive_momentum) { + return; + } // First matrix // ( I 0)