Skip to content

Commit

Permalink
HLLC Flux for Ideal MHD 1D (#1702)
Browse files Browse the repository at this point in the history 
* start HLLC MHD 1D Cartesian

* Continue HLLC

* tests and format

* remove ode diff eq and plots

* Update src/equations/ideal_glm_mhd_1d.jl

Co-authored-by: Hendrik Ranocha <[email protected]>

* Update src/equations/ideal_glm_mhd_1d.jl

Co-authored-by: Hendrik Ranocha <[email protected]>

* comments

* further 1d simplifications

* label variables in hll flux comp

* update test vals for updated ode

* fmt

---------

Co-authored-by: Hendrik Ranocha <[email protected]>
  • Loading branch information
@DanielDoehring @ranocha
DanielDoehring and ranocha authored Nov 15, 2023
1 parent 7bb3b46 commit 3f3ac9a
Show file tree
Hide file tree
Showing 3 changed files with 167 additions and 0 deletions.
142 changes: 142 additions & 0 deletions src/equations/ideal_glm_mhd_1d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
24 changes: 24 additions & 0 deletions test/test_tree_1d_mhd.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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=[
Expand Down
1 change: 1 addition & 0 deletions test/test_unit.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 =#
Expand Down

0 comments on commit 3f3ac9a

Please sign in to comment.