Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

HLLC Flux for Ideal MHD 1D #1702

Merged
Merged
Show file tree
Hide file tree
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
tests and format
  • Loading branch information
DanielDoehring committed Nov 3, 2023
commit 44915bc9eb6304a43dd73b4fb08d7a41e096addc
69 changes: 42 additions & 27 deletions src/equations/ideal_glm_mhd_1d.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
DanielDoehring marked this conversation as resolved.
Show resolved Hide resolved

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)
Expand Down
25 changes: 25 additions & 0 deletions test/test_tree_1d_mhd.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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=[
Expand Down
1 change: 1 addition & 0 deletions test/test_unit.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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 =#
Expand Down