HLLC Flux for Ideal MHD 1D

src/equations/ideal_glm_mhd_1d.jl

@@ -259,6 +259,128 @@ Hindenlang and Gassner (2019), extending [`flux_ranocha`](@ref) to the MHD equat
     return SVector(f1, f2, f3, f4, f5, f6, f7, f8)
 end
 
+"""
+    flux_hllc(u_ll, u_rr, orientation, equations::IdealGlmMhdEquations1D)
+
+- Li (2005)
+An HLLC Riemann solver for magneto-hydrodynamics
+[DOI: 10.1016/j.jcp.2004.08.020](https://doi.org/10.1016/j.jcp.2004.08.020).
+"""
+function flux_hllc(u_ll, u_rr, orientation::Integer,
+                   equations::IdealGlmMhdEquations1D)
+    # Unpack left and right states
+    rho_ll, v1_ll, v2_ll, v3_ll, p_ll, B1_ll, B2_ll, B3_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_rr, v3_rr, p_rr, B1_rr, B2_rr, B3_rr = cons2prim(u_rr, equations)
+
+    # Obtain left and right fluxes
+    f_ll = flux(u_ll, orientation, equations)
+    f_rr = flux(u_rr, orientation, equations)
+
+    SsL, SsR = min_max_speed_naive(u_ll, u_rr, orientation, equations)
+    sMu_L = SsL - v1_ll
+    sMu_R = SsR - v1_rr
+    if SsL >= 0
+        f1 = f_ll[1]
+        f2 = f_ll[2]
+        f3 = f_ll[3]
+        f4 = f_ll[4]
+        f5 = f_ll[5]
+        f6 = f_ll[6]
+        f7 = f_ll[7]
+        f8 = f_ll[8]
+    elseif SsR <= 0
+        f1 = f_rr[1]
+        f2 = f_rr[2]
+        f3 = f_rr[3]
+        f4 = f_rr[4]
+        f5 = f_rr[5]
+        f6 = f_rr[6]
+        f7 = f_rr[7]
+        f8 = f_rr[8]
+    else
+        #=
+        SStar = (rho_rr * v1_rr * sMu_R - rho_ll * v1_ll * sMu_L + p_ll - p_rr - B1_ll^2 + B1_rr^2 ) /
+                (rho_rr * sMu_R - rho_ll * sMu_L)
+        =#
+        # 1D
+
+        SStar = (rho_rr * v1_rr * sMu_R - rho_ll * v1_ll * sMu_L + p_ll - p_rr) /
+                (rho_rr * sMu_R - rho_ll * sMu_L)
+
+        Sdiff = SsR - SsL
+
+        # Compute HLL values for vStar, BStar
+        rho_HLL = (SsR * u_rr[1] - SsL * u_ll[1] - (f_rr[1] - f_ll[1])) / Sdiff
+
+        v1Star = (SsR * u_rr[2] - SsL * u_ll[2] - (f_rr[2] - f_ll[2])) /
+                 (Sdiff * rho_HLL)
+        v2Star = (SsR * u_rr[3] - SsL * u_ll[3] - (f_rr[3] - f_ll[3])) /
+                 (Sdiff * rho_HLL)
+        v3Star = (SsR * u_rr[4] - SsL * u_ll[4] - (f_rr[4] - f_ll[4])) /
+                 (Sdiff * rho_HLL)
+
+        B1Star = (SsR * u_rr[6] - SsL * u_ll[6] - (f_rr[6] - f_ll[6])) / Sdiff
+        B2Star = (SsR * u_rr[7] - SsL * u_ll[7] - (f_rr[7] - f_ll[7])) / Sdiff
+        B3Star = (SsR * u_rr[8] - SsL * u_ll[8] - (f_rr[8] - f_ll[8])) / Sdiff
+        if SsL <= SStar
+            SdiffStar = SsL - SStar
+
+            densStar = rho_ll * sMu_L / SdiffStar
+
+            mom_1_Star = densStar * SStar
+            mom_2_Star = densStar * v2_ll -
+                         (B1Star * B2Star - B1_ll * B2_ll) / SdiffStar
+            mom_3_Star = densStar * v3_ll -
+                         (B1Star * B3Star - B1_ll * B3_ll) / SdiffStar
+
+            pstar = rho_ll * sMu_L * (SStar - v1_ll) + p_ll - B1_ll^2 + B1Star^2
+
+            enerStar = u_ll[5] * sMu_L / SdiffStar +
+                       (pstar * SStar - p_ll * v1_ll - (B1Star *
+                         (B1Star * v1Star + B2Star * v2Star + B3Star * v3Star) -
+                         B1_ll * (B1_ll * v1_ll + B2_ll * v2_ll + B3_ll * v3_ll))) /
+                       SdiffStar
+
+            f1 = f_ll[1] + SsL * (densStar - u_ll[1])
+            f2 = f_ll[2] + SsL * (mom_1_Star - u_ll[2])
+            f3 = f_ll[3] + SsL * (mom_2_Star - u_ll[3])
+            f4 = f_ll[4] + SsL * (mom_3_Star - u_ll[4])
+            f5 = f_ll[5] + SsL * (enerStar - u_ll[5])
+            f6 = f_ll[6] + SsL * (B1Star - u_ll[6])
+            f7 = f_ll[7] + SsL * (B2Star - u_ll[7])
+            f8 = f_ll[8] + SsL * (B3Star - u_ll[8])
+        else # SStar <= Ssr
+            SdiffStar = SsR - SStar
+
+            densStar = rho_rr * sMu_R / SdiffStar
+
+            mom_1_Star = densStar * SStar
+            mom_2_Star = densStar * v2_rr -
+                         (B1Star * B2Star - B1_rr * B2_rr) / SdiffStar
+            mom_3_Star = densStar * v3_rr -
+                         (B1Star * B3Star - B1_rr * B3_rr) / SdiffStar
+
+            pstar = rho_rr * sMu_R * (SStar - v1_rr) + p_rr - B1_rr^2 + B1Star^2
+
+            enerStar = u_rr[5] * sMu_R / SdiffStar +
+                       (pstar * SStar - p_rr * v1_rr - (B1Star *
+                         (B1Star * v1Star + B2Star * v2Star + B3Star * v3Star) -
+                         B1_rr * (B1_rr * v1_rr + B2_rr * v2_rr + B3_rr * v3_rr))) /
+                       SdiffStar
+
+            f1 = f_rr[1] + SsR * (densStar - u_rr[1])
+            f2 = f_rr[2] + SsR * (mom_1_Star - u_rr[2])
+            f3 = f_rr[3] + SsR * (mom_2_Star - u_rr[3])
+            f4 = f_rr[4] + SsR * (mom_3_Star - u_rr[4])
+            f5 = f_rr[5] + SsR * (enerStar - u_rr[5])
+            f6 = f_rr[6] + SsR * (B1Star - u_rr[6])
+            f7 = f_rr[7] + SsR * (B2Star - u_rr[7])
+            f8 = f_rr[8] + SsR * (B3Star - u_rr[8])
+        end
+    end
+    return SVector(f1, f2, f3, f4, f5, f6, f7, f8)
+end
+
 # Calculate maximum wave speed for local Lax-Friedrichs-type dissipation
 @inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer,
                                      equations::IdealGlmMhdEquations1D)

test/test_tree_1d_mhd.jl

@@ -206,6 +206,31 @@ end
     end
 end
 
+@trixi_testset "elixir_mhd_torrilhon_shock_tube.jl (HLLC)" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_torrilhon_shock_tube.jl"),
+                        surface_flux=flux_hllc,
+                        l2=[
+                            0.4573896571821181, 0.4794022228630343, 0.34069729746971966,
+                            0.44795514335603914, 0.9206813325913942,
+                            1.7232531493769943e-16, 0.2889672868492896,
+                            0.25527942207797844,
+                        ],
+                        linf=[
+                            1.2181099854600026, 0.8869319941714484, 0.8763562906332127,
+                            0.9712221036087139, 1.6734231113527835,
+                            6.661338147750939e-16, 0.7035011427822784,
+                            0.6562884129650187,
+                        ])
+    # 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_mhd_ryujones_shock_tube.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ryujones_shock_tube.jl"),
                         l2=[

test/test_unit.jl

@@ -968,6 +968,7 @@ end
 
     for u in u_values
         @test flux_hll(u, u, 1, equations) ≈ flux(u, 1, equations)
+        @test flux_hllc(u, u, 1, equations) ≈ flux(u, 1, equations)
     end
 
     equations = IdealGlmMhdEquations2D(1.4, 5.0) #= c_h =#