HLLE CEE 2D3D NonCartesian Meshes

src/equations/compressible_euler_2d.jl

@@ -1341,6 +1341,96 @@ function flux_hlle(u_ll, u_rr, orientation::Integer,
     return SVector(f1, f2, f3, f4)
 end
 
+"""
+    flux_hlle(u_ll, u_rr, normal_direction, equations::CompressibleEulerEquations2D)
+
+Computes the HLLE (Harten-Lax-van Leer-Einfeldt) flux for the compressible Euler equations.
+Special estimates of the signal velocites and linearization of the Riemann problem developed
+by Einfeldt to ensure that the internal energy and density remain positive during the computation
+of the numerical flux.
+
+- Bernd Einfeldt (1988)
+  On Godunov-type methods for gas dynamics.
+  [DOI: 10.1137/0725021](https://doi.org/10.1137/0725021)
+- Bernd Einfeldt, Claus-Dieter Munz, Philip L. Roe and Björn Sjögreen (1991)
+  On Godunov-type methods near low densities.
+  [DOI: 10.1016/0021-9991(91)90211-3](https://doi.org/10.1016/0021-9991(91)90211-3)
+"""
+function flux_hlle(u_ll, u_rr, normal_direction::AbstractVector,
+                   equations::CompressibleEulerEquations2D)
+    # Calculate primitive variables, enthalpy 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)
+
+    # `u_ll[4]` is total energy `rho_e_ll` on the left
+    H_ll = (u_ll[4] + p_ll) / rho_ll
+    c_ll = sqrt(equations.gamma * p_ll / rho_ll) * norm_
+
+    # `u_rr[4]` is total energy `rho_e_rr` on the right
+    H_rr = (u_rr[4] + p_rr) / rho_rr
+    c_rr = sqrt(equations.gamma * p_rr / rho_rr) * norm_
+
+    # Compute Roe averages
+    sqrt_rho_ll = sqrt(rho_ll)
+    sqrt_rho_rr = sqrt(rho_rr)
+    inv_sum_sqrt_rho = inv(sqrt_rho_ll + sqrt_rho_rr)
+
+    v1_roe = (sqrt_rho_ll * v1_ll + sqrt_rho_rr * v1_rr) * inv_sum_sqrt_rho
+    v2_roe = (sqrt_rho_ll * v2_ll + sqrt_rho_rr * v2_rr) * inv_sum_sqrt_rho
+    v_roe = v1_roe * normal_direction[1] + v2_roe * normal_direction[2]
+    v_roe_mag = (v1_roe * normal_direction[1])^2 + (v2_roe * normal_direction[2])^2
+
+    H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) * inv_sum_sqrt_rho
+    c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5 * v_roe_mag)) * norm_
+
+    # Compute convenience constant for positivity preservation, see
+    # https://doi.org/10.1016/0021-9991(91)90211-3
+    beta = sqrt(0.5 * (equations.gamma - 1) / equations.gamma)
+
+    # Estimate the edges of the Riemann fan (with positivity conservation)
+    SsL = min(v_roe - c_roe, v_dot_n_ll - beta * c_ll, zero(v_roe))
+    SsR = max(v_roe + c_roe, v_dot_n_rr + beta * c_rr, zero(v_roe))
+
+    if SsL >= 0.0 && SsR > 0.0
+        # Positive supersonic speed
+        f_ll = flux(u_ll, normal_direction, equations)
+
+        f1 = f_ll[1]
+        f2 = f_ll[2]
+        f3 = f_ll[3]
+        f4 = f_ll[4]
+    elseif SsR <= 0.0 && SsL < 0.0
+        # Negative supersonic speed
+        f_rr = flux(u_rr, normal_direction, equations)
+
+        f1 = f_rr[1]
+        f2 = f_rr[2]
+        f3 = f_rr[3]
+        f4 = f_rr[4]
+    else
+        # Subsonic case
+        # Compute left and right fluxes
+        f_ll = flux(u_ll, normal_direction, equations)
+        f_rr = flux(u_rr, normal_direction, equations)
+
+        f1 = (SsR * f_ll[1] - SsL * f_rr[1] + SsL * SsR * (u_rr[1] - u_ll[1])) /
+             (SsR - SsL)
+        f2 = (SsR * f_ll[2] - SsL * f_rr[2] + SsL * SsR * (u_rr[2] - u_ll[2])) /
+             (SsR - SsL)
+        f3 = (SsR * f_ll[3] - SsL * f_rr[3] + SsL * SsR * (u_rr[3] - u_ll[3])) /
+             (SsR - SsL)
+        f4 = (SsR * f_ll[4] - SsL * f_rr[4] + SsL * SsR * (u_rr[4] - u_ll[4])) /
+             (SsR - SsL)
+    end
+
+    return SVector(f1, f2, f3, f4)
+end
+
 @inline function max_abs_speeds(u, equations::CompressibleEulerEquations2D)
     rho, v1, v2, p = cons2prim(u, equations)
     c = sqrt(equations.gamma * p / rho)

src/equations/compressible_euler_3d.jl

@@ -1418,6 +1418,105 @@ function flux_hlle(u_ll, u_rr, orientation::Integer,
     return SVector(f1, f2, f3, f4, f5)
 end
 
+"""
+    flux_hlle(u_ll, u_rr, normal_direction, equations::CompressibleEulerEquations3D)
+
+Computes the HLLE (Harten-Lax-van Leer-Einfeldt) flux for the compressible Euler equations.
+Special estimates of the signal velocites and linearization of the Riemann problem developed
+by Einfeldt to ensure that the internal energy and density remain positive during the computation
+of the numerical flux.
+
+- Bernd Einfeldt (1988)
+  On Godunov-type methods for gas dynamics.
+  [DOI: 10.1137/0725021](https://doi.org/10.1137/0725021)
+- Bernd Einfeldt, Claus-Dieter Munz, Philip L. Roe and Björn Sjögreen (1991)
+  On Godunov-type methods near low densities.
+  [DOI: 10.1016/0021-9991(91)90211-3](https://doi.org/10.1016/0021-9991(91)90211-3)
+"""
+function flux_hlle(u_ll, u_rr, normal_direction::AbstractVector,
+                   equations::CompressibleEulerEquations3D)
+    # Calculate primitive variables, enthalpy 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)
+
+    # `u_ll[5]` is total energy `rho_e_ll` on the left
+    H_ll = (u_ll[5] + p_ll) / rho_ll
+    c_ll = sqrt(equations.gamma * p_ll / rho_ll) * norm_
+
+    # `u_rr[5]` is total energy `rho_e_rr` on the right
+    H_rr = (u_rr[5] + p_rr) / rho_rr
+    c_rr = sqrt(equations.gamma * p_rr / rho_rr) * norm_
+
+    # Compute Roe averages
+    sqrt_rho_ll = sqrt(rho_ll)
+    sqrt_rho_rr = sqrt(rho_rr)
+    inv_sum_sqrt_rho = inv(sqrt_rho_ll + sqrt_rho_rr)
+
+    v1_roe = (sqrt_rho_ll * v1_ll + sqrt_rho_rr * v1_rr) * inv_sum_sqrt_rho
+    v2_roe = (sqrt_rho_ll * v2_ll + sqrt_rho_rr * v2_rr) * inv_sum_sqrt_rho
+    v3_roe = (sqrt_rho_ll * v3_ll + sqrt_rho_rr * v3_rr) * inv_sum_sqrt_rho
+    v_roe = v1_roe * normal_direction[1] + v2_roe * normal_direction[2] + 
+            v3_roe * normal_direction[3]
+    v_roe_mag = v1_roe^2 + v2_roe^2 + v3_roe^2
+
+    H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) * inv_sum_sqrt_rho
+    c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5 * v_roe_mag)) * norm_
+
+    # Compute convenience constant for positivity preservation, see
+    # https://doi.org/10.1016/0021-9991(91)90211-3
+    beta = sqrt(0.5 * (equations.gamma - 1) / equations.gamma)
+
+    # Estimate the edges of the Riemann fan (with positivity conservation)
+    SsL = min(v_roe - c_roe, v_dot_n_ll - beta * c_ll, zero(v_roe))
+    SsR = max(v_roe + c_roe, v_dot_n_rr + beta * c_rr, zero(v_roe))
+
+    if SsL >= 0.0 && SsR > 0.0
+        # Positive supersonic speed
+        f_ll = flux(u_ll, normal_direction, equations)
+
+        f1 = f_ll[1]
+        f2 = f_ll[2]
+        f3 = f_ll[3]
+        f4 = f_ll[4]
+        f5 = f_ll[5]
+    elseif SsR <= 0.0 && SsL < 0.0
+        # Negative supersonic speed
+        f_rr = flux(u_rr, normal_direction, equations)
+
+        f1 = f_rr[1]
+        f2 = f_rr[2]
+        f3 = f_rr[3]
+        f4 = f_rr[4]
+        f5 = f_rr[5]
+    else
+        # Subsonic case
+        # Compute left and right fluxes
+        f_ll = flux(u_ll, normal_direction, equations)
+        f_rr = flux(u_rr, normal_direction, equations)
+
+        f1 = (SsR * f_ll[1] - SsL * f_rr[1] + SsL * SsR * (u_rr[1] - u_ll[1])) /
+             (SsR - SsL)
+        f2 = (SsR * f_ll[2] - SsL * f_rr[2] + SsL * SsR * (u_rr[2] - u_ll[2])) /
+             (SsR - SsL)
+        f3 = (SsR * f_ll[3] - SsL * f_rr[3] + SsL * SsR * (u_rr[3] - u_ll[3])) /
+             (SsR - SsL)
+        f4 = (SsR * f_ll[4] - SsL * f_rr[4] + SsL * SsR * (u_rr[4] - u_ll[4])) /
+             (SsR - SsL)
+        f5 = (SsR * f_ll[5] - SsL * f_rr[5] + SsL * SsR * (u_rr[5] - u_ll[5])) /
+             (SsR - SsL)
+    end
+
+    return SVector(f1, f2, f3, f4, f5)
+end
+
+
 @inline function max_abs_speeds(u, equations::CompressibleEulerEquations3D)
     rho, v1, v2, v3, p = cons2prim(u, equations)
     c = sqrt(equations.gamma * p / rho)

test/test_p4est_2d.jl

@@ -160,6 +160,22 @@ isdir(outdir) && rm(outdir, recursive=true)
       end
   end
 
+  @trixi_testset "elixir_euler_sedov.jl HLLE" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov.jl"),
+      l2   = [0.411541263324004, 0.2558929632770186, 0.2558929632770193, 1.298715766843915],
+      linf = [1.3457201726152221, 1.3138961427140758, 1.313896142714079, 6.293305112638921],
+      surface_flux = flux_hlle,
+      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_blast_wave_amr.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave_amr.jl"),
       l2   = [6.32183914e-01, 3.86914231e-01, 3.86869171e-01, 1.06575688e+00],

test/test_p4est_3d.jl

@@ -207,6 +207,22 @@ isdir(outdir) && rm(outdir, recursive=true)
       end
   end
 
+  @trixi_testset "elixir_euler_sedov.jl HLLE" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov.jl"),
+      l2   = [0.09946224487902565, 0.04863386374672001, 0.048633863746720116, 0.04863386374672032, 0.3751015774232693],
+      linf = [0.789241521871487, 0.42046970270100276, 0.42046970270100276, 0.4204697027010028, 4.730877375538398],
+      tspan = (0.0, 0.3),
+      surface_flux = flux_hlle)
+      # 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_source_terms_nonconforming_earth.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_nonconforming_earth.jl"),
       l2 = [6.040180337738628e-6, 5.4254175153621895e-6, 5.677698851333843e-6, 5.8017136892469794e-6, 1.3637854615117974e-5],

test/test_unit.jl

@@ -875,7 +875,7 @@ isdir(outdir) && rm(outdir, recursive=true)
                           SVector(-1.2, 0.3)]
 
     for normal_direction in normal_directions
-      @test flux_hll(u, u, normal_direction, equations) ≈ flux(u, normal_direction, equations)
+      @test flux_hlle(u, u, normal_direction, equations) ≈ flux(u, normal_direction, equations)
     end
 
     equations = CompressibleEulerEquations3D(1.4)
@@ -893,7 +893,7 @@ isdir(outdir) && rm(outdir, recursive=true)
                         SVector(-1.2, 0.3, 1.4)]
 
     for normal_direction in normal_directions
-      @test flux_hll(u, u, normal_direction, equations) ≈ flux(u, normal_direction, equations)
+      @test flux_hlle(u, u, normal_direction, equations) ≈ flux(u, normal_direction, equations)
     end
   end