diff --git a/CMakeLists.txt b/CMakeLists.txt index c895bde2..12000d20 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -43,7 +43,8 @@ message(STATUS "Git revision: ${HERMES_REVISION}") option(HERMES_BUILD_BOUT "Build BOUT++ in external/BOUT-dev" ON) if(HERMES_BUILD_BOUT) set(BOUT_BUILD_EXAMPLES, OFF) # Don't create example makefiles - add_subdirectory(external/BOUT-dev) + set(HERMES_BOUT_SRC "external/BOUT-dev" CACHE STRING "Directory of BOUT++ dir") + add_subdirectory(${HERMES_BOUT_SRC} ${CMAKE_CURRENT_BINARY_DIR}/external/BOUT-dev) else() find_package(bout++ REQUIRED) endif() @@ -224,6 +225,7 @@ if(HERMES_TESTS) hermes_add_integrated_test(sod-shock) hermes_add_integrated_test(sod-shock-energy) hermes_add_integrated_test(drift-wave) + hermes_add_integrated_test(alfven-wave) # Unit tests option(HERMES_UNIT_TESTS "Build the unit tests" ON) @@ -264,5 +266,8 @@ endif() configure_file(include/hermes_build_config.hxx.in include/hermes_build_config.hxx) # Installation - -install(TARGETS hermes-3) +if (BUILD_SHARED_LIBS) + install(TARGETS hermes-3 hermes-3-lib) +else() + install(TARGETS hermes-3) +endif() diff --git a/docs/sphinx/code_structure.rst b/docs/sphinx/code_structure.rst index a0d49347..09dd8a21 100644 --- a/docs/sphinx/code_structure.rst +++ b/docs/sphinx/code_structure.rst @@ -54,7 +54,7 @@ for new variables which are added to the state: * `charge` Charge, in units of proton charge (i.e. electron=-1) * `density` - * `momentum` + * `momentum` Parallel momentum * `pressure` * `velocity` Parallel velocity * `temperature` @@ -74,12 +74,14 @@ for new variables which are added to the state: * `fields` * `vorticity` - * `phi` Electrostatic potential - * `DivJdia` Divergence of diamagnetic current - * `DivJcol` Divergence of collisional current - * `DivJextra` Divergence of current, including 2D parallel current - closures. Not including diamagnetic, parallel current due to - flows, or polarisation currents + * `phi` Electrostatic potential + * `Apar` Electromagnetic potential b dot A in induction terms + * `Apar_flutter` The electromagnetic potential (b dot A) in flutter terms + * `DivJdia` Divergence of diamagnetic current + * `DivJcol` Divergence of collisional current + * `DivJextra` Divergence of current, including 2D parallel current + closures. Not including diamagnetic, parallel current due to + flows, or polarisation currents For example to get the electron density:: diff --git a/docs/sphinx/components.rst b/docs/sphinx/components.rst index 7051575d..d9046a50 100644 --- a/docs/sphinx/components.rst +++ b/docs/sphinx/components.rst @@ -969,7 +969,8 @@ listed after all the species groups in the component list, so that all the species are present in the state. One of the most important is the `collisions`_ component. This sets collision -times for all species, which are then used +times for all species, which are then used in other components to calculate +quantities like heat diffusivities and viscosity closures. .. _sound_speed: @@ -1176,12 +1177,18 @@ species :math:`b` due to temperature differences, is given by: - Ion-neutral and electron-neutral collisions + *Note*: These are disabled by default. If enabled, care is needed to + avoid double-counting collisions in atomic reactions e.g charge-exchange + reactions. + The cross-section for elastic collisions between charged and neutral particles can vary significantly. Here for simplicity we just take a value of :math:`5\times 10^{-19}m^2` from the NRL formulary. - Neutral-neutral collisions + *Note* This is enabled by default. + The cross-section is given by .. math:: @@ -1916,6 +1923,20 @@ This functionality is not yet currently implemented for helium or neon reactions | R_multiplier | Impurity species | Fixed frac. impurity radiation rate | +-----------------------+------------------+---------------------------------------+ +The charge exchange reaction can also be modified so that the momentum transfer channel is disabled. This can be useful when +testing the impact of the full neutral momentum equation equation compared to purely diffusive neutrals. A diffusive only model +leads to all of the ion momentum being lost during charge exchange due to the lack of a neutral momentum equation. +Enabling neutral momentum introduces a more accurate transport model but also prevents CX momentum from being lost, which +can have a significant impact on the solution and may be difficult to analyse. +Disabling the momentum transfer channel allows you to study the impact of the improved transport only and is set as: + +.. code-block:: ini + + [hermes] + components = ..., c, ... + + [reactions] + no_neutral_cx_mom_gain = true Electromagnetic fields ---------------------- @@ -2017,9 +2038,26 @@ the potential in time as a diffusion equation. .. doxygenstruct:: RelaxPotential :members: +.. _electromagnetic: + electromagnetic ~~~~~~~~~~~~~~~ +**Notes**: When using this module, + +1. Set ``sound_speed:alfven_wave=true`` so that the shear Alfven wave + speed is included in the calculation of the fastest parallel wave + speed for numerical dissipation. +2. For tokamak simulations use Neumann boundary condition on the core + and Dirichlet on SOL and PF boundaries by setting + ``electromagnetic:apar_core_neumann=true`` (this is the default). +3. Set the potential core boundary to be constant in Y by setting + ``vorticity:phi_core_averagey = true`` +4. Magnetic flutter terms must be enabled to be active + (``electromagnetic:magnetic_flutter=true``). They use an + ``Apar_flutter`` field, not the ``Apar`` field that is calculated + from the induction terms. + This component modifies the definition of momentum of all species, to include the contribution from the electromagnetic potential :math:`A_{||}`. @@ -2032,8 +2070,17 @@ Assumes that "momentum" :math:`p_s` calculated for all species p_s = m_s n_s v_{||s} + Z_s e n_s A_{||} which arises once the electromagnetic contribution to the force on -each species is included in the momentum equation. This is normalised -so that in dimensionless quantities +each species is included in the momentum equation. This requires +an additional term in the momentum equation: + +.. math:: + + \frac{\partial p_s}{\partial t} = \cdots + Z_s e A_{||} \frac{\partial n_s}{\partial t} + +This is implemented so that the density time-derivative is calculated using the lowest order +terms (parallel flow, ExB drift and a low density numerical diffusion term). + +The above equations are normalised so that in dimensionless quantities: .. math:: @@ -2068,5 +2115,22 @@ To convert the species momenta into a current, we take the sum of - \frac{1}{\beta_{em}} \nabla_\perp^2 A_{||} + \sum_s \frac{Z^2 n_s}{A}A_{||} = \sum_s \frac{Z}{A} p_s +The toroidal variation of density :math:`n_s` must be kept in this +equation. By default the iterative "Naulin" solver is used to do +this: A fast FFT-based method is used in a fixed point iteration, +correcting for the density variation. + +Magnetic flutter terms are disabled by default, and can be enabled by setting + +.. code-block:: ini + + [electromagnetic] + magnetic_flutter = true + +This writes an ``Apar_flutter`` field to the state, which then enables perturbed +parallel derivative terms in the ``evolve_density``, ``evolve_pressure``, ``evolve_energy`` and +``evolve_momentum`` components. Parallel flow terms are modified, and parallel heat +conduction. + .. doxygenstruct:: Electromagnetic :members: diff --git a/docs/sphinx/figs/alfven-wave.png b/docs/sphinx/figs/alfven-wave.png new file mode 120000 index 00000000..3b8c8e8a --- /dev/null +++ b/docs/sphinx/figs/alfven-wave.png @@ -0,0 +1 @@ +../../../tests/integrated/alfven-wave/alfven-wave.png \ No newline at end of file diff --git a/docs/sphinx/figs/drift-wave.png b/docs/sphinx/figs/drift-wave.png new file mode 120000 index 00000000..8a4762d7 --- /dev/null +++ b/docs/sphinx/figs/drift-wave.png @@ -0,0 +1 @@ +../../../tests/integrated/drift-wave/drift-wave.png \ No newline at end of file diff --git a/docs/sphinx/tests.rst b/docs/sphinx/tests.rst index 4c8c731e..cd080576 100644 --- a/docs/sphinx/tests.rst +++ b/docs/sphinx/tests.rst @@ -203,3 +203,54 @@ where This is a cubic dispersion relation, so we find the three roots (using NumPy), and choose the root with the most positive growth rate (imaginary component of :math:`\omega`). + +.. figure:: figs/drift-wave.png + :name: drift-wave + :alt: Comparison of drift-wave growth rate (top) and frequency (bottom) + :width: 60% + +Alfven wave +----------- + +The equations solved are + +.. math:: + + \begin{aligned} + \frac{\partial}{\partial t}\nabla\cdot\left(\frac{n_0 m_i}{B^2}\nabla_\perp\phi\right) =& \nabla_{||}J_{||} = -\nabla_{||}\left(en_ev_{||e}\right) \\ + \frac{\partial}{\partial t}\left(m_en_ev_{||e} - en_eA_{||}\right) =& -\nabla\cdot\left(m_en_ev_{||e} \mathbf{b}v_{||e}\right) + en_e\partial_{||}\phi - 0.51\nu_{ei}n_im_ev_{||e} \\ + J_{||} =& \frac{1}{\mu_0}\nabla_\perp^2 A_{||} + \end{aligned} + +Linearising around a stationary background with constant density +:math:`n_0` and temperature :math:`T_0`, using +:math:`\frac{\partial}{\partial t}\rightarrow -i\omega` gives: + +.. math:: + + \begin{aligned} + \tilde{\phi} =& -\frac{k_{||}}{\omega k_\perp^2}\frac{eB^2}{m_i}\tilde{v_{||e}} \\ + \omega \left( m_e \tilde{v_{||e}} - e\tilde{A}_{||}\right) =& -ek_{||}\tilde{\phi} - i0.51\nu_{ei}m_e\tilde{v_{||e}} \\ + en_0\tilde{v_{||e}} =& -\frac{k_\perp^2}{\mu_0}\tilde{A}_{||} + \end{aligned} + +Rearranging results in a quadratic dispersion relation: + +.. math:: + + \omega^2\left(1 + \frac{k_\perp^2 c^2}{\omega_{pe}^2}\right) + i 0.51\nu_{ei}\frac{k_\perp^2 c^2}{\omega_{pe}^2}\omega - k_{||}^2V_A^2 = 0 + +where :math:`V_A = B / \sqrt{\mu_0 n_0 m_i}` is the Alfven speed, and +:math:`c / \omega_{pe} = \sqrt{m_e / \left(\mu_0 n_0 e^2\right)}` is +the electron skin depth. + +When collisions are neglected, we obtain the result + +.. math:: + + \omega^2 = V_A^2\frac{k_{||}^2}{1 + k_\perp^2 c^2 / \omega_{pe}^2} + +.. figure:: figs/alfven-wave.png + :name: alfven-wave + :alt: Alfven wave speed, as function of parallel and perpendicular wavenumbers + :width: 60% diff --git a/include/amjuel_reaction.hxx b/include/amjuel_reaction.hxx index a8a9f96e..92eb3724 100644 --- a/include/amjuel_reaction.hxx +++ b/include/amjuel_reaction.hxx @@ -164,7 +164,7 @@ protected: Ne.getRegion("RGN_NOBNDRY"))(Ne, N1, Te); // Loss is reduced by heating - energy_loss -= (electron_heating / Tnorm) * reaction_rate; + energy_loss -= (electron_heating / Tnorm) * reaction_rate * radiation_multiplier; subtract(electron["energy_source"], energy_loss); } diff --git a/include/binormal_stpm.hxx b/include/binormal_stpm.hxx index e68659be..8652c4fd 100644 --- a/include/binormal_stpm.hxx +++ b/include/binormal_stpm.hxx @@ -31,12 +31,15 @@ struct BinormalSTPM : public Component { /// - density correction /// void transform(Options& state) override; + void outputVars(Options &state) override; private: std::string name; ///< Short name of the species e.g. h+ + bool diagnose; ///< Output diagnostics? Field3D Theta, chi, D, nu; ///< Field line pitch, anomalous thermal, momentum diffusion Field3D nu_Theta, chi_Theta, D_Theta; ///< nu/Theta, chi/Theta, D/Theta, precalculated + Field3D Theta_inv; ///< Precalculate 1/Theta }; diff --git a/include/component.hxx b/include/component.hxx index 68726df4..61969da2 100644 --- a/include/component.hxx +++ b/include/component.hxx @@ -319,4 +319,15 @@ void set_with_attrs(Options& option, T value, std::initializer_list= 1 +template<> +inline void set_with_attrs(Options& option, Field3D value, std::initializer_list> attrs) { + if (!value.isAllocated()) { + throw BoutException("set_with_attrs: Field3D assigned to {:s} is not allocated", option.str()); + } + option.force(value); + option.setAttributes(attrs); +} +#endif + #endif // HERMES_COMPONENT_H diff --git a/include/div_ops.hxx b/include/div_ops.hxx index 09bcccee..a04fff93 100644 --- a/include/div_ops.hxx +++ b/include/div_ops.hxx @@ -23,8 +23,8 @@ */ -#ifndef HERMES_DIV_OPS_H -#define HERMES_DIV_OPS_H +#ifndef DIV_OPS_H +#define DIV_OPS_H #include #include @@ -44,6 +44,11 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D& n, const Field3D& f, bool bndry_flux = true, bool poloidal = false, bool positive = false); +/// This version has an extra coefficient 'g' that is linearly interpolated +/// onto cell faces +const Field3D Div_n_g_bxGrad_f_B_XZ(const Field3D &n, const Field3D &g, const Field3D &f, + bool bndry_flux = true, bool positive = false); + const Field3D Div_Perp_Lap_FV_Index(const Field3D& a, const Field3D& f, bool xflux); const Field3D Div_Z_FV_Index(const Field3D& a, const Field3D& f); @@ -210,9 +215,10 @@ const Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, flux = bndryval * vpar * vpar; } else { // Add flux due to difference in boundary values - flux = s.R * vpar * sv.R - + BOUTMAX(wave_speed(i, j, k), fabs(v(i, j, k)), fabs(v(i, j + 1, k))) - * (s.R * sv.R - bndryval * vpar); + flux = + s.R * sv.R * sv.R + + BOUTMAX(wave_speed(i, j, k), fabs(v(i, j, k)), fabs(v(i, j + 1, k))) + * (s.R * sv.R - bndryval * vpar); } } else { // Maximum wave speed in the two cells @@ -238,9 +244,10 @@ const Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, flux = bndryval * vpar * vpar; } else { // Add flux due to difference in boundary values - flux = s.L * vpar * sv.L - - BOUTMAX(wave_speed(i, j, k), fabs(v(i, j, k)), fabs(v(i, j - 1, k))) - * (s.L * sv.L - bndryval * vpar); + flux = + s.L * sv.L * sv.L + - BOUTMAX(wave_speed(i, j, k), fabs(v(i, j, k)), fabs(v(i, j - 1, k))) + * (s.L * sv.L - bndryval * vpar); } } else { // Maximum wave speed in the two cells @@ -683,4 +690,4 @@ const Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, const Fi } // namespace FV -#endif // HERMES_DIV_OPS_H +#endif // DIV_OPS_H diff --git a/include/electromagnetic.hxx b/include/electromagnetic.hxx index eb43b70e..15a3588e 100644 --- a/include/electromagnetic.hxx +++ b/include/electromagnetic.hxx @@ -64,6 +64,13 @@ private: std::unique_ptr aparSolver; // Laplacian solver in X-Z + bool const_gradient; // Set neumann boundaries by extrapolation + BoutReal apar_boundary_timescale; // Relaxation timescale + BoutReal last_time; // The last time the boundaries were updated + + bool magnetic_flutter; ///< Set the magnetic flutter term? + Field3D Apar_flutter; + bool diagnose; ///< Output additional diagnostics? }; diff --git a/include/evolve_momentum.hxx b/include/evolve_momentum.hxx index 00a90314..c13f9157 100644 --- a/include/evolve_momentum.hxx +++ b/include/evolve_momentum.hxx @@ -35,7 +35,8 @@ private: std::string name; ///< Short name of species e.g "e" Field3D NV; ///< Species parallel momentum (normalised, evolving) - Field3D NV_solver; ///< Momentum as input from solver + Field3D NV_err; ///< Difference from momentum as input from solver + Field3D NV_solver; ///< Momentum as calculated in the solver Field3D V; ///< Species parallel velocity Field3D momentum_source; ///< From other components. Stored for diagnostic output diff --git a/include/fixed_fraction_radiation.hxx b/include/fixed_fraction_radiation.hxx index d1b3b563..127d0404 100644 --- a/include/fixed_fraction_radiation.hxx +++ b/include/fixed_fraction_radiation.hxx @@ -302,7 +302,7 @@ struct FixedFractionRadiation : public Component { AUTO_TRACE(); if (diagnose) { - set_with_attrs(state[std::string("R") + name], radiation, + set_with_attrs(state[std::string("R") + name], -radiation, {{"time_dimension", "t"}, {"units", "W / m^3"}, {"conversion", SI::qe * Tnorm * Nnorm * FreqNorm}, diff --git a/include/hydrogen_charge_exchange.hxx b/include/hydrogen_charge_exchange.hxx index 503eb561..1abe89ba 100644 --- a/include/hydrogen_charge_exchange.hxx +++ b/include/hydrogen_charge_exchange.hxx @@ -72,7 +72,8 @@ protected: Field3D& R, Field3D& atom_mom, Field3D& ion_mom, Field3D& atom_energy, Field3D& ion_energy, Field3D& atom_rate, Field3D& ion_rate, - BoutReal& rate_multiplier); + BoutReal& rate_multiplier, + bool& no_neutral_cx_mom_gain); }; /// Hydrogen charge exchange @@ -119,13 +120,24 @@ struct HydrogenChargeExchangeIsotope : public HydrogenChargeExchange { HydrogenChargeExchangeIsotope(std::string name, Options& alloptions, Solver* solver) : HydrogenChargeExchange(name, alloptions, solver) { + //// Options under [reactions] diagnose = alloptions[name]["diagnose"] .doc("Output additional diagnostics?") .withDefault(false); + // This is useful for testing the impact of enabling the neutral momentum equation. + // When set to true, CX behaves as if using diffusive neutrals but the neutral transport + // still enjoys the full momentum equation treatment. + no_neutral_cx_mom_gain = alloptions[name]["no_neutral_cx_mom_gain"] + .doc("If true, ion momentum in CX is still lost but not given to the neutrals") + .withDefault(false); + + // Options under neutral species of isotope 1 (on LHS of reaction) rate_multiplier = alloptions[{Isotope1}]["K_cx_multiplier"] .doc("Scale the charge exchange rate by this factor") .withDefault(1.0); + + } void transform(Options& state) override { @@ -137,7 +149,8 @@ struct HydrogenChargeExchangeIsotope : public HydrogenChargeExchange { state["species"][{Isotope1, '+'}], // e.g. "h+" R, atom_mom, ion_mom, atom_energy, ion_energy, // Transfer channels atom_rate, ion_rate, // Collision rates in s^-1 - rate_multiplier); // Arbitrary user set multiplier + rate_multiplier, // Arbitrary user set multiplier + no_neutral_cx_mom_gain); // Make CX behave as in diffusive neutrals? if (diagnose) { // Calculate diagnostics to be written to dump file @@ -252,6 +265,7 @@ private: Field3D F, F2; ///< Momentum exchange Field3D E, E2; ///< Energy exchange Field3D atom_rate, ion_rate; ///< Collision rates in s^-1 + bool no_neutral_cx_mom_gain; ///< Make CX behave as in diffusive neutrals? }; namespace { diff --git a/include/neutral_parallel_diffusion.hxx b/include/neutral_parallel_diffusion.hxx index 561cd42c..5366e2a2 100644 --- a/include/neutral_parallel_diffusion.hxx +++ b/include/neutral_parallel_diffusion.hxx @@ -32,6 +32,18 @@ struct NeutralParallelDiffusion : public Component { diagnose = options["diagnose"] .doc("Output additional diagnostics?") .withDefault(false); + + equation_fix = options["equation_fix"] + .doc("Fix correcting pressure advection and conductivity factors?") + .withDefault(true); + + thermal_conduction = options["thermal_conducton"] + .doc("Enable conduction?") + .withDefault(true); + + viscosity = options["viscosity"] + .doc("Enable viscosity?") + .withDefault(true); } /// @@ -60,6 +72,9 @@ private: BoutReal dneut; ///< cross-field diffusion projection (B / Bpol)^2 bool diagnose; ///< Output diagnostics? + bool equation_fix; ///< Fix incorrect 3/2 factor in pressure advection? + bool thermal_conduction; ///< Enable conduction? + bool viscosity; ///< Enable viscosity? /// Per-species diagnostics struct Diagnostics { diff --git a/include/vorticity.hxx b/include/vorticity.hxx index 18fd8b94..acf2c281 100644 --- a/include/vorticity.hxx +++ b/include/vorticity.hxx @@ -37,6 +37,8 @@ struct Vorticity : public Component { /// Relax radial phi boundaries towards zero-gradient? /// - phi_boundary_timescale: float, 1e-4 /// Timescale for phi boundary relaxation [seconds] + /// - phi_core_averagey: bool, default false + /// Average phi core boundary in Y? (if phi_boundary_relax) /// - phi_dissipation: bool, default true /// Parallel dissipation of potential (Recommended) /// - poloidal_flows: bool, default true @@ -128,7 +130,8 @@ private: bool phi_boundary_relax; ///< Relax boundary to zero-gradient BoutReal phi_boundary_timescale; ///< Relaxation timescale [normalised] BoutReal phi_boundary_last_update; ///< Time when last updated - + bool phi_core_averagey; ///< Average phi core boundary in Y? + bool split_n0; // Split phi into n=0 and n!=0 components LaplaceXY* laplacexy; // Laplacian solver in X-Y (n=0) diff --git a/src/binormal_stpm.cxx b/src/binormal_stpm.cxx index f67eadd3..e96d6dba 100644 --- a/src/binormal_stpm.cxx +++ b/src/binormal_stpm.cxx @@ -39,7 +39,11 @@ BinormalSTPM::BinormalSTPM(std::string name, Options& alloptions, Solver* solver chi_Theta = chi/Theta; D_Theta = D/Theta; nu_Theta = nu/Theta; - + Theta_inv = 1./Theta; + + diagnose = options["diagnose"] + .doc("Output diagnostics?") + .withDefault(false); } void BinormalSTPM::transform(Options& state) { @@ -62,15 +66,51 @@ void BinormalSTPM::transform(Options& state) { ? GET_NOBOUNDARY(Field3D, species["momentum"]) : 0.0; - add(species["pressure_source"], - (2. / 3) * (1/Theta) * FV::Div_par_K_Grad_par(chi_Theta*N, T, false)); + add(species["energy_source"], + Theta_inv * FV::Div_par_K_Grad_par(chi_Theta*N, T, false)); add(species["momentum_source"], - (1/Theta) * FV::Div_par_K_Grad_par(AA*nu_Theta, NV, false)); + Theta_inv * FV::Div_par_K_Grad_par(AA*nu_Theta, NV, false)); add(species["density_source"], - (1/Theta) * FV::Div_par_K_Grad_par(D_Theta, N, false)); + Theta_inv * FV::Div_par_K_Grad_par(D_Theta, N, false)); } } +void BinormalSTPM::outputVars(Options& state) { + AUTO_TRACE(); + // Normalisations + auto Omega_ci = get(state["Omega_ci"]); + auto rho_s0 = get(state["rho_s0"]); + + if (diagnose) { + + AUTO_TRACE(); + // Save particle, momentum and energy channels + + set_with_attrs(state[{std::string("D_") + name}], D, + {{"time_dimension", "t"}, + {"units", "m^2 s^-1"}, + {"conversion", rho_s0 * rho_s0 * Omega_ci}, + {"standard_name", "anomalous density diffusion"}, + {"long_name", std::string("Binormal Stellarator 2pt model density diffusion of ") + name}, + {"source", "binormal_stpm"}}); + + set_with_attrs(state[{std::string("chi_") + name}], chi, + {{"time_dimension", "t"}, + {"units", "m^2 s^-1"}, + {"conversion", rho_s0 * rho_s0 * Omega_ci}, + {"standard_name", "anomalous thermal diffusion"}, + {"long_name", std::string("Binormal Stellarator 2pt model thermal diffusion of ") + name}, + {"source", "binormal_stpm"}}); + + set_with_attrs(state[{std::string("nu_") + name}], nu, + {{"time_dimension", "t"}, + {"units", "m^2 s^-1"}, + {"conversion", rho_s0 * rho_s0 * Omega_ci}, + {"standard_name", "anomalous momentum diffusion"}, + {"long_name", std::string("Binormal Stellarator 2pt model momentum diffusion of ") + name}, + {"source", "binormal_stpm"}}); + } +} diff --git a/src/collisions.cxx b/src/collisions.cxx index 92313e8e..d41ba097 100644 --- a/src/collisions.cxx +++ b/src/collisions.cxx @@ -33,13 +33,13 @@ Collisions::Collisions(std::string name, Options& alloptions, Solver*) { .withDefault(true); electron_neutral = options["electron_neutral"] .doc("Include electron-neutral elastic collisions?") - .withDefault(true); + .withDefault(false); ion_ion = options["ion_ion"] .doc("Include ion-ion elastic collisions?") .withDefault(true); ion_neutral = options["ion_neutral"] .doc("Include ion-neutral elastic collisions?") - .withDefault(true); + .withDefault(false); neutral_neutral = options["neutral_neutral"] .doc("Include neutral-neutral elastic collisions?") .withDefault(true); diff --git a/src/div_ops.cxx b/src/div_ops.cxx index af235251..64ff8143 100644 --- a/src/div_ops.cxx +++ b/src/div_ops.cxx @@ -1140,3 +1140,161 @@ const Field3D Div_a_Grad_perp_upwind_flows(const Field3D& a, const Field3D& f, return result; } + +const Field3D Div_n_g_bxGrad_f_B_XZ(const Field3D &n, const Field3D &g, const Field3D &f, + bool bndry_flux, bool positive) { + Field3D result{0.0}; + + Coordinates *coord = mesh->getCoordinates(); + + ////////////////////////////////////////// + // X-Z advection. + // + // Z + // | + // + // fmp --- vU --- fpp + // | nU | + // | | + // vL nL nR vR -> X + // | | + // | nD | + // fmm --- vD --- fpm + // + + int nz = mesh->LocalNz; + for (int i = mesh->xstart; i <= mesh->xend; i++) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + for (int k = 0; k < nz; k++) { + int kp = (k + 1) % nz; + int kpp = (kp + 1) % nz; + int km = (k - 1 + nz) % nz; + int kmm = (km - 1 + nz) % nz; + + // 1) Interpolate stream function f onto corners fmp, fpp, fpm + + BoutReal fmm = 0.25 * (f(i, j, k) + f(i - 1, j, k) + f(i, j, km) + + f(i - 1, j, km)); + BoutReal fmp = 0.25 * (f(i, j, k) + f(i, j, kp) + f(i - 1, j, k) + + f(i - 1, j, kp)); // 2nd order accurate + BoutReal fpp = 0.25 * (f(i, j, k) + f(i, j, kp) + f(i + 1, j, k) + + f(i + 1, j, kp)); + BoutReal fpm = 0.25 * (f(i, j, k) + f(i + 1, j, k) + f(i, j, km) + + f(i + 1, j, km)); + + // 2) Calculate velocities on cell faces + + BoutReal vU = coord->J(i, j) * (fmp - fpp) / coord->dx(i, j); // -J*df/dx + BoutReal vD = coord->J(i, j) * (fmm - fpm) / coord->dx(i, j); // -J*df/dx + + BoutReal vR = 0.5 * (coord->J(i, j) + coord->J(i + 1, j)) * (fpp - fpm) / + coord->dz(i, j); // J*df/dz + BoutReal vL = 0.5 * (coord->J(i, j) + coord->J(i - 1, j)) * (fmp - fmm) / + coord->dz(i, j); // J*df/dz + + // 3) Calculate g on cell faces + + BoutReal gU = 0.5 * (g(i, j, kp) + g(i, j, k)); + BoutReal gD = 0.5 * (g(i, j, km) + g(i, j, k)); + BoutReal gR = 0.5 * (g(i + 1, j, k) + g(i, j, k)); + BoutReal gL = 0.5 * (g(i - 1, j, k) + g(i, j, k)); + + // 4) Calculate n on the cell faces. The sign of the + // velocity determines which side is used. + + // X direction + Stencil1D s; + s.c = n(i, j, k); + s.m = n(i - 1, j, k); + s.mm = n(i - 2, j, k); + s.p = n(i + 1, j, k); + s.pp = n(i + 2, j, k); + + MC(s, coord->dx(i, j)); + + // Right side + if ((i == mesh->xend) && (mesh->lastX())) { + // At right boundary in X + + if (bndry_flux) { + BoutReal flux; + if (vR > 0.0) { + // Flux to boundary + flux = vR * s.R * gR; + } else { + // Flux in from boundary + flux = vR * 0.5 * (n(i + 1, j, k) + n(i, j, k)) * gR; + } + + result(i, j, k) += flux / (coord->dx(i, j) * coord->J(i, j)); + result(i + 1, j, k) -= + flux / (coord->dx(i + 1, j) * coord->J(i + 1, j)); + } + } else { + // Not at a boundary + if (vR > 0.0) { + // Flux out into next cell + BoutReal flux = vR * s.R * gR; + result(i, j, k) += flux / (coord->dx(i, j) * coord->J(i, j)); + result(i + 1, j, k) -= + flux / (coord->dx(i + 1, j) * coord->J(i + 1, j)); + } + } + + // Left side + + if ((i == mesh->xstart) && (mesh->firstX())) { + // At left boundary in X + + if (bndry_flux) { + BoutReal flux; + + if (vL < 0.0) { + // Flux to boundary + flux = vL * s.L * gL; + + } else { + // Flux in from boundary + flux = vL * 0.5 * (n(i - 1, j, k) + n(i, j, k)) * gL; + } + result(i, j, k) -= flux / (coord->dx(i, j) * coord->J(i, j)); + result(i - 1, j, k) += + flux / (coord->dx(i - 1, j) * coord->J(i - 1, j)); + } + } else { + // Not at a boundary + + if (vL < 0.0) { + BoutReal flux = vL * s.L * gL; + result(i, j, k) -= flux / (coord->dx(i, j) * coord->J(i, j)); + result(i - 1, j, k) += + flux / (coord->dx(i - 1, j) * coord->J(i - 1, j)); + } + } + + /// NOTE: Need to communicate fluxes + + // Z direction + s.m = n(i, j, km); + s.mm = n(i, j, kmm); + s.p = n(i, j, kp); + s.pp = n(i, j, kpp); + + MC(s, coord->dz(i, j)); + + if (vU > 0.0) { + BoutReal flux = vU * s.R * gU/ (coord->J(i, j) * coord->dz(i, j)); + result(i, j, k) += flux; + result(i, j, kp) -= flux; + } + if (vD < 0.0) { + BoutReal flux = vD * s.L * gD / (coord->J(i, j) * coord->dz(i, j)); + result(i, j, k) -= flux; + result(i, j, km) += flux; + } + } + } + } + FV::communicateFluxes(result); + return result; +} diff --git a/src/electromagnetic.cxx b/src/electromagnetic.cxx index 971f0c30..7048b47d 100644 --- a/src/electromagnetic.cxx +++ b/src/electromagnetic.cxx @@ -20,19 +20,59 @@ Electromagnetic::Electromagnetic(std::string name, Options &alloptions, Solver*) auto& options = alloptions[name]; + // Use the "Naulin" solver because we need to include toroidal + // variations of the density (A coefficient) + if (!options["laplacian"].isSet("type")) { + options["laplacian"]["type"] = "naulin"; + } aparSolver = Laplacian::create(&options["laplacian"]); - // Set zero-gradient (neumann) boundary conditions - aparSolver->setInnerBoundaryFlags(INVERT_DC_GRAD + INVERT_AC_GRAD); - aparSolver->setOuterBoundaryFlags(INVERT_DC_GRAD + INVERT_AC_GRAD); + + const_gradient = options["const_gradient"] + .doc("Extrapolate gradient of Apar into all radial boundaries?") + .withDefault(false); + + // Give Apar an initial value because we solve Apar by iteration + // starting from the previous solution + Apar = 0.0; + + if (const_gradient) { + // Set flags to take the gradient from the RHS + aparSolver->setInnerBoundaryFlags(INVERT_DC_GRAD + INVERT_AC_GRAD + INVERT_RHS); + aparSolver->setOuterBoundaryFlags(INVERT_DC_GRAD + INVERT_AC_GRAD + INVERT_RHS); + last_time = 0.0; + + apar_boundary_timescale = options["apar_boundary_timescale"] + .doc("Timescale for Apar boundary relaxation [seconds]") + .withDefault(1e-8) + / get(alloptions["units"]["seconds"]); + + } else if (options["apar_boundary_neumann"] + .doc("Neumann on all radial boundaries?") + .withDefault(false)) { + // Set zero-gradient (neumann) boundary condition DC on the core + aparSolver->setInnerBoundaryFlags(INVERT_DC_GRAD + INVERT_AC_GRAD); + aparSolver->setOuterBoundaryFlags(INVERT_DC_GRAD + INVERT_AC_GRAD); + + } else if (options["apar_core_neumann"] + .doc("Neumann radial boundary in the core? False => Dirichlet") + .withDefault(true) + and bout::globals::mesh->periodicY(bout::globals::mesh->xstart)) { + // Set zero-gradient (neumann) boundary condition DC on the core + aparSolver->setInnerBoundaryFlags(INVERT_DC_GRAD); + } diagnose = options["diagnose"] .doc("Output additional diagnostics?") .withDefault(false); + + magnetic_flutter = options["magnetic_flutter"] + .doc("Set magnetic flutter terms (Apar_flutter)?") + .withDefault(false); } void Electromagnetic::transform(Options &state) { AUTO_TRACE(); - + Options& allspecies = state["species"]; // Sum coefficients over species @@ -58,7 +98,7 @@ void Electromagnetic::transform(Options &state) { const BoutReal A = get(species["AA"]); // Coefficient in front of A_|| - alpha_em += N * (SQ(Z) / A); + alpha_em += floor(N, 1e-5) * (SQ(Z) / A); // Right hand side Ajpar += mom * (Z / A); @@ -66,7 +106,46 @@ void Electromagnetic::transform(Options &state) { // Invert Helmholtz equation for Apar aparSolver->setCoefA((-beta_em) * alpha_em); - Apar = aparSolver->solve((-beta_em) * Ajpar); + + if (const_gradient) { + // Set gradient boundary condition from gradient inside boundary + Field3D rhs = (-beta_em) * Ajpar; + + const auto* mesh = Apar.getMesh(); + const auto* coords = Apar.getCoordinates(); + + BoutReal time = get(state["time"]); + BoutReal weight = 1.0; + if (time > last_time) { + weight = exp((last_time - time) / apar_boundary_timescale); + } + last_time = time; + + if (mesh->firstX()) { + const int x = mesh->xstart - 1; + for (int y = mesh->ystart; y <= mesh->yend; y++) { + for (int z = mesh->zstart; z <= mesh->zend; z++) { + rhs(x, y, z) = (weight * (Apar(x + 1, y, z) - Apar(x, y, z)) + + (1 - weight) * (Apar(x + 2, y, z) - Apar(x + 1, y, z))) / + (sqrt(coords->g_11(x, y)) * coords->dx(x, y)); + } + } + } + if (mesh->lastX()) { + const int x = mesh->xend + 1; + for (int y = mesh->ystart; y <= mesh->yend; y++) { + for (int z = mesh->zstart; z <= mesh->zend; z++) { + rhs(x, y, z) = (weight * (Apar(x, y, z) - Apar(x - 1, y, z)) + + (1 - weight) * (Apar(x - 1, y, z) - Apar(x - 2, y, z))) / + sqrt(coords->g_11(x, y)) / coords->dx(x, y); + } + } + } + // Use previous value of Apar as initial guess + Apar = aparSolver->solve(rhs, Apar); + } else { + Apar = aparSolver->solve((-beta_em) * Ajpar, Apar); + } // Save in the state set(state["fields"]["Apar"], Apar); @@ -89,13 +168,45 @@ void Electromagnetic::transform(Options &state) { nv -= Z * N * Apar; // Note: velocity is momentum / (A * N) Field3D v = getNonFinal(species["velocity"]); - v -= (Z / A) * Apar; + v -= (Z / A) * N * Apar / floor(N, 1e-5); // Need to update the guard cells bout::globals::mesh->communicate(nv, v); + v.applyBoundary("dirichlet"); + nv.applyBoundary("dirichlet"); set(species["momentum"], nv); set(species["velocity"], v); } + + if (magnetic_flutter) { + // Magnetic flutter terms + Apar_flutter = Apar - DC(Apar); + + // Ensure that guard cells are communicated + Apar.getMesh()->communicate(Apar_flutter); + + set(state["fields"]["Apar_flutter"], Apar_flutter); + +#if 0 + // Create a vector A from covariant components + // (A^x, A^y, A^z) + // Note: b = e_y / (JB) + const auto* coords = Apar.getCoordinates(); + Vector3D A; + A.covariant = true; + A.x = A.z = 0.0; + A.y = Apar_flutter * (coords->J * coords->Bxy); + + // Perturbed magnetic field vector + // Note: Contravariant components (dB_x, dB_y, dB_z) + Vector3D delta_B = Curl(A); + + // Set components of the perturbed unit vector + // Note: Options can't (yet) contain vectors + set(state["fields"]["deltab_flutter_x"], delta_B.x / coords->Bxy); + set(state["fields"]["deltab_flutter_z"], delta_B.z / coords->Bxy); +#endif + } } void Electromagnetic::outputVars(Options &state) { @@ -104,8 +215,8 @@ void Electromagnetic::outputVars(Options &state) { auto rho_s0 = get(state["rho_s0"]); set_with_attrs(state["beta_em"], beta_em, { - {"long_name", "Helmholtz equation parameter"} - }); + {"long_name", "Helmholtz equation parameter"} + }); set_with_attrs(state["Apar"], Apar, { {"time_dimension", "t"}, @@ -115,6 +226,16 @@ void Electromagnetic::outputVars(Options &state) { {"long_name", "Parallel component of vector potential A"} }); + if (magnetic_flutter) { + set_with_attrs(state["Apar_flutter"], Apar_flutter, { + {"time_dimension", "t"}, + {"units", "T m"}, + {"conversion", Bnorm * rho_s0}, + {"standard_name", "b dot A"}, + {"long_name", "Vector potential A|| used in flutter terms"} + }); + } + if (diagnose) { set_with_attrs(state["Ajpar"], Ajpar, { {"time_dimension", "t"}, diff --git a/src/evolve_density.cxx b/src/evolve_density.cxx index 892627b4..33ed120b 100644 --- a/src/evolve_density.cxx +++ b/src/evolve_density.cxx @@ -124,8 +124,6 @@ EvolveDensity::EvolveDensity(std::string name, Options& alloptions, Solver* solv } } - - neumann_boundary_average_z = alloptions[std::string("N") + name]["neumann_boundary_average_z"] .doc("Apply neumann boundary with Z average?") .withDefault(false); @@ -232,6 +230,14 @@ void EvolveDensity::finally(const Options& state) { } ddt(N) -= FV::Div_par_mod(N, V, fastest_wave); + + if (state.isSection("fields") and state["fields"].isSet("Apar_flutter")) { + // Magnetic flutter term + const Field3D Apar_flutter = get(state["fields"]["Apar_flutter"]); + // Note: Using -Apar_flutter rather than reversing sign in front, + // so that upwinding is handled correctly + ddt(N) -= Div_n_g_bxGrad_f_B_XZ(N, V, -Apar_flutter); + } } if (low_n_diffuse) { diff --git a/src/evolve_energy.cxx b/src/evolve_energy.cxx index 75b8f806..a6ce0e64 100644 --- a/src/evolve_energy.cxx +++ b/src/evolve_energy.cxx @@ -257,6 +257,12 @@ void EvolveEnergy::finally(const Options& state) { } ddt(E) -= FV::Div_par_mod(E + P, V, fastest_wave); + + if (state.isSection("fields") and state["fields"].isSet("Apar_flutter")) { + // Magnetic flutter term + const Field3D Apar_flutter = get(state["fields"]["Apar_flutter"]); + ddt(E) -= Div_n_g_bxGrad_f_B_XZ(E + P, V, -Apar_flutter); + } } if (species.isSet("low_n_coeff")) { @@ -322,6 +328,20 @@ void EvolveEnergy::finally(const Options& state) { // Note: Flux through boundary turned off, because sheath heat flux // is calculated and removed separately ddt(E) += FV::Div_par_K_Grad_par(kappa_par, T, false); + + if (state.isSection("fields") and state["fields"].isSet("Apar_flutter")) { + // Magnetic flutter term. The operator splits into 4 pieces: + // Div(k b b.Grad(T)) = Div(k b0 b0.Grad(T)) + Div(k d0 db.Grad(T)) + // + Div(k db b0.Grad(T)) + Div(k db db.Grad(T)) + // The first term is already calculated above. + // Here we add the terms containing db + const Field3D Apar_flutter = get(state["fields"]["Apar_flutter"]); + Field3D db_dot_T = bracket(T, Apar_flutter, BRACKET_ARAKAWA); + Field3D b0_dot_T = Grad_par(T); + mesh->communicate(db_dot_T, b0_dot_T); + ddt(E) += Div_par(kappa_par * db_dot_T) + - Div_n_g_bxGrad_f_B_XZ(kappa_par, db_dot_T + b0_dot_T, Apar_flutter); + } } if (hyper_z > 0.) { @@ -336,6 +356,11 @@ void EvolveEnergy::finally(const Options& state) { if (species.isSet("energy_source")) { Se += get(species["energy_source"]); // For diagnostic output } +#if CHECKLEVEL >= 1 + if (species.isSet("pressure_source")) { + throw BoutException("Components must evolve `energy_source` rather then `pressure_source`"); + } +#endif if (species.isSet("momentum_source")) { Se += V * get(species["momentum_source"]); } diff --git a/src/evolve_momentum.cxx b/src/evolve_momentum.cxx index 38631e91..9aa1289a 100644 --- a/src/evolve_momentum.cxx +++ b/src/evolve_momentum.cxx @@ -58,6 +58,9 @@ EvolveMomentum::EvolveMomentum(std::string name, Options &alloptions, Solver *so fix_momentum_boundary_flux = options["fix_momentum_boundary_flux"] .doc("Fix Y boundary momentum flux to boundary midpoint value?") .withDefault(false); + + // Set to zero so set for output + momentum_source = 0.0; } void EvolveMomentum::transform(Options &state) { @@ -78,6 +81,7 @@ void EvolveMomentum::transform(Options &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); } @@ -88,6 +92,7 @@ void EvolveMomentum::finally(const Options &state) { BoutReal AA = get(species["AA"]); // Get updated momentum with boundary conditions + // Note: This momentum may be modified by electromagnetic terms NV = get(species["momentum"]); // Get the species density @@ -95,15 +100,27 @@ void EvolveMomentum::finally(const Options &state) { // Apply a floor to the density Field3D Nlim = floor(N, density_floor); + // Typical wave speed used for numerical diffusion + Field3D fastest_wave; + if (state.isSet("fastest_wave")) { + fastest_wave = get(state["fastest_wave"]); + } else { + Field3D T = get(species["temperature"]); + fastest_wave = sqrt(T / AA); + } + + // Parallel flow + V = get(species["velocity"]); + if (state.isSection("fields") and state["fields"].isSet("phi") and species.isSet("charge")) { - BoutReal Z = get(species["charge"]); + const BoutReal Z = get(species["charge"]); if (Z != 0.0) { // Electrostatic potential set and species has charge // -> include ExB flow and parallel force - Field3D phi = get(state["fields"]["phi"]); + const Field3D phi = get(state["fields"]["phi"]); ddt(NV) = -Div_n_bxGrad_f_B_XPPM(NV, phi, bndry_flux, poloidal_flows, true); // ExB drift @@ -111,35 +128,56 @@ void EvolveMomentum::finally(const Options &state) { // Parallel electric field // Force density = - Z N ∇ϕ ddt(NV) -= Z * N * Grad_par(phi); + + if (state["fields"].isSet("Apar")) { + // Include a correction term for electromagnetic simulations + const Field3D Apar = get(state["fields"]["Apar"]); + + const Field3D density_source = species.isSet("density_source") ? + get(species["density_source"]) + : zeroFrom(Apar); + + // This is Z * Apar * dn/dt, keeping just leading order terms + Field3D dndt = density_source + - FV::Div_par_mod(N, V, fastest_wave) + - Div_n_bxGrad_f_B_XPPM(N, phi, bndry_flux, poloidal_flows, true) + ; + if (low_n_diffuse_perp) { + dndt += Div_Perp_Lap_FV_Index(density_floor / floor(N, 1e-3 * density_floor), N, + bndry_flux); + } + ddt(NV) += Z * Apar * dndt; + } + } else { + ddt(NV) = 0.0; } } else { ddt(NV) = 0.0; } - // Parallel flow - V = get(species["velocity"]); - - // Typical wave speed used for numerical diffusion - Field3D fastest_wave; - if (state.isSet("fastest_wave")) { - fastest_wave = get(state["fastest_wave"]); - } else { - Field3D T = get(species["temperature"]); - fastest_wave = sqrt(T / AA); - } - // Note: // - Density floor should be consistent with calculation of V // 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"]); ddt(NV) -= Grad_par(P); } + if (state.isSection("fields") and state["fields"].isSet("Apar_flutter")) { + // Magnetic flutter term + const Field3D Apar_flutter = get(state["fields"]["Apar_flutter"]); + ddt(NV) -= Div_n_g_bxGrad_f_B_XZ(NV, V, -Apar_flutter); + + if (species.isSet("pressure")) { + Field3D P = get(species["pressure"]); + ddt(NV) -= bracket(P, Apar_flutter, BRACKET_ARAKAWA); + } + } + if (species.isSet("low_n_coeff")) { // Low density parallel diffusion Field3D low_n_coeff = get(species["low_n_coeff"]); @@ -168,7 +206,9 @@ void EvolveMomentum::finally(const Options &state) { // If N < density_floor then NV and NV_solver may differ // -> Add term to force NV_solver towards NV - ddt(NV) += NV - NV_solver; + // 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")) { @@ -193,6 +233,11 @@ void EvolveMomentum::finally(const Options &state) { flow_ylow = get(species["momentum_flow_ylow"]); } } + // Restore NV to the value returned by the solver + // so that restart files contain the correct values + // Note: Copy boundary condition so dump file has correct boundary. + NV_solver.setBoundaryTo(NV); + NV = NV_solver; } void EvolveMomentum::outputVars(Options &state) { diff --git a/src/evolve_pressure.cxx b/src/evolve_pressure.cxx index 97289538..85b41a95 100644 --- a/src/evolve_pressure.cxx +++ b/src/evolve_pressure.cxx @@ -261,6 +261,13 @@ void EvolvePressure::finally(const Options& state) { ddt(P) += (2. / 3) * V * Grad_par(P); } + + 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) += (2. / 3) * V * bracket(P, Apar_flutter, BRACKET_ARAKAWA); + } } if (species.isSet("low_n_coeff")) { @@ -339,6 +346,20 @@ void EvolvePressure::finally(const Options& state) { // Note: Flux through boundary turned off, because sheath heat flux // is calculated and removed separately ddt(P) += (2. / 3) * FV::Div_par_K_Grad_par(kappa_par, T, false); + + if (state.isSection("fields") and state["fields"].isSet("Apar_flutter")) { + // Magnetic flutter term. The operator splits into 4 pieces: + // Div(k b b.Grad(T)) = Div(k b0 b0.Grad(T)) + Div(k d0 db.Grad(T)) + // + Div(k db b0.Grad(T)) + Div(k db db.Grad(T)) + // The first term is already calculated above. + // Here we add the terms containing db + const Field3D Apar_flutter = get(state["fields"]["Apar_flutter"]); + Field3D db_dot_T = bracket(T, Apar_flutter, BRACKET_ARAKAWA); + Field3D b0_dot_T = Grad_par(T); + mesh->communicate(db_dot_T, b0_dot_T); + ddt(P) += (2. / 3) * (Div_par(kappa_par * db_dot_T) - + Div_n_g_bxGrad_f_B_XZ(kappa_par, db_dot_T + b0_dot_T, Apar_flutter)); + } } if (hyper_z > 0.) { @@ -365,6 +386,11 @@ void EvolvePressure::finally(const Options& state) { if (species.isSet("energy_source")) { Sp += (2. / 3) * get(species["energy_source"]); // For diagnostic output } +#if CHECKLEVEL >= 1 + if (species.isSet("pressure_source")) { + throw BoutException("Components must evolve `energy_source` rather then `pressure_source`"); + } +#endif ddt(P) += Sp; // Term to force evolved P towards N * T @@ -427,7 +453,7 @@ void EvolvePressure::outputVars(Options& state) { set_with_attrs(state[std::string("kappa_par_") + name], kappa_par, {{"time_dimension", "t"}, {"units", "W / m / eV"}, - {"conversion", Pnorm * Omega_ci * SQ(rho_s0)}, + {"conversion", (Pnorm * Omega_ci * SQ(rho_s0) )/ Tnorm}, {"long_name", name + " heat conduction coefficient"}, {"species", name}, {"source", "evolve_pressure"}}); diff --git a/src/hydrogen_charge_exchange.cxx b/src/hydrogen_charge_exchange.cxx index b1781a14..39551b28 100644 --- a/src/hydrogen_charge_exchange.cxx +++ b/src/hydrogen_charge_exchange.cxx @@ -6,7 +6,8 @@ void HydrogenChargeExchange::calculate_rates(Options& atom1, Options& ion1, Field3D &atom_mom, Field3D &ion_mom, Field3D &atom_energy, Field3D &ion_energy, Field3D &atom_rate, Field3D &ion_rate, - BoutReal &rate_multiplier) { + BoutReal &rate_multiplier, + bool &no_neutral_cx_mom_gain) { // Temperatures and masses of initial atom and ion const Field3D Tatom = get(atom1["temperature"]); @@ -63,7 +64,9 @@ void HydrogenChargeExchange::calculate_rates(Options& atom1, Options& ion1, // Transfer fom atom1 to ion2 atom_mom = R * Aatom * atom1_velocity; subtract(atom1["momentum_source"], atom_mom); - add(ion2["momentum_source"], atom_mom); + if (no_neutral_cx_mom_gain == false) { + add(ion2["momentum_source"], atom_mom); + } // Transfer from ion1 to atom2 ion_mom = R * Aion * ion1_velocity; diff --git a/src/neutral_parallel_diffusion.cxx b/src/neutral_parallel_diffusion.cxx index f533544f..744a2653 100644 --- a/src/neutral_parallel_diffusion.cxx +++ b/src/neutral_parallel_diffusion.cxx @@ -26,29 +26,44 @@ void NeutralParallelDiffusion::transform(Options& state) { const Field3D Pn = IS_SET(species["pressure"]) ? GET_VALUE(Field3D, species["pressure"]) : Nn * Tn; - // Diffusion coefficient + BoutReal advection_factor = 0; + BoutReal kappa_factor = 0; + + if (equation_fix) { + advection_factor = (5. / 2); // This is equivalent to 5/3 if on pressure basis + kappa_factor = (5. / 2); + } else { + advection_factor = (3. / 2); + kappa_factor = 1; + } + + // Pressure-diffusion coefficient Field3D Dn = dneut * Tn / (AA * nu); Dn.applyBoundary("dirichlet_o2"); mesh->communicate(Dn); // Cross-field diffusion calculated from pressure gradient + // This is the pressure-diffusion approximation Field3D logPn = log(floor(Pn, 1e-7)); logPn.applyBoundary("neumann"); - // Particle diffusion + // Particle advection Field3D S = FV::Div_par_K_Grad_par(Dn * Nn, logPn); add(species["density_source"], S); - Field3D kappa_n = Nn * Dn; + Field3D kappa_n = kappa_factor * Nn * Dn; kappa_n.applyBoundary("neumann"); - // Heat conduction - Field3D E = FV::Div_par_K_Grad_par(kappa_n, Tn) // Parallel - + FV::Div_par_K_Grad_par(Dn * (3. / 2) * Pn, logPn); // Perpendicular diffusion + // Heat transfer + Field3D E = + FV::Div_par_K_Grad_par( + Dn * advection_factor * Pn, logPn); // Pressure advection + if (thermal_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"])) { + if (IS_SET(species["velocity"]) and 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/sheath_boundary.cxx b/src/sheath_boundary.cxx index bb3d9b08..9d416c61 100644 --- a/src/sheath_boundary.cxx +++ b/src/sheath_boundary.cxx @@ -200,8 +200,9 @@ void SheathBoundary::transform(Options &state) { // Note: Needed to get past initial conditions, perhaps transients // but this shouldn't happen in steady state - if (fabs(grad_ni) < 1e-3) { - grad_ni = grad_ne = 1e-3; // Remove kinetic correction term + // Note: Factor of 2 to be consistent with later calculation + if (fabs(grad_ni) < 2e-3) { + grad_ni = grad_ne = 2e-3; // Remove kinetic correction term } BoutReal C_i_sq = clip( @@ -234,12 +235,12 @@ void SheathBoundary::transform(Options &state) { BoutReal ti = Ti[i]; // Equation (9) in Tskhakaya 2005 - + BoutReal grad_ne = Ne[i] - Ne[im]; BoutReal grad_ni = Ni[i] - Ni[im]; - if (fabs(grad_ni) < 1e-3) { - grad_ni = grad_ne = 1e-3; // Remove kinetic correction term + if (fabs(grad_ni) < 2e-3) { + grad_ni = grad_ne = 2e-3; // Remove kinetic correction term } BoutReal C_i_sq = clip( @@ -247,8 +248,7 @@ void SheathBoundary::transform(Options &state) { / Mi, 0, 100); // Limit for e.g. Ni zero gradient - ion_sum(r.ind, mesh->yend, jz) += s_i * Zi * sqrt(C_i_sq); - + ion_sum[i] += s_i * Zi * sin_alpha * sqrt(C_i_sq); } } } @@ -348,6 +348,10 @@ void SheathBoundary::transform(Options &state) { BoutReal q = ((gamma_e - 1 - 1 / (electron_adiabatic - 1)) * tesheath - 0.5 * Me * SQ(vesheath)) * nesheath * vesheath; + if (q > 0.0) { + // Note: This could happen if tesheath > 2*phisheath + q = 0.0; + } // Multiply by cell area to get power BoutReal flux = q * (coord->J[i] + coord->J[im]) @@ -410,7 +414,10 @@ void SheathBoundary::transform(Options &state) { BoutReal q = ((gamma_e - 1 - 1 / (electron_adiabatic - 1)) * tesheath - 0.5 * Me * SQ(vesheath)) * nesheath * vesheath; - + if (q < 0.0) { + // Note: This could happen if tesheath > 2*phisheath + q = 0.0; + } // Multiply by cell area to get power BoutReal flux = q * (coord->J[i] + coord->J[ip]) / (sqrt(coord->g_22[i]) + sqrt(coord->g_22[ip])); @@ -429,6 +436,10 @@ void SheathBoundary::transform(Options &state) { } // Set electron density and temperature, now with boundary conditions + // Note: Clear parallel slices because they do not contain correct boundary conditions + Ne.clearParallelSlices(); + Te.clearParallelSlices(); + Pe.clearParallelSlices(); setBoundary(electrons["density"], fromFieldAligned(Ne)); setBoundary(electrons["temperature"], fromFieldAligned(Te)); setBoundary(electrons["pressure"], fromFieldAligned(Pe)); @@ -438,17 +449,18 @@ void SheathBoundary::transform(Options &state) { set(electrons["energy_source"], fromFieldAligned(electron_energy_source)); if (IS_SET_NOBOUNDARY(electrons["velocity"])) { + Ve.clearParallelSlices(); setBoundary(electrons["velocity"], fromFieldAligned(Ve)); } if (IS_SET_NOBOUNDARY(electrons["momentum"])) { + NVe.clearParallelSlices(); setBoundary(electrons["momentum"], fromFieldAligned(NVe)); } if (always_set_phi or IS_SET_NOBOUNDARY(state["fields"]["phi"])) { // Set the potential, including boundary conditions - phi = fromFieldAligned(phi); - //output.write("-> phi {}\n", phi(10, mesh->yend+1, 0)); - setBoundary(state["fields"]["phi"], phi); + phi.clearParallelSlices(); + setBoundary(state["fields"]["phi"], fromFieldAligned(phi)); } ////////////////////////////////////////////////////////////////// @@ -525,7 +537,7 @@ void SheathBoundary::transform(Options &state) { // 1 / (1 + ∂_{ln n_e} ln s_i = s_i ∂_z n_e / ∂_z n_i // (from comparing C_i^2 in eq. 9 with eq. 20 // - // + // BoutReal s_i = clip(nisheath / floor(nesheath, 1e-10), 0, 1); // Concentration BoutReal grad_ne = Ne[i] - nesheath; BoutReal grad_ni = Ni[i] - nisheath; @@ -600,7 +612,7 @@ void SheathBoundary::transform(Options &state) { // Ion sheath heat transmission coefficient // // 1 / (1 + ∂_{ln n_e} ln s_i = s_i * ∂n_e / (s_i * ∂n_e + ∂ n_i) - BoutReal s_i = (nesheath > 1e-5) ? nisheath / nesheath : 0.0; // Concentration + BoutReal s_i = clip(nisheath / floor(nesheath, 1e-10), 0, 1); // Concentration BoutReal grad_ne = Ne[i] - nesheath; BoutReal grad_ni = Ni[i] - nisheath; @@ -650,15 +662,20 @@ void SheathBoundary::transform(Options &state) { // Finished boundary conditions for this species // Put the modified fields back into the state. + Ni.clearParallelSlices(); + Ti.clearParallelSlices(); + Pi.clearParallelSlices(); setBoundary(species["density"], fromFieldAligned(Ni)); setBoundary(species["temperature"], fromFieldAligned(Ti)); setBoundary(species["pressure"], fromFieldAligned(Pi)); if (species.isSet("velocity")) { + Vi.clearParallelSlices(); setBoundary(species["velocity"], fromFieldAligned(Vi)); } if (species.isSet("momentum")) { + NVi.clearParallelSlices(); setBoundary(species["momentum"], fromFieldAligned(NVi)); } diff --git a/src/sheath_boundary_simple.cxx b/src/sheath_boundary_simple.cxx index f987d787..fde1ce56 100644 --- a/src/sheath_boundary_simple.cxx +++ b/src/sheath_boundary_simple.cxx @@ -409,6 +409,10 @@ void SheathBoundarySimple::transform(Options& state) { } // Set electron density and temperature, now with boundary conditions + // Note: Clear parallel slices because they do not contain boundary conditions. + Ne.clearParallelSlices(); + Te.clearParallelSlices(); + Pe.clearParallelSlices(); setBoundary(electrons["density"], fromFieldAligned(Ne)); setBoundary(electrons["temperature"], fromFieldAligned(Te)); setBoundary(electrons["pressure"], fromFieldAligned(Pe)); @@ -421,14 +425,17 @@ void SheathBoundarySimple::transform(Options& state) { add(electrons["energy_flow_ylow"], fromFieldAligned(electron_sheath_power_ylow)); if (IS_SET_NOBOUNDARY(electrons["velocity"])) { + Ve.clearParallelSlices(); setBoundary(electrons["velocity"], fromFieldAligned(Ve)); } if (IS_SET_NOBOUNDARY(electrons["momentum"])) { + NVe.clearParallelSlices(); setBoundary(electrons["momentum"], fromFieldAligned(NVe)); } if (always_set_phi or (state.isSection("fields") and state["fields"].isSet("phi"))) { // Set the potential, including boundary conditions + phi.clearParallelSlices(); setBoundary(state["fields"]["phi"], fromFieldAligned(phi)); } @@ -605,15 +612,20 @@ void SheathBoundarySimple::transform(Options& state) { } // Finished boundary conditions for this species // Put the modified fields back into the state. + Ni.clearParallelSlices(); + Ti.clearParallelSlices(); + Pi.clearParallelSlices(); setBoundary(species["density"], fromFieldAligned(Ni)); setBoundary(species["temperature"], fromFieldAligned(Ti)); setBoundary(species["pressure"], fromFieldAligned(Pi)); if (species.isSet("velocity")) { + Vi.clearParallelSlices(); setBoundary(species["velocity"], fromFieldAligned(Vi)); } if (species.isSet("momentum")) { + NVi.clearParallelSlices(); setBoundary(species["momentum"], fromFieldAligned(NVi)); } diff --git a/src/vorticity.cxx b/src/vorticity.cxx index 936f71e1..496f7296 100644 --- a/src/vorticity.cxx +++ b/src/vorticity.cxx @@ -161,6 +161,10 @@ Vorticity::Vorticity(std::string name, Options& alloptions, Solver* solver) { // the first time RHS function is called phi_boundary_last_update = -1.; + phi_core_averagey = options["phi_core_averagey"] + .doc("Average phi core boundary in Y?") + .withDefault(false) and mesh->periodicY(mesh->xstart); + phi_boundary_timescale = options["phi_boundary_timescale"] .doc("Timescale for phi boundary relaxation [seconds]") .withDefault(1e-4) @@ -217,6 +221,7 @@ Vorticity::Vorticity(std::string name, Options& alloptions, Solver* solver) { void Vorticity::transform(Options& state) { AUTO_TRACE(); + phi.name = "phi"; auto& fields = state["fields"]; // Set the boundary of phi. Both 2D and 3D fields are kept, though the 3D field @@ -270,12 +275,32 @@ void Vorticity::transform(Options& state) { phi_boundary_last_update = time; if (mesh->firstX()) { + BoutReal phivalue = 0.0; + if (phi_core_averagey) { + BoutReal philocal = 0.0; + for (int j = mesh->ystart; j <= mesh->yend; j++) { + for (int k = 0; k < mesh->LocalNz; k++) { + philocal += phi(mesh->xstart, j, k); + } + } + MPI_Comm comm_inner = mesh->getYcomm(0); + int np; + MPI_Comm_size(comm_inner, &np); + MPI_Allreduce(&philocal, + &phivalue, + 1, MPI_DOUBLE, + MPI_SUM, comm_inner); + phivalue /= (np * mesh->LocalNz * mesh->LocalNy); + } + for (int j = mesh->ystart; j <= mesh->yend; j++) { - BoutReal phivalue = 0.0; - for (int k = 0; k < mesh->LocalNz; k++) { - phivalue += phi(mesh->xstart, j, k); + if (!phi_core_averagey) { + phivalue = 0.0; // Calculate phi boundary for each Y index separately + for (int k = 0; k < mesh->LocalNz; k++) { + phivalue += phi(mesh->xstart, j, k); + } + phivalue /= mesh->LocalNz; // Average in Z of point next to boundary } - phivalue /= mesh->LocalNz; // Average in Z of point next to boundary // Old value of phi at boundary BoutReal oldvalue = diff --git a/tests/integrated/1D-recycling/runtest b/tests/integrated/1D-recycling/runtest index a8e9ab1b..6664f3ad 100755 --- a/tests/integrated/1D-recycling/runtest +++ b/tests/integrated/1D-recycling/runtest @@ -44,9 +44,9 @@ if Te_up < 50 or Te_up > 70: Ti = collect("Td+", tind=-1, path=path) Ti_up = Ti[-1,0,0,0] * Tnorm # Upstream ion temperature should be about 140eV -if Ti_up < 130 or Ti_up > 150: +if Ti_up < 140 or Ti_up > 160: success = False - print("Ion temperature failed: {}eV. Expecting about 140eV".format(Ti_up)) + print("Ion temperature failed: {}eV. Expecting about 150eV".format(Ti_up)) if success: print(" => Test passed") diff --git a/tests/integrated/alfven-wave/README.md b/tests/integrated/alfven-wave/README.md new file mode 100644 index 00000000..8dc61a6d --- /dev/null +++ b/tests/integrated/alfven-wave/README.md @@ -0,0 +1,74 @@ +# Electromagnetic Alfven wave with finite electron mass + +![Alfven wave speed](alfven-wave.png) + +## Equations + +The ion density is fixed, with zero velocity and fixed temperature: +``` +[i] +type = fixed_density, fixed_velocity, fixed_temperature + +charge = 1 +AA = 1 + +density = 1e19 # [m^-3] +velocity = 0 +temperature = 100 # eV +``` + +The electron density is set to the ion density by quasi-neutrality, +and only the parallel momentum is evolved: +``` +[e] +type = quasineutral, evolve_momentum, fixed_temperature + +charge = -1 +AA = 1./1836 + +temperature = 100 # eV +``` +Note that the momentum is the canonical momentum because the +`electromagnetic` component calculates and corrects for the vector +potential term in the parallel momentum of all (charged) species. + +Finally, the potential is evolved by a vorticity equation. + +## Slab domain + +The domain is a thin slab, with one cell in the X (radial) direction, +and periodic in both Y (parallel) and Z (binormal) directions. +The Y direction is used to set k_parallel, and the Z direction sets k_perp. + +The number of cells in each dimension are specified by `nx`, `ny` and +`nz`. Note that `nx` includes 2 boundary cells on either side, so `nx += 5` is one cell in the middle, and 4 boundary cells. + +The size of the domain in each dimension in meters is set by `Lx`, `Ly` and `Lz`. +These are not directly used by Hermes-3 or BOUT++, but are used to calculate the +metric tensor components and cell sizes. + +The Y and Z cell sizes are `dy = Ly / ny` and `dz = Lz / nz`. +The X cell size is more complicated because the field-aligned +operators assume a Clebsch coordinate system, in which B = nabla z +cross nabla x. See the [coordinates manual page](https://bout-dev.readthedocs.io/en/stable/user_docs/coordinates.html#magnetic-field) +for details. The `dx` cell size therefore includes a factor of the +magnetic field strength `B`. With this choice the metric tensor is +diagonal, with elements +``` +g11 = B^2 +g22 = 1 +g33 = 1 +``` +and the Jacobian is `J = 1 / B` + +## Boundary conditions + +The radial boundary conditions on potential are Neumann for the +non-constant (AC) components, and zero for the constant (DC) component: +``` +[vorticity:laplacian] +inner_boundary_flags = 2 +outer_boundary_flags = 2 +``` +All other boundary conditions are Neumann. diff --git a/tests/integrated/alfven-wave/alfven-wave.png b/tests/integrated/alfven-wave/alfven-wave.png new file mode 100644 index 00000000..3d91cd0a Binary files /dev/null and b/tests/integrated/alfven-wave/alfven-wave.png differ diff --git a/tests/integrated/alfven-wave/data/BOUT.inp b/tests/integrated/alfven-wave/data/BOUT.inp new file mode 100644 index 00000000..db8580a7 --- /dev/null +++ b/tests/integrated/alfven-wave/data/BOUT.inp @@ -0,0 +1,83 @@ +nout = 50 +timestep = 10 + +[mesh] +nx = 5 # Radial. Includes 4 boundary cells +ny = 32 # Parallel direction +nz = 27 # Perpendicular direction + +Lx = 0.01 # X domain size [meters] +Ly = 2 # Y domain size [meters] +Lz = 0.01 # Z domain size [meters] + +B = 0.2 # Magnetic field [Tesla] + +ixseps1 = nx # Domain is periodic in Y and Z +ixseps2 = nx + +# Define an orthogonal slab mesh in Clebsch coordinates +# In terms of BOUT++ coordinates +# https://bout-dev.readthedocs.io/en/stable/user_docs/coordinates.html#id1 +# this corresponds to R = 1, hthe = 1, Bpol = B, Btor = 0, I = 0 + +dr = Lx / (nx - 4) +dx = dr * B +dy = Ly / ny +dz = Lz / nz + +g11 = B^2 +g22 = 1 +g33 = 1 + +J = 1 / B # Note: Volume = J * dx * dy * dx + +[mesh:paralleltransform] +type = identity + +[solver] +mxstep = 10000 + +[hermes] +components = i, e, electromagnetic, vorticity, sound_speed + +loadmetric = false # Load metric from Bxy, Rxy etc? + +[electromagnetic] +diagnose = true +apar_boundary_neumann = true # Neumann radial boundary conditions + +[vorticity] + +diamagnetic = false # Include diamagnetic current? +diamagnetic_polarisation = false # Include diamagnetic drift in polarisation current? +average_atomic_mass = 1 # Weighted average atomic mass, for polarisaion current +bndry_flux = false # Allow flows through radial boundaries +poloidal_flows = false # Include poloidal ExB flow + +[vorticity:laplacian] +inner_boundary_flags = 2 +outer_boundary_flags = 2 + +[Vort] +function = 1e-3 * sin(z - y) + +### Ions +[i] +type = fixed_density, fixed_velocity, fixed_temperature + +charge = 1 +AA = 1 + +density = 1e19 # [m^-3] +velocity = 0 +temperature = 100 # eV + +### Electrons +[e] +type = quasineutral, evolve_momentum, fixed_temperature + +charge = -1 +AA = 1/1836 + +temperature = 100 # eV + diff --git a/tests/integrated/alfven-wave/runtest b/tests/integrated/alfven-wave/runtest new file mode 100755 index 00000000..cf568a4e --- /dev/null +++ b/tests/integrated/alfven-wave/runtest @@ -0,0 +1,111 @@ +#!/usr/bin/env python3 +# +# Alfven wave test +# +# References: +# - M.V.Umansky et al. CPP 2008 https://doi.org/10.1002/ctpp.200810004 + +plotting = False +fail_early = True # Stop on the first failure? + +import numpy as np +from boututils.run_wrapper import shell, launch, getmpirun +import sys +from boutdata import collect + +# Link to the executable +shell("ln -s ../../../hermes-3 hermes-3") + +B = 0.2 # T +n0 = 1e19 # m^-3 +AA = 1 # amu + +mu0 = 4e-7 * np.pi +m_e = 9.11e-31 +qe = 1.602e-19 + +V_A = B / np.sqrt(mu0 * n0 * 1.67e-27 * AA) # Alfven speed +delta_e = np.sqrt(m_e / (mu0 * n0 * qe**2)) # Electron skin depth + +# List of +if plotting: + Lys = [2, 3, 4] + Lzs = [1, 0.1, 0.05, 0.02, 0.01] +else: + Lys = [2] + Lzs = [0.1, 0.01] + +if plotting: + import matplotlib.pyplot as plt + +passed = True + +for Ly in Lys: + kdeltas = [] + expected_results = [] + simulation_results = [] + + for Lz in Lzs: + # Run hermes-3 using MPI + print(f"Running Ly = {Ly} Lz = {Lz}") + shell("rm data/BOUT.*.nc") # Delete data files + s, out = launch(f"./hermes-3 mesh:B={B} i:density={n0} i:AA={AA} mesh:Ly={Ly} mesh:Lz={Lz}", nproc=1, pipe=True) + + # Save output to log file + f = open(f"run.log_{Ly}_{Lz}", "w") + f.write(out) + f.close() + + Omega_ci = collect("Omega_ci", path="data") + t_array = collect("t_array", path="data") + dt = (t_array[1] - t_array[0]) / Omega_ci # Seconds + + phi = collect("phi", xind=2, path="data").squeeze() + phi_ms = np.mean(phi**2, axis=(1,2)) + dphidt = np.gradient(phi_ms[1:]) + # Find zero crossings + inds = np.where(dphidt[1:] * dphidt[:-1] < 0)[0] + crossings = (inds * abs(dphidt[inds + 1]) + + (inds + 1) * abs(dphidt[inds])) / abs(dphidt[inds + 1] - dphidt[inds]) + + # Mean^2 value doubles frequency; two crossings per period + period = 4 * np.mean(crossings[1:] - crossings[:-1]) * dt # Seconds + + omega = 2.*np.pi / period + k_par = 2.*np.pi / Ly + k_perp = 2.*np.pi / Lz + + result = omega / k_par + expected = V_A / np.sqrt(1 + (k_perp * delta_e)**2) + + err = abs(result - expected) / expected + if err > 0.02: + # Failed! + print(f"Ly = {Ly}, Lz = {Lz} (k_par = {k_par}, k_perp = {k_perp}): Expected {expected} but got {result}. Error: {err}") + if fail_early: + sys.exit(1) + passed = False + + simulation_results.append(result) + expected_results.append(expected) + kdeltas.append(k_perp * delta_e) + + if plotting: + plt.plot(kdeltas, simulation_results, 'o', label=f'L = {Ly}m') + +if plotting: + kds = np.linspace(0, max(kdeltas), 100) + plt.plot(kds, V_A / np.sqrt(1 + kds**2), 'k', label='Analytic') + plt.legend() + plt.xlabel(r"$k_\perp \delta_e$") + plt.ylabel(r"$\omega / k_{||}$") + plt.savefig("alfven-wave.pdf") + plt.savefig("alfven-wave.png") + plt.show() + +if passed: + print(" => Test passed") + sys.exit(0) + +print(" => Test failed") +sys.exit(1)