From 43bc74d77ebd0e3ec0eea85b53555b5f8d2df6e1 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 1 May 2024 16:07:51 +0200 Subject: [PATCH 1/2] FluxChandrashekar_NormalDir --- src/equations/compressible_euler_2d.jl | 34 +++++++++++++++++++++- src/equations/compressible_euler_3d.jl | 40 +++++++++++++++++++++++++- test/test_p4est_2d.jl | 26 +++++++++++++++++ test/test_p4est_3d.jl | 29 +++++++++++++++++++ test/test_unit.jl | 6 ++-- 5 files changed, 131 insertions(+), 4 deletions(-) diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index 43f15a3cfb9..2bfee52904a 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -552,7 +552,7 @@ end end """ - flux_chandrashekar(u_ll, u_rr, orientation, equations::CompressibleEulerEquations2D) +flux_chandrashekar(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquations2D) Entropy conserving two-point flux by - Chandrashekar (2013) @@ -598,6 +598,38 @@ Entropy conserving two-point flux by return SVector(f1, f2, f3, f4) end +@inline function flux_chandrashekar(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerEquations2D) + # Unpack left and right state + 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] + beta_ll = 0.5 * rho_ll / p_ll + beta_rr = 0.5 * rho_rr / p_rr + specific_kin_ll = 0.5 * (v1_ll^2 + v2_ll^2) + specific_kin_rr = 0.5 * (v1_rr^2 + v2_rr^2) + + # Compute the necessary mean values + rho_avg = 0.5 * (rho_ll + rho_rr) + rho_mean = ln_mean(rho_ll, rho_rr) + beta_mean = ln_mean(beta_ll, beta_rr) + beta_avg = 0.5 * (beta_ll + beta_rr) + v1_avg = 0.5 * (v1_ll + v1_rr) + v2_avg = 0.5 * (v2_ll + v2_rr) + p_mean = 0.5 * rho_avg / beta_avg + velocity_square_avg = specific_kin_ll + specific_kin_rr + + # Multiply with average of normal velocities + f1 = rho_mean * 0.5 * (v_dot_n_ll + v_dot_n_rr) + f2 = f1 * v1_avg + p_mean * normal_direction[1] + f3 = f1 * v2_avg + p_mean * normal_direction[2] + f4 = f1 * 0.5 * (1 / (equations.gamma - 1) / beta_mean - velocity_square_avg) + + f2 * v1_avg + f3 * v2_avg + + return SVector(f1, f2, f3, f4) +end + """ flux_ranocha(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquations2D) diff --git a/src/equations/compressible_euler_3d.jl b/src/equations/compressible_euler_3d.jl index 292b912f009..f156aa29689 100644 --- a/src/equations/compressible_euler_3d.jl +++ b/src/equations/compressible_euler_3d.jl @@ -602,7 +602,7 @@ end end """ - flux_chandrashekar(u_ll, u_rr, orientation, equations::CompressibleEulerEquations3D) + flux_chandrashekar(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquations3D) Entropy conserving two-point flux by - Chandrashekar (2013) @@ -659,6 +659,44 @@ Entropy conserving two-point flux by return SVector(f1, f2, f3, f4, f5) end +@inline function flux_chandrashekar(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerEquations3D) + # Unpack left and right state + 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] + + beta_ll = 0.5 * rho_ll / p_ll + beta_rr = 0.5 * rho_rr / p_rr + specific_kin_ll = 0.5 * (v1_ll^2 + v2_ll^2 + v3_ll^2) + specific_kin_rr = 0.5 * (v1_rr^2 + v2_rr^2 + v3_rr^2) + + # Compute the necessary mean values + rho_avg = 0.5 * (rho_ll + rho_rr) + rho_mean = ln_mean(rho_ll, rho_rr) + beta_mean = ln_mean(beta_ll, beta_rr) + beta_avg = 0.5 * (beta_ll + beta_rr) + v1_avg = 0.5 * (v1_ll + v1_rr) + v2_avg = 0.5 * (v2_ll + v2_rr) + v3_avg = 0.5 * (v3_ll + v3_rr) + p_mean = 0.5 * rho_avg / beta_avg + velocity_square_avg = specific_kin_ll + specific_kin_rr + + # Multiply with average of normal velocities + f1 = rho_mean * 0.5 * (v_dot_n_ll + v_dot_n_rr) + f2 = f1 * v1_avg + p_mean * normal_direction[1] + f3 = f1 * v2_avg + p_mean * normal_direction[2] + f4 = f1 * v3_avg + p_mean * normal_direction[3] + f5 = f1 * 0.5 * (1 / (equations.gamma - 1) / beta_mean - velocity_square_avg) + + f2 * v1_avg + f3 * v2_avg + f4 * v3_avg + + return SVector(f1, f2, f3, f4, f5) +end + """ flux_ranocha(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquations3D) diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index 1bacccbf812..fbc94fdfd6d 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -178,6 +178,32 @@ end end end +@trixi_testset "elixir_euler_shockcapturing_ec.jl (flux_chandrashekar)" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_shockcapturing_ec.jl"), + l2=[ + 0.09527896382082567, + 0.10557894830184737, + 0.10559379376154387, + 0.3503791205165925, + ], + linf=[ + 0.2733486454092644, + 0.3877283966722886, + 0.38650482703821426, + 1.0053712251056308, + ], + tspan=(0.0, 1.0), + volume_flux=flux_chandrashekar) + # 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" 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 8b684b22d09..7483cde2752 100644 --- a/test/test_p4est_3d.jl +++ b/test/test_p4est_3d.jl @@ -313,6 +313,35 @@ end end end +@trixi_testset "elixir_euler_ec.jl (flux_chandrashekar)" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"), + l2=[ + 0.010368548525287055, + 0.006216054794583285, + 0.006020401857347216, + 0.006019175682769779, + 0.026228080232814154, + ], + linf=[ + 0.3169376449662026, + 0.28950510175646726, + 0.4402523227566396, + 0.4869168122387365, + 0.7999141641954051, + ], + tspan=(0.0, 0.2), + volume_flux=flux_chandrashekar, + coverage_override=(polydeg = 3,)) # Prevent long compile time in CI + # 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" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov.jl"), l2=[ diff --git a/test/test_unit.jl b/test/test_unit.jl index b23d4aa60e2..c72b51113f5 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -1431,7 +1431,8 @@ 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, - FluxLMARS(340), flux_hll, FluxHLL(min_max_speed_davis), flux_hlle, flux_hllc, + FluxLMARS(340), flux_hll, FluxHLL(min_max_speed_davis), flux_hlle, + flux_hllc, flux_chandrashekar, ] for f_std in fluxes @@ -1455,7 +1456,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, flux_hllc, + FluxLMARS(340), flux_hll, FluxHLL(min_max_speed_davis), flux_hlle, + flux_hllc, flux_chandrashekar, ] for f_std in fluxes From f04132224ba104330c1e3b800e2306b7c32cf136 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 1 May 2024 16:09:50 +0200 Subject: [PATCH 2/2] fmt --- src/equations/compressible_euler_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index 2bfee52904a..fa7af285aa4 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -552,7 +552,7 @@ end end """ -flux_chandrashekar(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquations2D) + flux_chandrashekar(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquations2D) Entropy conserving two-point flux by - Chandrashekar (2013)