diff --git a/examples/tree_2d_dgsem/elixir_mhdmultiion_collisions.jl b/examples/tree_2d_dgsem/elixir_mhdmultiion_collisions.jl new file mode 100644 index 0000000000..17e7c0aed1 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_mhdmultiion_collisions.jl @@ -0,0 +1,152 @@ + +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 +# 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 = 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 + +# Return the constant electron temperature Te = 1 keV +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] + 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) + v11 = 0.65508770000000 + v21 = 0.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.0), + 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) +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 + save_everystep = false, callback = callbacks); + +summary_callback() # print the timer summary diff --git a/src/Trixi.jl b/src/Trixi.jl index f8cc3a1a4d..6cce44d995 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -217,7 +217,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 diff --git a/src/equations/ideal_glm_mhd_multiion.jl b/src/equations/ideal_glm_mhd_multiion.jl index 557cb8a12d..86d882cc65 100644 --- a/src/equations/ideal_glm_mhd_multiion.jl +++ b/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 + v_mag = sqrt(v1^2 + v2^2 + v3^2) + gamma = equations.gammas[k] + p[k] = (gamma - 1) * + (rho_e - 0.5f0 * rho * v_mag^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,172 @@ 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_{k'}\bar{\nu}_{kk'}(\vec{v}_{k'} - \vec{v}_k),\\ + s_{E_k} =& + 3 \sum_{k'} \left( + \bar{\nu}_{kk'} \frac{\rho_k M_1}{M_{k'} + M_k} R_1 (T_{k'} - T_k) + \right) + + \sum_{k'} \left( + \bar{\nu}_{kk'} \rho_k \frac{M_{k'}}{M_{k'} + M_k} \|\vec{v}_{k'} - \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}_{kk'}`` is the effective collision frequency of species `k` with species `k'`, which is computed as +```math +\begin{aligned} + \bar{\nu}_{kk'} = \bar{\nu}^1_{kk'} \tilde{B}_{kk'} \frac{\rho_{k'}}{T_{k k'}^{3/2}}, +\end{aligned} +``` +with the so-called reduced temperature ``T_{k k'}`` and the ion-ion collision constants ``\tilde{B}_{kk'}`` provided +in `equations.ion_electron_collision_constants` (see [`IdealGlmMhdMultiIonEquations2D`](@ref)). + +The additional coefficient ``\bar{\nu}^1_{kk'}`` 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. +- Schunk, R. W., Nagy, A. F. (2000). Ionospheres: Physics, plasma physics, and chemistry. + Cambridge university press. +""" +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}_{kk'}`` 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. +- Schunk, R. W., Nagy, A. F. (2000). Ionospheres: Physics, plasma physics, and chemistry. + Cambridge university press. +""" +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_e32 = 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_e32 + + 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 diff --git a/src/equations/ideal_glm_mhd_multiion_2d.jl b/src/equations/ideal_glm_mhd_multiion_2d.jl index fe452665ed..af14b84468 100644 --- a/src/equations/ideal_glm_mhd_multiion_2d.jl +++ b/src/equations/ideal_glm_mhd_multiion_2d.jl @@ -7,7 +7,17 @@ @doc raw""" IdealGlmMhdMultiIonEquations2D(; gammas, charge_to_mass, + gas_constants = zero(SVector{length(gammas), + eltype(gammas)}), + molar_masses = zero(SVector{length(gammas), + eltype(gammas)}), + ion_ion_collision_constants = zeros(eltype(gammas), + length(gammas), + length(gammas)), + ion_electron_collision_constants = zero(SVector{length(gammas), + eltype(gammas)}), electron_pressure = electron_pressure_zero, + electron_temperature = electron_pressure_zero, initial_c_h = NaN) The ideal compressible multi-ion MHD equations in two space dimensions augmented with a @@ -19,9 +29,41 @@ assumes that the equations are non-dimensionalized such, that the vacuum permeab In case of more than one ion species, the specific heat capacity ratios `gammas` and the charge-to-mass ratios `charge_to_mass` should be passed as tuples, e.g., `gammas=(1.4, 1.667)`. +The ion-ion and ion-electron collision source terms can be computed using the functions +[`source_terms_collision_ion_ion`](@ref) and [`source_terms_collision_ion_electron`](@ref), respectively. + +For ion-ion collision terms, the optional arguments `gas_constants`, `molar_masses`, and `ion_ion_collision_constants` +must be provided. For ion-electron collision terms, the optional arguments `gas_constants`, `molar_masses`, +`ion_electron_collision_constants`, and `electron_temperature` are required. + +- **`gas_constants`** and **`molar_masses`** are tuples containing the gas constant and molar mass of each + ion species, respectively. The **molar masses** can be provided in any unit system, as they are only used to + compute ratios and are independent of the other arguments. + +- **`ion_ion_collision_constants`** is a symmetric matrix that contains coefficients to compute the collision + frequencies between pairs of ion species. For example, `ion_ion_collision_constants[2, 3]` contains the collision + coefficient for collisions between the ion species 2 and the ion species 3. These constants are derived using the kinetic + theory of gases (see, e.g., *Schunk & Nagy, 2000*). They are related to the collision coefficients ``B_{st}`` listed + in Table 4.3 of *Schunk & Nagy (2000)*, but are scaled by the molecular mass of ion species ``t`` (i.e., + `ion_ion_collision_constants[2, 3] = ` ``B_{st}/m_{t}``) and must be provided in consistent physical units + (Schunk & Nagy use ``cm^3 K^{3/2} / s``). + See [`source_terms_collision_ion_ion`](@ref) for more details on how these constants are used to compute the collision + frequencies. + +- **`ion_electron_collision_constants`** is a tuple containing coefficients to compute the ion-electron collision frequency + for each ion species. They correspond to the collision coefficients `B_{se}` divided by the elementary charge. + The ion-electron collision frequencies can also be computed using the kinetic theory + of gases (see, e.g., *Schunk & Nagy, 2000*). See [`source_terms_collision_ion_electron`](@ref) for more details on how these + constants are used to compute the collision frequencies. + +- **`electron_temperature`** is a function with the signature `electron_temperature(u, equations)` that can be used + compute the electron temperature as a function of the state `u`. The electron temperature is relevant for the computation + of the ion-electron collision source terms. + The argument `electron_pressure` can be used to pass a function that computes the electron -pressure as a function of the state `u` with the signature `electron_pressure(u, equations::IdealGlmMhdMultiIonEquations2D)`. -By default, the electron pressure is zero. +pressure as a function of the state `u` with the signature `electron_pressure(u, equations)`. +The gradient of the electron pressure is relevant for the computation of the Lorentz flux +and non-conservative term. By default, the electron pressure is zero. The argument `initial_c_h` can be used to set the GLM divergence-cleaning speed. Note that `initial_c_h = 0` deactivates the divergence cleaning. The callback [`GlmSpeedCallback`](@ref) @@ -33,58 +75,121 @@ References: - A. Rueda-Ramírez, A. Sikstel, G. Gassner, An Entropy-Stable Discontinuous Galerkin Discretization of the Ideal Multi-Ion Magnetohydrodynamics System (2024). Journal of Computational Physics. [DOI: 10.1016/j.jcp.2024.113655](https://doi.org/10.1016/j.jcp.2024.113655). +- Schunk, R. W., & Nagy, A. F. (2000). Ionospheres: Physics, plasma physics, and chemistry. + Cambridge university press. !!! info "The multi-ion GLM-MHD equations require source terms" In case of more than one ion species, the multi-ion GLM-MHD equations should ALWAYS be used with [`source_terms_lorentz`](@ref). """ mutable struct IdealGlmMhdMultiIonEquations2D{NVARS, NCOMP, RealT <: Real, - ElectronPressure} <: + ElectronPressure, ElectronTemperature} <: AbstractIdealGlmMhdMultiIonEquations{2, NVARS, NCOMP} gammas::SVector{NCOMP, RealT} # Heat capacity ratios charge_to_mass::SVector{NCOMP, RealT} # Charge to mass ratios + gas_constants::SVector{NCOMP, RealT} # Specific gas constants + molar_masses::SVector{NCOMP, RealT} # Molar masses (can be provided in any units as they are only used to compute ratios) + ion_ion_collision_constants::Array{RealT, 2} # Symmetric matrix of collision frequency coefficients + ion_electron_collision_constants::SVector{NCOMP, RealT} # Constants for the ion-electron collision frequencies. The collision frequency is obtained as constant * (e * n_e) / T_e^1.5 electron_pressure::ElectronPressure # Function to compute the electron pressure + electron_temperature::ElectronTemperature # Function to compute the electron temperature c_h::RealT # GLM cleaning speed - function IdealGlmMhdMultiIonEquations2D{NVARS, NCOMP, RealT, - ElectronPressure}(gammas - ::SVector{NCOMP, RealT}, - charge_to_mass - ::SVector{NCOMP, RealT}, - electron_pressure - ::ElectronPressure, - c_h::RealT) where - {NVARS, NCOMP, RealT <: Real, ElectronPressure} + + function IdealGlmMhdMultiIonEquations2D{NVARS, NCOMP, RealT, ElectronPressure, + ElectronTemperature}(gammas + ::SVector{NCOMP, + RealT}, + charge_to_mass + ::SVector{NCOMP, + RealT}, + gas_constants + ::SVector{NCOMP, + RealT}, + molar_masses + ::SVector{NCOMP, + RealT}, + ion_ion_collision_constants + ::Array{RealT, 2}, + ion_electron_collision_constants + ::SVector{NCOMP, + RealT}, + electron_pressure + ::ElectronPressure, + electron_temperature + ::ElectronTemperature, + c_h::RealT) where + {NVARS, NCOMP, RealT <: Real, ElectronPressure, ElectronTemperature} NCOMP >= 1 || throw(DimensionMismatch("`gammas` and `charge_to_mass` have to be filled with at least one value")) - new(gammas, charge_to_mass, electron_pressure, c_h) + new(gammas, charge_to_mass, gas_constants, molar_masses, + ion_ion_collision_constants, + ion_electron_collision_constants, electron_pressure, electron_temperature, + c_h) end end function IdealGlmMhdMultiIonEquations2D(; gammas, charge_to_mass, + gas_constants = zero(SVector{length(gammas), + eltype(gammas)}), + molar_masses = zero(SVector{length(gammas), + eltype(gammas)}), + ion_ion_collision_constants = zeros(eltype(gammas), + length(gammas), + length(gammas)), + ion_electron_collision_constants = zero(SVector{length(gammas), + eltype(gammas)}), electron_pressure = electron_pressure_zero, + electron_temperature = electron_pressure_zero, initial_c_h = convert(eltype(gammas), NaN)) _gammas = promote(gammas...) _charge_to_mass = promote(charge_to_mass...) - RealT = promote_type(eltype(_gammas), eltype(_charge_to_mass)) + _gas_constants = promote(gas_constants...) + _molar_masses = promote(molar_masses...) + _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(ion_ion_collision_constants), + eltype(_ion_electron_collision_constants)) __gammas = SVector(map(RealT, _gammas)) __charge_to_mass = SVector(map(RealT, _charge_to_mass)) + __gas_constants = SVector(map(RealT, _gas_constants)) + __molar_masses = SVector(map(RealT, _molar_masses)) + __ion_ion_collision_constants = map(RealT, ion_ion_collision_constants) + __ion_electron_collision_constants = SVector(map(RealT, + _ion_electron_collision_constants)) NVARS = length(_gammas) * 5 + 4 NCOMP = length(_gammas) return IdealGlmMhdMultiIonEquations2D{NVARS, NCOMP, RealT, - typeof(electron_pressure)}(__gammas, - __charge_to_mass, - electron_pressure, - initial_c_h) + typeof(electron_pressure), + typeof(electron_temperature)}(__gammas, + __charge_to_mass, + __gas_constants, + __molar_masses, + __ion_ion_collision_constants, + __ion_electron_collision_constants, + electron_pressure, + electron_temperature, + initial_c_h) end # Outer constructor for `@reset` works correctly -function IdealGlmMhdMultiIonEquations2D(gammas, charge_to_mass, electron_pressure, c_h) +function IdealGlmMhdMultiIonEquations2D(gammas, charge_to_mass, gas_constants, + molar_masses, ion_ion_collision_constants, + ion_electron_collision_constants, + electron_pressure, + electron_temperature, + c_h) return IdealGlmMhdMultiIonEquations2D(gammas = gammas, charge_to_mass = charge_to_mass, + gas_constants = gas_constants, + molar_masses = molar_masses, + ion_ion_collision_constants = ion_ion_collision_constants, + ion_electron_collision_constants = ion_electron_collision_constants, electron_pressure = electron_pressure, + electron_temperature = electron_temperature, initial_c_h = c_h) end @@ -214,7 +319,7 @@ end orientation::Integer, equations::IdealGlmMhdMultiIonEquations2D) -Entropy-conserving non-conservative two-point "flux"" as described in +Entropy-conserving non-conservative two-point "flux" as described in - A. Rueda-Ramírez, A. Sikstel, G. Gassner, An Entropy-Stable Discontinuous Galerkin Discretization of the Ideal Multi-Ion Magnetohydrodynamics System (2024). Journal of Computational Physics. [DOI: 10.1016/j.jcp.2024.113655](https://doi.org/10.1016/j.jcp.2024.113655). diff --git a/test/test_tree_2d_mhdmultiion.jl b/test/test_tree_2d_mhdmultiion.jl index 8a9a261c30..fdf3c58658 100644 --- a/test/test_tree_2d_mhdmultiion.jl +++ b/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, + 0.0, + 0.0595534208484378, + 0.0, + 0.0, + 0.019718485574500753, + 0.0, + 0.059553420848437816, + 0.0, + 0.0, + 0.01738507024352939, + 0.0 + ], + linf=[ + 0.0, + 0.0, + 0.0, + 0.0, + 0.059553420848437816, + 0.0, + 0.0, + 0.019718485574500757, + 0.0, + 0.05955342084843786, + 0.0, + 0.0, + 0.017385070243529404, + 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 diff --git a/test/test_type.jl b/test/test_type.jl index 6f42b7600a..4efb3b2da1 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -1533,6 +1533,12 @@ isdir(outdir) && rm(outdir, recursive = true) @test eltype(@inferred source_terms_lorentz(u, x, t, equations)) == RealT + @test eltype(@inferred source_terms_collision_ion_ion(u, x, t, equations)) == + RealT + + @test eltype(@inferred source_terms_collision_ion_electron(u, x, t, equations)) == + RealT + for orientation in orientations @test eltype(@inferred flux(u, orientation, equations)) == RealT @test eltype(@inferred flux_nonconservative_ruedaramirez_etal(u_ll, u_rr,