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.

Sign up for GitHub

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

Add HLLC flux for non-cartesian meshes to CompressibleEulerEquations{2,3}D #1790

Merged
merged 16 commits into from
Dec 30, 2023
Merged

Add HLLC flux for non-cartesian meshes to CompressibleEulerEquations{2,3}D #1790

Show file tree
Hide file tree
Changes from 15 commits
Commits
Show all changes
16 commits
Select commit Hold shift + click to select a range
03eb725
add HLLC flux for non-cartesian meshes
benegee Dec 21, 2023
6dbb5e3
add tests for HLLC flux
benegee Dec 21, 2023
5a7462a
Merge branch 'main' into add-hllc-non-cartesian
DanielDoehring Dec 25, 2023
bb43eda
Add 2D test with HLLC
DanielDoehring Dec 28, 2023
e84598b
Update test_p4est_3d.jl
DanielDoehring Dec 28, 2023
2294cd3
Update test_p4est_2d.jl
DanielDoehring Dec 28, 2023
434c3f0
Update test_p4est_3d.jl
DanielDoehring Dec 28, 2023
408a50c
Update src/equations/compressible_euler_3d.jl
DanielDoehring Dec 28, 2023
a58b799
Update src/equations/compressible_euler_2d.jl
DanielDoehring Dec 28, 2023
d5ab42e
Update compressible_euler_2d.jl
DanielDoehring Dec 28, 2023
52792d0
Update compressible_euler_3d.jl
DanielDoehring Dec 28, 2023
fc3fea6
Update test_p4est_2d.jl
DanielDoehring Dec 28, 2023
19f710a
Update test_p4est_3d.jl
DanielDoehring Dec 28, 2023
dde444d
Update compressible_euler_2d.jl
DanielDoehring Dec 30, 2023
1fb202a
Update compressible_euler_2d.jl
DanielDoehring Dec 30, 2023
7a62f33
Update NEWS.md
DanielDoehring Dec 30, 2023
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
102 changes: 97 additions & 5 deletions src/equations/compressible_euler_2d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1150,7 +1150,7 @@ end
end

"""
flux_hllc(u_ll, u_rr, orientation, equations::CompressibleEulerEquations2D)
flux_hllc(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquations2D)

Computes the HLLC flux (HLL with Contact) for compressible Euler equations developed by E.F. Toro
[Lecture slides](http://www.prague-sum.com/download/2012/Toro_2-HLLC-RiemannSolver.pdf)
Expand Down Expand Up @@ -1185,18 +1185,18 @@ function flux_hllc(u_ll, u_rr, orientation::Integer,
if orientation == 1 # x-direction
vel_L = v1_ll
vel_R = v1_rr
ekin_roe = (sqrt_rho_ll * v2_ll + sqrt_rho_rr * v2_rr)^2
elseif orientation == 2 # y-direction
vel_L = v2_ll
vel_R = v2_rr
ekin_roe = (sqrt_rho_ll * v1_ll + sqrt_rho_rr * v1_rr)^2
end
vel_roe = (sqrt_rho_ll * vel_L + sqrt_rho_rr * vel_R) / sum_sqrt_rho
ekin_roe = 0.5 * (vel_roe^2 + ekin_roe / sum_sqrt_rho^2)
v1_roe = sqrt_rho_ll * v1_ll + sqrt_rho_rr * v1_rr
v2_roe = sqrt_rho_ll * v2_ll + sqrt_rho_rr * v2_rr
vel_roe_mag = (v1_roe^2 + v2_roe^2) / sum_sqrt_rho^2
H_ll = (rho_e_ll + p_ll) / rho_ll
H_rr = (rho_e_rr + p_rr) / rho_rr
H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) / sum_sqrt_rho
c_roe = sqrt((equations.gamma - 1) * (H_roe - ekin_roe))
c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5 * vel_roe_mag))
DanielDoehring marked this conversation as resolved.
Show resolved Hide resolved
Ssl = min(vel_L - c_ll, vel_roe - c_roe)
Ssr = max(vel_R + c_rr, vel_roe + c_roe)
sMu_L = Ssl - vel_L
Expand Down Expand Up @@ -1252,6 +1252,98 @@ function flux_hllc(u_ll, u_rr, orientation::Integer,
return SVector(f1, f2, f3, f4)
end

function flux_hllc(u_ll, u_rr, normal_direction::AbstractVector,
DanielDoehring marked this conversation as resolved.
Show resolved Hide resolved
equations::CompressibleEulerEquations2D)
# Calculate primitive variables and speed of sound
rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations)
rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations)

v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2]
v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2]

norm_ = norm(normal_direction)
norm_sq = norm_ * norm_
inv_norm_sq = inv(norm_sq)

c_ll = sqrt(equations.gamma * p_ll / rho_ll) * norm_
c_rr = sqrt(equations.gamma * p_rr / rho_rr) * norm_

# Obtain left and right fluxes
f_ll = flux(u_ll, normal_direction, equations)
f_rr = flux(u_rr, normal_direction, equations)

# Compute Roe averages
sqrt_rho_ll = sqrt(rho_ll)
sqrt_rho_rr = sqrt(rho_rr)
sum_sqrt_rho = sqrt_rho_ll + sqrt_rho_rr

v1_roe = (sqrt_rho_ll * v1_ll + sqrt_rho_rr * v1_rr) / sum_sqrt_rho
v2_roe = (sqrt_rho_ll * v2_ll + sqrt_rho_rr * v2_rr) / sum_sqrt_rho
vel_roe = v1_roe * normal_direction[1] + v2_roe * normal_direction[2]
vel_roe_mag = v1_roe^2 + v2_roe^2

e_ll = u_ll[4] / rho_ll
e_rr = u_rr[4] / rho_rr

H_ll = (u_ll[4] + p_ll) / rho_ll
H_rr = (u_rr[4] + p_rr) / rho_rr

H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) / sum_sqrt_rho
c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5 * vel_roe_mag)) * norm_

Ssl = min(v_dot_n_ll - c_ll, vel_roe - c_roe)
Ssr = max(v_dot_n_rr + c_rr, vel_roe + c_roe)
sMu_L = Ssl - v_dot_n_ll
sMu_R = Ssr - v_dot_n_rr

if Ssl >= 0.0
f1 = f_ll[1]
f2 = f_ll[2]
f3 = f_ll[3]
f4 = f_ll[4]
elseif Ssr <= 0.0
f1 = f_rr[1]
f2 = f_rr[2]
f3 = f_rr[3]
f4 = f_rr[4]
else
SStar = (rho_ll * v_dot_n_ll * sMu_L - rho_rr * v_dot_n_rr * sMu_R +
(p_rr - p_ll) * norm_sq) / (rho_ll * sMu_L - rho_rr * sMu_R)
if Ssl <= 0.0 <= SStar
densStar = rho_ll * sMu_L / (Ssl - SStar)
enerStar = e_ll +
(SStar - v_dot_n_ll) *
(SStar * inv_norm_sq + p_ll / (rho_ll * sMu_L))
UStar1 = densStar
UStar2 = densStar *
(v1_ll + (SStar - v_dot_n_ll) * normal_direction[1] * inv_norm_sq)
UStar3 = densStar *
(v2_ll + (SStar - v_dot_n_ll) * normal_direction[2] * inv_norm_sq)
UStar4 = densStar * enerStar
f1 = f_ll[1] + Ssl * (UStar1 - u_ll[1])
f2 = f_ll[2] + Ssl * (UStar2 - u_ll[2])
f3 = f_ll[3] + Ssl * (UStar3 - u_ll[3])
f4 = f_ll[4] + Ssl * (UStar4 - u_ll[4])
else
densStar = rho_rr * sMu_R / (Ssr - SStar)
enerStar = e_rr +
(SStar - v_dot_n_rr) *
(SStar * inv_norm_sq + p_rr / (rho_rr * sMu_R))
UStar1 = densStar
UStar2 = densStar *
(v1_rr + (SStar - v_dot_n_rr) * normal_direction[1] * inv_norm_sq)
UStar3 = densStar *
(v2_rr + (SStar - v_dot_n_rr) * normal_direction[2] * inv_norm_sq)
UStar4 = densStar * enerStar
f1 = f_rr[1] + Ssr * (UStar1 - u_rr[1])
f2 = f_rr[2] + Ssr * (UStar2 - u_rr[2])
f3 = f_rr[3] + Ssr * (UStar3 - u_rr[3])
f4 = f_rr[4] + Ssr * (UStar4 - u_rr[4])
end
end
return SVector(f1, f2, f3, f4)
end

"""
min_max_speed_einfeldt(u_ll, u_rr, orientation, equations::CompressibleEulerEquations2D)

Expand Down
119 changes: 110 additions & 9 deletions src/equations/compressible_euler_3d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1192,7 +1192,7 @@ end
end

"""
flux_hllc(u_ll, u_rr, orientation, equations::CompressibleEulerEquations3D)
flux_hllc(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquations3D)

Computes the HLLC flux (HLL with Contact) for compressible Euler equations developed by E.F. Toro
[Lecture slides](http://www.prague-sum.com/download/2012/Toro_2-HLLC-RiemannSolver.pdf)
Expand Down Expand Up @@ -1231,25 +1231,22 @@ function flux_hllc(u_ll, u_rr, orientation::Integer,
if orientation == 1 # x-direction
vel_L = v1_ll
vel_R = v1_rr
ekin_roe = (sqrt_rho_ll * v2_ll + sqrt_rho_rr * v2_rr)^2 +
(sqrt_rho_ll * v3_ll + sqrt_rho_rr * v3_rr)^2
elseif orientation == 2 # y-direction
vel_L = v2_ll
vel_R = v2_rr
ekin_roe = (sqrt_rho_ll * v1_ll + sqrt_rho_rr * v1_rr)^2 +
(sqrt_rho_ll * v3_ll + sqrt_rho_rr * v3_rr)^2
else # z-direction
vel_L = v3_ll
vel_R = v3_rr
ekin_roe = (sqrt_rho_ll * v1_ll + sqrt_rho_rr * v1_rr)^2 +
(sqrt_rho_ll * v2_ll + sqrt_rho_rr * v2_rr)^2
end
vel_roe = (sqrt_rho_ll * vel_L + sqrt_rho_rr * vel_R) / sum_sqrt_rho
ekin_roe = 0.5 * (vel_roe^2 + ekin_roe / sum_sqrt_rho^2)
v1_roe = sqrt_rho_ll * v1_ll + sqrt_rho_rr * v1_rr
v2_roe = sqrt_rho_ll * v2_ll + sqrt_rho_rr * v2_rr
v3_roe = sqrt_rho_ll * v3_ll + sqrt_rho_rr * v3_rr
vel_roe_mag = (v1_roe^2 + v2_roe^2 + v3_roe^2) / sum_sqrt_rho^2
H_ll = (rho_e_ll + p_ll) / rho_ll
H_rr = (rho_e_rr + p_rr) / rho_rr
H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) / sum_sqrt_rho
c_roe = sqrt((equations.gamma - 1) * (H_roe - ekin_roe))
c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5 * vel_roe_mag))
DanielDoehring marked this conversation as resolved.
Show resolved Hide resolved
Ssl = min(vel_L - c_ll, vel_roe - c_roe)
Ssr = max(vel_R + c_rr, vel_roe + c_roe)
sMu_L = Ssl - vel_L
Expand Down Expand Up @@ -1321,6 +1318,110 @@ function flux_hllc(u_ll, u_rr, orientation::Integer,
return SVector(f1, f2, f3, f4, f5)
end

function flux_hllc(u_ll, u_rr, normal_direction::AbstractVector,
equations::CompressibleEulerEquations3D)
DanielDoehring marked this conversation as resolved.
Show resolved Hide resolved
# Calculate primitive variables and speed of sound
rho_ll, v1_ll, v2_ll, v3_ll, p_ll = cons2prim(u_ll, equations)
rho_rr, v1_rr, v2_rr, v3_rr, p_rr = cons2prim(u_rr, equations)

v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] +
v3_ll * normal_direction[3]
v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] +
v3_rr * normal_direction[3]

norm_ = norm(normal_direction)
norm_sq = norm_ * norm_
inv_norm_sq = inv(norm_sq)

c_ll = sqrt(equations.gamma * p_ll / rho_ll) * norm_
c_rr = sqrt(equations.gamma * p_rr / rho_rr) * norm_

# Obtain left and right fluxes
f_ll = flux(u_ll, normal_direction, equations)
f_rr = flux(u_rr, normal_direction, equations)

# Compute Roe averages
sqrt_rho_ll = sqrt(rho_ll)
sqrt_rho_rr = sqrt(rho_rr)
sum_sqrt_rho = sqrt_rho_ll + sqrt_rho_rr

v1_roe = (sqrt_rho_ll * v1_ll + sqrt_rho_rr * v1_rr) / sum_sqrt_rho
v2_roe = (sqrt_rho_ll * v2_ll + sqrt_rho_rr * v2_rr) / sum_sqrt_rho
v3_roe = (sqrt_rho_ll * v3_ll + sqrt_rho_rr * v3_rr) / sum_sqrt_rho
vel_roe = v1_roe * normal_direction[1] + v2_roe * normal_direction[2] +
v3_roe * normal_direction[3]
vel_roe_mag = v1_roe^2 + v2_roe^2 + v3_roe^2

e_ll = u_ll[5] / rho_ll
e_rr = u_rr[5] / rho_rr

H_ll = (u_ll[5] + p_ll) / rho_ll
H_rr = (u_rr[5] + p_rr) / rho_rr

H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) / sum_sqrt_rho
c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5 * vel_roe_mag)) * norm_

Ssl = min(v_dot_n_ll - c_ll, vel_roe - c_roe)
Ssr = max(v_dot_n_rr + c_rr, vel_roe + c_roe)
sMu_L = Ssl - v_dot_n_ll
sMu_R = Ssr - v_dot_n_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]
elseif Ssr <= 0.0
f1 = f_rr[1]
f2 = f_rr[2]
f3 = f_rr[3]
f4 = f_rr[4]
f5 = f_rr[5]
else
SStar = (rho_ll * v_dot_n_ll * sMu_L - rho_rr * v_dot_n_rr * sMu_R +
(p_rr - p_ll) * norm_sq) / (rho_ll * sMu_L - rho_rr * sMu_R)
if Ssl <= 0.0 <= SStar
densStar = rho_ll * sMu_L / (Ssl - SStar)
enerStar = e_ll +
(SStar - v_dot_n_ll) *
(SStar * inv_norm_sq + p_ll / (rho_ll * sMu_L))
UStar1 = densStar
UStar2 = densStar *
(v1_ll + (SStar - v_dot_n_ll) * normal_direction[1] * inv_norm_sq)
UStar3 = densStar *
(v2_ll + (SStar - v_dot_n_ll) * normal_direction[2] * inv_norm_sq)
UStar4 = densStar *
(v3_ll + (SStar - v_dot_n_ll) * normal_direction[3] * inv_norm_sq)
UStar5 = densStar * enerStar
f1 = f_ll[1] + Ssl * (UStar1 - u_ll[1])
f2 = f_ll[2] + Ssl * (UStar2 - u_ll[2])
f3 = f_ll[3] + Ssl * (UStar3 - u_ll[3])
f4 = f_ll[4] + Ssl * (UStar4 - u_ll[4])
f5 = f_ll[5] + Ssl * (UStar5 - u_ll[5])
else
densStar = rho_rr * sMu_R / (Ssr - SStar)
enerStar = e_rr +
(SStar - v_dot_n_rr) *
(SStar * inv_norm_sq + p_rr / (rho_rr * sMu_R))
UStar1 = densStar
UStar2 = densStar *
(v1_rr + (SStar - v_dot_n_rr) * normal_direction[1] * inv_norm_sq)
UStar3 = densStar *
(v2_rr + (SStar - v_dot_n_rr) * normal_direction[2] * inv_norm_sq)
UStar4 = densStar *
(v3_rr + (SStar - v_dot_n_rr) * normal_direction[3] * inv_norm_sq)
UStar5 = densStar * enerStar
f1 = f_rr[1] + Ssr * (UStar1 - u_rr[1])
f2 = f_rr[2] + Ssr * (UStar2 - u_rr[2])
f3 = f_rr[3] + Ssr * (UStar3 - u_rr[3])
f4 = f_rr[4] + Ssr * (UStar4 - u_rr[4])
f5 = f_rr[5] + Ssr * (UStar5 - u_rr[5])
end
end
return SVector(f1, f2, f3, f4, f5)
end

"""
min_max_speed_einfeldt(u_ll, u_rr, orientation, equations::CompressibleEulerEquations3D)

Expand Down
26 changes: 26 additions & 0 deletions test/test_p4est_2d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -203,6 +203,32 @@ end
end
end

@trixi_testset "elixir_euler_sedov.jl with HLLC Flux" begin
@test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov.jl"),
l2=[
0.4229948321239887,
0.2559038337457483,
0.2559038337457484,
1.2990046683564136,
],
linf=[
1.4989357969730492,
1.325456585141623,
1.3254565851416251,
6.331283015053501,
],
surface_flux=flux_hllc,
tspan=(0.0, 0.3))
# 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_euler_sedov.jl (HLLE)" begin
@test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov.jl"),
l2=[
Expand Down
28 changes: 28 additions & 0 deletions test/test_p4est_3d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -234,6 +234,34 @@ end
end
end

@trixi_testset "elixir_euler_free_stream_extruded.jl with HLLC FLux" begin
@test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream_extruded.jl"),
l2=[
8.444868392439035e-16,
4.889826056731442e-15,
2.2921260987087585e-15,
4.268460455702414e-15,
1.1356712092620279e-14,
],
linf=[
7.749356711883593e-14,
4.513472928735496e-13,
2.9790059308254513e-13,
1.057154364048074e-12,
1.6271428648906294e-12,
],
tspan=(0.0, 0.1),
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_euler_ec.jl" begin
@test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
l2=[
Expand Down
Loading
Loading