From 19781af753e8f1f74f626755a2e91114c2087753 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 3 Nov 2023 14:51:01 +0100 Subject: [PATCH 01/11] start HLLC MHD 1D Cartesian --- Project.toml | 2 + src/equations/ideal_glm_mhd_1d.jl | 80 +++++++++++++++++++++++++++++++ 2 files changed, 82 insertions(+) diff --git a/Project.toml b/Project.toml index f19f7fdecc3..c1ec26e144e 100644 --- a/Project.toml +++ b/Project.toml @@ -19,7 +19,9 @@ MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" Octavian = "6fd5a793-0b7e-452c-907f-f8bfe9c57db4" OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" +OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" P4est = "7d669430-f675-4ae7-b43e-fab78ec5a902" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" diff --git a/src/equations/ideal_glm_mhd_1d.jl b/src/equations/ideal_glm_mhd_1d.jl index 7e5c94c7bc3..6cd2d24eff8 100644 --- a/src/equations/ideal_glm_mhd_1d.jl +++ b/src/equations/ideal_glm_mhd_1d.jl @@ -259,6 +259,86 @@ 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) + + v1_ll = rho_v1_ll / rho_ll + e_ll = rho_e_ll / rho_ll + p_ll = (equations.gamma - 1) * (rho_e_ll - 1 / 2 * rho_ll * v1_ll^2) + c_ll = sqrt(equations.gamma * p_ll / rho_ll) + + v1_rr = rho_v1_rr / rho_rr + e_rr = rho_e_rr / rho_rr + p_rr = (equations.gamma - 1) * (rho_e_rr - 1 / 2 * rho_rr * v1_rr^2) + c_rr = sqrt(equations.gamma * p_rr / rho_rr) + + # 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.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.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 * sMu_R - rho_ll * sMu_L + p_ll - p_rr - B1_ll^2 + B1_rr^2 ) / + (rho_rr * sMu_R - rho_ll * sMu_L) + + # Compute HLL values for BStar + if Ssl <= SStar + densStar = rho_ll * sMu_L / (Ssl - SStar) + mom_x_Star = densStar * SStar + + Mom_y_Star = rho_ll*v2_ll*(SsL - u_ll)/sMu_L - () + Mom_z_Star = densStar * SStar + enerStar = e_ll + (SStar - v1_ll) * (SStar + p_ll / (rho_ll * sMu_L)) + UStar3 = densStar * enerStar + + f1 = f_ll[1] + Ssl * (densStar - rho_ll) + f2 = f_ll[2] + Ssl * (UStar2 - rho_v1_ll) + f3 = f_ll[3] + Ssl * (UStar3 - rho_e_ll) + else # SStar <= Ssr + densStar = rho_rr * sMu_R / (Ssr - SStar) + enerStar = e_rr + (SStar - vel_R) * (SStar + p_rr / (rho_rr * sMu_R)) + UStar2 = densStar * SStar + UStar3 = densStar * enerStar + + #end + f1 = f_rr[1] + Ssr * (densStar - rho_rr) + f2 = f_rr[2] + Ssr * (UStar2 - rho_v1_rr) + f3 = f_rr[3] + Ssr * (UStar3 - rho_e_rr) + end + end + return SVector(f1, f2, f3) +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) From 47ab3d9a3a8b88ca3a0eacdca0c0c62257acc61f Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 3 Nov 2023 16:07:47 +0100 Subject: [PATCH 02/11] Continue HLLC --- src/equations/ideal_glm_mhd_1d.jl | 101 +++++++++++++++++++----------- 1 file changed, 64 insertions(+), 37 deletions(-) diff --git a/src/equations/ideal_glm_mhd_1d.jl b/src/equations/ideal_glm_mhd_1d.jl index 6cd2d24eff8..7c77203f666 100644 --- a/src/equations/ideal_glm_mhd_1d.jl +++ b/src/equations/ideal_glm_mhd_1d.jl @@ -272,24 +272,14 @@ 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) - v1_ll = rho_v1_ll / rho_ll - e_ll = rho_e_ll / rho_ll - p_ll = (equations.gamma - 1) * (rho_e_ll - 1 / 2 * rho_ll * v1_ll^2) - c_ll = sqrt(equations.gamma * p_ll / rho_ll) - - v1_rr = rho_v1_rr / rho_rr - e_rr = rho_e_rr / rho_rr - p_rr = (equations.gamma - 1) * (rho_e_rr - 1 / 2 * rho_rr * v1_rr^2) - c_rr = sqrt(equations.gamma * p_rr / rho_rr) - # 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.0 + 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.0 f1 = f_ll[1] f2 = f_ll[2] f3 = f_ll[3] @@ -298,7 +288,7 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, f6 = f_ll[6] f7 = f_ll[7] f8 = f_ll[8] - elseif Ssr <= 0.0 + elseif SsR <= 0.0 f1 = f_rr[1] f2 = f_rr[2] f3 = f_rr[3] @@ -308,35 +298,72 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, f7 = f_rr[7] f8 = f_rr[8] else + Sdiff = SsR - SsL + #= SStar = (rho_rr * sMu_R - rho_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) - # Compute HLL values for BStar - if Ssl <= SStar - densStar = rho_ll * sMu_L / (Ssl - SStar) - mom_x_Star = densStar * SStar + # Compute HLL values for vStar, BStar + v1Star = (SsR * v1_rr - SsL * v1_ll - (f_rr[2] - f_ll[2]))/Sdiff + v2Star = (SsR * v2_rr - SsL * v2_ll - (f_rr[3] - f_ll[3]))/Sdiff + v3Star = (SsR * v3_rr - SsL * v3_ll - (f_rr[4] - f_ll[4]))/Sdiff + + B1Star = (SsR * B1_rr - SsL * B1_ll - (f_rr[5] - f_ll[5]))/Sdiff + B2Star = (SsR * B2_rr - SsL * B2_ll - (f_rr[6] - f_ll[6]))/Sdiff + B3Star = (SsR * B3_rr - SsL * B3_ll - (f_rr[7] - f_ll[7]))/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[8]*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 * (B1Star - u_ll[5]) + f6 = f_ll[6] + SsL * (B2Star - u_ll[6]) + f7 = f_ll[7] + SsL * (B3Star - u_ll[7]) + f8 = f_ll[8] + SsL * (enerStar - u_ll[8]) + else # SStar <= Ssr + SdiffStar = SsR - SStar - Mom_y_Star = rho_ll*v2_ll*(SsL - u_ll)/sMu_L - () - Mom_z_Star = densStar * SStar - enerStar = e_ll + (SStar - v1_ll) * (SStar + p_ll / (rho_ll * sMu_L)) - UStar3 = densStar * enerStar + densStar = rho_rr * sMu_R / SdiffStar - f1 = f_ll[1] + Ssl * (densStar - rho_ll) - f2 = f_ll[2] + Ssl * (UStar2 - rho_v1_ll) - f3 = f_ll[3] + Ssl * (UStar3 - rho_e_ll) - else # SStar <= Ssr - densStar = rho_rr * sMu_R / (Ssr - SStar) - enerStar = e_rr + (SStar - vel_R) * (SStar + p_rr / (rho_rr * sMu_R)) - UStar2 = densStar * SStar - UStar3 = densStar * enerStar - - #end - f1 = f_rr[1] + Ssr * (densStar - rho_rr) - f2 = f_rr[2] + Ssr * (UStar2 - rho_v1_rr) - f3 = f_rr[3] + Ssr * (UStar3 - rho_e_rr) + 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[8]*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 * (B1Star - u_rr[5]) + f6 = f_rr[6] + SsR * (B2Star - u_rr[6]) + f7 = f_rr[7] + SsR * (B3Star - u_rr[7]) + f8 = f_rr[8] + SsR * (enerStar - u_rr[8]) end end - return SVector(f1, f2, f3) + return SVector(f1, f2, f3, f4, f5, f6, f7, f8) end # Calculate maximum wave speed for local Lax-Friedrichs-type dissipation From 44915bc9eb6304a43dd73b4fb08d7a41e096addc Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 3 Nov 2023 16:53:17 +0100 Subject: [PATCH 03/11] tests and format --- src/equations/ideal_glm_mhd_1d.jl | 69 +++++++++++++++++++------------ test/test_tree_1d_mhd.jl | 25 +++++++++++ test/test_unit.jl | 1 + 3 files changed, 68 insertions(+), 27 deletions(-) diff --git a/src/equations/ideal_glm_mhd_1d.jl b/src/equations/ideal_glm_mhd_1d.jl index 7c77203f666..2d513d559f7 100644 --- a/src/equations/ideal_glm_mhd_1d.jl +++ b/src/equations/ideal_glm_mhd_1d.jl @@ -298,69 +298,84 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, f7 = f_rr[7] f8 = f_rr[8] else - Sdiff = SsR - SsL #= - SStar = (rho_rr * sMu_R - rho_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_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) - # Compute HLL values for vStar, BStar - v1Star = (SsR * v1_rr - SsL * v1_ll - (f_rr[2] - f_ll[2]))/Sdiff - v2Star = (SsR * v2_rr - SsL * v2_ll - (f_rr[3] - f_ll[3]))/Sdiff - v3Star = (SsR * v3_rr - SsL * v3_ll - (f_rr[4] - f_ll[4]))/Sdiff + Sdiff = SsR - SsL - B1Star = (SsR * B1_rr - SsL * B1_ll - (f_rr[5] - f_ll[5]))/Sdiff - B2Star = (SsR * B2_rr - SsL * B2_ll - (f_rr[6] - f_ll[6]))/Sdiff - B3Star = (SsR * B3_rr - SsL * B3_ll - (f_rr[7] - f_ll[7]))/Sdiff + # 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 + 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[8]*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 + 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 * (B1Star - u_ll[5]) - f6 = f_ll[6] + SsL * (B2Star - u_ll[6]) - f7 = f_ll[7] + SsL * (B3Star - u_ll[7]) - f8 = f_ll[8] + SsL * (enerStar - u_ll[8]) + 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 + 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[8]*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 + 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 * (B1Star - u_rr[5]) - f6 = f_rr[6] + SsR * (B2Star - u_rr[6]) - f7 = f_rr[7] + SsR * (B3Star - u_rr[7]) - f8 = f_rr[8] + SsR * (enerStar - u_rr[8]) + 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) diff --git a/test/test_tree_1d_mhd.jl b/test/test_tree_1d_mhd.jl index b34bdf0660c..a73166a2316 100644 --- a/test/test_tree_1d_mhd.jl +++ b/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=[ diff --git a/test/test_unit.jl b/test/test_unit.jl index a73dfab5504..cb2ccffa98c 100644 --- a/test/test_unit.jl +++ b/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 =# From 255e84ce7beafc51e84bf53a21d4b64bcaafedb9 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 3 Nov 2023 16:54:23 +0100 Subject: [PATCH 04/11] remove ode diff eq and plots --- Project.toml | 2 -- 1 file changed, 2 deletions(-) diff --git a/Project.toml b/Project.toml index c1ec26e144e..f19f7fdecc3 100644 --- a/Project.toml +++ b/Project.toml @@ -19,9 +19,7 @@ MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" Octavian = "6fd5a793-0b7e-452c-907f-f8bfe9c57db4" OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" -OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" P4est = "7d669430-f675-4ae7-b43e-fab78ec5a902" -Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" From fe6004a425af57db7523b2b69f5da18f067d927d Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Mon, 6 Nov 2023 12:10:19 +0100 Subject: [PATCH 05/11] Update src/equations/ideal_glm_mhd_1d.jl Co-authored-by: Hendrik Ranocha --- 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 2d513d559f7..62fd0fb7e33 100644 --- a/src/equations/ideal_glm_mhd_1d.jl +++ b/src/equations/ideal_glm_mhd_1d.jl @@ -288,7 +288,7 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, f6 = f_ll[6] f7 = f_ll[7] f8 = f_ll[8] - elseif SsR <= 0.0 + elseif SsR <= 0 f1 = f_rr[1] f2 = f_rr[2] f3 = f_rr[3] From 7c67dad11ae1a21d2fefabb8749d82cc9ab7bffa Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Mon, 6 Nov 2023 12:10:26 +0100 Subject: [PATCH 06/11] Update src/equations/ideal_glm_mhd_1d.jl Co-authored-by: Hendrik Ranocha --- 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 62fd0fb7e33..bbecef7a646 100644 --- a/src/equations/ideal_glm_mhd_1d.jl +++ b/src/equations/ideal_glm_mhd_1d.jl @@ -279,7 +279,7 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, 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.0 + if SsL >= 0 f1 = f_ll[1] f2 = f_ll[2] f3 = f_ll[3] From 76c263c041d7bbd97967d07556b64b8cf7f52f2c Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sat, 11 Nov 2023 14:16:42 +0100 Subject: [PATCH 07/11] comments --- src/equations/ideal_glm_mhd_1d.jl | 32 +++++++++++++++++-------------- 1 file changed, 18 insertions(+), 14 deletions(-) diff --git a/src/equations/ideal_glm_mhd_1d.jl b/src/equations/ideal_glm_mhd_1d.jl index bbecef7a646..31f62f01b19 100644 --- a/src/equations/ideal_glm_mhd_1d.jl +++ b/src/equations/ideal_glm_mhd_1d.jl @@ -298,18 +298,20 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, 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) =# - # 1D - + # 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 * 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])) / @@ -325,22 +327,23 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, if SsL <= SStar SdiffStar = SsL - SStar - densStar = rho_ll * sMu_L / SdiffStar + densStar = rho_ll * sMu_L / SdiffStar # (19) - mom_1_Star = densStar * SStar + mom_1_Star = densStar * SStar # (20) mom_2_Star = densStar * v2_ll - - (B1Star * B2Star - B1_ll * B2_ll) / SdiffStar + (B1Star * B2Star - B1_ll * B2_ll) / SdiffStar # (21) mom_3_Star = densStar * v3_ll - - (B1Star * B3Star - B1_ll * B3_ll) / SdiffStar + (B1Star * B3Star - B1_ll * B3_ll) / SdiffStar # (22) - pstar = rho_ll * sMu_L * (SStar - v1_ll) + p_ll - B1_ll^2 + B1Star^2 + pstar = rho_ll * sMu_L * (SStar - v1_ll) + p_ll - B1_ll^2 + B1Star^2 # (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 + 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]) @@ -352,22 +355,23 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, else # SStar <= Ssr SdiffStar = SsR - SStar - densStar = rho_rr * sMu_R / SdiffStar + densStar = rho_rr * sMu_R / SdiffStar # (19) - mom_1_Star = densStar * SStar + mom_1_Star = densStar * SStar # (20) mom_2_Star = densStar * v2_rr - - (B1Star * B2Star - B1_rr * B2_rr) / SdiffStar + (B1Star * B2Star - B1_rr * B2_rr) / SdiffStar # (21) mom_3_Star = densStar * v3_rr - - (B1Star * B3Star - B1_rr * B3_rr) / SdiffStar + (B1Star * B3Star - B1_rr * B3_rr) / SdiffStar # (22) - pstar = rho_rr * sMu_R * (SStar - v1_rr) + p_rr - B1_rr^2 + B1Star^2 + pstar = rho_rr * sMu_R * (SStar - v1_rr) + p_rr - B1_rr^2 + B1Star^2 # (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 + 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]) From bcf0d0ce6524128d3ed2a0ab3cf5196d388a921a Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sat, 11 Nov 2023 14:21:46 +0100 Subject: [PATCH 08/11] further 1d simplifications --- src/equations/ideal_glm_mhd_1d.jl | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/src/equations/ideal_glm_mhd_1d.jl b/src/equations/ideal_glm_mhd_1d.jl index 31f62f01b19..1c24162c5e5 100644 --- a/src/equations/ideal_glm_mhd_1d.jl +++ b/src/equations/ideal_glm_mhd_1d.jl @@ -321,7 +321,10 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, 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 + #B1Star = (SsR * u_rr[6] - SsL * u_ll[6] - (f_rr[6] - f_ll[6])) / Sdiff + # 1D B1 = constant => B1_ll = B1_rr = B1Star + B1Star = B1_ll + 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 @@ -335,7 +338,9 @@ 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) + #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 * @@ -363,7 +368,9 @@ 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) + #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 * From 8a9afcb9f170f6721f5808ed2418e28fa3375617 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sat, 11 Nov 2023 14:33:31 +0100 Subject: [PATCH 09/11] label variables in hll flux comp --- src/equations/ideal_glm_mhd_1d.jl | 23 ++++++++++++++++------- 1 file changed, 16 insertions(+), 7 deletions(-) diff --git a/src/equations/ideal_glm_mhd_1d.jl b/src/equations/ideal_glm_mhd_1d.jl index 1c24162c5e5..f4f8f9e92fc 100644 --- a/src/equations/ideal_glm_mhd_1d.jl +++ b/src/equations/ideal_glm_mhd_1d.jl @@ -272,6 +272,15 @@ 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) + # 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) @@ -312,21 +321,21 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, # 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 * u_rr[1] - SsL * u_ll[1] - (f_rr[1] - f_ll[1])) / Sdiff + rho_HLL = (SsR * rho_rr - SsL * rho_ll - (f_rr[1] - f_ll[1])) / Sdiff - v1Star = (SsR * u_rr[2] - SsL * u_ll[2] - (f_rr[2] - f_ll[2])) / + v1Star = (SsR * rho_v1_rr - SsL * rho_v1_ll - (f_rr[2] - f_ll[2])) / (Sdiff * rho_HLL) - v2Star = (SsR * u_rr[3] - SsL * u_ll[3] - (f_rr[3] - f_ll[3])) / + v2Star = (SsR * rho_v2_rr - SsL * rho_v2_ll - (f_rr[3] - f_ll[3])) / (Sdiff * rho_HLL) - v3Star = (SsR * u_rr[4] - SsL * u_ll[4] - (f_rr[4] - f_ll[4])) / + v3Star = (SsR * rho_v3_rr - SsL * rho_v3_ll - (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 + #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 * 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 + 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 From a52303bdce8046d69256ed1c1073797ad7564a58 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 12 Nov 2023 14:56:37 +0100 Subject: [PATCH 10/11] update test vals for updated ode --- test/test_tree_1d_mhd.jl | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/test/test_tree_1d_mhd.jl b/test/test_tree_1d_mhd.jl index a73166a2316..79e10de2e74 100644 --- a/test/test_tree_1d_mhd.jl +++ b/test/test_tree_1d_mhd.jl @@ -210,16 +210,10 @@ end @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, + 0.4574266553239646, 0.4794143154876439, 0.3407079689595056, 0.44797768430829343, 0.9206916204424165, 1.3216517820475193e-16, 0.2889748702415378, 0.25529778018020927 ], linf=[ - 1.2181099854600026, 0.8869319941714484, 0.8763562906332127, - 0.9712221036087139, 1.6734231113527835, - 6.661338147750939e-16, 0.7035011427822784, - 0.6562884129650187, + 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) From 65bdcc9bc735de2ede3d0790a551df1c71664000 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 12 Nov 2023 14:56:59 +0100 Subject: [PATCH 11/11] fmt --- test/test_tree_1d_mhd.jl | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/test/test_tree_1d_mhd.jl b/test/test_tree_1d_mhd.jl index 79e10de2e74..447572eee88 100644 --- a/test/test_tree_1d_mhd.jl +++ b/test/test_tree_1d_mhd.jl @@ -210,10 +210,15 @@ 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.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 + 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)