Implement collision source terms for multi-ion MHD (#2213)

* Add collision source terms for multi-ion MHD

* Fixed docstring issue

* Use zero(eltype(u)) instead of 0

* Removed unexisting function from export

* Fixed more docstring issues

* Skip computation of collision sources of an ion species with itself

* Added missing elixir

* Improved documentation

* Double precision variables in initial condition

* Removed extra backtick

* Use CarpenterKennedy2N54 for the collision sources' test

* Rewrote temperature functions in elixir to increase coverage

* Apply suggestions from code review

Co-authored-by: Daniel Doehring <daniel.doehring@rwth-aachen.de>

* Removed unnecessary sqrt(v)^2

* Added DOIs and bib codes

* Added note about mesh

* Apply suggestions from code review

Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com>
Co-authored-by: Daniel Doehring <daniel.doehring@rwth-aachen.de>

* Update examples/tree_2d_dgsem/elixir_mhdmultiion_collisions.jl

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

* Changed index of species k'->l to avoid confusion

* Changed variable name

* Update examples/tree_2d_dgsem/elixir_mhdmultiion_collisions.jl

Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com>

* Added missing DOI

* Made example compatible with other real types

* psi defined separately in initial condition

---------

Co-authored-by: Daniel Doehring <daniel.doehring@rwth-aachen.de>
Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com>
Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>
Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>

examples/tree_2d_dgsem/elixir_mhdmultiion_collisions.jl

@@ -0,0 +1,157 @@
+
+using OrdinaryDiffEqLowStorageRK
+using Trixi
+
+###############################################################################
+# This elixir describes the frictional slowing of an ionized carbon fluid (C6+) with respect to another species 
+# of a background ionized carbon fluid with an initially nonzero relative velocity. It is the second slow-down
+# test (fluids with different densities) described in:
+# - Ghosh, D., Chapman, T. D., Berger, R. L., Dimits, A., & Banks, J. W. (2019). A 
+#   multispecies, multifluid model for laser–induced counterstreaming plasma simulations. 
+#   Computers & Fluids, 186, 38-57. [DOI: 10.1016/j.compfluid.2019.04.012](https://doi.org/10.1016/j.compfluid.2019.04.012).
+#
+# This is effectively a zero-dimensional case because the spatial gradients are zero, and we use it to test the
+# collision source terms.
+#
+# To run this physically relevant test, we use the following characteristic quantities to non-dimensionalize
+# the equations:
+# Characteristic length: L_inf = 1.00E-03 m (domain size)
+# Characteristic density: rho_inf = 1.99E+00 kg/m^3 (corresponds to a number density of 1e20 cm^{-3})
+# Characteristic vacuum permeability: mu0_inf = 1.26E-06 N/A^2 (for equations with mu0 = 1)
+# Characteristic gas constant: R_inf = 6.92237E+02 J/kg/K (specific gas constant for a Carbon fluid)
+# Characteristic velocity: V_inf = 1.00E+06 m/s
+#
+# The results of the paper can be reproduced using `source_terms = source_terms_collision_ion_ion` (i.e., only
+# taking into account ion-ion collisions). However, we include ion-electron collisions assuming a constant 
+# electron temperature of 1 keV in this elixir to test the function `source_terms_collision_ion_electron`
+
+# Return the electron pressure for a constant electron temperature Te = 1 keV
+function electron_pressure_constantTe(u, equations::IdealGlmMhdMultiIonEquations2D)
+    @unpack charge_to_mass = equations
+    Te = electron_temperature_constantTe(u, equations)
+    total_electron_charge = zero(eltype(u))
+    for k in eachcomponent(equations)
+        rho_k = u[3 + (k - 1) * 5 + 1]
+        total_electron_charge += rho_k * charge_to_mass[k]
+    end
+
+    # Boltzmann constant divided by elementary charge
+    RealT = eltype(u)
+    kB_e = convert(RealT, 7.86319034E-02) #[nondimensional]
+
+    return total_electron_charge * kB_e * Te
+end
+
+# Return the constant electron temperature Te = 1 keV
+function electron_temperature_constantTe(u, equations::IdealGlmMhdMultiIonEquations2D)
+    RealT = eltype(u)
+    return convert(RealT, 0.008029953773) # [nondimensional] = 1 [keV]
+end
+
+# semidiscretization of the ideal MHD equations
+equations = IdealGlmMhdMultiIonEquations2D(gammas = (5 / 3, 5 / 3),
+                                           charge_to_mass = (76.3049060157692000,
+                                                             76.3049060157692000), # [nondimensional]
+                                           gas_constants = (1.0, 1.0), # [nondimensional]
+                                           molar_masses = (1.0, 1.0), # [nondimensional]
+                                           ion_ion_collision_constants = [0.0 0.4079382480442680;
+                                                                          0.4079382480442680 0.0], # [nondimensional] (computed with eq (4.142) of Schunk & Nagy (2009))
+                                           ion_electron_collision_constants = (8.56368379833E-06,
+                                                                               8.56368379833E-06), # [nondimensional] (computed with eq (9) of Ghosh et al. (2019))
+                                           electron_pressure = electron_pressure_constantTe,
+                                           electron_temperature = electron_temperature_constantTe,
+                                           initial_c_h = 0.0) # Deactivate GLM divergence cleaning
+
+# Frictional slowing of an ionized carbon fluid with respect to another background carbon fluid in motion
+function initial_condition_slow_down(x, t, equations::IdealGlmMhdMultiIonEquations2D)
+    RealT = eltype(x)
+
+    v11 = convert(RealT, 0.6550877)
+    v21 = zero(RealT)
+    v2 = v3 = zero(RealT)
+    B1 = B2 = B3 = zero(RealT)
+    rho1 = convert(RealT, 0.1)
+    rho2 = one(RealT)
+
+    p1 = convert(RealT, 0.00040170535986)
+    p2 = convert(RealT, 0.00401705359856)
+
+    psi = zero(RealT)
+
+    return prim2cons(SVector(B1, B2, B3, rho1, v11, v2, v3, p1, rho2, v21, v2, v3, p2, psi),
+                     equations)
+end
+
+# Temperature of ion 1
+function temperature1(u, equations::IdealGlmMhdMultiIonEquations2D)
+    rho_1, _ = Trixi.get_component(1, u, equations)
+    p = pressure(u, equations)
+
+    return p[1] / (rho_1 * equations.gas_constants[1])
+end
+
+# Temperature of ion 2
+function temperature2(u, equations::IdealGlmMhdMultiIonEquations2D)
+    rho_2, _ = Trixi.get_component(2, u, equations)
+    p = pressure(u, equations)
+
+    return p[2] / (rho_2 * equations.gas_constants[2])
+end
+
+initial_condition = initial_condition_slow_down
+tspan = (0.0, 0.1) # 100 [ps]
+
+# Entropy conservative volume numerical fluxes with standard LLF dissipation at interfaces
+volume_flux = (flux_ruedaramirez_etal, flux_nonconservative_ruedaramirez_etal)
+surface_flux = (flux_lax_friedrichs, flux_nonconservative_central)
+
+solver = DGSEM(polydeg = 3, surface_flux = surface_flux,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+coordinates_min = (0.0, 0.0)
+coordinates_max = (1.0, 1.0)
+# We use a very coarse mesh because this is a 0-dimensional case
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 1,
+                n_cells_max = 1_000_000)
+
+# Ion-ion and ion-electron collision source terms
+# In this particular case, we can omit source_terms_lorentz because the magnetic field is zero!
+function source_terms(u, x, t, equations::IdealGlmMhdMultiIonEquations2D)
+    source_terms_collision_ion_ion(u, x, t, equations) +
+    source_terms_collision_ion_electron(u, x, t, equations)
+end
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms = source_terms)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1
+analysis_callback = AnalysisCallback(semi,
+                                     save_analysis = true,
+                                     interval = analysis_interval,
+                                     extra_analysis_integrals = (temperature1,
+                                                                 temperature2))
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+stepsize_callback = StepsizeCallback(cfl = 0.01) # Very small CFL due to the stiff source terms 
+
+save_restart = SaveRestartCallback(interval = 100,
+                                   save_final_restart = true)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_restart,
+                        stepsize_callback)
+
+###############################################################################
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false);
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            ode_default_options()..., callback = callbacks);

src/Trixi.jl

@@ -219,7 +219,8 @@ export boundary_condition_do_nothing,
        BoundaryConditionCoupled
 
 export initial_condition_convergence_test, source_terms_convergence_test,
-       source_terms_lorentz
+       source_terms_lorentz, source_terms_collision_ion_electron,
+       source_terms_collision_ion_ion
 export source_terms_harmonic
 export initial_condition_poisson_nonperiodic, source_terms_poisson_nonperiodic,
        boundary_condition_poisson_nonperiodic

src/equations/ideal_glm_mhd_multiion.jl

@@ -183,6 +183,22 @@ divergence_cleaning_field(u, equations::AbstractIdealGlmMhdMultiIonEquations) =
     return rho
 end
 
+@inline function pressure(u, equations::AbstractIdealGlmMhdMultiIonEquations)
+    B1, B2, B3, _ = u
+    p = zero(MVector{ncomponents(equations), real(equations)})
+    for k in eachcomponent(equations)
+        rho, rho_v1, rho_v2, rho_v3, rho_e = get_component(k, u, equations)
+        v1 = rho_v1 / rho
+        v2 = rho_v2 / rho
+        v3 = rho_v3 / rho
+        gamma = equations.gammas[k]
+        p[k] = (gamma - 1) *
+               (rho_e - 0.5f0 * rho * (v1^2 + v2^2 + v3^2) -
+                0.5f0 * (B1^2 + B2^2 + B3^2))
+    end
+    return SVector{ncomponents(equations), real(equations)}(p)
+end
+
 #Convert conservative variables to primitive
 function cons2prim(u, equations::AbstractIdealGlmMhdMultiIonEquations)
     @unpack gammas = equations
@@ -471,4 +487,173 @@ end
 
     return dissipation
 end
+
+@doc raw"""
+    source_terms_collision_ion_ion(u, x, t,
+                                   equations::AbstractIdealGlmMhdMultiIonEquations)
+
+Compute the ion-ion collision source terms for the momentum and energy equations of each ion species as
+```math
+\begin{aligned}
+  \vec{s}_{\rho_k \vec{v}_k} =&  \rho_k\sum_{l}\bar{\nu}_{kl}(\vec{v}_{l} - \vec{v}_k),\\
+  s_{E_k}  =& 
+    3 \sum_{l} \left(
+    \bar{\nu}_{kl} \frac{\rho_k M_1}{M_{l} + M_k} R_1 (T_{l} - T_k)
+    \right) + 
+    \sum_{l} \left(
+        \bar{\nu}_{kl} \rho_k \frac{M_{l}}{M_{l} + M_k} \|\vec{v}_{l} - \vec{v}_k\|^2
+        \right)
+        +
+        \vec{v}_k \cdot \vec{s}_{\rho_k \vec{v}_k},
+\end{aligned}
+```
+where ``M_k`` is the molar mass of ion species `k` provided in `equations.molar_masses`, 
+``R_k`` is the specific gas constant of ion species `k` provided in `equations.gas_constants`, and
+ ``\bar{\nu}_{kl}`` is the effective collision frequency of species `k` with species `l`, which is computed as
+```math
+\begin{aligned}
+  \bar{\nu}_{kl} = \bar{\nu}^1_{kl} \tilde{B}_{kl} \frac{\rho_{l}}{T_{k l}^{3/2}},
+\end{aligned}
+```
+with the so-called reduced temperature ``T_{k l}`` and the ion-ion collision constants ``\tilde{B}_{kl}`` provided
+in `equations.ion_electron_collision_constants` (see [`IdealGlmMhdMultiIonEquations2D`](@ref)).
+
+The additional coefficient ``\bar{\nu}^1_{kl}`` is a non-dimensional drift correction factor proposed by Rambo and Denavit.
+
+References:
+- P. Rambo, J. Denavit, Interpenetration and ion separation in colliding plasmas, Physics of Plasmas 1 (1994) 4050–4060.
+  [DOI: 10.1063/1.870875](https://doi.org/10.1063/1.870875).
+- Schunk, R. W., Nagy, A. F. (2000). Ionospheres: Physics, plasma physics, and chemistry. 
+  Cambridge university press. [DOI: 10.1017/CBO9780511635342](https://doi.org/10.1017/CBO9780511635342).
+"""
+function source_terms_collision_ion_ion(u, x, t,
+                                        equations::AbstractIdealGlmMhdMultiIonEquations)
+    s = zero(MVector{nvariables(equations), eltype(u)})
+    @unpack gas_constants, molar_masses, ion_ion_collision_constants = equations
+
+    prim = cons2prim(u, equations)
+
+    for k in eachcomponent(equations)
+        rho_k, v1_k, v2_k, v3_k, p_k = get_component(k, prim, equations)
+        T_k = p_k / (rho_k * gas_constants[k])
+
+        S_q1 = zero(eltype(u))
+        S_q2 = zero(eltype(u))
+        S_q3 = zero(eltype(u))
+        S_E = zero(eltype(u))
+        for l in eachcomponent(equations)
+            # Do not compute collisions of an ion species with itself
+            k == l && continue
+
+            rho_l, v1_l, v2_l, v3_l, p_l = get_component(l, prim, equations)
+            T_l = p_l / (rho_l * gas_constants[l])
+
+            # Reduced temperature
+            T_kl = (molar_masses[l] * T_k + molar_masses[k] * T_l) /
+                   (molar_masses[k] + molar_masses[l])
+
+            delta_v2 = (v1_l - v1_k)^2 + (v2_l - v2_k)^2 + (v3_l - v3_k)^2
+
+            # Compute collision frequency without drifting correction
+            v_kl = ion_ion_collision_constants[k, l] * rho_l / T_kl^(3 / 2)
+
+            # Correct the collision frequency with the drifting effect
+            z2 = delta_v2 / (p_l / rho_l + p_k / rho_k)
+            v_kl /= (1 + (2 / (9 * pi))^(1 / 3) * z2)^(3 / 2)
+
+            S_q1 += rho_k * v_kl * (v1_l - v1_k)
+            S_q2 += rho_k * v_kl * (v2_l - v2_k)
+            S_q3 += rho_k * v_kl * (v3_l - v3_k)
+
+            S_E += (3 * molar_masses[1] * gas_constants[1] * (T_l - T_k)
+                    +
+                    molar_masses[l] * delta_v2) * v_kl * rho_k /
+                   (molar_masses[k] + molar_masses[l])
+        end
+
+        S_E += (v1_k * S_q1 + v2_k * S_q2 + v3_k * S_q3)
+
+        set_component!(s, k, 0, S_q1, S_q2, S_q3, S_E, equations)
+    end
+    return SVector{nvariables(equations), real(equations)}(s)
+end
+
+@doc raw"""
+    source_terms_collision_ion_electron(u, x, t,
+                                        equations::AbstractIdealGlmMhdMultiIonEquations)
+
+Compute the ion-electron collision source terms for the momentum and energy equations of each ion species. We assume ``v_e = v^+`` 
+(no effect of currents on the electron velocity).
+
+The collision sources read as
+```math
+\begin{aligned}
+    \vec{s}_{\rho_k \vec{v}_k} =&  \rho_k \bar{\nu}_{ke} (\vec{v}_{e} - \vec{v}_k),
+    \\
+    s_{E_k}  =& 
+    3  \left(
+    \bar{\nu}_{ke} \frac{\rho_k M_{1}}{M_k} R_1 (T_{e} - T_k)
+    \right) 
+        +
+        \vec{v}_k \cdot \vec{s}_{\rho_k \vec{v}_k},
+\end{aligned}
+```
+where ``T_e`` is the electron temperature computed with the function `equations.electron_temperature`, 
+``M_k`` is the molar mass of ion species `k` provided in `equations.molar_masses`, 
+``R_k`` is the specific gas constant of ion species `k` provided in `equations.gas_constants`, and
+``\bar{\nu}_{ke}`` is the collision frequency of species `k` with the electrons, which is computed as
+```math
+\begin{aligned}
+  \bar{\nu}_{ke} = \tilde{B}_{ke} \frac{e n_e}{T_e^{3/2}},
+\end{aligned}
+```
+with the total electron charge ``e n_e`` (computed assuming quasi-neutrality), and the
+ion-electron collision coefficient ``\tilde{B}_{ke}`` provided in `equations.ion_electron_collision_constants`,
+which is scaled with the elementary charge (see [`IdealGlmMhdMultiIonEquations2D`](@ref)).
+
+References:
+- P. Rambo, J. Denavit, Interpenetration and ion separation in colliding plasmas, Physics of Plasmas 1 (1994) 4050–4060.
+  [DOI: 10.1063/1.870875](https://doi.org/10.1063/1.870875).
+- Schunk, R. W., Nagy, A. F. (2000). Ionospheres: Physics, plasma physics, and chemistry. 
+  Cambridge university press. [DOI: 10.1017/CBO9780511635342](https://doi.org/10.1017/CBO9780511635342).
+"""
+function source_terms_collision_ion_electron(u, x, t,
+                                             equations::AbstractIdealGlmMhdMultiIonEquations)
+    s = zero(MVector{nvariables(equations), eltype(u)})
+    @unpack gas_constants, molar_masses, ion_electron_collision_constants, electron_temperature = equations
+
+    prim = cons2prim(u, equations)
+    T_e = electron_temperature(u, equations)
+    T_e_power32 = T_e^(3 / 2)
+
+    v1_plus, v2_plus, v3_plus, vk1_plus, vk2_plus, vk3_plus = charge_averaged_velocities(u,
+                                                                                         equations)
+
+    # Compute total electron charge
+    total_electron_charge = zero(real(equations))
+    for k in eachcomponent(equations)
+        rho, _ = get_component(k, u, equations)
+        total_electron_charge += rho * equations.charge_to_mass[k]
+    end
+
+    for k in eachcomponent(equations)
+        rho_k, v1_k, v2_k, v3_k, p_k = get_component(k, prim, equations)
+        T_k = p_k / (rho_k * gas_constants[k])
+
+        # Compute effective collision frequency
+        v_ke = ion_electron_collision_constants[k] * total_electron_charge / T_e_power32
+
+        S_q1 = rho_k * v_ke * (v1_plus - v1_k)
+        S_q2 = rho_k * v_ke * (v2_plus - v2_k)
+        S_q3 = rho_k * v_ke * (v3_plus - v3_k)
+
+        S_E = 3 * molar_masses[1] * gas_constants[1] * (T_e - T_k) * v_ke * rho_k /
+              molar_masses[k]
+
+        S_E += (v1_k * S_q1 + v2_k * S_q2 + v3_k * S_q3)
+
+        set_component!(s, k, 0, S_q1, S_q2, S_q3, S_E, equations)
+    end
+    return SVector{nvariables(equations), real(equations)}(s)
+end
 end