diff --git a/src/equations/ideal_glm_mhd_1d.jl b/src/equations/ideal_glm_mhd_1d.jl index eca06c8c203..5a523daf3f6 100644 --- a/src/equations/ideal_glm_mhd_1d.jl +++ b/src/equations/ideal_glm_mhd_1d.jl @@ -272,6 +272,10 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, 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) + # Total pressure, i.e., thermal + magnetic pressures (eq. (12)) + p_tot_ll = p_ll + 0.5 * (B1_ll^2 + B2_ll^2 + B3_ll^2) + p_tot_rr = p_rr + 0.5 * (B1_rr^2 + B2_rr^2 + B3_rr^2) + # Conserved variables rho_v1_ll = u_ll[2] rho_v2_ll = u_ll[3] @@ -309,11 +313,11 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, 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 ) / + SStar = (rho_rr * v1_rr * sMu_R - rho_ll * v1_ll * sMu_L + p_tot_ll - p_tot_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) / + SStar = (rho_rr * v1_rr * sMu_R - rho_ll * v1_ll * sMu_L + p_tot_ll - p_tot_rr) / (rho_rr * sMu_R - rho_ll * sMu_L) Sdiff = SsR - SsL @@ -347,12 +351,12 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, 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) + #p_tot_Star = rho_ll * sMu_L * (SStar - v1_ll) + p_tot_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) + p_tot_Star = rho_ll * sMu_L * (SStar - v1_ll) + p_tot_ll # (17) enerStar = u_ll[5] * sMu_L / SdiffStar + - (pstar * SStar - p_ll * v1_ll - (B1Star * + (p_tot_Star * SStar - p_tot_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) @@ -377,12 +381,12 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, 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) + #p_tot_Star = rho_rr * sMu_R * (SStar - v1_rr) + p_tot_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) + p_tot_Star = rho_rr * sMu_R * (SStar - v1_rr) + p_tot_rr # (17) enerStar = u_rr[5] * sMu_R / SdiffStar + - (pstar * SStar - p_rr * v1_rr - (B1Star * + (p_tot_Star * SStar - p_tot_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) diff --git a/test/test_tree_1d_mhd.jl b/test/test_tree_1d_mhd.jl index 8895fe30e8b..2150ddfd074 100644 --- a/test/test_tree_1d_mhd.jl +++ b/test/test_tree_1d_mhd.jl @@ -109,6 +109,31 @@ end end end +@trixi_testset "elixir_mhd_alfven_wave.jl with flux_hllc" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"), + l2=[ + 1.036850596986597e-5, 1.965192583650368e-6, + 3.5882124656715505e-5, 3.5882124656638764e-5, + 5.270975504780837e-6, 1.1963224165731992e-16, + 3.595811808912869e-5, 3.5958118089159453e-5, + ], + linf=[ + 2.887280521446378e-5, 7.310580790352001e-6, + 0.00012390046377899755, 0.00012390046377787345, + 1.5102711136583125e-5, 2.220446049250313e-16, + 0.0001261935452181312, 0.0001261935452182006, + ], + surface_flux=flux_hllc) + # 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_ec.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec.jl"), l2=[ @@ -210,16 +235,16 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_torrilhon_shock_tube.jl"), surface_flux=flux_hllc, l2=[ - 0.45738965718253993, 0.479402222862685, 0.34069729746967664, - 0.44795514335568865, 0.9206813325913135, - 1.3216517820475193e-16, 0.2889672868491632, - 0.2552794220777942, + 0.4573799618744708, 0.4792633358230866, 0.34064852506872795, + 0.4479668434955162, 0.9203891782415092, + 1.3216517820475193e-16, 0.28887826520860815, + 0.255281629265771, ], linf=[ - 1.2181099854251536, 0.8869319941747589, 0.8763562906332134, - 0.9712221036087284, 1.6734231113527818, - 2.220446049250313e-16, 0.7035011427822779, - 0.6562884129650286, + 1.2382842201671505, 0.8929169308132259, 0.871298623806198, + 0.9822415614542821, 1.6726170732132717, + 2.220446049250313e-16, 0.7016155888023747, + 0.6556091522071984, ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities)