Hll 2 wave improvements

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)

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)

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)

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)

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)

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)

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)

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)

examples/p4est_2d_dgsem/elixir_mhd_alfven_wave.jl

@@ -13,7 +13,7 @@ initial_condition = initial_condition_convergence_test
 
 # Get the DG approximation space
 volume_flux = (flux_central, flux_nonconservative_powell)
-solver = DGSEM(polydeg=4, surface_flux=(flux_hll, flux_nonconservative_powell),
+solver = DGSEM(polydeg=4, surface_flux=(FluxHLL(min_max_speed_einfeldt), flux_nonconservative_powell),
                volume_integral=VolumeIntegralFluxDifferencing(volume_flux))
 
 coordinates_min = (0.0      , 0.0      )

examples/p4est_3d_dgsem/elixir_mhd_alfven_wave_nonconforming.jl

@@ -10,7 +10,7 @@ equations = IdealGlmMhdEquations3D(5/3)
 initial_condition = initial_condition_convergence_test
 
 volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell)
-solver = DGSEM(polydeg=3, surface_flux=(flux_hll, flux_nonconservative_powell),
+solver = DGSEM(polydeg=3, surface_flux=(FluxHLL(min_max_speed_einfeldt), flux_nonconservative_powell),
                volume_integral=VolumeIntegralFluxDifferencing(volume_flux))
 
 coordinates_min = (-1.0, -1.0, -1.0)

examples/structured_1d_dgsem/elixir_euler_sedov_hll_Davis.jl

@@ -0,0 +1,93 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations1D(1.4)
+
+"""
+    initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations1D)
+
+The Sedov blast wave setup based on Flash
+- http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000
+"""
+function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations1D)
+  # Set up polar coordinates
+  inicenter = SVector(0.0)
+  x_norm = x[1] - inicenter[1]
+  r = abs(x_norm)
+
+  # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000
+  r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6)
+  # r0 = 0.5 # = more reasonable setup
+  E = 1.0
+  p0_inner = 6 * (equations.gamma - 1) * E / (3 * pi * r0)
+  p0_outer = 1.0e-5 # = true Sedov setup
+
+  # Calculate primitive variables
+  rho = 1.0
+  v1  = 0.0
+  p   = r > r0 ? p0_outer : p0_inner
+
+  return prim2cons(SVector(rho, v1, p), equations)
+end
+
+initial_condition = initial_condition_sedov_blast_wave
+surface_flux = FluxHLL(min_max_speed_davis)
+volume_flux  = flux_ranocha
+basis = LobattoLegendreBasis(3)
+shock_indicator_variable = density_pressure
+indicator_sc = IndicatorHennemannGassner(equations, basis,
+                                         alpha_max=1.0,
+                                         alpha_min=0.001,
+                                         alpha_smooth=true,
+                                         variable=shock_indicator_variable)
+volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+                                                 volume_flux_dg=volume_flux,
+                                                 volume_flux_fv=surface_flux)
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+cells_per_dimension = (64,)
+coordinates_min = (-2.0,)
+coordinates_max = ( 2.0,)
+mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max, periodicity=true)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 12.5)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+save_solution = SaveSolutionCallback(interval=100,
+                                     save_initial_solution=true,
+                                     save_final_solution=true,
+                                     solution_variables=cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl=0.5)
+
+callbacks = CallbackSet(summary_callback, 
+                        analysis_callback, 
+                        alive_callback,
+                        save_solution,
+                        stepsize_callback)
+
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt=stepsize_callback(ode), # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep=false, callback=callbacks);
+summary_callback() # print the timer summary

examples/structured_3d_dgsem/elixir_mhd_alfven_wave.jl

@@ -10,7 +10,7 @@ equations = IdealGlmMhdEquations3D(5/3)
 initial_condition = initial_condition_convergence_test
 
 volume_flux = (flux_central, flux_nonconservative_powell)
-solver = DGSEM(polydeg=5, surface_flux=(flux_hll, flux_nonconservative_powell),
+solver = DGSEM(polydeg=5, surface_flux=(FluxHLL(min_max_speed_einfeldt), flux_nonconservative_powell),
                volume_integral=VolumeIntegralFluxDifferencing(volume_flux))
 
 # Create the mesh

examples/tree_1d_dgsem/elixir_mhd_briowu_shock_tube.jl

@@ -32,7 +32,7 @@ initial_condition = initial_condition_briowu_shock_tube
 
 boundary_conditions = BoundaryConditionDirichlet(initial_condition)
 
-surface_flux = flux_hll
+surface_flux = FluxHLL(min_max_speed_einfeldt)
 volume_flux  = flux_derigs_etal
 basis = LobattoLegendreBasis(4)
 

examples/tree_1d_dgsem/elixir_mhd_ryujones_shock_tube.jl

@@ -39,7 +39,7 @@ initial_condition = initial_condition_ryujones_shock_tube
 
 boundary_conditions = BoundaryConditionDirichlet(initial_condition)
 
-surface_flux = flux_hll
+surface_flux = FluxHLL(min_max_speed_einfeldt)
 volume_flux  = flux_hindenlang_gassner
 basis = LobattoLegendreBasis(3)
 

examples/tree_1d_dgsem/elixir_mhd_shu_osher_shock_tube.jl

@@ -60,7 +60,7 @@ initial_condition = initial_condition_shu_osher_shock_tube
 
 boundary_conditions = BoundaryConditionDirichlet(initial_condition)
 
-surface_flux = flux_hll
+surface_flux = FluxHLL(min_max_speed_einfeldt)
 volume_flux  = flux_hindenlang_gassner
 basis = LobattoLegendreBasis(4)
 

examples/tree_2d_dgsem/elixir_mhd_alfven_wave_mortar.jl

@@ -10,7 +10,7 @@ equations = IdealGlmMhdEquations2D(gamma)
 initial_condition = initial_condition_convergence_test
 
 volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell)
-solver = DGSEM(polydeg=3, surface_flux=(flux_hll, flux_nonconservative_powell),
+solver = DGSEM(polydeg=3, surface_flux=(FluxHLL(min_max_speed_einfeldt), flux_nonconservative_powell),
                volume_integral=VolumeIntegralFluxDifferencing(volume_flux))
 
 coordinates_min = (0.0, 0.0)

examples/tree_3d_dgsem/elixir_mhd_alfven_wave_mortar.jl

@@ -10,7 +10,7 @@ equations = IdealGlmMhdEquations3D(5/3)
 initial_condition = initial_condition_convergence_test
 
 volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell)
-solver = DGSEM(polydeg=3, surface_flux=(flux_hll, flux_nonconservative_powell),
+solver = DGSEM(polydeg=3, surface_flux=(FluxHLL(min_max_speed_einfeldt), flux_nonconservative_powell),
                volume_integral=VolumeIntegralFluxDifferencing(volume_flux))
 
 coordinates_min = (-1.0, -1.0, -1.0)

examples/unstructured_2d_dgsem/elixir_mhd_alfven_wave.jl

@@ -11,7 +11,7 @@ equations = IdealGlmMhdEquations2D(gamma)
 initial_condition = initial_condition_convergence_test
 
 volume_flux = (flux_central, flux_nonconservative_powell)
-solver = DGSEM(polydeg=7, surface_flux=(flux_hll, flux_nonconservative_powell),
+solver = DGSEM(polydeg=7, surface_flux=(FluxHLL(min_max_speed_einfeldt), flux_nonconservative_powell),
                volume_integral=VolumeIntegralFluxDifferencing(volume_flux))
 
 # Get the unstructured quad mesh from a file (downloads the file if not available locally)

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,

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)
 

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)