diff --git a/NEWS.md b/NEWS.md index e7bdd14eab2..cf695912ed7 100644 --- a/NEWS.md +++ b/NEWS.md @@ -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 diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index a992f99eaf4..b0fd5c53f45 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -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) @@ -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 @@ -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) diff --git a/src/equations/compressible_euler_3d.jl b/src/equations/compressible_euler_3d.jl index fc56f58025b..82c4a7efa32 100644 --- a/src/equations/compressible_euler_3d.jl +++ b/src/equations/compressible_euler_3d.jl @@ -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) @@ -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 @@ -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) diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index db34aecc168..cebc2917d52 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -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=[ diff --git a/test/test_p4est_3d.jl b/test/test_p4est_3d.jl index f2467f30204..4a2d2112c99 100644 --- a/test/test_p4est_3d.jl +++ b/test/test_p4est_3d.jl @@ -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=[ diff --git a/test/test_unit.jl b/test/test_unit.jl index b3ed29d38e3..817b4cd550d 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -1076,6 +1076,46 @@ end end end +@timed_testset "Consistency check for HLLC flux: CEE" begin + # Set up equations and dummy conservative variables state + equations = CompressibleEulerEquations2D(1.4) + u = SVector(1.1, -0.5, 2.34, 5.5) + + orientations = [1, 2] + for orientation in orientations + @test flux_hllc(u, u, orientation, equations) ≈ flux(u, orientation, equations) + end + + normal_directions = [SVector(1.0, 0.0), + SVector(0.0, 1.0), + SVector(0.5, -0.5), + SVector(-1.2, 0.3)] + + for normal_direction in normal_directions + @test flux_hllc(u, u, normal_direction, equations) ≈ + flux(u, normal_direction, equations) + end + + equations = CompressibleEulerEquations3D(1.4) + u = SVector(1.1, -0.5, 2.34, 2.4, 5.5) + + orientations = [1, 2, 3] + for orientation in orientations + @test flux_hllc(u, u, orientation, equations) ≈ flux(u, orientation, equations) + end + + normal_directions = [SVector(1.0, 0.0, 0.0), + SVector(0.0, 1.0, 0.0), + SVector(0.0, 0.0, 1.0), + SVector(0.5, -0.5, 0.2), + SVector(-1.2, 0.3, 1.4)] + + for normal_direction in normal_directions + @test flux_hllc(u, u, normal_direction, equations) ≈ + flux(u, normal_direction, equations) + end +end + @timed_testset "Consistency check for Godunov flux" begin # Set up equations and dummy conservative variables state # Burgers' Equation @@ -1280,7 +1320,7 @@ end u_values = [SVector(1.0, 0.5, -0.7, 1.0), SVector(1.5, -0.2, 0.1, 5.0)] fluxes = [flux_central, flux_ranocha, flux_shima_etal, flux_kennedy_gruber, - flux_hll, FluxHLL(min_max_speed_davis), flux_hlle] + flux_hll, FluxHLL(min_max_speed_davis), flux_hlle, flux_hllc] for f_std in fluxes f_rot = FluxRotated(f_std) @@ -1303,8 +1343,8 @@ end u_values = [SVector(1.0, 0.5, -0.7, 0.1, 1.0), SVector(1.5, -0.2, 0.1, 0.2, 5.0)] fluxes = [flux_central, flux_ranocha, flux_shima_etal, flux_kennedy_gruber, - FluxLMARS(340), - flux_hll, FluxHLL(min_max_speed_davis), flux_hlle] + FluxLMARS(340), flux_hll, FluxHLL(min_max_speed_davis), flux_hlle, flux_hllc, + ] for f_std in fluxes f_rot = FluxRotated(f_std)