diff --git a/src/equations/ideal_glm_mhd_1d.jl b/src/equations/ideal_glm_mhd_1d.jl index 85030e8a5ad..a465571989b 100644 --- a/src/equations/ideal_glm_mhd_1d.jl +++ b/src/equations/ideal_glm_mhd_1d.jl @@ -259,6 +259,148 @@ 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) + + # Conserved variables + rho_v1_ll = u_ll[2] + rho_v2_ll = u_ll[3] + rho_v3_ll = u_ll[4] + + rho_v1_rr = u_rr[2] + rho_v2_rr = u_rr[3] + rho_v3_rr = u_rr[4] + + # 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 + # Compute the "HLLC-speed", eq. (14) from paper mentioned above + #= + 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) + =# + # Simplification for 1D: B1 is constant + 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 + # These correspond to eq. (28) and (30) from the referenced paper + # and the classic HLL intermediate state given by (2) + rho_HLL = (SsR * rho_rr - SsL * rho_ll - (f_rr[1] - f_ll[1])) / Sdiff + + v1Star = (SsR * rho_v1_rr - SsL * rho_v1_ll - (f_rr[2] - f_ll[2])) / + (Sdiff * rho_HLL) + v2Star = (SsR * rho_v2_rr - SsL * rho_v2_ll - (f_rr[3] - f_ll[3])) / + (Sdiff * rho_HLL) + v3Star = (SsR * rho_v3_rr - SsL * rho_v3_ll - (f_rr[4] - f_ll[4])) / + (Sdiff * rho_HLL) + + #B1Star = (SsR * B1_rr - SsL * B1_ll - (f_rr[6] - f_ll[6])) / Sdiff + # 1D B1 = constant => B1_ll = B1_rr = B1Star + B1Star = B1_ll + + B2Star = (SsR * B2_rr - SsL * B2_ll - (f_rr[7] - f_ll[7])) / Sdiff + B3Star = (SsR * B3_rr - SsL * B3_ll - (f_rr[8] - f_ll[8])) / Sdiff + if SsL <= SStar + SdiffStar = SsL - SStar + + densStar = rho_ll * sMu_L / SdiffStar # (19) + + mom_1_Star = densStar * SStar # (20) + mom_2_Star = densStar * v2_ll - + (B1Star * B2Star - B1_ll * B2_ll) / SdiffStar # (21) + mom_3_Star = densStar * v3_ll - + (B1Star * B3Star - B1_ll * B3_ll) / SdiffStar # (22) + + #pstar = rho_ll * sMu_L * (SStar - v1_ll) + p_ll - B1_ll^2 + B1Star^2 # (17) + # 1D B1 = constant => B1_ll = B1_rr = B1Star + pstar = rho_ll * sMu_L * (SStar - v1_ll) + p_ll # (17) + + 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 # (23) + + # Classic HLLC update (32) + 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 # (19) + + mom_1_Star = densStar * SStar # (20) + mom_2_Star = densStar * v2_rr - + (B1Star * B2Star - B1_rr * B2_rr) / SdiffStar # (21) + mom_3_Star = densStar * v3_rr - + (B1Star * B3Star - B1_rr * B3_rr) / SdiffStar # (22) + + #pstar = rho_rr * sMu_R * (SStar - v1_rr) + p_rr - B1_rr^2 + B1Star^2 # (17) + # 1D B1 = constant => B1_ll = B1_rr = B1Star + pstar = rho_rr * sMu_R * (SStar - v1_rr) + p_rr # (17) + + 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 # (23) + + # Classic HLLC update (32) + 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) diff --git a/test/Project.toml b/test/Project.toml index fcb96b9e48f..ecae0ac0900 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -12,7 +12,7 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] -Aqua = "0.7" +Aqua = "0.8" CairoMakie = "0.6, 0.7, 0.8, 0.9, 0.10" Downloads = "1" ForwardDiff = "0.10" diff --git a/test/test_aqua.jl b/test/test_aqua.jl index 9f57791406f..93457caba28 100644 --- a/test/test_aqua.jl +++ b/test/test_aqua.jl @@ -11,8 +11,8 @@ include("test_trixi.jl") ambiguities = false, # exceptions necessary for adding a new method `StartUpDG.estimate_h` # in src/solvers/dgmulti/sbp.jl - piracy = (treat_as_own = [Trixi.StartUpDG.RefElemData, - Trixi.StartUpDG.MeshData],)) + piracies = (treat_as_own = [Trixi.StartUpDG.RefElemData, + Trixi.StartUpDG.MeshData],)) end end #module diff --git a/test/test_tree_1d_mhd.jl b/test/test_tree_1d_mhd.jl index b34bdf0660c..447572eee88 100644 --- a/test/test_tree_1d_mhd.jl +++ b/test/test_tree_1d_mhd.jl @@ -206,6 +206,30 @@ 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.4574266553239646, 0.4794143154876439, 0.3407079689595056, + 0.44797768430829343, 0.9206916204424165, + 1.3216517820475193e-16, 0.2889748702415378, + 0.25529778018020927, + ], + linf=[ + 1.217943947570543, 0.8868438459815245, 0.878215340656725, + 0.9710882819266371, 1.6742759645320984, + 2.220446049250313e-16, 0.704710220504591, 0.6562122176458641, + ]) + # 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=[ diff --git a/test/test_unit.jl b/test/test_unit.jl index 3d8f0c45c4d..077b8272d5c 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -958,6 +958,7 @@ end for u in u_values @test flux_hlle(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 =#