Add 2D multi-ion GLM-MHD equations for the TreeMesh solver

examples/tree_2d_dgsem/elixir_mhdmultiion_ec.jl

@@ -3,7 +3,7 @@ using OrdinaryDiffEq
 using Trixi
 
 ###############################################################################
-# semidiscretization of the ideal MHD equations
+# semidiscretization of the ideal multi-ion MHD equations
 equations = IdealGlmMhdMultiIonEquations2D(gammas = (1.4, 1.667),
                                            charge_to_mass = (1.0, 2.0))
 

src/Trixi.jl

@@ -232,7 +232,6 @@ export entropy, energy_total, energy_kinetic, energy_internal, energy_magnetic,
        enstrophy, magnetic_field, divergence_cleaning_field
 export lake_at_rest_error
 export ncomponents, eachcomponent
-export get_component
 
 export TreeMesh, StructuredMesh, StructuredMeshView, UnstructuredMesh2D, P4estMesh,
        T8codeMesh

src/equations/ideal_glm_mhd_multiion.jl

@@ -38,11 +38,12 @@ function default_analysis_integrals(::AbstractIdealGlmMhdMultiIonEquations)
 end
 
 """
-         source_terms_lorentz(u, x, t, equations::AbstractIdealGlmMhdMultiIonEquations)
+    source_terms_lorentz(u, x, t, equations::AbstractIdealGlmMhdMultiIonEquations)
 
 Source terms due to the Lorentz' force for plasmas with more than one ion species. These source 
 terms are a fundamental, inseparable part of the multi-ion GLM-MHD equations, and vanish for 
-a single-species plasma.
+a single-species plasma. In particular, they have to be used for every
+simulation of [`IdealGlmMhdMultiIonEquations2D`](@ref).
 """
 function source_terms_lorentz(u, x, t, equations::AbstractIdealGlmMhdMultiIonEquations)
     @unpack charge_to_mass = equations
@@ -54,9 +55,10 @@ function source_terms_lorentz(u, x, t, equations::AbstractIdealGlmMhdMultiIonEqu
 
     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
+        rho_inv = 1 / rho
+        v1 = rho_v1 * rho_inv
+        v2 = rho_v2 * rho_inv
+        v3 = rho_v3 * rho_inv
         v1_diff = v1_plus - v1
         v2_diff = v2_plus - v2
         v3_diff = v3_plus - v3
@@ -73,7 +75,7 @@ function source_terms_lorentz(u, x, t, equations::AbstractIdealGlmMhdMultiIonEqu
 end
 
 """
-         electron_pressure_zero(u, equations::AbstractIdealGlmMhdMultiIonEquations)
+    electron_pressure_zero(u, equations::AbstractIdealGlmMhdMultiIonEquations)
 
 Returns the value of zero for the electron pressure. Needed for consistency with the 
 single-fluid MHD equations in the limit of one ion species.
@@ -83,13 +85,16 @@ function electron_pressure_zero(u, equations::AbstractIdealGlmMhdMultiIonEquatio
 end
 
 """
-v1, v2, v3, vk1, vk2, vk3 = charge_averaged_velocities(u,
+    v1, v2, v3, vk1, vk2, vk3 = charge_averaged_velocities(u,
                                                        equations::AbstractIdealGlmMhdMultiIonEquations)
 
 
-Compute the charge-averaged velocities (v1, v2, and v3) and each ion species' contribution
-to the charge-averaged velocities (vk1, vk2, and vk3). The output variables vk1, vk2, and vk3
-are SVectors of size ncomponents(equations).
+Compute the charge-averaged velocities (`v1`, `v2`, and `v3`) and each ion species' contribution
+to the charge-averaged velocities (`vk1`, `vk2`, and `vk3`). The output variables `vk1`, `vk2`, and `vk3`
+are `SVectors` of size `ncomponents(equations)`.
+
+!!! warning "Experimental implementation"
+    This is an experimental feature and may change in future releases.
 """
 @inline function charge_averaged_velocities(u,
                                             equations::AbstractIdealGlmMhdMultiIonEquations)
@@ -121,7 +126,10 @@ end
 """
     get_component(k, u, equations::AbstractIdealGlmMhdMultiIonEquations)
 
-Get the hydrodynamic variables of component (ion species) k.
+Get the hydrodynamic variables of component (ion species) `k`.
+
+!!! warning "Experimental implementation"
+    This is an experimental feature and may change in future releases.
 """
 @inline function get_component(k, u, equations::AbstractIdealGlmMhdMultiIonEquations)
     return SVector(u[3 + (k - 1) * 5 + 1],
@@ -135,7 +143,10 @@ end
     set_component!(u, k, u1, u2, u3, u4, u5,
                    equations::AbstractIdealGlmMhdMultiIonEquations)
 
-Set the hydrodynamic variables of component (ion species) k.
+Set the hydrodynamic variables (`u1` to `u5`) of component (ion species) `k`.
+
+!!! warning "Experimental implementation"
+    This is an experimental feature and may change in future releases.
 """
 @inline function set_component!(u, k, u1, u2, u3, u4, u5,
                                 equations::AbstractIdealGlmMhdMultiIonEquations)
@@ -176,10 +187,10 @@ function cons2prim(u, equations::AbstractIdealGlmMhdMultiIonEquations)
     prim[3] = B3
     for k in eachcomponent(equations)
         rho, rho_v1, rho_v2, rho_v3, rho_e = get_component(k, u, equations)
-        srho = 1 / rho
-        v1 = srho * rho_v1
-        v2 = srho * rho_v2
-        v3 = srho * rho_v3
+        rho_inv = 1 / rho
+        v1 = rho_inv * rho_v1
+        v2 = rho_inv * rho_v2
+        v3 = rho_inv * rho_v3
 
         p = (gammas[k] - 1) * (rho_e -
              0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3
@@ -241,7 +252,7 @@ end
         rho_v2 = rho * v2
         rho_v3 = rho * v3
 
-        rho_e = p / (gammas[k] - 1.0) +
+        rho_e = p / (gammas[k] - 1) +
                 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) +
                 0.5f0 * (B1^2 + B2^2 + B3^2) +
                 0.5f0 * psi^2

src/equations/ideal_glm_mhd_multiion_2d.jl

@@ -34,8 +34,9 @@ References:
   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).
 
-!!! ATTENTION: In case of more than one ion species, these equations should ALWAYS be used
-    with `source_terms_lorentz`.
+!!! 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} <:
@@ -66,13 +67,12 @@ function IdealGlmMhdMultiIonEquations2D(; gammas, charge_to_mass,
     _gammas = promote(gammas...)
     _charge_to_mass = promote(charge_to_mass...)
     RealT = promote_type(eltype(_gammas), eltype(_charge_to_mass))
+    __gammas = SVector(map(RealT, _gammas))
+    __charge_to_mass = SVector(map(RealT, _charge_to_mass))
 
     NVARS = length(_gammas) * 5 + 4
     NCOMP = length(_gammas)
 
-    __gammas = SVector(map(RealT, _gammas))
-    __charge_to_mass = SVector(map(RealT, _charge_to_mass))
-
     return IdealGlmMhdMultiIonEquations2D{NVARS, NCOMP, RealT,
                                           typeof(electron_pressure)}(__gammas,
                                                                      __charge_to_mass,
@@ -110,27 +110,23 @@ function initial_condition_weak_blast_wave(x, t,
     # Adapted MHD version of the weak blast wave from Hennemann & Gassner JCP paper 2020 (Sec. 6.3)
     # Same discontinuity in the velocities but with magnetic fields
     # Set up polar coordinates
+    RealT = eltype(x)
     inicenter = (0, 0)
     x_norm = x[1] - inicenter[1]
     y_norm = x[2] - inicenter[2]
     r = sqrt(x_norm^2 + y_norm^2)
     phi = atan(y_norm, x_norm)
 
     # Calculate primitive variables
-    rho = zero(real(equations))
-    if r > 0.5
-        rho = 1.0
-    else
-        rho = 1.1691
-    end
-    v1 = r > 0.5 ? 0.0 : 0.1882 * cos(phi)
-    v2 = r > 0.5 ? 0.0 : 0.1882 * sin(phi)
-    p = r > 0.5 ? 1.0 : 1.245
+    rho = r > 0.5f0 ? one(RealT) : convert(RealT, 1.1691)
+    v1 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos(phi)
+    v2 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * sin(phi)
+    p = r > 0.5f0 ? one(RealT) : convert(RealT, 1.245)
 
     prim = zero(MVector{nvariables(equations), real(equations)})
-    prim[1] = 1.0
-    prim[2] = 1.0
-    prim[3] = 1.0
+    prim[1] = 1
+    prim[2] = 1
+    prim[3] = 1
     for k in eachcomponent(equations)
         set_component!(prim, k,
                        2^(k - 1) * (1 - 2) / (1 - 2^ncomponents(equations)) * rho, v1,
@@ -161,9 +157,10 @@ end
 
         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
+            rho_inv = 1 / rho
+            v1 = rho_v1 * rho_inv
+            v2 = rho_v2 * rho_inv
+            v3 = rho_v3 * rho_inv
             kin_en = 0.5f0 * rho * (v1^2 + v2^2 + v3^2)
 
             gamma = equations.gammas[k]
@@ -187,9 +184,10 @@ end
 
         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
+            rho_inv = 1 / rho
+            v1 = rho_v1 * rho_inv
+            v2 = rho_v2 * rho_inv
+            v3 = rho_v3 * rho_inv
             kin_en = 0.5f0 * rho * (v1^2 + v2^2 + v3^2)
 
             gamma = equations.gammas[k]
@@ -212,7 +210,7 @@ end
 end
 
 """
-        flux_nonconservative_ruedaramirez_etal(u_ll, u_rr,
+    flux_nonconservative_ruedaramirez_etal(u_ll, u_rr,
                                                orientation::Integer,
                                                equations::IdealGlmMhdMultiIonEquations2D)
 
@@ -221,11 +219,12 @@ Entropy-conserving non-conservative two-point "flux"" as described in
   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).
 
-ATTENTION: The non-conservative fluxes derived in the reference above are written as the product
-           of local and symmetric parts and are meant to be used in the same way as the conservative
-           fluxes (i.e., flux + flux_noncons in both volume and surface integrals). In this routine, 
-           the fluxes are multiplied by 2 because the non-conservative fluxes are always multiplied 
-           by 0.5 whenever they are used in the Trixi code
+!!! info "Usage and Scaling of Non-Conservative Fluxes in Trixi.jl"
+    The non-conservative fluxes derived in the reference above are written as the product
+    of local and symmetric parts and are meant to be used in the same way as the conservative
+    fluxes (i.e., flux + flux_noncons in both volume and surface integrals). In this routine, 
+    the fluxes are multiplied by 2 because the non-conservative fluxes are always multiplied 
+    by 0.5 whenever they are used in the Trixi code.
 
 The term is composed of four individual non-conservative terms:
 1. The Godunov-Powell term, which arises for plasmas with non-vanishing magnetic field divergence, and
@@ -375,19 +374,21 @@ The term is composed of four individual non-conservative terms:
 end
 
 """
-        flux_nonconservative_central(u_ll, u_rr, orientation::Integer,
-                                     equations::IdealGlmMhdMultiIonEquations2D)
+    flux_nonconservative_central(u_ll, u_rr, orientation::Integer,
+                                 equations::IdealGlmMhdMultiIonEquations2D)
 
 Central non-conservative two-point "flux", where the symmetric parts are computed with standard averages.
-The use of this term together with flux_central with VolumeIntegralFluxDifferencing yields a "standard"
+The use of this term together with [`flux_central`](@ref) 
+with [`VolumeIntegralFluxDifferencing`](@ref) yields a "standard"
 (weak-form) DGSEM discretization of the multi-ion GLM-MHD system. This flux can also be used to construct a
 standard local Lax-Friedrichs flux using `surface_flux = (flux_lax_friedrichs, flux_nonconservative_central)`.
 
-ATTENTION: The central non-conservative fluxes are written as the product
-           of local and symmetric parts and are meant to be used in the same way as the conservative
-           fluxes (i.e., flux + flux_noncons in both volume and surface integrals). In this routine, 
-           we omit the 0.5 when computing averages because the non-conservative flux is multiplied by 
-           0.5 whenever it's used in the Trixi code
+!!! info "Usage and Scaling of Non-Conservative Fluxes in Trixi"
+    The central non-conservative fluxes implemented in this function are written as the product
+    of local and symmetric parts, where the symmetric part is a standard average. These fluxes
+    are meant to be used in the same way as the conservative fluxes (i.e., flux + flux_noncons 
+    in both volume and surface integrals). In this routine, the fluxes are multiplied by 2 because 
+    the non-conservative fluxes are always multiplied by 0.5 whenever they are used in the Trixi code.
 
 The term is composed of four individual non-conservative terms:
 1. The Godunov-Powell term, which arises for plasmas with non-vanishing magnetic field divergence, and
@@ -520,7 +521,7 @@ The term is composed of four individual non-conservative terms:
 end
 
 """
-flux_ruedaramirez_etal(u_ll, u_rr, orientation, equations::IdealGlmMhdMultiIonEquations2D)
+    flux_ruedaramirez_etal(u_ll, u_rr, orientation, equations::IdealGlmMhdMultiIonEquations2D)
 
 Entropy conserving two-point flux for the multi-ion GLM-MHD equations from
 - A. Rueda-Ramírez, A. Sikstel, G. Gassner, An Entropy-Stable Discontinuous Galerkin Discretization
@@ -582,13 +583,14 @@ function flux_ruedaramirez_etal(u_ll, u_rr, orientation::Integer,
                                                                               equations)
             rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr = get_component(k, u_rr,
                                                                               equations)
-
-            v1_ll = rho_v1_ll / rho_ll
-            v2_ll = rho_v2_ll / rho_ll
-            v3_ll = rho_v3_ll / rho_ll
-            v1_rr = rho_v1_rr / rho_rr
-            v2_rr = rho_v2_rr / rho_rr
-            v3_rr = rho_v3_rr / rho_rr
+            rho_inv_ll = 1 / rho_ll
+            v1_ll = rho_v1_ll * rho_inv_ll
+            v2_ll = rho_v2_ll * rho_inv_ll
+            v3_ll = rho_v3_ll * rho_inv_ll
+            rho_inv_rr = 1 / rho_rr
+            v1_rr = rho_v1_rr * rho_inv_rr
+            v2_rr = rho_v2_rr * rho_inv_rr
+            v3_rr = rho_v3_rr * rho_inv_rr
             vel_norm_ll = v1_ll^2 + v2_ll^2 + v3_ll^2
             vel_norm_rr = v1_rr^2 + v2_rr^2 + v3_rr^2
 
@@ -680,12 +682,14 @@ function flux_ruedaramirez_etal(u_ll, u_rr, orientation::Integer,
             rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr = get_component(k, u_rr,
                                                                               equations)
 
-            v1_ll = rho_v1_ll / rho_ll
-            v2_ll = rho_v2_ll / rho_ll
-            v3_ll = rho_v3_ll / rho_ll
-            v1_rr = rho_v1_rr / rho_rr
-            v2_rr = rho_v2_rr / rho_rr
-            v3_rr = rho_v3_rr / rho_rr
+            rho_inv_ll = 1 / rho_ll
+            v1_ll = rho_v1_ll * rho_inv_ll
+            v2_ll = rho_v2_ll * rho_inv_ll
+            v3_ll = rho_v3_ll * rho_inv_ll
+            rho_inv_rr = 1 / rho_rr
+            v1_rr = rho_v1_rr * rho_inv_rr
+            v2_rr = rho_v2_rr * rho_inv_rr
+            v3_rr = rho_v3_rr * rho_inv_rr
             vel_norm_ll = v1_ll^2 + v2_ll^2 + v3_ll^2
             vel_norm_rr = v1_rr^2 + v2_rr^2 + v3_rr^2
 
@@ -820,20 +824,21 @@ end
     for k in eachcomponent(equations)
         rho, rho_v1, rho_v2, rho_v3, rho_e = get_component(k, cons, equations)
 
-        v1 = rho_v1 / rho
-        v2 = rho_v2 / rho
-        v3 = rho_v3 / rho
+        rho_inv = 1 / rho
+        v1 = rho_v1 * rho_inv
+        v2 = rho_v2 * rho_inv
+        v3 = rho_v3 * rho_inv
         v_mag = sqrt(v1^2 + v2^2 + v3^2)
         gamma = equations.gammas[k]
         p = (gamma - 1) *
             (rho_e - 0.5f0 * rho * v_mag^2 - 0.5f0 * (B1^2 + B2^2 + B3^2) -
              0.5f0 * psi^2)
-        a_square = gamma * p / rho
-        sqrt_rho = sqrt(rho)
+        a_square = gamma * p * rho_inv
+        inv_sqrt_rho = 1 / sqrt(rho)
 
-        b1 = B1 / sqrt_rho
-        b2 = B2 / sqrt_rho
-        b3 = B3 / sqrt_rho
+        b1 = B1 * inv_sqrt_rho
+        b2 = B2 * inv_sqrt_rho
+        b3 = B3 * inv_sqrt_rho
         b_square = b1^2 + b2^2 + b3^2
 
         if orientation == 1

test/test_tree_2d_mhdmultiion.jl

@@ -1,4 +1,4 @@
-module TestExamples2DEulerMulticomponent
+module TestExamples2DIdealGlmMhdMultiIon
 
 using Test
 using Trixi