Fixed problem with ion-electron collisions and added test for collisi…

…on terms (we still have some allocations problems)

examples/tree_2d_dgsem/elixir_mhdmultiion_collisions.jl

@@ -0,0 +1,149 @@
+
+using OrdinaryDiffEq
+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.
+#
+# This is effectively a zero-dimensional case because the spatial gradients are zero, and we use it to test the
+# collission 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
+# Characteristic density:	rho_inf = 1.99E+00	kg/m^3
+# Characteristic velocity: V_inf = 1.00E+06	m/s
+# Characteristic vacuum permeability: mu0_inf = 1.26E-06	N/A^2
+# Characteristic gas constant: R_inf = 6.92237E+02	J/kg/K
+#
+# 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 = 0.5 keV
+function electron_pressure_constantTe(u, equations::IdealGlmMhdMultiIonEquations2D)
+    @unpack charge_to_mass = equations
+    Te = 0.008029953773 # [nondimensional] = 1 [keV]
+    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
+    kB_e = 7.86319034E-02 #[nondimensional]
+
+    return total_electron_charge * kB_e * Te
+end
+
+function electron_temperature_constantTe(u, equations::IdealGlmMhdMultiIonEquations2D)
+    return 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]
+                                           collision_frequency = 0.4079382480442680, # [nondimensional]
+                                           ion_electron_collision_constants = (8.56368379833E-06,
+                                                                               8.56368379833E-06), # [nondimensional]
+                                           electron_pressure = electron_pressure_constantTe,
+                                           electron_temperature = electron_temperature_constantTe,
+                                           initial_c_h = 0.0) # Deactivate GLM divergence cleaning
+
+"""
+
+"""
+function initial_condition_thermal_equilibration(x, t,
+                                                 equations::IdealGlmMhdMultiIonEquations2D)
+    v11 = 0.65508770000000
+    v21 = 0
+    v2 = v3 = 0.0
+    B1 = B2 = B3 = 0.0
+    rho1 = 0.1
+    rho2 = 1.0
+
+    p1 = 0.00040170535986
+    p2 = 0.00401705359856
+
+    return prim2cons(SVector(B1, B2, B3, rho1, v11, v2, v3, p1, rho2, v21, v2, v3, p2, 0),
+                     equations)
+end
+
+function temperature1(cons, equations::IdealGlmMhdMultiIonEquations2D)
+    prim = cons2prim(cons, equations)
+    rho, _, _, _, p = Trixi.get_component(1, prim, equations)
+
+    return p / rho / equations.gas_constants[1]
+end
+function temperature2(cons, equations::IdealGlmMhdMultiIonEquations2D)
+    prim = cons2prim(cons, equations)
+    rho, _, _, _, p = Trixi.get_component(2, prim, equations)
+
+    return p / rho / equations.gas_constants[2]
+end
+
+initial_condition = initial_condition_thermal_equilibration
+tspan = (0.0, 0.1) # 100 [ps]
+
+# Entropy conservative numerical fluxes
+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)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 1,
+                n_cells_max = 1_000_000)
+
+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
+
+# In this particular case, we can omit source_terms_lorentz because the magnetic field is zero!
+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 = Trixi.solve(ode, Trixi.SimpleSSPRK33();
+                  dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+                  save_everystep = false, callback = callbacks);
+
+summary_callback() # print the timer summary

src/equations/ideal_glm_mhd_multiion.jl

@@ -478,8 +478,6 @@ Ion-ion collision source terms, cf. Rueda-Ramirez et al. (2023) and Rubin et al.
 """
 function source_terms_collision_ion_ion(u, x, t,
                                         equations::AbstractIdealGlmMhdMultiIonEquations)
-    S_std = source_terms_standard(u, x, t, equations)
-
     s = zero(MVector{nvariables(equations), eltype(u)})
     @unpack gammas, gas_constants, molar_masses, collision_frequency = equations
 
@@ -530,7 +528,7 @@ function source_terms_collision_ion_ion(u, x, t,
 
         set_component!(s, k, 0, S_q1, S_q2, S_q3, S_E, equations)
     end
-    return SVector{nvariables(equations), real(equations)}(S_std .+ s)
+    return SVector{nvariables(equations), real(equations)}(s)
 end
 
 """
@@ -549,8 +547,15 @@ function source_terms_collision_ion_electron(u, x, t,
     T_e = electron_temperature(u, equations)
     T_e32 = T_e^(3 / 2)
 
-    v1_plus, v2_plus, v3_plus, vk1_plus, vk2_plus, vk3_plus, total_electron_charge = charge_averaged_velocities(u,
-                                                                                                                equations)
+    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)
@@ -589,8 +594,15 @@ function source_terms_collision_ion_electron_ohm(u, x, t,
     T_e = electron_temperature(u, equations)
     T_e32 = T_e^(3 / 2)
 
-    v1_plus, v2_plus, v3_plus, vk1_plus, vk2_plus, vk3_plus, total_electron_charge = charge_averaged_velocities(u,
-                                                                                                                equations)
+    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
 
     Se_q1 = 0
     Se_q2 = 0

src/equations/ideal_glm_mhd_multiion_2d.jl

@@ -99,7 +99,7 @@ function IdealGlmMhdMultiIonEquations2D(; gammas, charge_to_mass,
     _ion_electron_collision_constants = promote(ion_electron_collision_constants...)
     RealT = promote_type(eltype(_gammas), eltype(_charge_to_mass),
                          eltype(_gas_constants), eltype(_molar_masses),
-                         eltype(collision_frequency),
+                         typeof(collision_frequency),
                          eltype(_ion_electron_collision_constants))
     __gammas = SVector(map(RealT, _gammas))
     __charge_to_mass = SVector(map(RealT, _charge_to_mass))

test/test_tree_2d_mhdmultiion.jl

@@ -146,6 +146,50 @@ end
         @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
     end
 end
+
+@trixi_testset "Multi-ion GLM-MHD collision source terms" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmultiion_collisions.jl"),
+                        l2=[
+                            0.0,
+                            0.0,
+                            0.0,
+                            1.6829845497998036e-17,
+                            0.059553423755556875,
+                            5.041277811626498e-19,
+                            0.0,
+                            0.01971848646448756,
+                            7.144325530681256e-17,
+                            0.059553423755556924,
+                            5.410518863721139e-18,
+                            0.0,
+                            0.017385071119051767,
+                            0.0
+                        ],
+                        linf=[
+                            0.0,
+                            0.0,
+                            0.0,
+                            4.163336342344337e-17,
+                            0.05955342375555689,
+                            1.609474550734496e-18,
+                            0.0,
+                            0.019718486464487567,
+                            2.220446049250313e-16,
+                            0.059553423755556965,
+                            1.742701984446815e-17,
+                            0.0,
+                            0.01738507111905178,
+                            0.0
+                        ])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+        t = sol.t[end]
+        u_ode = sol.u[end]
+        du_ode = similar(u_ode)
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+    end
+end
 end
 
 end # module