Hll 2 wave improvements non breaking (#1561)

* Add Classical and Naive HLL 2 Wave solver to classic Hyperbolic PDEs * Format Code * HLLE wave speeds for SWE * Fix typos * Update tests for HLL * Unit test 1D MHD HLL, HLLE * Add example for classical HLL 2 wave * remove plots * Use lowercase for flux * Use einfeldt for mhd * Use hlle for mhd tets * Missing comma causes failing tests * Correct bug in SWE 2D Roe eigval comp, unit tests * format * Revert "format" This reverts commit 047a5e7. * format equations * Add unit tests for HLL naive * Revert default hll flux * Rename min_max_speed to min_max_speed_davis and reduce documentation * Update src/equations/shallow_water_1d.jl: Comments Co-authored-by: Hendrik Ranocha <[email protected]> * Add published resource for Roe averages for SWE * Add tests for rotation * Remove breaking portionv from PR * fix copy paste error * Lowercase davis * Update src/equations/numerical_fluxes.jl Co-authored-by: Hendrik Ranocha <[email protected]> * Update src/equations/numerical_fluxes.jl Co-authored-by: Hendrik Ranocha <[email protected]> * Update src/equations/numerical_fluxes.jl Co-authored-by: Hendrik Ranocha <[email protected]> * Update src/equations/numerical_fluxes.jl Co-authored-by: Hendrik Ranocha <[email protected]> * Update src/equations/numerical_fluxes.jl Co-authored-by: Hendrik Ranocha <[email protected]> * Update src/equations/numerical_fluxes.jl Co-authored-by: Hendrik Ranocha <[email protected]> * Update test/test_tree_2d_mhd.jl Co-authored-by: Hendrik Ranocha <[email protected]> * Update src/equations/ideal_glm_mhd_1d.jl Co-authored-by: Hendrik Ranocha <[email protected]> * Update src/equations/ideal_glm_mhd_2d.jl Co-authored-by: Hendrik Ranocha <[email protected]> * Update src/equations/ideal_glm_mhd_3d.jl Co-authored-by: Hendrik Ranocha <[email protected]> * Update test/test_tree_3d_mhd.jl Co-authored-by: Hendrik Ranocha <[email protected]> * Remove hll_davis test * Split consistency checks * Try to resolve conflict with 5ff677c * Add tests * More tests --------- Co-authored-by: Hendrik Ranocha <[email protected]>
trixi-framework · Jul 13, 2023 · a191e39 · a191e39
1 parent 42732db
commit a191e39
Show file tree

Hide file tree

Showing 21 changed files with 946 additions and 31 deletions.
diff --git a/examples/dgmulti_2d/elixir_euler_bilinear.jl b/examples/dgmulti_2d/elixir_euler_bilinear.jl
@@ -2,7 +2,7 @@
 using Trixi, OrdinaryDiffEq
 
 dg = DGMulti(polydeg = 3, element_type = Quad(), approximation_type = SBP(),
-             surface_integral = SurfaceIntegralWeakForm(FluxHLL()),
+             surface_integral = SurfaceIntegralWeakForm(flux_hll),
              volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha))
 
 equations = CompressibleEulerEquations2D(1.4)

diff --git a/examples/dgmulti_2d/elixir_euler_curved.jl b/examples/dgmulti_2d/elixir_euler_curved.jl
@@ -2,7 +2,7 @@
 using Trixi, OrdinaryDiffEq
 
 dg = DGMulti(polydeg = 3, element_type = Quad(), approximation_type = SBP(),
-             surface_integral = SurfaceIntegralWeakForm(FluxHLL()),
+             surface_integral = SurfaceIntegralWeakForm(flux_hll),
              volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha))
 
 equations = CompressibleEulerEquations2D(1.4)

diff --git a/examples/dgmulti_2d/elixir_euler_triangulate_pkg_mesh.jl b/examples/dgmulti_2d/elixir_euler_triangulate_pkg_mesh.jl
@@ -1,7 +1,7 @@
 using Trixi, OrdinaryDiffEq
 
 dg = DGMulti(polydeg = 3, element_type = Tri(),
-             surface_integral = SurfaceIntegralWeakForm(FluxHLL()),
+             surface_integral = SurfaceIntegralWeakForm(flux_hll),
              volume_integral = VolumeIntegralWeakForm())
 
 equations = CompressibleEulerEquations2D(1.4)

diff --git a/examples/dgmulti_2d/elixir_euler_weakform.jl b/examples/dgmulti_2d/elixir_euler_weakform.jl
@@ -2,7 +2,7 @@
 using Trixi, OrdinaryDiffEq
 
 dg = DGMulti(polydeg = 3, element_type = Tri(), approximation_type = Polynomial(),
-             surface_integral = SurfaceIntegralWeakForm(FluxHLL()),
+             surface_integral = SurfaceIntegralWeakForm(flux_hll),
              volume_integral = VolumeIntegralWeakForm())
 
 equations = CompressibleEulerEquations2D(1.4)

diff --git a/examples/dgmulti_2d/elixir_euler_weakform_periodic.jl b/examples/dgmulti_2d/elixir_euler_weakform_periodic.jl
@@ -2,7 +2,7 @@
 using Trixi, OrdinaryDiffEq
 
 dg = DGMulti(polydeg = 3, element_type = Tri(), approximation_type = Polynomial(),
-             surface_integral = SurfaceIntegralWeakForm(FluxHLL()),
+             surface_integral = SurfaceIntegralWeakForm(flux_hll),
              volume_integral = VolumeIntegralWeakForm())
 
 equations = CompressibleEulerEquations2D(1.4)

diff --git a/examples/dgmulti_3d/elixir_euler_curved.jl b/examples/dgmulti_3d/elixir_euler_curved.jl
@@ -2,7 +2,7 @@
 using Trixi, OrdinaryDiffEq
 
 dg = DGMulti(polydeg = 3, element_type = Hex(), approximation_type=SBP(),
-             surface_integral = SurfaceIntegralWeakForm(FluxHLL()),
+             surface_integral = SurfaceIntegralWeakForm(flux_hll),
              volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha))
 
 equations = CompressibleEulerEquations3D(1.4)

diff --git a/examples/dgmulti_3d/elixir_euler_weakform.jl b/examples/dgmulti_3d/elixir_euler_weakform.jl
@@ -2,7 +2,7 @@
 using Trixi, OrdinaryDiffEq
 
 dg = DGMulti(polydeg = 3, element_type = Tet(),
-             surface_integral = SurfaceIntegralWeakForm(FluxHLL()),
+             surface_integral = SurfaceIntegralWeakForm(flux_hll),
              volume_integral = VolumeIntegralWeakForm())
 
 equations = CompressibleEulerEquations3D(1.4)

diff --git a/examples/dgmulti_3d/elixir_euler_weakform_periodic.jl b/examples/dgmulti_3d/elixir_euler_weakform_periodic.jl
@@ -2,7 +2,7 @@
 using Trixi, OrdinaryDiffEq
 
 dg = DGMulti(polydeg = 3, element_type = Tet(), approximation_type = Polynomial(),
-             surface_integral = SurfaceIntegralWeakForm(FluxHLL()),
+             surface_integral = SurfaceIntegralWeakForm(flux_hll),
              volume_integral = VolumeIntegralWeakForm())
 
 equations = CompressibleEulerEquations3D(1.4)

diff --git a/src/Trixi.jl b/src/Trixi.jl
@@ -164,7 +164,7 @@ export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle,
        hydrostatic_reconstruction_audusse_etal, flux_nonconservative_audusse_etal,
        FluxPlusDissipation, DissipationGlobalLaxFriedrichs, DissipationLocalLaxFriedrichs,
        FluxLaxFriedrichs, max_abs_speed_naive,
-       FluxHLL, min_max_speed_naive,
+       FluxHLL, min_max_speed_naive, min_max_speed_davis, min_max_speed_einfeldt,
        FluxLMARS,
        FluxRotated,
        flux_shima_etal_turbo, flux_ranocha_turbo,

diff --git a/src/equations/compressible_euler_1d.jl b/src/equations/compressible_euler_1d.jl
@@ -628,7 +628,7 @@ end
     return SVector(f1m, f2m, f3m)
 end
 
-# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation as the
+# Calculate estimates for maximum wave speed for local Lax-Friedrichs-type dissipation as the
 # maximum velocity magnitude plus the maximum speed of sound
 @inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer,
                                      equations::CompressibleEulerEquations1D)
@@ -648,7 +648,7 @@ end
     λ_max = max(v_mag_ll, v_mag_rr) + max(c_ll, c_rr)
 end
 
-# Calculate minimum and maximum wave speeds for HLL-type fluxes
+# Calculate estimates for minimum and maximum wave speeds for HLL-type fluxes
 @inline function min_max_speed_naive(u_ll, u_rr, orientation::Integer,
                                      equations::CompressibleEulerEquations1D)
     rho_ll, v1_ll, p_ll = cons2prim(u_ll, equations)
@@ -660,6 +660,21 @@ end
     return λ_min, λ_max
 end
 
+# More refined estimates for minimum and maximum wave speeds for HLL-type fluxes
+@inline function min_max_speed_davis(u_ll, u_rr, orientation::Integer,
+                                     equations::CompressibleEulerEquations1D)
+    rho_ll, v1_ll, p_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, p_rr = cons2prim(u_rr, equations)
+
+    c_ll = sqrt(equations.gamma * p_ll / rho_ll)
+    c_rr = sqrt(equations.gamma * p_rr / rho_rr)
+
+    λ_min = min(v1_ll - c_ll, v1_rr - c_rr)
+    λ_max = max(v1_ll + c_ll, v1_rr + c_rr)
+
+    return λ_min, λ_max
+end
+
 """
     flux_hllc(u_ll, u_rr, orientation, equations::CompressibleEulerEquations1D)
 

diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl
@@ -1032,7 +1032,7 @@ end
     return max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) * norm(normal_direction)
 end
 
-# Calculate minimum and maximum wave speeds for HLL-type fluxes
+# Calculate estimate for minimum and maximum wave speeds for HLL-type fluxes
 @inline function min_max_speed_naive(u_ll, u_rr, orientation::Integer,
                                      equations::CompressibleEulerEquations2D)
     rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations)
@@ -1065,6 +1065,47 @@ end
     return λ_min, λ_max
 end
 
+# More refined estimates for minimum and maximum wave speeds for HLL-type fluxes
+@inline function min_max_speed_davis(u_ll, u_rr, orientation::Integer,
+                                     equations::CompressibleEulerEquations2D)
+    rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations)
+
+    c_ll = sqrt(equations.gamma * p_ll / rho_ll)
+    c_rr = sqrt(equations.gamma * p_rr / rho_rr)
+
+    if orientation == 1 # x-direction
+        λ_min = min(v1_ll - c_ll, v1_rr - c_rr)
+        λ_max = max(v1_ll + c_ll, v1_rr + c_rr)
+    else # y-direction
+        λ_min = min(v2_ll - c_ll, v2_rr - c_rr)
+        λ_max = max(v2_ll + c_ll, v2_rr + c_rr)
+    end
+
+    return λ_min, λ_max
+end
+
+# More refined estimates for minimum and maximum wave speeds for HLL-type fluxes
+@inline function min_max_speed_davis(u_ll, u_rr, normal_direction::AbstractVector,
+                                     equations::CompressibleEulerEquations2D)
+    rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations)
+
+    norm_ = norm(normal_direction)
+
+    c_ll = sqrt(equations.gamma * p_ll / rho_ll) * norm_
+    c_rr = sqrt(equations.gamma * p_rr / rho_rr) * norm_
+
+    v_normal_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2]
+    v_normal_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2]
+
+    # The v_normals are already scaled by the norm
+    λ_min = min(v_normal_ll - c_ll, v_normal_rr - c_rr)
+    λ_max = max(v_normal_ll + c_ll, v_normal_rr + c_rr)
+
+    return λ_min, λ_max
+end
+
 # Called inside `FluxRotated` in `numerical_fluxes.jl` so the direction
 # has been normalized prior to this rotation of the state vector
 @inline function rotate_to_x(u, normal_vector, equations::CompressibleEulerEquations2D)

diff --git a/src/equations/compressible_euler_3d.jl b/src/equations/compressible_euler_3d.jl
@@ -1070,7 +1070,7 @@ end
     return max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) * norm(normal_direction)
 end
 
-# Calculate minimum and maximum wave speeds for HLL-type fluxes
+# Calculate estimates for minimum and maximum wave speeds for HLL-type fluxes
 @inline function min_max_speed_naive(u_ll, u_rr, orientation::Integer,
                                      equations::CompressibleEulerEquations3D)
     rho_ll, v1_ll, v2_ll, v3_ll, p_ll = cons2prim(u_ll, equations)
@@ -1108,6 +1108,54 @@ end
     return λ_min, λ_max
 end
 
+# More refined estimates for minimum and maximum wave speeds for HLL-type fluxes
+@inline function min_max_speed_davis(u_ll, u_rr, orientation::Integer,
+                                     equations::CompressibleEulerEquations3D)
+    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)
+
+    c_ll = sqrt(equations.gamma * p_ll / rho_ll)
+    c_rr = sqrt(equations.gamma * p_rr / rho_rr)
+
+    if orientation == 1 # x-direction
+        λ_min = min(v1_ll - c_ll, v1_rr - c_rr)
+        λ_max = max(v1_ll + c_ll, v1_rr + c_rr)
+    elseif orientation == 2 # y-direction
+        λ_min = min(v2_ll - c_ll, v2_rr - c_rr)
+        λ_max = max(v2_ll + c_ll, v2_rr + c_rr)
+    else # z-direction
+        λ_min = min(v3_ll - c_ll, v3_rr - c_rr)
+        λ_max = max(v3_ll + c_ll, v3_rr + c_rr)
+    end
+
+    return λ_min, λ_max
+end
+
+# More refined estimates for minimum and maximum wave speeds for HLL-type fluxes
+@inline function min_max_speed_davis(u_ll, u_rr, normal_direction::AbstractVector,
+                                     equations::CompressibleEulerEquations3D)
+    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)
+
+    norm_ = norm(normal_direction)
+
+    c_ll = sqrt(equations.gamma * p_ll / rho_ll) * norm_
+    c_rr = sqrt(equations.gamma * p_rr / rho_rr) * norm_
+
+    v_normal_ll = v1_ll * normal_direction[1] +
+                  v2_ll * normal_direction[2] +
+                  v3_ll * normal_direction[3]
+    v_normal_rr = v1_rr * normal_direction[1] +
+                  v2_rr * normal_direction[2] +
+                  v3_rr * normal_direction[3]
+
+    # The v_normals are already scaled by the norm
+    λ_min = min(v_normal_ll - c_ll, v_normal_rr - c_rr)
+    λ_max = max(v_normal_ll + c_ll, v_normal_rr + c_rr)
+
+    return λ_min, λ_max
+end
+
 # Rotate normal vector to x-axis; normal, tangent1 and tangent2 need to be orthonormal
 # Called inside `FluxRotated` in `numerical_fluxes.jl` so the directions
 # has been normalized prior to this rotation of the state vector

diff --git a/src/equations/ideal_glm_mhd_1d.jl b/src/equations/ideal_glm_mhd_1d.jl
@@ -277,13 +277,33 @@ end
     λ_max = max(abs(v_ll), abs(v_rr)) + max(cf_ll, cf_rr)
 end
 
+# More refined estimates for minimum and maximum wave speeds for HLL-type fluxes
+@inline function min_max_speed_davis(u_ll, u_rr, orientation::Integer,
+                                     equations::IdealGlmMhdEquations1D)
+    rho_ll, rho_v1_ll, _ = u_ll
+    rho_rr, rho_v1_rr, _ = u_rr
+
+    # Calculate primitive variables
+    v1_ll = rho_v1_ll / rho_ll
+    v1_rr = rho_v1_rr / rho_rr
+
+    # Approximate the left-most and right-most eigenvalues in the Riemann fan
+    c_f_ll = calc_fast_wavespeed(u_ll, orientation, equations)
+    c_f_rr = calc_fast_wavespeed(u_rr, orientation, equations)
+
+    λ_min = min(v1_ll - c_f_ll, v1_rr - c_f_rr)
+    λ_max = max(v1_ll + c_f_ll, v1_rr + c_f_rr)
+
+    return λ_min, λ_max
+end
+
 """
-    min_max_speed_naive(u_ll, u_rr, orientation, equations::IdealGlmMhdEquations1D)
+    min_max_speed_naive(u_ll, u_rr, orientation::Integer, equations::IdealGlmMhdEquations1D)
 
 Calculate minimum and maximum wave speeds for HLL-type fluxes as in
 - Li (2005)
   An HLLC Riemann solver for magneto-hydrodynamics
-  [DOI: 10.1016/j.jcp.2004.08.020](https://doi.org/10.1016/j.jcp.2004.08.020)
+  [DOI: 10.1016/j.jcp.2004.08.020](https://doi.org/10.1016/j.jcp.2004.08.020).
 """
 @inline function min_max_speed_naive(u_ll, u_rr, orientation::Integer,
                                      equations::IdealGlmMhdEquations1D)

diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl
@@ -585,13 +585,70 @@ end
     return max(abs(v_ll), abs(v_rr)) + max(cf_ll, cf_rr)
 end
 
+# More refined estimates for minimum and maximum wave speeds for HLL-type fluxes
+@inline function min_max_speed_davis(u_ll, u_rr, orientation::Integer,
+                                     equations::IdealGlmMhdEquations2D)
+    rho_ll, rho_v1_ll, rho_v2_ll, _ = u_ll
+    rho_rr, rho_v1_rr, rho_v2_rr, _ = u_rr
+
+    # Calculate primitive velocity variables
+    v1_ll = rho_v1_ll / rho_ll
+    v2_ll = rho_v2_ll / rho_ll
+
+    v1_rr = rho_v1_rr / rho_rr
+    v2_rr = rho_v2_rr / rho_rr
+
+    # Approximate the left-most and right-most eigenvalues in the Riemann fan
+    if orientation == 1 # x-direction
+        c_f_ll = calc_fast_wavespeed(u_ll, orientation, equations)
+        c_f_rr = calc_fast_wavespeed(u_rr, orientation, equations)
+
+        λ_min = min(v1_ll - c_f_ll, v1_rr - c_f_rr)
+        λ_max = max(v1_ll + c_f_ll, v1_rr + c_f_rr)
+    else # y-direction
+        c_f_ll = calc_fast_wavespeed(u_ll, orientation, equations)
+        c_f_rr = calc_fast_wavespeed(u_rr, orientation, equations)
+
+        λ_min = min(v2_ll - c_f_ll, v2_rr - c_f_rr)
+        λ_max = max(v2_ll + c_f_ll, v1_rr + c_f_rr)
+    end
+
+    return λ_min, λ_max
+end
+
+# More refined estimates for minimum and maximum wave speeds for HLL-type fluxes
+@inline function min_max_speed_davis(u_ll, u_rr, normal_direction::AbstractVector,
+                                     equations::IdealGlmMhdEquations2D)
+    rho_ll, rho_v1_ll, rho_v2_ll, _ = u_ll
+    rho_rr, rho_v1_rr, rho_v2_rr, _ = u_rr
+
+    # Calculate primitive velocity variables
+    v1_ll = rho_v1_ll / rho_ll
+    v2_ll = rho_v2_ll / rho_ll
+
+    v1_rr = rho_v1_rr / rho_rr
+    v2_rr = rho_v2_rr / rho_rr
+
+    v_normal_ll = (v1_ll * normal_direction[1] + v2_ll * normal_direction[2])
+    v_normal_rr = (v1_rr * normal_direction[1] + v2_rr * normal_direction[2])
+
+    c_f_ll = calc_fast_wavespeed(u_ll, normal_direction, equations)
+    c_f_rr = calc_fast_wavespeed(u_rr, normal_direction, equations)
+
+    # Estimate the min/max eigenvalues in the normal direction
+    λ_min = min(v_normal_ll - c_f_ll, v_normal_rr - c_f_rr)
+    λ_max = max(v_normal_ll + c_f_ll, v_normal_rr + c_f_rr)
+
+    return λ_min, λ_max
+end
+
 """
-    min_max_speed_naive(u_ll, u_rr, orientation, equations::IdealGlmMhdEquations2D)
+    min_max_speed_naive(u_ll, u_rr, orientation::Integer, equations::IdealGlmMhdEquations2D)
 
 Calculate minimum and maximum wave speeds for HLL-type fluxes as in
 - Li (2005)
   An HLLC Riemann solver for magneto-hydrodynamics
-  [DOI: 10.1016/j.jcp.2004.08.020](https://doi.org/10.1016/j.jcp.2004.08.020)
+  [DOI: 10.1016/j.jcp.2004.08.020](https://doi.org/10.1016/j.jcp.2004.08.020).
 """
 @inline function min_max_speed_naive(u_ll, u_rr, orientation::Integer,
                                      equations::IdealGlmMhdEquations2D)
@@ -635,10 +692,8 @@ end
     v1_rr = rho_v1_rr / rho_rr
     v2_rr = rho_v2_rr / rho_rr
 
-    v_normal_ll = (v1_ll * normal_direction[1] +
-                   v2_ll * normal_direction[2])
-    v_normal_rr = (v1_rr * normal_direction[1] +
-                   v2_rr * normal_direction[2])
+    v_normal_ll = (v1_ll * normal_direction[1] + v2_ll * normal_direction[2])
+    v_normal_rr = (v1_rr * normal_direction[1] + v2_rr * normal_direction[2])
 
     c_f_ll = calc_fast_wavespeed(u_ll, normal_direction, equations)
     c_f_rr = calc_fast_wavespeed(u_rr, normal_direction, equations)