diff --git a/CMakeLists.txt b/CMakeLists.txt index 741bba680..b00a8bee4 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -89,6 +89,7 @@ set(HERMES_SOURCES src/noflow_boundary.cxx src/neutral_parallel_diffusion.cxx src/neutral_boundary.cxx + src/parallel_inertia_curvature_drift.cxx src/polarisation_drift.cxx src/solkit_neutral_parallel_diffusion.cxx src/hydrogen_charge_exchange.cxx @@ -137,6 +138,7 @@ set(HERMES_SOURCES include/neutral_parallel_diffusion.hxx include/solkit_neutral_parallel_diffusion.hxx include/noflow_boundary.hxx + include/parallel_inertia_curvature_drift.hxx include/polarisation_drift.hxx include/quasineutral.hxx include/radiation.hxx diff --git a/docs/sphinx/components.rst b/docs/sphinx/components.rst index 2b02e4bb2..87b22dfaf 100644 --- a/docs/sphinx/components.rst +++ b/docs/sphinx/components.rst @@ -210,6 +210,12 @@ By default parallel thermal conduction is included, which requires a collision time. If collisions are not calculated, then thermal conduction should be turned off by setting `thermal_conduction = false` in the input options. +The choice of collision frequency is set by the flag `conduction_collisions_mode`: `legacy` uses +all available collision frequencies involving the chosen species, while `braginskii` uses only +self-collisions .The default is `legacy` and it is recommended for use in multispecies. If you are solving +for a single ion and want to recover Braginskii, use the `braginskii` mode. + + If the component option ``diagnose = true`` then additional fields will be saved to the dump files: The species temperature ``T + name`` (e.g. ``Td+`` or ``Te``), the time derivative ``ddt(P + name)`` @@ -405,6 +411,11 @@ must be calculated after collisions: [hermes] components = ..., collisions, ion_viscosity +The choice of collision frequency is set by the flag `viscosity_collisions_mode`: `legacy` uses +all available collision frequencies involving the chosen species, while `braginskii` uses only +ii collisions. The default is `legacy` and it is recommended for use in multispecies. If you are solving +for a single ion and want to recover Braginskii, use the `braginskii` mode. + By default only the parallel diffusion of momentum is included, adding a force to each ion's momentum equation: @@ -639,7 +650,7 @@ with velocity :math:`v_D = -D \nabla_\perp N / N`. .. doxygenstruct:: AnomalousDiffusion :members: -Neutral gas models +2D Neutral transport ------------------ In 1D, neutral transport is currently done through the same components as for plasma, i.e. `evolve_density`, @@ -722,6 +733,14 @@ In an additional effort to limit the diffusivitiy to more physical values, a flu This formulation is equivalent to defining a :math:`D_n` with a free streaming velocity while accounting for the pseudo collisionality due to the maximum vessel mean free path :math:`l_{max}`. The flux limiter :math:`f_l` is set to 1.0 by default. + \mathbf{v}_{\perp n} = -D_{nn}\frac{1}{p_n}\nabla_\perp p_n + +The choice of collision frequency is set by the flag `diffusion_collisions_mode`: `legacy` uses +all available collision frequencies involving the chosen species, while `afn` uses only +CX and IZ rates. The default is `afn` and corresponds to the choice in UEDGE and +the SOLPS-ITER AFN (Advanced Fluid Neutral) model. + + Sources ------------------- Applying sources using the input file @@ -1239,6 +1258,12 @@ direction on the parallel transport, and is the `dneut` input setting. Currently use case for this component is to represent the neutrals diffusing orthogonal to the target wall, and it is recommended to set `dneut` according to the field line pitch at the target. +The choice of collision frequency is set by the flag `diffusion_collisions_mode`: `legacy` uses +all available collision frequencies involving the chosen species, while `afn` uses only +CX and IZ rates. The default is `afn` and corresponds to the choice in UEDGE +and the SOLPS-ITER AFN (Advanced Fluid Neutral) model. + + .. doxygenstruct:: NeutralParallelDiffusion :members: @@ -1296,7 +1321,7 @@ The frequency of charged species `a` colliding with charged species `b` is .. math:: - \nu_{ab} = \frac{1}{3\pi^{3/2}\epsilon_0^2}\frac{Z_a^2 Z_b^2 n_b \ln\Lambda}{\left(v_a^2 + v_b^2\right)^{3/2}}\frac{\left(1 + m_a / m_b\right)}{m_a^2} + \nu_{ab} = \frac{1}{3\pi^{3/2}\varepsilon_0^2}\frac{Z_a^2 Z_b^2 n_b \ln\Lambda}{\left(v_a^2 + v_b^2\right)^{3/2}}\frac{\left(1 + m_a / m_b\right)}{m_a^2} Note that the cgs expression in Hinton is divided by :math:`\left(4\pi\epsilon_0\right)^2` to get @@ -1310,6 +1335,15 @@ Note that with this definition we recover the `Braginskii expressions `_ for e-i and i-i collision times. +The electron-electron collision time definition follows Braginskii (note that Fitzpatrick uses +a different definition in his `notes `_, +these are not consistent with Braginskii): + +.. math:: + + \nu_{ee} = \frac{ln \Lambda e^4 n_e} { 12 \pi^{3/2} \varepsilon_0^2 m_{e}^{1/2} T_{e}^{3/2} } + + The electron-electron collision time definition follows Braginskii (note that Fitzpatrick uses a different definition in his `notes `_, these are not consistent with Braginskii): diff --git a/hermes-3.cxx b/hermes-3.cxx index ebf725505..bb810f40e 100644 --- a/hermes-3.cxx +++ b/hermes-3.cxx @@ -54,6 +54,7 @@ #include "include/neutral_mixed.hxx" #include "include/neutral_parallel_diffusion.hxx" #include "include/noflow_boundary.hxx" +#include "include/parallel_inertia_curvature_drift.hxx" #include "include/polarisation_drift.hxx" #include "include/quasineutral.hxx" #include "include/recycling.hxx" diff --git a/include/amjuel_helium.hxx b/include/amjuel_helium.hxx index d4b26900a..73e531f8e 100644 --- a/include/amjuel_helium.hxx +++ b/include/amjuel_helium.hxx @@ -23,11 +23,12 @@ struct AmjuelHeIonisation01 : public AmjuelReaction { void calculate_rates(Options& state, Field3D &reaction_rate, Field3D &momentum_exchange, + Field3D &heavy_particle_frequency, Field3D &electron_frequency, Field3D &energy_exchange, Field3D &energy_loss, BoutReal &rate_multiplier, BoutReal &radiation_multiplier); void transform(Options& state) override{ - Field3D reaction_rate, momentum_exchange, energy_exchange, energy_loss; - calculate_rates(state, reaction_rate, momentum_exchange, energy_exchange, energy_loss, rate_multiplier, radiation_multiplier); + Field3D heavy_particle_frequency, electron_frequency, reaction_rate, momentum_exchange, energy_exchange, energy_loss; + calculate_rates(state, heavy_particle_frequency, electron_frequency, reaction_rate, momentum_exchange, energy_exchange, energy_loss, rate_multiplier, radiation_multiplier); }; private: @@ -56,11 +57,12 @@ struct AmjuelHeRecombination10 : public AmjuelReaction { void calculate_rates(Options& state, Field3D &reaction_rate, Field3D &momentum_exchange, + Field3D &heavy_particle_frequency, Field3D &electron_frequency, Field3D &energy_exchange, Field3D &energy_loss, BoutReal &rate_multiplier, BoutReal &radiation_multiplier); void transform(Options& state) override{ - Field3D reaction_rate, momentum_exchange, energy_exchange, energy_loss; - calculate_rates(state, reaction_rate, momentum_exchange, energy_exchange, energy_loss, rate_multiplier, radiation_multiplier); + Field3D heavy_particle_frequency, electron_frequency, reaction_rate, momentum_exchange, energy_exchange, energy_loss; + calculate_rates(state, heavy_particle_frequency, electron_frequency, reaction_rate, momentum_exchange, energy_exchange, energy_loss, rate_multiplier, radiation_multiplier); }; private: diff --git a/include/amjuel_hyd_ionisation.hxx b/include/amjuel_hyd_ionisation.hxx index 7540937f7..f1af1040c 100644 --- a/include/amjuel_hyd_ionisation.hxx +++ b/include/amjuel_hyd_ionisation.hxx @@ -12,6 +12,7 @@ struct AmjuelHydIonisation : public AmjuelReaction { : AmjuelReaction(name, alloptions, solver) {} void calculate_rates(Options& electron, Options& atom, Options& ion, + Field3D &heavy_particle_frequency, Field3D &electron_frequency, Field3D& reaction_rate, Field3D& momentum_exchange, Field3D& energy_exchange, Field3D& energy_loss, BoutReal& rate_multiplier, BoutReal& radiation_multiplier); }; @@ -40,9 +41,9 @@ struct AmjuelHydIonisationIsotope : public AmjuelHydIonisation { Options& electron = state["species"]["e"]; Options& atom = state["species"][{Isotope}]; // e.g. "h" Options& ion = state["species"][{Isotope, '+'}]; // e.g. "h+" - Field3D reaction_rate, momentum_exchange, energy_exchange, energy_loss; + Field3D heavy_particle_frequency, electron_frequency, reaction_rate, momentum_exchange, energy_exchange, energy_loss; - calculate_rates(electron, atom, ion, reaction_rate, momentum_exchange, + calculate_rates(electron, atom, ion, heavy_particle_frequency, electron_frequency, reaction_rate, momentum_exchange, energy_exchange, energy_loss, rate_multiplier, radiation_multiplier); if (diagnose) { diff --git a/include/amjuel_hyd_recombination.hxx b/include/amjuel_hyd_recombination.hxx index d58df4e57..4a62ebe10 100644 --- a/include/amjuel_hyd_recombination.hxx +++ b/include/amjuel_hyd_recombination.hxx @@ -14,6 +14,7 @@ struct AmjuelHydRecombination : public AmjuelReaction { : AmjuelReaction(name, alloptions, solver) {} void calculate_rates(Options& electron, Options& atom, Options& ion, + Field3D &heavy_particle_frequency, Field3D &electron_frequency, Field3D& reaction_rate, Field3D& momentum_exchange, Field3D& energy_exchange, Field3D& energy_loss, BoutReal& rate_multiplier, BoutReal& radiation_multiplier); }; @@ -42,9 +43,9 @@ struct AmjuelHydRecombinationIsotope : public AmjuelHydRecombination { Options& electron = state["species"]["e"]; Options& atom = state["species"][{Isotope}]; // e.g. "h" Options& ion = state["species"][{Isotope, '+'}]; // e.g. "h+" - Field3D reaction_rate, momentum_exchange, energy_exchange, energy_loss; + Field3D heavy_particle_frequency, electron_frequency, reaction_rate, momentum_exchange, energy_exchange, energy_loss; - calculate_rates(electron, atom, ion, reaction_rate, momentum_exchange, + calculate_rates(electron, atom, ion, heavy_particle_frequency, electron_frequency, reaction_rate, momentum_exchange, energy_exchange, energy_loss, rate_multiplier, radiation_multiplier); if (diagnose) { diff --git a/include/amjuel_reaction.hxx b/include/amjuel_reaction.hxx index 92eb37241..b81d2190e 100644 --- a/include/amjuel_reaction.hxx +++ b/include/amjuel_reaction.hxx @@ -63,12 +63,22 @@ protected: /// Coefficients from Amjuel: /// - rate_coefs Double-polynomial log fit [T][n] for <σv> /// - radiation_coefs Double-polynomial log fit [T][n] for electron loss - /// electron_heating Energy added to electrons per reaction [eV] + /// - electron_heating Energy added to electrons per reaction [eV] + /// - heavy_particle_frequency Collisionality for the heavy particle (ion/neutral) [s^-1] + /// - electron_frequency Collisionality for the electron [s^-1] + /// - reaction_rate Density source rate in [m^-3 s^-1] + /// - momentum_exchange Momentum source rate [kg m^-2 s^-1] + /// - energy_exchange Energy source rate [W m^-3] + /// - energy_loss Electron energy loss (radiation or ionisation potential) [W m^-3] + /// - rate_multiplier User-set arbitrary rate multiplier [-] + /// - energy_loss User-set arbitrary multiplier on the radiation rate associated with the reaction [-] template void electron_reaction(Options& electron, Options& from_ion, Options& to_ion, const BoutReal (&rate_coefs)[rows][cols], const BoutReal (&radiation_coefs)[rows][cols], BoutReal electron_heating, + Field3D &heavy_particle_frequency, + Field3D &electron_frequency, Field3D &reaction_rate, Field3D &momentum_exchange, Field3D &energy_exchange, @@ -93,7 +103,15 @@ protected: const BoutReal to_charge = to_ion.isSet("charge") ? get(to_ion["charge"]) : 0.0; - // Calculate reaction rate using cell averaging. Optionally scale by multiplier + // Check for reaction type + std::string reaction_type; + if (from_ion.name().find("+") != std::string::npos) { + reaction_type = "rec"; + } else { + reaction_type = "iz"; + } + + // Calculate reaction rate in [m^3 s^-1] using cell averaging. Optionally scale by multiplier reaction_rate = cellAverage( [&](BoutReal ne, BoutReal n1, BoutReal te) { return ne * n1 * evaluate(rate_coefs, te * Tnorm, ne * Nnorm) * Nnorm @@ -101,6 +119,29 @@ protected: }, Ne.getRegion("RGN_NOBNDRY"))(Ne, N1, Te); + // Same as reaction_rate but without the n1 factor: returns atom/ion collisionality in [s^-1] + heavy_particle_frequency = cellAverage( + [&](BoutReal ne, BoutReal n1, BoutReal te) { + return ne * evaluate(rate_coefs, te * Tnorm, ne * Nnorm) * Nnorm + / FreqNorm * rate_multiplier; + }, + Ne.getRegion("RGN_NOBNDRY"))(Ne, N1, Te); + + // Same as reaction_rate but without the n1 factor: returns electron collisionality in [s^-1] + electron_frequency = cellAverage( + [&](BoutReal ne, BoutReal n1, BoutReal te) { + return n1 * evaluate(rate_coefs, te * Tnorm, ne * Nnorm) * Nnorm + / FreqNorm * rate_multiplier; + }, + Ne.getRegion("RGN_NOBNDRY"))(Ne, N1, Te); + + // Add individual reaction collision frequency to each species + if (reaction_type == "iz") { + set(from_ion["collision_frequencies"][from_ion.name() + std::string("_") + to_ion.name() + std::string("_iz")], heavy_particle_frequency); + } else if (reaction_type == "rec") { + set(to_ion["collision_frequencies"][from_ion.name() + std::string("_") + to_ion.name() + std::string("_rec")], heavy_particle_frequency); + } + // Particles // For ionisation, "from_ion" is the neutral and "to_ion" is the ion subtract(from_ion["density_source"], reaction_rate); diff --git a/include/evolve_energy.hxx b/include/evolve_energy.hxx index 331aad029..9eb0c903e 100644 --- a/include/evolve_energy.hxx +++ b/include/evolve_energy.hxx @@ -86,6 +86,9 @@ private: bool neumann_boundary_average_z; ///< Apply neumann boundary with Z average? bool poloidal_flows; bool thermal_conduction; ///< Include thermal conduction? + std::vector collision_names; ///< Collisions used for collisionality + std::string conduction_collisions_mode; ///< Collision selection, either multispecies or braginskii + Field3D nu; ///< Collision frequency for conduction BoutReal kappa_coefficient; ///< Leading numerical coefficient in parallel heat flux ///< calculation BoutReal kappa_limit_alpha; ///< Flux limit if >0 diff --git a/include/evolve_pressure.hxx b/include/evolve_pressure.hxx index a1b510710..2d81b3cae 100644 --- a/include/evolve_pressure.hxx +++ b/include/evolve_pressure.hxx @@ -79,6 +79,9 @@ private: bool neumann_boundary_average_z; ///< Apply neumann boundary with Z average? bool poloidal_flows; bool thermal_conduction; ///< Include thermal conduction? + std::vector collision_names; ///< Collisions used for collisionality + std::string conduction_collisions_mode; ///< Collision selection, either multispecies or braginskii + Field3D nu; ///< Collision frequency for conduction BoutReal kappa_coefficient; ///< Leading numerical coefficient in parallel heat flux calculation BoutReal kappa_limit_alpha; ///< Flux limit if >0 @@ -96,6 +99,8 @@ private: Field3D kappa_par; ///< Parallel heat conduction coefficient + Field3D conduction_div; ///< Divergence of heat conduction flux + Field3D source, final_source; ///< External pressure source Field3D Sp; ///< Total pressure source FieldGeneratorPtr source_prefactor_function; diff --git a/include/hermes_utils.hxx b/include/hermes_utils.hxx index 1efc086f0..b87b04173 100644 --- a/include/hermes_utils.hxx +++ b/include/hermes_utils.hxx @@ -26,11 +26,14 @@ inline T clamp(const T& var, BoutReal lo, BoutReal hi, const std::string& rgn = return result; } +// TODO: Replace the later SpeciesType with this one. This one +// doesn't work in elec + /// Enum that identifies the type of a species: electron, ion, neutral BOUT_ENUM_CLASS(SpeciesType, electron, ion, neutral); /// Identify species name string as electron, ion or neutral -inline SpeciesType identifySpeciesType(const std::string& species) { +inline SpeciesType identifySpeciesTypeEnum(const std::string& species) { if (species == "e") { return SpeciesType::electron; } else if ((species == "i") or @@ -49,3 +52,100 @@ Ind3D indexAt(const T& f, int x, int y, int z) { } #endif // HERMES_UTILS_H + +/// Function which returns true if any of a list of subtstrings is contained within string +inline bool containsAnySubstring(const std::string& mainString, const std::vector& substrings) { + for (const auto& subString : substrings) { + if (mainString.find(subString) != std::string::npos) { + return true; // Found at least one substring + } + } + return false; // None of the substrings found +} + +/// Identify species name string as electron, ion or neutral +inline std::string identifySpeciesType(const std::string& species) { + + std::string type = ""; + + if (species == "e") { + type = "electron"; + } else if (species.find(std::string("+")) != std::string::npos) { + type = "ion"; + } else { + type = "neutral"; + } + + return type; +} + + +/// Takes a string representing a collision, e.g. d+_e_coll +/// Splits it using underscores and finds species, e.g. d+ and e +/// Does partial match against these, e.g. True if species1 = d+ or + and species2 = e +/// Used across all processes that require collisions for identifying the right ones. +inline bool collisionSpeciesMatch(std::string input, const std::string& species1, const std::string& species2, const std::string& reaction, const std::string& mode) { + // Split the input string into substrings using underscore as delimiter + std::vector substrings; + size_t pos = 0; + std::string token; + while ((pos = input.find('_')) != std::string::npos) { + token = input.substr(0, pos); + substrings.push_back(token); + input.erase(0, pos + 1); // Erase string up to the current underscore + } + substrings.push_back(input); // Add the last substring after the last underscore + + bool species1_found = false; + bool species2_found = false; + bool reaction_found = false; + + // output << std::string("\n############################\n"); + + if (mode == "partial") { + if (substrings[0].find(species1) != std::string::npos) { + species1_found = true; + } + + if (substrings[1].find(species2) != std::string::npos) { + species2_found = true; + } + + if (substrings[2].find(reaction) != std::string::npos) { + reaction_found = true; + } + + + } else if (mode == "exact") { + // output << substrings[0] << std::string(" and ") << substrings[1] << std::string(" -- Testing for ") << species1 << std::string(" and ") << species2; + if (substrings[0] == species1) { + species1_found = true; + // output << std::string(" <-- ") << species1 << std::string(" FOUND"); + } + + if (substrings[1] == species2) { + species2_found = true; + // output << std::string(" <-- ") << species2 << std::string(" FOUND"); + } + + if (substrings[2] == reaction) { + reaction_found = true; + // output << std::string(" <-- ") << species2 << std::string(" FOUND"); + } + + // else if (mode == "species") { + + // } + + + // output << std::endl; + } else { + throw BoutException("Collision species match mode must be 'exact' or 'partial'"); + } + + // output << std::string("############################\n"); + + // Check if the first species matches species1 and the second species matches species2 + return (species1_found && species2_found && reaction_found); +} + diff --git a/include/hydrogen_charge_exchange.hxx b/include/hydrogen_charge_exchange.hxx index 1abe89ba1..1da468145 100644 --- a/include/hydrogen_charge_exchange.hxx +++ b/include/hydrogen_charge_exchange.hxx @@ -215,8 +215,8 @@ struct HydrogenChargeExchangeIsotope : public HydrogenChargeExchange { {"units", "s^-1"}, {"conversion", Omega_ci}, {"standard_name", "collision frequency"}, - {"long_name", (std::string("CX collision frequency between") + atom1 + " and " - + ion1 + " producing" + ion2 + " and" + atom2 + ". Note Kab != Kba")}, + {"long_name", (std::string("Collision frequency of CX of ") + atom1 + " and " + + ion1 + " producing " + ion2 + " and " + atom2 + ". Note Kab != Kba")}, {"source", "hydrogen_charge_exchange"}}); if (Isotope1 != Isotope2) { diff --git a/include/ion_viscosity.hxx b/include/ion_viscosity.hxx index 118caf423..f47536edc 100644 --- a/include/ion_viscosity.hxx +++ b/include/ion_viscosity.hxx @@ -53,6 +53,9 @@ struct IonViscosity : public Component { private: BoutReal eta_limit_alpha; ///< Flux limit coefficient bool perpendicular; ///< Include perpendicular flow? (Requires phi) + std::vector collision_names; ///< Collisions used for collisionality + std::string viscosity_collisions_mode; ///< Collision selection, either multispecies or braginskii + Field3D nu; ///< Collision frequency for conduction Vector2D Curlb_B; ///< Curvature vector Curl(b/B) bool diagnose; ///< Output additional diagnostics? diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index 652c5b00e..f7a2a7389 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -41,16 +41,20 @@ private: BoutReal AA; ///< Atomic mass (proton = 1) + std::vector collision_names; ///< Collisions used for collisionality + std::string diffusion_collisions_mode; ///< Collision selection, either afn or multispecies + Field3D nu; ///< Collisionality to use for diffusion Field3D Dnn; ///< Diffusion coefficient - Field3D DnnNn, DnnPn, DnnNVn; + Field3D DnnNn, DnnPn, DnnTn, DnnNVn; ///< Used for operators + BoutReal flux_limit; ///< Diffusive flux limit + BoutReal diffusion_limit; ///< Maximum diffusion coefficient bool sheath_ydown, sheath_yup; BoutReal density_floor; ///< Minimum Nn used when dividing NVn by Nn to get Vn. BoutReal pressure_floor; ///< Minimum Pn used when dividing Pn by Nn to get Tn. - BoutReal flux_limit; ///< Diffusive flux limit - BoutReal diffusion_limit; ///< Maximum diffusion coefficient + bool neutral_viscosity; ///< include viscosity? bool neutral_conduction; ///< Include heat conduction? diff --git a/include/neutral_parallel_diffusion.hxx b/include/neutral_parallel_diffusion.hxx index 5366e2a20..aaa182c12 100644 --- a/include/neutral_parallel_diffusion.hxx +++ b/include/neutral_parallel_diffusion.hxx @@ -33,16 +33,20 @@ struct NeutralParallelDiffusion : public Component { .doc("Output additional diagnostics?") .withDefault(false); + diffusion_collisions_mode = options["diffusion_collisions_mode"] + .doc("Can be afn: CX and IZ, or multispecies: CX and nn, ni, ne (if enabled)") + .withDefault("afn"); + equation_fix = options["equation_fix"] .doc("Fix correcting pressure advection and conductivity factors?") .withDefault(true); - thermal_conduction = options["thermal_conducton"] - .doc("Enable conduction?") + perpendicular_conduction = options["perpendicular_conduction"] + .doc("Enable parallel projection of perpendicular conduction?") .withDefault(true); - viscosity = options["viscosity"] - .doc("Enable viscosity?") + perpendicular_viscosity = options["perpendicular_viscosity"] + .doc("Enable parallel projection of perpendicular viscosity?") .withDefault(true); } @@ -72,9 +76,12 @@ private: BoutReal dneut; ///< cross-field diffusion projection (B / Bpol)^2 bool diagnose; ///< Output diagnostics? + std::vector collision_names; ///< Collisions used for collisionality + std::string diffusion_collisions_mode; ///< Collision selection, either afn or multispecies + Field3D nu; ///< Collision frequency for conduction bool equation_fix; ///< Fix incorrect 3/2 factor in pressure advection? - bool thermal_conduction; ///< Enable conduction? - bool viscosity; ///< Enable viscosity? + bool perpendicular_conduction; ///< Enable conduction? + bool perpendicular_viscosity; ///< Enable viscosity? /// Per-species diagnostics struct Diagnostics { diff --git a/include/parallel_inertia_curvature_drift.hxx b/include/parallel_inertia_curvature_drift.hxx new file mode 100644 index 000000000..4daa5afa2 --- /dev/null +++ b/include/parallel_inertia_curvature_drift.hxx @@ -0,0 +1,48 @@ +#pragma once +#ifndef PARALLEL_INERTIA_CURVATURE_DRIFT_H +#define PARALLEL_INERTIA_CURVATURE_DRIFT_H + +#include "component.hxx" + +/// Calculate parallel inertia curvature flow + +/// This term comes from the parallel inertia term in the momentum equation. +/// It is due to parallel inertia in a curved magnetic field (centrifugal-like force). +/// See ~Vpar^2 term in the paper by in A.V. Chankin, +/// Journal of Nuclear Materials 241-243 (1997) 199-213 + +struct ParallelInertiaCurvatureDrift : public Component { + ParallelInertiaCurvatureDrift(std::string name, Options &options, Solver *UNUSED(solver)); + + /// For every species, if it has: + /// - parallel velocity + /// - charge + /// + /// Modifies: + /// - density_source + /// - energy_source + /// - momentum_source + /// - currents + void transform(Options &state) override; + + /// Save variables to the output + void outputVars(Options& state) override; + +private: + Vector2D Curlb_B; + bool bndry_flux; + + // Diagnostic outputs + Field3D DivJicurv; // Divergence of parallel inertia curvature current + + bool diagnose; ///< Output additional diagnostics? + +}; + +namespace { +RegisterComponent registercomponentparallelinertiacurvature("parallel_inertia_curvature_drift"); +} + +#endif // PARALLEL_INERTIA_CURVATURE_DRIFT_H + + diff --git a/include/sheath_boundary_simple.hxx b/include/sheath_boundary_simple.hxx index bfb82bd23..6f4cf975c 100644 --- a/include/sheath_boundary_simple.hxx +++ b/include/sheath_boundary_simple.hxx @@ -73,6 +73,7 @@ struct SheathBoundarySimple : public Component { /// /// void transform(Options &state) override; + void outputVars(Options &state) override; private: BoutReal Ge; // Secondary electron emission coefficient BoutReal sin_alpha; // sin of angle between magnetic field and wall. @@ -88,6 +89,13 @@ private: Field3D wall_potential; ///< Voltage of the wall. Normalised units. + Field3D hflux_e; // Electron heat flux through sheath + Field3D phi; // Phi at sheath + Field3D ion_sum; // Sum of ion current at sheath + + bool diagnose; // Save diagnostic variables? + Options diagnostics; // Options object to store diagnostic fields like a dict + bool no_flow; ///< No flow speed, only remove energy BoutReal density_boundary_mode, pressure_boundary_mode, temperature_boundary_mode; ///< BC mode: 0=LimitFree, 1=ExponentialFree, 2=LinearFree diff --git a/src/amjuel_helium.cxx b/src/amjuel_helium.cxx index 7a2a03b6c..7be384644 100644 --- a/src/amjuel_helium.cxx +++ b/src/amjuel_helium.cxx @@ -67,6 +67,7 @@ static constexpr const BoutReal he01_radiation_coefs[9][9] = { -1.713225271579e-14}}; void AmjuelHeIonisation01::calculate_rates(Options& state, + Field3D &heavy_particle_frequency, Field3D &electron_frequency, Field3D &reaction_rate, Field3D &momentum_exchange, Field3D &energy_exchange, Field3D &energy_loss, BoutReal &rate_multiplier, BoutReal &radiation_multiplier) { electron_reaction(state["species"]["e"], @@ -74,7 +75,7 @@ void AmjuelHeIonisation01::calculate_rates(Options& state, state["species"]["he+"], // To helium ions he01_rate_coefs, he01_radiation_coefs, 0.0, // Note: Ionisation potential included in radiation_coefs, - reaction_rate, momentum_exchange, energy_exchange, energy_loss, rate_multiplier, radiation_multiplier + heavy_particle_frequency, electron_frequency, reaction_rate, momentum_exchange, energy_exchange, energy_loss, rate_multiplier, radiation_multiplier ); } @@ -146,6 +147,7 @@ static constexpr const BoutReal he10_radiation_coefs[9][9] = { -1.678705755876e-14}}; void AmjuelHeRecombination10::calculate_rates(Options& state, + Field3D &heavy_particle_frequency, Field3D &electron_frequency, Field3D &reaction_rate, Field3D &momentum_exchange, Field3D &energy_exchange, Field3D &energy_loss, BoutReal &rate_multiplier, BoutReal &radiation_multiplier) { electron_reaction(state["species"]["e"], @@ -153,6 +155,6 @@ void AmjuelHeRecombination10::calculate_rates(Options& state, state["species"]["he"], // To helium atoms he10_rate_coefs, he10_radiation_coefs, 24.586, // Ionisation potential heating of electrons - reaction_rate, momentum_exchange, energy_exchange, energy_loss, rate_multiplier, radiation_multiplier + heavy_particle_frequency, electron_frequency, reaction_rate, momentum_exchange, energy_exchange, energy_loss, rate_multiplier, radiation_multiplier ); } diff --git a/src/amjuel_hyd_ionisation.cxx b/src/amjuel_hyd_ionisation.cxx index 143b01254..7ff84f03b 100644 --- a/src/amjuel_hyd_ionisation.cxx +++ b/src/amjuel_hyd_ionisation.cxx @@ -66,10 +66,11 @@ static constexpr const BoutReal radiation_coefs[9][9] = { void AmjuelHydIonisation::calculate_rates( Options& electron, Options& atom, Options& ion, + Field3D &heavy_particle_frequency, Field3D &electron_frequency, Field3D &reaction_rate, Field3D &momentum_exchange, Field3D &energy_exchange, Field3D &energy_loss, BoutReal &rate_multiplier, BoutReal &radiation_multiplier) { electron_reaction(electron, atom, ion, rate_coefs, radiation_coefs, 0.0, // Note: Ionisation potential included in radiation_coefs - reaction_rate, momentum_exchange, energy_exchange, energy_loss, rate_multiplier, radiation_multiplier + heavy_particle_frequency, electron_frequency, reaction_rate, momentum_exchange, energy_exchange, energy_loss, rate_multiplier, radiation_multiplier ); } diff --git a/src/amjuel_hyd_recombination.cxx b/src/amjuel_hyd_recombination.cxx index 3808758a1..e36ba9223 100644 --- a/src/amjuel_hyd_recombination.cxx +++ b/src/amjuel_hyd_recombination.cxx @@ -66,10 +66,11 @@ static constexpr const BoutReal radiation_coefs[9][9] = { void AmjuelHydRecombination::calculate_rates( Options& electron, Options& atom, Options& ion, + Field3D &heavy_particle_frequency, Field3D &electron_frequency, Field3D &reaction_rate, Field3D &momentum_exchange, Field3D &energy_exchange, Field3D &energy_loss, BoutReal &rate_multiplier, BoutReal &radiation_multiplier) { electron_reaction(electron, ion, atom, rate_coefs, radiation_coefs, 13.6, // Potential energy loss [eV] heats electrons - reaction_rate, momentum_exchange, energy_exchange, energy_loss, rate_multiplier, radiation_multiplier + heavy_particle_frequency, electron_frequency, reaction_rate, momentum_exchange, energy_exchange, energy_loss, rate_multiplier, radiation_multiplier ); } diff --git a/src/collisions.cxx b/src/collisions.cxx index 241f12553..8dbf5ad4a 100644 --- a/src/collisions.cxx +++ b/src/collisions.cxx @@ -71,8 +71,10 @@ Collisions::Collisions(std::string name, Options& alloptions, Solver*) { void Collisions::collide(Options& species1, Options& species2, const Field3D& nu_12, BoutReal momentum_coefficient) { AUTO_TRACE(); - add(species1["collision_frequency"], nu_12); - set(collision_rates[species1.name()][species2.name()], nu_12); + add(species1["collision_frequency"], nu_12); // Total collision frequency + std::string coll_name = species1.name() + std::string("_") + species2.name() + std::string("_coll"); + set(species1["collision_frequencies"][coll_name], nu_12); // Collision frequency for individual reaction + set(collision_rates[species1.name()][species2.name()], nu_12); // Individual collision frequency used for diagnostics if (&species1 != &species2) { // For collisions between different species @@ -88,8 +90,10 @@ void Collisions::collide(Options& species1, Options& species2, const Field3D& nu return nu_12[i] * (A1 / A2) * density1[i] / floor(density2[i], 1e-5); }); - add(species2["collision_frequency"], nu); - set(collision_rates[species2.name()][species1.name()], nu); + add(species2["collision_frequency"], nu); // Total collision frequency + std::string coll_name = species2.name() + std::string("_") + species1.name() + std::string("_coll"); + set(species2["collision_frequencies"][coll_name], nu); // Collision frequency for individual reaction + set(collision_rates[species2.name()][species1.name()], nu); // Individual collision frequency used for diagnostics // Momentum exchange if (isSetFinalNoBoundary(species1["velocity"]) or diff --git a/src/evolve_energy.cxx b/src/evolve_energy.cxx index 9e4a058c9..c49e78656 100644 --- a/src/evolve_energy.cxx +++ b/src/evolve_energy.cxx @@ -31,6 +31,11 @@ EvolveEnergy::EvolveEnergy(std::string name, Options& alloptions, Solver* solver .withDefault(false); density_floor = options["density_floor"].doc("Minimum density floor").withDefault(1e-7); + + conduction_collisions_mode = options["conduction_collisions_mode"] + .doc("Can be multispecies: all collisions, or braginskii: self collisions and ie") + .withDefault("multispecies"); + if (evolve_log) { // Evolve logarithm of energy solver->add(logE, std::string("logE") + name); @@ -275,8 +280,84 @@ void EvolveEnergy::finally(const Options& state) { // Parallel heat conduction if (thermal_conduction) { + // Collisionality + // Braginskii mode: plasma - self collisions and ei, neutrals - CX, IZ + if (collision_names.empty()) { /// Calculate only once - at the beginning + + if (conduction_collisions_mode == "braginskii") { + for (const auto& collision : species["collision_frequencies"].getChildren()) { + + std::string collision_name = collision.second.name(); + + if (identifySpeciesType(species.name()) == "neutral") { + if (/// Charge exchange + (collisionSpeciesMatch( + collision_name, species.name(), "+", "cx", "partial")) or + /// Ionisation + (collisionSpeciesMatch( + collision_name, species.name(), "+", "iz", "partial"))) { + collision_names.push_back(collision_name); + } + + } else if (identifySpeciesType(species.name()) == "electron") { + if (/// Electron-electron collisions + (collisionSpeciesMatch( + collision_name, species.name(), "e", "coll", "exact"))) { + collision_names.push_back(collision_name); + } + + } else if (identifySpeciesType(species.name()) == "ion") { + if (/// Self-collisions + (collisionSpeciesMatch( + collision_name, species.name(), species.name(), "coll", "exact"))) { + collision_names.push_back(collision_name); + } + } + + } + // Multispecies mode: all collisions and CX are included + } else if (conduction_collisions_mode == "multispecies") { + for (const auto& collision : species["collision_frequencies"].getChildren()) { + + std::string collision_name = collision.second.name(); + + if (/// Charge exchange + (collisionSpeciesMatch( + collision_name, species.name(), "", "cx", "partial")) or + /// Any collision (en, in, ee, ii, nn) + (collisionSpeciesMatch( + collision_name, species.name(), "", "coll", "partial"))) { + collision_names.push_back(collision_name); + } + } + + } else { + throw BoutException("\tconduction_collisions_mode for {:s} must be either multispecies or braginskii", species.name()); + } + + if (collision_names.empty()) { + throw BoutException("\tNo collisions found for {:s} in evolve_pressure for selected collisions mode", species.name()); + } + + /// Write chosen collisions to log file + output_info.write("\t{:s} conduction collisionality mode: '{:s}' using ", + species.name(), conduction_collisions_mode); + for (const auto& collision : collision_names) { + output_info.write("{:s} ", collision); + } + + output_info.write("\n"); + + } + + /// Collect the collisionalities based on list of names + nu = 0; + for (const auto& collision_name : collision_names) { + nu += GET_VALUE(Field3D, species["collision_frequencies"][collision_name]); + } + // Calculate ion collision times - const Field3D tau = 1. / floor(get(species["collision_frequency"]), 1e-10); + const Field3D tau = 1. / floor(nu, 1e-10); const BoutReal AA = get(species["AA"]); // Atomic mass // Parallel heat conduction diff --git a/src/evolve_pressure.cxx b/src/evolve_pressure.cxx index 6f0607b4d..56e663ebe 100644 --- a/src/evolve_pressure.cxx +++ b/src/evolve_pressure.cxx @@ -41,6 +41,10 @@ EvolvePressure::EvolvePressure(std::string name, Options& alloptions, Solver* so .doc("Perpendicular diffusion at low pressure") .withDefault(false); + conduction_collisions_mode = options["conduction_collisions_mode"] + .doc("Can be multispecies: all collisions, or braginskii: self collisions and ie") + .withDefault("multispecies"); + if (evolve_log) { // Evolve logarithm of pressure solver->add(logP, std::string("logP") + name); @@ -120,6 +124,7 @@ EvolvePressure::EvolvePressure(std::string name, Options& alloptions, Solver* so source_prefactor_function = FieldFactory::get()->parse(str, &p_options); } + if (p_options["source_only_in_core"] .doc("Zero the source outside the closed field-line region?") .withDefault(false)) { @@ -155,7 +160,7 @@ EvolvePressure::EvolvePressure(std::string name, Options& alloptions, Solver* so BoutReal default_kappa; // default conductivity, changes depending on species - switch(identifySpeciesType(name)) { + switch(identifySpeciesTypeEnum(name)) { case SpeciesType::ion: default_kappa = 3.9; break; @@ -252,6 +257,8 @@ void EvolvePressure::finally(const Options& state) { T = get(species["temperature"]); N = get(species["density"]); + + Coordinates* coord = mesh->getCoordinates(); if (species.isSet("charge") and (fabs(get(species["charge"])) > 1e-5) and state.isSection("fields") and state["fields"].isSet("phi")) { @@ -337,8 +344,95 @@ void EvolvePressure::finally(const Options& state) { // Parallel heat conduction if (thermal_conduction) { - // Calculate ion collision times - const Field3D tau = 1. / floor(get(species["collision_frequency"]), 1e-10); + // Collisionality + // Braginskii mode: plasma - self collisions and ei, neutrals - CX, IZ + if (collision_names.empty()) { /// Calculate only once - at the beginning + + if (conduction_collisions_mode == "braginskii") { + for (const auto& collision : species["collision_frequencies"].getChildren()) { + + std::string collision_name = collision.second.name(); + + if (identifySpeciesType(species.name()) == "neutral") { + throw BoutException("\tBraginskii conduction collisions mode not available for neutrals, choose multispecies or afn"); + } else if (identifySpeciesType(species.name()) == "electron") { + if (/// Electron-electron collisions + (collisionSpeciesMatch( + collision_name, species.name(), "e", "coll", "exact"))) { + collision_names.push_back(collision_name); + } + + } else if (identifySpeciesType(species.name()) == "ion") { + if (/// Self-collisions + (collisionSpeciesMatch( + collision_name, species.name(), species.name(), "coll", "exact"))) { + collision_names.push_back(collision_name); + } + } + } + + // Multispecies mode: all collisions and CX are included + } else if (conduction_collisions_mode == "multispecies") { + for (const auto& collision : species["collision_frequencies"].getChildren()) { + + std::string collision_name = collision.second.name(); + + if (/// Charge exchange + (collisionSpeciesMatch( + collision_name, species.name(), "", "cx", "partial")) or + /// Any collision (en, in, ee, ii, nn) + (collisionSpeciesMatch( + collision_name, species.name(), "", "coll", "partial"))) { + collision_names.push_back(collision_name); + } + } + + } else if (conduction_collisions_mode == "afn") { + for (const auto& collision : species["collision_frequencies"].getChildren()) { + + std::string collision_name = collision.second.name(); + + if (identifySpeciesType(species.name()) != "neutral") { + throw BoutException("\tAFN conduction collisions mode not available for ions or electrons, choose braginskii or multispecies"); + } + if (/// Charge exchange + (collisionSpeciesMatch( + collision_name, species.name(), "+", "cx", "partial")) or + /// Ionisation + (collisionSpeciesMatch( + collision_name, species.name(), "+", "iz", "partial"))) { + collision_names.push_back(collision_name); + } + } + + } else { + throw BoutException("\tConduction_collisions_mode incorrect", species.name()); + } + + if (collision_names.empty()) { + throw BoutException("\tNo collisions found for {:s} in evolve_pressure for selected collisions mode", species.name()); + } + + /// Write chosen collisions to log file + output_info.write("\t{:s} conduction collisionality mode: '{:s}' using ", + species.name(), conduction_collisions_mode); + for (const auto& collision : collision_names) { + output_info.write("{:s} ", collision); + } + + output_info.write("\n"); + + } + + /// Collect the collisionalities based on list of names + nu = 0; + for (const auto& collision_name : collision_names) { + nu += GET_VALUE(Field3D, species["collision_frequencies"][collision_name]); + } + + + // Calculate ion collision times + const Field3D tau = 1. / floor(nu, 1e-10); const BoutReal AA = get(species["AA"]); // Atomic mass // Parallel heat conduction @@ -390,7 +484,8 @@ void EvolvePressure::finally(const Options& state) { // Note: Flux through boundary turned off, because sheath heat flux // is calculated and removed separately flow_ylow_conduction; - ddt(P) += (2. / 3) * Div_par_K_Grad_par_mod(kappa_par, T, flow_ylow_conduction, false); + conduction_div = (2. / 3) * Div_par_K_Grad_par_mod(kappa_par, T, flow_ylow_conduction, false); + ddt(P) += conduction_div; flow_ylow += flow_ylow_conduction; if (state.isSection("fields") and state["fields"].isSet("Apar_flutter")) { @@ -505,6 +600,23 @@ void EvolvePressure::outputVars(Options& state) { {"long_name", name + " heat conduction coefficient"}, {"species", name}, {"source", "evolve_pressure"}}); + + set_with_attrs(state[std::string("K") + name + std::string("_cond")], nu, + {{"time_dimension", "t"}, + {"units", "s^-1"}, + {"conversion", Omega_ci}, + {"long_name", "collision frequency for conduction"}, + {"species", name}, + {"source", "evolve_pressure"}}); + + set_with_attrs(state[std::string("div_cond_par_") + name], conduction_div * 3/2, + {{"time_dimension", "t"}, + {"units", "W m^-3"}, + {"conversion", Pnorm * Omega_ci}, + {"long_name", name + " parallel energy flow divergence due to heat conduction"}, + {"species", name}, + {"source", "evolve_pressure"}}); + } set_with_attrs(state[std::string("T") + name], T, {{"time_dimension", "t"}, diff --git a/src/hydrogen_charge_exchange.cxx b/src/hydrogen_charge_exchange.cxx index 39551b288..39c89827a 100644 --- a/src/hydrogen_charge_exchange.cxx +++ b/src/hydrogen_charge_exchange.cxx @@ -47,7 +47,7 @@ void HydrogenChargeExchange::calculate_rates(Options& atom1, Options& ion1, const Field3D Natom = floor(get(atom1["density"]), 1e-5); const Field3D Nion = floor(get(ion1["density"]), 1e-5); - R = Natom * Nion * sigmav; // Rate coefficient. This is an output parameter. + R = Natom * Nion * sigmav; // Rate coefficient in [m^-3 s^-1] if ((&atom1 != &atom2) or (&ion1 != &ion2)) { // Transfer particles atom1 -> ion2, ion1 -> atom2 @@ -94,9 +94,15 @@ void HydrogenChargeExchange::calculate_rates(Options& atom1, Options& ion1, subtract(ion1["energy_source"], ion_energy); add(atom2["energy_source"], ion_energy); - // Update collision frequency for the two colliding species in s^-1 - atom_rate = Nion * sigmav; - ion_rate = Natom * sigmav; + // Update collision frequency for the two colliding species + atom_rate = Nion * sigmav; // [s^-1] + ion_rate = Natom * sigmav; // [s^-1] + + // Add to total collision frequency add(atom1["collision_frequency"], atom_rate); add(ion1["collision_frequency"], ion_rate); + + // Set individual collision frequencies + set(atom1["collision_frequencies"][atom1.name() + std::string("_") + ion1.name() + std::string("_cx")], atom_rate); + set(ion1["collision_frequencies"][ion1.name() + std::string("_") + atom1.name() + std::string("_cx")], ion_rate); } diff --git a/src/ion_viscosity.cxx b/src/ion_viscosity.cxx index 12be7e6ce..8e182ebfc 100644 --- a/src/ion_viscosity.cxx +++ b/src/ion_viscosity.cxx @@ -5,6 +5,7 @@ #include #include #include +#include "../include/hermes_utils.hxx" #include "../include/ion_viscosity.hxx" #include "../include/div_ops.hxx" @@ -25,6 +26,10 @@ IonViscosity::IonViscosity(std::string name, Options& alloptions, Solver*) { perpendicular = options["perpendicular"] .doc("Include perpendicular flow? (Requires phi)") .withDefault(false); + + viscosity_collisions_mode = options["viscosity_collisions_mode"] + .doc("Can be multispecies: all collisions, or braginskii: self collisions") + .withDefault("multispecies"); if (perpendicular) { // Read curvature vector @@ -92,7 +97,65 @@ void IonViscosity::transform(Options &state) { continue; } - const Field3D tau = 1. / get(species["collision_frequency"]); + if (collision_names.empty()) { /// Calculate only once - at the beginning + + if (viscosity_collisions_mode == "braginskii") { + for (const auto& collision : species["collision_frequencies"].getChildren()) { + + std::string collision_name = collision.second.name(); + + if (/// Self-collisions + (collisionSpeciesMatch( + collision_name, species.name(), species.name(), "coll", "exact")) or + /// Ion-electron collisions + (collisionSpeciesMatch( + collision_name, species.name(), "+", "coll", "partial"))) { + + collision_names.push_back(collision_name); + } + } + // Multispecies mode: all collisions and CX are included + } else if (viscosity_collisions_mode == "multispecies") { + for (const auto& collision : species["collision_frequencies"].getChildren()) { + + std::string collision_name = collision.second.name(); + + if (/// Charge exchange + (collisionSpeciesMatch( + collision_name, species.name(), "", "cx", "partial")) or + /// Any collision (en, in, ee, ii, nn) + (collisionSpeciesMatch( + collision_name, species.name(), "", "coll", "partial"))) { + + collision_names.push_back(collision_name); + } + } + } else { + throw BoutException("\tviscosity_collisions_mode for {:s} must be either multispecies or braginskii", species.name()); + } + + if (collision_names.empty()) { + throw BoutException("\tNo collisions found for {:s} in ion_viscosity for selected collisions mode", species.name()); + } + + /// Write chosen collisions to log file + output_info.write("\t{:s} viscosity collisionality mode: '{:s}' using ", + species.name(), viscosity_collisions_mode); + for (const auto& collision : collision_names) { + output_info.write("{:s} ", collision); + } + + output_info.write("\n"); + + } + + /// Collect the collisionalities based on list of names + nu = 0; + for (const auto& collision_name : collision_names) { + nu += GET_VALUE(Field3D, species["collision_frequencies"][collision_name]); + } + + const Field3D tau = 1. / nu; const Field3D P = get(species["pressure"]); const Field3D V = get(species["velocity"]); diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index e670f9341..ee3f3b5a5 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -5,6 +5,7 @@ #include #include +#include "../include/hermes_utils.hxx" #include "../include/div_ops.hxx" #include "../include/hermes_build_config.hxx" #include "../include/neutral_mixed.hxx" @@ -88,6 +89,10 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* .doc("Include neutral gas heat conduction?") .withDefault(true); + diffusion_collisions_mode = options["diffusion_collisions_mode"] + .doc("Can be multispecies: all enabled collisions excl. IZ, or afn: CX, IZ and NN collisions") + .withDefault("multispecies"); + if (precondition) { inv = std::unique_ptr(Laplacian::create(&options["precon_laplace"])); @@ -285,8 +290,70 @@ void NeutralMixed::finally(const Options& state) { sqrt(floor(Tn, 1e-5) / AA) / neutral_lmax; // Neutral-neutral collisions [normalised frequency] if (localstate.isSet("collision_frequency")) { + + // Collisionality + // Braginskii mode: plasma - self collisions and ei, neutrals - CX, IZ + if (collision_names.empty()) { /// Calculate only once - at the beginning + + if (diffusion_collisions_mode == "afn") { + for (const auto& collision : localstate["collision_frequencies"].getChildren()) { + + std::string collision_name = collision.second.name(); + + if (/// Charge exchange + (collisionSpeciesMatch( + collision_name, name, "+", "cx", "partial")) or + /// Ionisation + (collisionSpeciesMatch( + collision_name, name, "+", "iz", "partial")) or + /// Neutral-neutral collisions + (collisionSpeciesMatch( + collision_name, name, name, "coll", "exact"))) { + collision_names.push_back(collision_name); + } + } + // Multispecies mode: all collisions and CX are included + } else if (diffusion_collisions_mode == "multispecies") { + for (const auto& collision : localstate["collision_frequencies"].getChildren()) { + + std::string collision_name = collision.second.name(); + + if (/// Charge exchange + (collisionSpeciesMatch( + collision_name, name, "", "cx", "partial")) or + /// Any collision (en, in, ee, ii, nn) + (collisionSpeciesMatch( + collision_name, name, "", "coll", "partial"))) { + collision_names.push_back(collision_name); + } + } + + } else { + throw BoutException("\ndiffusion_collisions_mode for {:s} must be either multispecies or braginskii", name); + } + + if (collision_names.empty()) { + throw BoutException("\tNo collisions found for {:s} in neutral_mixed for selected collisions mode", name); + } + + /// Write chosen collisions to log file + output_info.write("\t{:s} neutral collisionality mode: '{:s}' using ", + name, diffusion_collisions_mode); + for (const auto& collision : collision_names) { + output_info.write("{:s} ", collision); + } + output_info.write("\n"); + } + + /// Collect the collisionalities based on list of names + nu = 0; + for (const auto& collision_name : collision_names) { + nu += GET_VALUE(Field3D, localstate["collision_frequencies"][collision_name]); + } + + // Dnn = Vth^2 / sigma - Dnn = (Tn / AA) / (get(localstate["collision_frequency"]) + Rnn); + Dnn = (Tn / AA) / (nu + Rnn); } else { Dnn = (Tn / AA) / Rnn; } diff --git a/src/neutral_parallel_diffusion.cxx b/src/neutral_parallel_diffusion.cxx index 744a2653d..74a62e75e 100644 --- a/src/neutral_parallel_diffusion.cxx +++ b/src/neutral_parallel_diffusion.cxx @@ -1,4 +1,5 @@ #include "../include/neutral_parallel_diffusion.hxx" +#include "../include/hermes_utils.hxx" #include #include @@ -19,13 +20,74 @@ void NeutralParallelDiffusion::transform(Options& state) { continue; } - auto nu = GET_VALUE(Field3D, species["collision_frequency"]); const BoutReal AA = GET_VALUE(BoutReal, species["AA"]); // Atomic mass const Field3D Nn = GET_VALUE(Field3D, species["density"]); const Field3D Tn = GET_VALUE(Field3D, species["temperature"]); const Field3D Pn = IS_SET(species["pressure"]) ? GET_VALUE(Field3D, species["pressure"]) : Nn * Tn; + // Collisionality + // Multispecies mode: in, en, nn, cx + // New mode: cx, iz (in line with SOLPS AFN, Horsten 2017) + if (collision_names.empty()) { /// Calculate only once - at the beginning + + if (diffusion_collisions_mode == "afn") { + for (const auto& collision : species["collision_frequencies"].getChildren()) { + + std::string collision_name = collision.second.name(); + + if (/// Charge exchange + (collisionSpeciesMatch( + collision_name, species.name(), "+", "cx", "partial")) or + /// Ionisation + (collisionSpeciesMatch( + collision_name, species.name(), "+", "iz", "partial"))) { + + collision_names.push_back(collision_name); + } + } + // Multispecies mode: all collisions and CX are included + } else if (diffusion_collisions_mode == "multispecies") { + for (const auto& collision : species["collision_frequencies"].getChildren()) { + + std::string collision_name = collision.second.name(); + + if (/// Charge exchange + (collisionSpeciesMatch( + collision_name, species.name(), "", "cx", "partial")) or + /// Any collision (ne, ni, nn) + (collisionSpeciesMatch( + collision_name, species.name(), "", "coll", "partial"))) { + + collision_names.push_back(collision_name); + } + } + + } else { + throw BoutException("\tdiffusion_collisions_mode for {:s} must be either multispecies or afn", species.name()); + } + + if (collision_names.empty()) { + throw BoutException("\tNo collisions found for {:s} in neutral_parallel_diffusion for selected collisions mode", species.name()); + } + + /// Write chosen collisions to log file + output_info.write("\t{:s} neutral diffusion collisionality mode: '{:s}' using ", + species.name(), diffusion_collisions_mode); + for (const auto& collision : collision_names) { + output_info.write("{:s} ", collision); + } + output_info.write("\n"); + + } + + /// Collect the collisionalities based on list of names + nu = 0; + for (const auto& collision_name : collision_names) { + nu += GET_VALUE(Field3D, species["collision_frequencies"][collision_name]); + } + + // Diffusion coefficient BoutReal advection_factor = 0; BoutReal kappa_factor = 0; @@ -57,13 +119,13 @@ void NeutralParallelDiffusion::transform(Options& state) { // Heat transfer Field3D E = + FV::Div_par_K_Grad_par( Dn * advection_factor * Pn, logPn); // Pressure advection - if (thermal_conduction) { + if (perpendicular_conduction) { E += FV::Div_par_K_Grad_par(kappa_n, Tn); // Conduction } add(species["energy_source"], E); Field3D F = 0.0; - if (IS_SET(species["velocity"]) and viscosity) { + if (IS_SET(species["velocity"]) and perpendicular_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 diff --git a/src/parallel_inertia_curvature_drift.cxx b/src/parallel_inertia_curvature_drift.cxx new file mode 100644 index 000000000..0c3ee3d3e --- /dev/null +++ b/src/parallel_inertia_curvature_drift.cxx @@ -0,0 +1,141 @@ +#include +#include +#include + +#include "../include/parallel_inertia_curvature_drift.hxx" + +using bout::globals::mesh; + +ParallelInertiaCurvatureDrift::ParallelInertiaCurvatureDrift(std::string name, Options& alloptions, + Solver* UNUSED(solver)) { + + // Get options for this component + auto& options = alloptions[name]; + + bndry_flux = + options["bndry_flux"].doc("Allow fluxes through boundary?").withDefault(true); + + diagnose = + options["diagnose"].doc("Output additional diagnostics?").withDefault(false); + + // Read curvature vector + Curlb_B.covariant = false; // Contravariant + if (mesh->get(Curlb_B, "bxcv")) { + Curlb_B.x = Curlb_B.y = Curlb_B.z = 0.0; + } + + Options& paralleltransform = Options::root()["mesh"]["paralleltransform"]; + if (paralleltransform.isSet("type") and + paralleltransform["type"].as() == "shifted") { + Field2D I; + if (mesh->get(I, "sinty")) { + I = 0.0; + } + Curlb_B.z += I * Curlb_B.x; + } + + // Normalise + + // Get the units + const auto& units = alloptions["units"]; + BoutReal Bnorm = get(units["Tesla"]); + BoutReal Lnorm = get(units["meters"]); + + Curlb_B.x /= Bnorm; + Curlb_B.y *= SQ(Lnorm); + Curlb_B.z *= SQ(Lnorm); + + Curlb_B *= 2. / mesh->getCoordinates()->Bxy; + + // Set drift to zero through sheath boundaries. + // Flux through those cell faces should be set by sheath. + for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { + Curlb_B.y(r.ind, mesh->ystart - 1) = -Curlb_B.y(r.ind, mesh->ystart); + } + for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { + Curlb_B.y(r.ind, mesh->yend + 1) = -Curlb_B.y(r.ind, mesh->yend); + } + + // Initialise current field + DivJicurv = 0.0; +} + +void ParallelInertiaCurvatureDrift::transform(Options& state) { + // Iterate through all subsections + Options& allspecies = state["species"]; + + for (auto& kv : allspecies.getChildren()) { + Options& species = allspecies[kv.first]; // Note: Need non-const + + //NOTE: Should we skip electrons? + // if (kv.first == "e") { + // continue; // Skip electrons -> only ions + // } + + if (!(species.isSet("charge") and species.isSet("velocity") and species.isSet("AA"))) + continue; // Skip, go to next species + + // Calculate parallel inertia curvature drift velocity for this species + auto q = get(species["charge"]); + if (fabs(q) < 1e-5) { + continue; + } + + auto Vpar = GET_VALUE(Field3D, species["velocity"]); // Parallel velocity + const auto AA = get(species["AA"]); + + // Parallel inertia curvature drift velocity + Vector3D vI = 0.5 * AA / q * SQ(Vpar) * Curlb_B; + + if (IS_SET(species["density"])) { + auto N = GET_VALUE(Field3D, species["density"]); + + Field3D density_source = FV::Div_f_v(N, vI, bndry_flux); + + subtract(species["density_source"], density_source); + + + Vector3D Jicurv = N * vI * q; // Parallel inertia curvature current for this species + + // Divergence of current in vorticity equation + DivJicurv += Div(Jicurv); + + } + + if (IS_SET(species["pressure"])) { + auto P = get(species["pressure"]); + + Field3D energy_source = (5. / 2) * FV::Div_f_v(P, vI, bndry_flux); + subtract(species["energy_source"], energy_source ); + } + + if (IS_SET(species["momentum"])) { + auto NV = get(species["momentum"]); + Field3D momentum_source = FV::Div_f_v(NV, vI, bndry_flux); + subtract(species["momentum_source"], momentum_source); + } + + } + + add(state["fields"]["DivJextra"], DivJicurv); + +} + +void ParallelInertiaCurvatureDrift::outputVars(Options& state) { + AUTO_TRACE(); + + if (diagnose) { + + // Normalisations + auto Nnorm = get(state["Nnorm"]); + auto Omega_ci = get(state["Omega_ci"]); + + set_with_attrs(state["DivJicurv"], DivJicurv, + {{"time_dimension", "t"}, + {"units", "A m^-3"}, + {"conversion", SI::qe * Nnorm * Omega_ci}, + {"long_name", "Divergence of parallel inertia curvature current"}, + {"source", "parallel_inertia_curvature_drift"}}); + + } +} \ No newline at end of file diff --git a/src/sheath_boundary_simple.cxx b/src/sheath_boundary_simple.cxx index d46e40fb2..ce4fc5a24 100644 --- a/src/sheath_boundary_simple.cxx +++ b/src/sheath_boundary_simple.cxx @@ -138,6 +138,9 @@ SheathBoundarySimple::SheathBoundarySimple(std::string name, Options& alloptions .doc("BC mode: 0=LimitFree, 1=ExponentialFree, 2=LinearFree") .withDefault(1); + diagnose = options["diagnose"] + .doc("Save additional output diagnostics") + .withDefault(false); } @@ -184,7 +187,7 @@ void SheathBoundarySimple::transform(Options& state) { // // To avoid looking up species for every grid point, this // loops over the boundaries once per species. - Field3D ion_sum = 0.0; + ion_sum = 0.0; // Iterate through charged ion species for (auto& kv : allspecies.getChildren()) { @@ -336,6 +339,8 @@ void SheathBoundarySimple::transform(Options& state) { ? toFieldAligned(getNonFinal(electrons["energy_source"])) : zeroFrom(Ne); + hflux_e = zeroFrom(electron_energy_source); // sheath heat flux for diagnostics + if (lower_y) { for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { for (int jz = 0; jz < mesh->LocalNz; jz++) { @@ -362,6 +367,7 @@ void SheathBoundarySimple::transform(Options& state) { floor(0.5 * (phi[im] + phi[i]), phi_wall); // Electron saturation at phi = phi_wall // Electron velocity into sheath (< 0) + // Equal to Bohm for single ions and no currents BoutReal vesheath = -sqrt(tesheath / (TWOPI * Me)) * (1. - Ge) * exp(-(phisheath - phi_wall) / floor(tesheath, 1e-5)); @@ -379,16 +385,21 @@ void SheathBoundarySimple::transform(Options& state) { // This is additional energy flux through the sheath q -= (2.5 * tesheath + 0.5 * Me * SQ(vesheath)) * nesheath * vesheath; - // Multiply by cell area to get power - BoutReal heatflow = q * (coord->J[i] + coord->J[im]) - / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[im])); // This omits dx*dz because we divide by dx*dz next - - // Divide by volume of cell to get energy loss rate (< 0) - BoutReal power = heatflow / (coord->dy[i] * coord->J[i]); + // Cross-sectional area in XZ plane and cell volume + BoutReal da = (coord->J[i] + coord->J[im]) / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[im])) + * 0.5*(coord->dx[i] + coord->dx[im]) * 0.5*(coord->dz[i] + coord->dz[im]); // [m^2] + BoutReal dv = (coord->dx[i] * coord->dy[i] * coord->dz[i] * coord->J[i]); // [m^3] + // Get power and energy source + BoutReal heatflow = q * da; // [W] + BoutReal power = heatflow / dv; // [Wm^-3] electron_energy_source[i] += power; - electron_sheath_power_ylow[i] += heatflow * coord->dx[i] * coord->dz[i]; // lower Y, so power placed in final domain cell + // Total heat flux for diagnostic purposes + q = gamma_e * tesheath * nesheath * vesheath; // [Wm^-2] + hflux_e[i] += q * da / dv; // [Wm^-3] + electron_sheath_power_ylow[i] += heatflow; // [W], lower Y, so sheath boundary power placed in final domain cell + } } } @@ -421,8 +432,9 @@ void SheathBoundarySimple::transform(Options& state) { // Electron velocity into sheath (> 0) BoutReal vesheath = - sqrt(tesheath / (TWOPI * Me)) * (1. - Ge) * exp(-(phisheath - phi_wall) / floor(tesheath, 1e-5)); + sqrt(tesheath / (TWOPI * Me)) * (1. - Ge) * exp(-(phisheath - phi_wall) / floor(tesheath, 1e-5)); + // Heat flux. Note: Here this is positive because vesheath > 0 BoutReal q = gamma_e * tesheath * nesheath * vesheath; if (no_flow) { @@ -431,26 +443,32 @@ void SheathBoundarySimple::transform(Options& state) { Ve[ip] = 2 * vesheath - Ve[i]; NVe[ip] = 2. * Me * nesheath * vesheath - NVe[i]; + // Take into account the flow of energy due to fluid flow // This is additional energy flux through the sheath - // Note: Here this is positive because vesheath > 0 q -= (2.5 * tesheath + 0.5 * Me * SQ(vesheath)) * nesheath * vesheath; - // Multiply by cell area to get power - BoutReal heatflow = q * (coord->J[i] + coord->J[ip]) - / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[ip])); // This omits dx*dz because we divide by dx*dz next - - // Divide by volume of cell to get energy loss rate (> 0) - BoutReal power = heatflow / (coord->dy[i] * coord->J[i]); + // Cross-sectional area in XZ plane and cell volume + BoutReal da = (coord->J[i] + coord->J[ip]) / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[ip])) + * 0.5*(coord->dx[i] + coord->dx[ip]) * 0.5*(coord->dz[i] + coord->dz[ip]); // [m^2] + BoutReal dv = (coord->dx[i] * coord->dy[i] * coord->dz[i] * coord->J[i]); // [m^3] + // Get power and energy source + BoutReal heatflow = q * da; // [W] + BoutReal power = heatflow / dv; // [Wm^-3] electron_energy_source[i] -= power; - // Diagnostic contains energy removed in the sheath - electron_sheath_power_ylow[ip] += heatflow * coord->dx[i] * coord->dz[i]; // upper Y, so power placed in first guard cell + // Total heat flux for diagnostic purposes + q = gamma_e * tesheath * nesheath * vesheath; // [Wm^-2] + hflux_e[i] -= q * da / dv; // [Wm^-3] + electron_sheath_power_ylow[ip] += heatflow; // [W] Upper Y, so sheath boundary power on ylow side of inner guard cell + } } } + set(diagnostics["e"]["energy_source"], hflux_e); + // Set electron density and temperature, now with boundary conditions // Note: Clear parallel slices because they do not contain boundary conditions. Ne.clearParallelSlices(); @@ -482,6 +500,8 @@ void SheathBoundarySimple::transform(Options& state) { setBoundary(state["fields"]["phi"], fromFieldAligned(phi)); } + + ////////////////////////////////////////////////////////////////// // Iterate through all ions for (auto& kv : allspecies.getChildren()) { @@ -523,6 +543,10 @@ void SheathBoundarySimple::transform(Options& state) { ? toFieldAligned(getNonFinal(species["energy_source"])) : zeroFrom(Ni); + // Initialise sheath ion heat flux. This will be created for each species + // saved in diagnostics struct and then destroyed and re-created for next species + Field3D hflux_i = zeroFrom(energy_source); + Field3D particle_source = zeroFrom(energy_source); // Field to capture total sheath heat flux for diagnostics Field3D ion_sheath_power_ylow = zeroFrom(Ne); @@ -573,16 +597,22 @@ void SheathBoundarySimple::transform(Options& state) { // This is additional energy flux through the sheath q -= (2.5 * tisheath + 0.5 * Mi * SQ(visheath)) * nisheath * visheath; - // Multiply by cell area to get power - BoutReal heatflow = q * (coord->J[i] + coord->J[im]) - / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[im])); // This omits dx*dz because we divide by dx*dz next - - // Divide by volume of cell to get energy loss rate (< 0) - BoutReal power = heatflow / (coord->dy[i] * coord->J[i]); + // Cross-sectional area in XZ plane and cell volume + BoutReal da = (coord->J[i] + coord->J[ip]) / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[ip])) + * 0.5*(coord->dx[i] + coord->dx[ip]) * 0.5*(coord->dz[i] + coord->dz[ip]); // [m^2] + BoutReal dv = (coord->dx[i] * coord->dy[i] * coord->dz[i] * coord->J[i]); // [m^3] - energy_source[i] += power; + // Get power and energy source + BoutReal heatflow = q * da; // [W] + BoutReal power = heatflow / dv; // [Wm^-3] + ASSERT2(std::isfinite(power)); + energy_source[i] += power; // Note: Sign negative because power > 0 + particle_source[i] -= nisheath * visheath * da / dv; // [m^-3s^-1] Diagnostics only - ion_sheath_power_ylow[i] += heatflow * coord->dx[i] * coord->dz[i]; // lower Y, so power placed in final domain cell + // Total heat flux for diagnostic purposes + q = gamma_i * tisheath * nisheath * visheath; // [Wm^-2] + hflux_i[i] += q * da / dv; // [Wm^-3] + ion_sheath_power_ylow[i] += heatflow; // [W] lower Y, so power placed in final domain cell } } } @@ -636,20 +666,28 @@ void SheathBoundarySimple::transform(Options& state) { // Note: Here this is positive because visheath > 0 q -= (2.5 * tisheath + 0.5 * Mi * SQ(visheath)) * nisheath * visheath; - // Multiply by cell area to get power - BoutReal heatflow = q * (coord->J[i] + coord->J[ip]) - / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[ip])); // This omits dx*dz because we divide by dx*dz next + // Cross-sectional area in XZ plane and cell volume + BoutReal da = (coord->J[i] + coord->J[ip]) / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[ip])) + * 0.5*(coord->dx[i] + coord->dx[ip]) * 0.5*(coord->dz[i] + coord->dz[ip]); // [m^2] + BoutReal dv = (coord->dx[i] * coord->dy[i] * coord->dz[i] * coord->J[i]); // [m^3] - // Divide by volume of cell to get energy loss rate (> 0) - BoutReal power = heatflow / (coord->dy[i] * coord->J[i]); + // Get power and energy source + BoutReal heatflow = q * da; // [W] + BoutReal power = heatflow / dv; // [Wm^-3] ASSERT2(std::isfinite(power)); - energy_source[i] -= power; // Note: Sign negative because power > 0 + particle_source[i] -= nisheath * visheath * da / dv; // [m^-3s^-1] Diagnostics only + + // Total heat flux for diagnostic purposes + q = gamma_i * tisheath * nisheath * visheath; // [Wm^-2] + hflux_i[i] -= q * da / dv; // [Wm^-3] + ion_sheath_power_ylow[ip] += heatflow; // [W] Upper Y, so sheath boundary power on ylow side of inner guard cell - ion_sheath_power_ylow[ip] += heatflow * coord->dx[i] * coord->dz[i]; // Upper Y, so power placed in first guard cell } } } + + // Finished boundary conditions for this species // Put the modified fields back into the state. @@ -676,5 +714,48 @@ void SheathBoundarySimple::transform(Options& state) { // Add the total sheath power flux to the tracker of y power flows add(species["energy_flow_ylow"], fromFieldAligned(ion_sheath_power_ylow)); + + set(diagnostics[species.name()]["energy_source"], hflux_i); + set(diagnostics[species.name()]["particle_source"], particle_source); + } } + +void SheathBoundarySimple::outputVars(Options& state) { + AUTO_TRACE(); + // Normalisations + auto Nnorm = get(state["Nnorm"]); + auto rho_s0 = get(state["rho_s0"]); + auto Omega_ci = get(state["Omega_ci"]); + auto Tnorm = get(state["Tnorm"]); + BoutReal Pnorm = SI::qe * Tnorm * Nnorm; // Pressure normalisation + + if (diagnose) { + /// Iterate through the first species in each collision pair + const std::map& level1 = diagnostics.getChildren(); + for (auto s1 = std::begin(level1); s1 != std::end(level1); ++s1) { + auto species_name = s1->first; + const Options& section = diagnostics[species_name]; + + set_with_attrs(state[{"E" + species_name + "_sheath"}], getNonFinal(section["energy_source"]), + {{"time_dimension", "t"}, + {"units", "W / m^3"}, + {"conversion", Pnorm * Omega_ci}, + {"standard_name", "energy source"}, + {"long_name", species_name + " sheath energy source"}, + {"source", "sheath_boundary_simple"}}); + + if (species_name != "e") { + set_with_attrs(state[{"S" + species_name + "_sheath"}], getNonFinal(section["particle_source"]), + {{"time_dimension", "t"}, + {"units", "m^-3 s^-1"}, + {"conversion", Nnorm * Omega_ci}, + {"standard_name", "energy source"}, + {"long_name", species_name + " sheath energy source"}, + {"source", "sheath_boundary_simple"}}); + } + + } + } + + }