Added new ES numerical fluxes

src/Trixi.jl

@@ -178,6 +178,7 @@ export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle,
        hydrostatic_reconstruction_chen_noelle, flux_nonconservative_chen_noelle,
        flux_hll_chen_noelle,
        FluxPlusDissipation, DissipationGlobalLaxFriedrichs, DissipationLocalLaxFriedrichs,
+       DissipationEntropyStableLump, DissipationEntropyStable,
        FluxLaxFriedrichs, max_abs_speed_naive,
        FluxHLL, min_max_speed_naive, min_max_speed_davis, min_max_speed_einfeldt,
        min_max_speed_chen_noelle,

src/equations/ideal_mhd_multiion_2d.jl

@@ -1043,4 +1043,272 @@ Computes the sum of the densities times the sum of the pressures
     end
     return rho_total * p_total
 end
+
+"""
+    DissipationEntropyStableLump(max_abs_speed=max_abs_speed_naive)
+
+Create a local Lax-Friedrichs dissipation operator where the maximum absolute wave speed
+is estimated as
+`max_abs_speed(u_ll, u_rr, orientation_or_normal_direction, equations)`,
+defaulting to [`max_abs_speed_naive`](@ref).
+"""
+struct DissipationEntropyStableLump{MaxAbsSpeed}
+    max_abs_speed::MaxAbsSpeed
+end
+
+DissipationEntropyStableLump() = DissipationEntropyStableLump(max_abs_speed_naive)
+
+@inline function (dissipation::DissipationEntropyStableLump)(u_ll, u_rr,
+                                                              orientation_or_normal_direction,
+                                                              equations::IdealMhdMultiIonEquations2D)
+    @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.5 * (B1_ll + B1_rr)
+    B2_avg = 0.5 * (B2_ll + B2_rr)
+    B3_avg = 0.5 * (B3_ll + B3_rr)
+    psi_avg = 0.5 * (psi_ll + psi_rr)
+
+    dissipation = zero(MVector{nvariables(equations), eltype(u_ll)})
+
+    beta_plus_ll = 0
+    beta_plus_rr = 0
+    # Get the lumped dissipation for 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.5 * rho_ll / p_ll
+        beta_rr = 0.5 * 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.5 * (rho_ll + rho_rr)
+        v1_avg = 0.5 * (v1_ll + v1_rr)
+        v2_avg = 0.5 * (v2_ll + v2_rr)
+        v3_avg = 0.5 * (v3_ll + v3_rr)
+        beta_avg = 0.5 * (beta_ll + beta_rr)
+        tau = 1 / (beta_ll + beta_rr) 
+        p_mean = 0.5 * rho_avg / beta_avg
+        p_star = 0.5 * rho_ln / beta_ln
+        vel_norm_avg = 0.5 * (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.5 * rho_ln * (2 * vel_avg_norm - vel_norm_avg)
+
+        h11 = rho_ln
+        h22 = rho_ln * v1_avg^2 + p_mean
+        h33 = rho_ln * v2_avg^2 + p_mean
+        h44 = rho_ln * v3_avg^2 + p_mean
+        h55 = ((p_star^2 / (gammas[k] - 1) + E_bar * E_bar) / rho_ln 
+              + vel_avg_norm * p_mean 
+              + tau * ( B1_avg^2 + B2_avg^2 + B3_avg^2 + psi_avg^2 ))
+
+        d1 = -0.5 * λ * h11 * (w1_rr - w1_ll)
+        d2 = -0.5 * λ * h22 * (w2_rr - w2_ll)
+        d3 = -0.5 * λ * h33 * (w3_rr - w3_ll)
+        d4 = -0.5 * λ * h44 * (w4_rr - w4_ll)
+        d5 = -0.5 * λ * 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
+
+    # Set the magnetic field and psi terms
+    h_B_psi = 1 / (beta_plus_ll + beta_plus_rr) 
+    dissipation[1] = -0.5 * λ * h_B_psi * (w_rr[1] - w_ll[1])
+    dissipation[2] = -0.5 * λ * h_B_psi * (w_rr[2] - w_ll[2])
+    dissipation[3] = -0.5 * λ * h_B_psi * (w_rr[3] - w_ll[3])
+    dissipation[end] = -0.5 * λ * h_B_psi * (w_rr[end] - w_ll[end])
+
+    return dissipation
+end
+
+function Base.show(io::IO, d::DissipationEntropyStableLump)
+    print(io, "DissipationEntropyStableLump(", d.max_abs_speed, ")")
+end
+
+"""
+DissipationEntropyStable(max_abs_speed=max_abs_speed_naive)
+
+Create a local Lax-Friedrichs dissipation operator where the maximum absolute wave speed
+is estimated as
+`max_abs_speed(u_ll, u_rr, orientation_or_normal_direction, equations)`,
+defaulting to [`max_abs_speed_naive`](@ref).
+"""
+struct DissipationEntropyStable{MaxAbsSpeed}
+    max_abs_speed::MaxAbsSpeed
+end
+
+DissipationEntropyStable() = DissipationEntropyStable(max_abs_speed_naive)
+
+@inline function (dissipation::DissipationEntropyStable)(u_ll, u_rr,
+                                                              orientation_or_normal_direction,
+                                                              equations::IdealMhdMultiIonEquations2D)
+    @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.5 * (B1_ll + B1_rr)
+    B2_avg = 0.5 * (B2_ll + B2_rr)
+    B3_avg = 0.5 * (B3_ll + B3_rr)
+    psi_avg = 0.5 * (psi_ll + psi_rr)
+
+    dissipation = zero(MVector{nvariables(equations), eltype(u_ll)})
+
+    beta_plus_ll = 0
+    beta_plus_rr = 0
+    # Get the lumped dissipation for 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.5 * rho_ll / p_ll
+        beta_rr = 0.5 * 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.5 * (rho_ll + rho_rr)
+        v1_avg = 0.5 * (v1_ll + v1_rr)
+        v2_avg = 0.5 * (v2_ll + v2_rr)
+        v3_avg = 0.5 * (v3_ll + v3_rr)
+        beta_avg = 0.5 * (beta_ll + beta_rr)
+        tau = 1 / (beta_ll + beta_rr) 
+        p_mean = 0.5 * rho_avg / beta_avg
+        p_star = 0.5 * rho_ln / beta_ln
+        vel_norm_avg = 0.5 * (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.5 * 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.5 * λ * (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.5 * λ * (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.5 * λ * (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.5 * λ * (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 
+              + tau * ( B1_avg^2 + B2_avg^2 + B3_avg^2 + psi_avg^2 ))
+        d5 = -0.5 * λ * (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
+
+    # Set the magnetic field and psi terms
+    h_B_psi = 1 / (beta_plus_ll + beta_plus_rr) 
+
+    # diagonal entries
+    dissipation[1] = -0.5 * λ * h_B_psi * (w_rr[1] - w_ll[1])
+    dissipation[2] = -0.5 * λ * h_B_psi * (w_rr[2] - w_ll[2])
+    dissipation[3] = -0.5 * λ * h_B_psi * (w_rr[3] - w_ll[3])
+    dissipation[end] = -0.5 * λ * h_B_psi * (w_rr[end] - w_ll[end])
+    # Off-diagonal entries
+    for k in eachcomponent(equations)
+        _, _, _, _, w5_ll = get_component(k, w_ll, equations)
+        _, _, _, _, w5_rr = get_component(k, w_rr, equations)
+
+        dissipation[1] -= 0.5 * λ * h_B_psi * B1_avg * (w5_rr - w5_ll)
+        dissipation[2] -= 0.5 * λ * h_B_psi * B2_avg * (w5_rr - w5_ll)
+        dissipation[3] -= 0.5 * λ * h_B_psi * B3_avg * (w5_rr - w5_ll)
+        dissipation[end] -= 0.5 * λ * h_B_psi * psi_avg * (w5_rr - w5_ll)
+
+        ind_E = 3 + (k - 1) * 5 + 5
+        dissipation[ind_E] -= 0.5 * λ * h_B_psi * B1_avg * (w_rr[1] - w_ll[1])
+        dissipation[ind_E] -= 0.5 * λ * h_B_psi * B2_avg * (w_rr[2] - w_ll[2])
+        dissipation[ind_E] -= 0.5 * λ * h_B_psi * B3_avg * (w_rr[3] - w_ll[3])
+        dissipation[ind_E] -= 0.5 * λ * h_B_psi * psi_avg * (w_rr[end] - w_ll[end])
+    end
+
+    return dissipation
+end
+
+function Base.show(io::IO, d::DissipationEntropyStable)
+    print(io, "DissipationEntropyStable(", d.max_abs_speed, ")")
+end
+
 end # @muladd