From cb3060ca8cff3956bec79808dc079f613c05c4ca Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 22 Nov 2023 15:52:20 +0100 Subject: [PATCH 1/2] Use total pressure --- src/equations/ideal_glm_mhd_1d.jl | 22 ++++++++++------- test/test_tree_1d_mhd.jl | 40 +++++++++++++++++++++++++------ 2 files changed, 46 insertions(+), 16 deletions(-) diff --git a/src/equations/ideal_glm_mhd_1d.jl b/src/equations/ideal_glm_mhd_1d.jl index a465571989b..0088bb7e5a2 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 (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] @@ -285,7 +289,7 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, 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) + SsL, SsR = min_max_speed_einfeldt(u_ll, u_rr, orientation, equations) sMu_L = SsL - v1_ll sMu_R = SsR - v1_rr if SsL >= 0 @@ -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 447572eee88..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,15 +235,16 @@ end @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, + 0.4573799618744708, 0.4792633358230866, 0.34064852506872795, + 0.4479668434955162, 0.9203891782415092, + 1.3216517820475193e-16, 0.28887826520860815, + 0.255281629265771, ], linf=[ - 1.217943947570543, 0.8868438459815245, 0.878215340656725, - 0.9710882819266371, 1.6742759645320984, - 2.220446049250313e-16, 0.704710220504591, 0.6562122176458641, + 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) From 26503aba8791e6f4a2feb8706a3c8f79a48e486c Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Fri, 24 Nov 2023 15:15:23 +0100 Subject: [PATCH 2/2] Update src/equations/ideal_glm_mhd_1d.jl Co-authored-by: Andrew Winters --- src/equations/ideal_glm_mhd_1d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/equations/ideal_glm_mhd_1d.jl b/src/equations/ideal_glm_mhd_1d.jl index 0088bb7e5a2..5a523daf3f6 100644 --- a/src/equations/ideal_glm_mhd_1d.jl +++ b/src/equations/ideal_glm_mhd_1d.jl @@ -272,7 +272,7 @@ 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 (eq. (12)) + # 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)