Skip to content

Commit

Permalink
Add HLLC flux for non-cartesian meshes to CompressibleEulerEquations{…
Browse files Browse the repository at this point in the history 
…2,3}D (#1790)

* add HLLC flux for non-cartesian meshes

* add tests for HLLC flux

* Add 2D test with HLLC

* Update test_p4est_3d.jl

* Update test_p4est_2d.jl

* Update test_p4est_3d.jl

* Update src/equations/compressible_euler_3d.jl

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

* Update src/equations/compressible_euler_2d.jl

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

* Update compressible_euler_2d.jl

* Update compressible_euler_3d.jl

* Update test_p4est_2d.jl

* Update test_p4est_3d.jl

* Update compressible_euler_2d.jl

* Update compressible_euler_2d.jl

---------

Co-authored-by: Daniel Doehring <[email protected]>
Co-authored-by: Hendrik Ranocha <[email protected]>
  • Loading branch information
@benegee @DanielDoehring @ranocha
3 people authored Dec 30, 2023
1 parent dc8e927 commit a633f9c
Show file tree
Hide file tree
Showing 6 changed files with 305 additions and 17 deletions.
1 change: 1 addition & 0 deletions NEWS.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@ for human readability.

#### Added
- AMR for hyperbolic-parabolic equations on 3D `P4estMesh`
- `flux_hllc` on non-cartesian meshes for `CompressibleEulerEquations{2,3}D`

## Changes when updating to v0.6 from v0.5.x

Expand Down
102 changes: 97 additions & 5 deletions src/equations/compressible_euler_2d.jl
View file
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))
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,
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
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))
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)
# 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
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
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

0 comments on commit a633f9c

Please sign in to comment.