Add a Lax-Friedrichs-type entropy-stable dissipation operator and its specialization for multi-ion MHD

examples/tree_2d_dgsem/elixir_mhdmultiion_ec.jl

@@ -9,8 +9,15 @@ equations = IdealGlmMhdMultiIonEquations2D(gammas = (1.4, 1.667),
 
 initial_condition = initial_condition_weak_blast_wave
 
+# Entropy conservative numerical fluxes
 volume_flux = (flux_ruedaramirez_etal, flux_nonconservative_ruedaramirez_etal)
 surface_flux = (flux_ruedaramirez_etal, flux_nonconservative_ruedaramirez_etal)
+# For provably entropy-stable surface fluxes, use
+# surface_flux = (FluxPlusDissipation(flux_ruedaramirez_etal, DissipationEntropyStable()), 
+#                 flux_nonconservative_ruedaramirez_etal)
+# For a standard local Lax-Friedrichs surface flux, use
+# surface_flux = (flux_lax_friedrichs, flux_nonconservative_central)
+
 solver = DGSEM(polydeg = 3, surface_flux = surface_flux,
                volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
 

src/Trixi.jl

@@ -188,6 +188,7 @@ export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle,
        flux_chan_etal, flux_nonconservative_chan_etal, flux_winters_etal,
        hydrostatic_reconstruction_audusse_etal, flux_nonconservative_audusse_etal,
        FluxPlusDissipation, DissipationGlobalLaxFriedrichs, DissipationLocalLaxFriedrichs,
+       DissipationEntropyStable,
        FluxLaxFriedrichs, max_abs_speed_naive,
        FluxHLL, min_max_speed_naive, min_max_speed_davis, min_max_speed_einfeldt,
        FluxLMARS,

src/equations/ideal_glm_mhd_multiion.jl

@@ -86,7 +86,7 @@ end
 
 """
     v1, v2, v3, vk1, vk2, vk3 = charge_averaged_velocities(u,
-                                                       equations::AbstractIdealGlmMhdMultiIonEquations)
+                                                           equations::AbstractIdealGlmMhdMultiIonEquations)
 
 
 Compute the charge-averaged velocities (`v1`, `v2`, and `v3`) and each ion species' contribution
@@ -263,4 +263,171 @@ end
 
     return SVector(cons)
 end
+
+# Specialization of DissipationEntropyStable for the multi-ion GLM-MHD equations
+# See
+# - 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).
+@inline function (dissipation::DissipationEntropyStable)(u_ll, u_rr,
+                                                         orientation_or_normal_direction,
+                                                         equations::AbstractIdealGlmMhdMultiIonEquations)
+    @unpack gammas = equations
+    λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction,
+                                  equations)
+
+    w_ll = cons2entropy(u_ll, equations)
+    w_rr = cons2entropy(u_rr, equations)
+    prim_ll = cons2prim(u_ll, equations)
+    prim_rr = cons2prim(u_rr, equations)
+    B1_ll, B2_ll, B3_ll = magnetic_field(u_ll, equations)
+    B1_rr, B2_rr, B3_rr = magnetic_field(u_rr, equations)
+    psi_ll = divergence_cleaning_field(u_ll, equations)
+    psi_rr = divergence_cleaning_field(u_rr, equations)
+
+    # Some global averages
+    B1_avg = 0.5f0 * (B1_ll + B1_rr)
+    B2_avg = 0.5f0 * (B2_ll + B2_rr)
+    B3_avg = 0.5f0 * (B3_ll + B3_rr)
+    psi_avg = 0.5f0 * (psi_ll + psi_rr)
+
+    dissipation = zero(MVector{nvariables(equations), eltype(u_ll)})
+
+    beta_plus_ll = 0
+    beta_plus_rr = 0
+    # Compute the dissipation for the hydrodynamic quantities of all components
+    for k in eachcomponent(equations)
+        rho_ll, v1_ll, v2_ll, v3_ll, p_ll = get_component(k, prim_ll, equations)
+        rho_rr, v1_rr, v2_rr, v3_rr, p_rr = get_component(k, prim_rr, equations)
+
+        w1_ll, w2_ll, w3_ll, w4_ll, w5_ll = get_component(k, w_ll, equations)
+        w1_rr, w2_rr, w3_rr, w4_rr, w5_rr = get_component(k, w_rr, equations)
+
+        # Auxiliary variables
+        beta_ll = 0.5f0 * rho_ll / p_ll
+        beta_rr = 0.5f0 * rho_rr / p_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
+
+        # Mean variables
+        rho_ln = ln_mean(rho_ll, rho_rr)
+        beta_ln = ln_mean(beta_ll, beta_rr)
+        rho_avg = 0.5f0 * (rho_ll + rho_rr)
+        v1_avg = 0.5f0 * (v1_ll + v1_rr)
+        v2_avg = 0.5f0 * (v2_ll + v2_rr)
+        v3_avg = 0.5f0 * (v3_ll + v3_rr)
+        beta_avg = 0.5f0 * (beta_ll + beta_rr)
+        tau = 1 / (beta_ll + beta_rr)
+        p_mean = 0.5f0 * rho_avg / beta_avg
+        p_star = 0.5f0 * rho_ln / beta_ln
+        vel_norm_avg = 0.5f0 * (vel_norm_ll + vel_norm_rr)
+        vel_avg_norm = v1_avg^2 + v2_avg^2 + v3_avg^2
+        E_bar = p_star / (gammas[k] - 1) +
+                0.5f0 * rho_ln * (2 * vel_avg_norm - vel_norm_avg)
+
+        h11 = rho_ln
+        h12 = rho_ln * v1_avg
+        h13 = rho_ln * v2_avg
+        h14 = rho_ln * v3_avg
+        h15 = E_bar
+        d1 = -0.5f0 * λ *
+             (h11 * (w1_rr - w1_ll) +
+              h12 * (w2_rr - w2_ll) +
+              h13 * (w3_rr - w3_ll) +
+              h14 * (w4_rr - w4_ll) +
+              h15 * (w5_rr - w5_ll))
+
+        h21 = h12
+        h22 = rho_ln * v1_avg^2 + p_mean
+        h23 = h21 * v2_avg
+        h24 = h21 * v3_avg
+        h25 = (E_bar + p_mean) * v1_avg
+        d2 = -0.5f0 * λ *
+             (h21 * (w1_rr - w1_ll) +
+              h22 * (w2_rr - w2_ll) +
+              h23 * (w3_rr - w3_ll) +
+              h24 * (w4_rr - w4_ll) +
+              h25 * (w5_rr - w5_ll))
+
+        h31 = h13
+        h32 = h23
+        h33 = rho_ln * v2_avg^2 + p_mean
+        h34 = h31 * v3_avg
+        h35 = (E_bar + p_mean) * v2_avg
+        d3 = -0.5f0 * λ *
+             (h31 * (w1_rr - w1_ll) +
+              h32 * (w2_rr - w2_ll) +
+              h33 * (w3_rr - w3_ll) +
+              h34 * (w4_rr - w4_ll) +
+              h35 * (w5_rr - w5_ll))
+
+        h41 = h14
+        h42 = h24
+        h43 = h34
+        h44 = rho_ln * v3_avg^2 + p_mean
+        h45 = (E_bar + p_mean) * v3_avg
+        d4 = -0.5f0 * λ *
+             (h41 * (w1_rr - w1_ll) +
+              h42 * (w2_rr - w2_ll) +
+              h43 * (w3_rr - w3_ll) +
+              h44 * (w4_rr - w4_ll) +
+              h45 * (w5_rr - w5_ll))
+
+        h51 = h15
+        h52 = h25
+        h53 = h35
+        h54 = h45
+        h55 = ((p_star^2 / (gammas[k] - 1) + E_bar * E_bar) / rho_ln
+               +
+               vel_avg_norm * p_mean)
+        d5 = -0.5f0 * λ *
+             (h51 * (w1_rr - w1_ll) +
+              h52 * (w2_rr - w2_ll) +
+              h53 * (w3_rr - w3_ll) +
+              h54 * (w4_rr - w4_ll) +
+              h55 * (w5_rr - w5_ll))
+
+        beta_plus_ll += beta_ll
+        beta_plus_rr += beta_rr
+
+        set_component!(dissipation, k, d1, d2, d3, d4, d5, equations)
+    end
+
+    # Compute the dissipation related to the magnetic and divergence-cleaning fields
+    h_B_psi = 1 / (beta_plus_ll + beta_plus_rr)
+
+    # Dissipation due to the diagonal entries of the dissipation matrix H
+    dissipation[1] = -0.5f0 * λ * h_B_psi * (w_rr[1] - w_ll[1])
+    dissipation[2] = -0.5f0 * λ * h_B_psi * (w_rr[2] - w_ll[2])
+    dissipation[3] = -0.5f0 * λ * h_B_psi * (w_rr[3] - w_ll[3])
+    dissipation[end] = -0.5f0 * λ * h_B_psi * (w_rr[end] - w_ll[end])
+    # Dissipation due to the off-diagonal entries of the dissipation matrix H
+    for k in eachcomponent(equations)
+        _, _, _, _, w5_ll = get_component(k, w_ll, equations)
+        _, _, _, _, w5_rr = get_component(k, w_rr, equations)
+
+        dissipation[1] -= 0.5f0 * λ * h_B_psi * B1_avg * (w5_rr - w5_ll)
+        dissipation[2] -= 0.5f0 * λ * h_B_psi * B2_avg * (w5_rr - w5_ll)
+        dissipation[3] -= 0.5f0 * λ * h_B_psi * B3_avg * (w5_rr - w5_ll)
+        dissipation[end] -= 0.5f0 * λ * h_B_psi * psi_avg * (w5_rr - w5_ll)
+
+        # Dissipation for the energy equation of species k depending on w_1, w_2, w_3 and w_end
+        ind_E = 3 + (k - 1) * 5 + 5
+        dissipation[ind_E] -= 0.5f0 * λ * h_B_psi * B1_avg * (w_rr[1] - w_ll[1])
+        dissipation[ind_E] -= 0.5f0 * λ * h_B_psi * B2_avg * (w_rr[2] - w_ll[2])
+        dissipation[ind_E] -= 0.5f0 * λ * h_B_psi * B3_avg * (w_rr[3] - w_ll[3])
+        dissipation[ind_E] -= 0.5f0 * λ * h_B_psi * psi_avg * (w_rr[end] - w_ll[end])
+
+        # Dissipation for the energy equation of all ion species depending on w_5
+        for kk in eachcomponent(equations)
+            ind_E = 3 + (kk - 1) * 5 + 5
+            dissipation[ind_E] -= 0.5f0 * λ *
+                                  (h_B_psi *
+                                   (B1_avg^2 + B2_avg^2 + B3_avg^2 + psi_avg^2)) *
+                                  (w5_rr - w5_ll)
+        end
+    end
+
+    return dissipation
+end
 end

src/equations/ideal_glm_mhd_multiion_2d.jl

@@ -43,8 +43,8 @@ mutable struct IdealGlmMhdMultiIonEquations2D{NVARS, NCOMP, RealT <: Real,
                AbstractIdealGlmMhdMultiIonEquations{2, NVARS, NCOMP}
     gammas::SVector{NCOMP, RealT} # Heat capacity ratios
     charge_to_mass::SVector{NCOMP, RealT} # Charge to mass ratios
-    electron_pressure::ElectronPressure       # Function to compute the electron pressure
-    c_h::RealT                 # GLM cleaning speed
+    electron_pressure::ElectronPressure # Function to compute the electron pressure
+    c_h::RealT # GLM cleaning speed
     function IdealGlmMhdMultiIonEquations2D{NVARS, NCOMP, RealT,
                                             ElectronPressure}(gammas
                                                               ::SVector{NCOMP, RealT},
@@ -110,7 +110,7 @@ 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 = real(equations)
+    RealT = eltype(x)
     inicenter = (0, 0)
     x_norm = x[1] - inicenter[1]
     y_norm = x[2] - inicenter[2]

src/equations/numerical_fluxes.jl

@@ -221,6 +221,46 @@ See [`FluxLaxFriedrichs`](@ref).
 """
 const flux_lax_friedrichs = FluxLaxFriedrichs()
 
+"""
+    DissipationEntropyStable(max_abs_speed=max_abs_speed_naive)
+
+Create a local Lax-Friedrichs-type dissipation operator that is provably entropy stable. This operator
+must be used together with an entropy-conservative two-point flux function (e.g., `flux_ec`) to yield 
+an entropy-stable surface flux. The surface flux function can be initialized as:
+```
+flux_es = FluxPlusDissipation(flux_ec, DissipationEntropyStable())
+```
+
+In particular, the numerical flux has the form
+```math
+f^{ES} = f^{EC} + \frac{1}{2} λ_{max} H (w_r - w_l),
+````
+where ``f^{EC}`` is the entropy-conservative two-point flux function (computed with, e.g., `flux_ec`), ``λ_{max}`` 
+is the maximum wave speed estimated as `max_abs_speed(u_l, u_r, orientation_or_normal_direction, equations)`,
+defaulting to [`max_abs_speed_naive`](@ref), ``H`` is a symmetric positive-definite dissipation matrix that
+depends on the left and right states `u_l` and `u_r`, and ``(w_r - w_l)`` is the jump in entropy variables.
+Ideally, ``H (w_r - w_l) = (u_r - u_l)``, such that the dissipation operator is consistent with the local
+Lax-Friedrichs dissipation.
+
+The entropy-stable dissipation operator is computed with the function
+`function (dissipation::DissipationEntropyStable)(u_l, u_r, orientation_or_normal_direction, equations)`,
+which must be specialized for each equation.
+
+For the derivation of the dissipation matrix for the multi-ion GLM-MHD equations, see:
+- 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).
+"""
+struct DissipationEntropyStable{MaxAbsSpeed}
+    max_abs_speed::MaxAbsSpeed
+end
+
+DissipationEntropyStable() = DissipationEntropyStable(max_abs_speed_naive)
+
+function Base.show(io::IO, d::DissipationEntropyStable)
+    print(io, "DissipationEntropyStable(", d.max_abs_speed, ")")
+end
+
 """
     FluxHLL(min_max_speed=min_max_speed_davis)
 

test/test_tree_2d_mhdmultiion.jl

@@ -6,7 +6,7 @@ using Trixi
 include("test_trixi.jl")
 
 # pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
-EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_2d_dgsem")
+EXAMPLES_DIR = joinpath(examples_dir(), "tree_2d_dgsem")
 
 @testset "MHD Multi-ion" begin
 #! format: noindent
@@ -55,6 +55,53 @@ EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_2
     end
 end
 
+@trixi_testset "Provably entropy-stable LLF-type fluxes for multi-ion GLM-MHD" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmultiion_ec.jl"),
+                        l2=[
+                            0.017668017558288736,
+                            0.01779783612885502,
+                            0.027841673842076285,
+                            0.015603429086471762,
+                            0.017849042999817964,
+                            0.01814196379994667,
+                            0.005478212889809162,
+                            0.20585517887094282,
+                            0.021301245733548135,
+                            0.03018506565829777,
+                            0.02938517728342881,
+                            0.01837279433780041,
+                            0.11810307914710033,
+                            0.0002962677911603057
+                        ],
+                        linf=[
+                            0.06594754030722516,
+                            0.06587779693691242,
+                            0.09451464686853495,
+                            0.06787230638663028,
+                            0.08910065803824378,
+                            0.08828064474448032,
+                            0.023647579422062297,
+                            0.8059383650828509,
+                            0.1224367642558366,
+                            0.15930418161523857,
+                            0.15382860284948224,
+                            0.08695364286964764,
+                            0.4949375933243716,
+                            0.003287251595115295
+                        ],
+                        surface_flux=(FluxPlusDissipation(flux_ruedaramirez_etal,
+                                                          DissipationEntropyStable()),
+                                      flux_nonconservative_ruedaramirez_etal))
+    # 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
+
 @trixi_testset "elixir_mhdmultiion_ec.jl with local Lax-Friedrichs at the surface" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmultiion_ec.jl"),
                         l2=[

test/test_type.jl

@@ -1524,6 +1524,7 @@ isdir(outdir) && rm(outdir, recursive = true)
                                              one(RealT),
                                              one(RealT),
                                              one(RealT))
+            dissipation_es = DissipationEntropyStable()
             orientations = [1, 2]
 
             @test eltype(@inferred initial_condition_weak_blast_wave(x, t, equations)) ==
@@ -1547,6 +1548,8 @@ isdir(outdir) && rm(outdir, recursive = true)
                 @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation,
                                                            equations)) ==
                       RealT
+                @test eltype(@inferred dissipation_es(u_ll, u_rr, orientation, equations)) ==
+                      RealT
             end
 
             @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT