FluxChandrashekar_NormalDir

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)

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)

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=[

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=[

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