HLLE CEE 2D3D NonCartesian Meshes (#1692)

* HLLE CEE 2D3D NonCartesian Meshes * format * hlle via hll * format test * format test * format * do not export hlle * Correct test vals * test values CI * Update src/equations/compressible_euler_2d.jl Co-authored-by: Hendrik Ranocha <[email protected]> * Update src/equations/compressible_euler_1d.jl Co-authored-by: Hendrik Ranocha <[email protected]> * Update src/equations/compressible_euler_2d.jl Co-authored-by: Hendrik Ranocha <[email protected]> * Update src/equations/compressible_euler_3d.jl Co-authored-by: Hendrik Ranocha <[email protected]> * Update src/equations/compressible_euler_3d.jl Co-authored-by: Hendrik Ranocha <[email protected]> * apply suggestions * additional sentence * Fix typo * typos * correct test vals --------- Co-authored-by: Hendrik Ranocha <[email protected]>
trixi-framework · Nov 1, 2023 · c438d34 · c438d34
1 parent abd97b1
commit c438d34
Show file tree

Hide file tree

Showing 10 changed files with 265 additions and 113 deletions.
diff --git a/src/equations/compressible_euler_1d.jl b/src/equations/compressible_euler_1d.jl
@@ -808,7 +808,7 @@ function flux_hllc(u_ll, u_rr, orientation::Integer,
 end
 
 """
-    flux_hlle(u_ll, u_rr, orientation, equations::CompressibleEulerEquations1D)
+    min_max_speed_einfeldt(u_ll, u_rr, orientation, equations::CompressibleEulerEquations1D)
 
 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
@@ -825,8 +825,8 @@ Compactly summarized:
   Numerical methods for conservation laws and related equations.
   [Link](https://metaphor.ethz.ch/x/2019/hs/401-4671-00L/literature/mishra_hyperbolic_pdes.pdf)
 """
-function flux_hlle(u_ll, u_rr, orientation::Integer,
-                   equations::CompressibleEulerEquations1D)
+@inline function min_max_speed_einfeldt(u_ll, u_rr, orientation::Integer,
+                                        equations::CompressibleEulerEquations1D)
     # Calculate primitive variables, enthalpy and speed of sound
     rho_ll, v_ll, p_ll = cons2prim(u_ll, equations)
     rho_rr, v_rr, p_rr = cons2prim(u_rr, equations)
@@ -858,35 +858,7 @@ function flux_hlle(u_ll, u_rr, orientation::Integer,
     SsL = min(v_roe - c_roe, v_ll - beta * c_ll, zero(v_roe))
     SsR = max(v_roe + c_roe, v_rr + beta * c_rr, zero(v_roe))
 
-    if SsL >= 0.0 && SsR > 0.0
-        # Positive supersonic speed
-        f_ll = flux(u_ll, orientation, equations)
-
-        f1 = f_ll[1]
-        f2 = f_ll[2]
-        f3 = f_ll[3]
-    elseif SsR <= 0.0 && SsL < 0.0
-        # Negative supersonic speed
-        f_rr = flux(u_rr, orientation, equations)
-
-        f1 = f_rr[1]
-        f2 = f_rr[2]
-        f3 = f_rr[3]
-    else
-        # Subsonic case
-        # Compute left and right fluxes
-        f_ll = flux(u_ll, orientation, equations)
-        f_rr = flux(u_rr, orientation, 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)
-    end
-
-    return SVector(f1, f2, f3)
+    return SsL, SsR
 end
 
 @inline function max_abs_speeds(u, equations::CompressibleEulerEquations1D)

diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl
@@ -1253,7 +1253,7 @@ function flux_hllc(u_ll, u_rr, orientation::Integer,
 end
 
 """
-    flux_hlle(u_ll, u_rr, orientation, equations::CompressibleEulerEquations2D)
+    min_max_speed_einfeldt(u_ll, u_rr, orientation, 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
@@ -1267,8 +1267,8 @@ of the numerical flux.
   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, orientation::Integer,
-                   equations::CompressibleEulerEquations2D)
+@inline function min_max_speed_einfeldt(u_ll, u_rr, orientation::Integer,
+                                        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)
@@ -1306,39 +1306,65 @@ function flux_hlle(u_ll, u_rr, orientation::Integer,
         SsR = max(v2_roe + c_roe, v2_rr + beta * c_rr, zero(v2_roe))
     end
 
-    if SsL >= 0.0 && SsR > 0.0
-        # Positive supersonic speed
-        f_ll = flux(u_ll, orientation, equations)
+    return SsL, SsR
+end
 
-        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, orientation, equations)
+"""
+    min_max_speed_einfeldt(u_ll, u_rr, normal_direction, equations::CompressibleEulerEquations2D)
 
-        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, orientation, equations)
-        f_rr = flux(u_rr, orientation, 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
+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.
 
-    return SVector(f1, f2, f3, f4)
+- 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)
+"""
+@inline function min_max_speed_einfeldt(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))
+
+    return SsL, SsR
 end
 
 @inline function max_abs_speeds(u, equations::CompressibleEulerEquations2D)

diff --git a/src/equations/compressible_euler_3d.jl b/src/equations/compressible_euler_3d.jl
@@ -1322,7 +1322,7 @@ function flux_hllc(u_ll, u_rr, orientation::Integer,
 end
 
 """
-    flux_hlle(u_ll, u_rr, orientation, equations::CompressibleEulerEquations3D)
+    min_max_speed_einfeldt(u_ll, u_rr, orientation, 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
@@ -1336,8 +1336,8 @@ of the numerical flux.
   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, orientation::Integer,
-                   equations::CompressibleEulerEquations3D)
+@inline function min_max_speed_einfeldt(u_ll, u_rr, orientation::Integer,
+                                        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)
@@ -1379,43 +1379,69 @@ function flux_hlle(u_ll, u_rr, orientation::Integer,
         SsR = max(v3_roe + c_roe, v3_rr + beta * c_rr, zero(v3_roe))
     end
 
-    if SsL >= 0.0 && SsR > 0.0
-        # Positive supersonic speed
-        f_ll = flux(u_ll, orientation, equations)
+    return SsL, SsR
+end
 
-        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, orientation, equations)
+"""
+    min_max_speed_einfeldt(u_ll, u_rr, normal_direction, equations::CompressibleEulerEquations3D)
 
-        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, orientation, equations)
-        f_rr = flux(u_rr, orientation, 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
+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.
 
-    return SVector(f1, f2, f3, f4, f5)
+- 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)
+"""
+@inline function min_max_speed_einfeldt(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))
+
+    return SsL, SsR
 end
 
 @inline function max_abs_speeds(u, equations::CompressibleEulerEquations3D)

diff --git a/src/equations/numerical_fluxes.jl b/src/equations/numerical_fluxes.jl
@@ -304,6 +304,14 @@ See [`FluxHLL`](@ref).
 """
 const flux_hll = FluxHLL()
 
+"""
+    flux_hlle
+
+See [`min_max_speed_einfeldt`](@ref).
+This is a [`FluxHLL`](@ref)-type two-wave solver with special estimates of the wave speeds.
+"""
+const flux_hlle = FluxHLL(min_max_speed_einfeldt)
+
 # TODO: TrixiShallowWater: move the chen_noelle flux structure to the new package
 
 # An empty version of the `min_max_speed_chen_noelle` function is declared here

diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl
@@ -200,6 +200,32 @@ end
     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=[

diff --git a/test/test_p4est_3d.jl b/test/test_p4est_3d.jl
@@ -287,6 +287,34 @@ end
     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"),