diff --git a/NEWS.md b/NEWS.md
index ecc91581e9a..d70504d8c85 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -10,6 +10,9 @@ for human readability.
 
 #### Changed
 
+- The default wave speed estimate used within `flux_hll` is now `min_max_speed_davis` 
+  instead of `min_max_speed_naive`.
+
 #### Deprecated
 
 #### Removed
@@ -17,6 +20,7 @@ for human readability.
   Trixi.jl, but are moved to a dedicated repository: [TrixiShallowWater.jl](https://github.com/trixi-framework/TrixiShallowWater.jl). This includes all features related to wetting and drying, as well as the `ShallowWaterTwoLayerEquations1D` and `ShallowWaterTwoLayerEquations2D`.
   However, the basic shallow water equations are still part of Trixi.jl. We'll also be updating the TrixiShallowWater.jl documentation with instructions on how to use these relocated features in the future.
 
+
 ## Changes in the v0.6 lifecycle
 
 #### Added
@@ -27,6 +31,7 @@ for human readability.
 - Subcell (positivity) limiting support for nonlinear variables in 2D for `TreeMesh`
 - Added Lighthill-Whitham-Richards (LWR) traffic model
 
+
 ## Changes when updating to v0.6 from v0.5.x
 
 #### Added
diff --git a/examples/p4est_2d_dgsem/elixir_eulergravity_convergence.jl b/examples/p4est_2d_dgsem/elixir_eulergravity_convergence.jl
index d55a59ca5ce..974466e3b3b 100644
--- a/examples/p4est_2d_dgsem/elixir_eulergravity_convergence.jl
+++ b/examples/p4est_2d_dgsem/elixir_eulergravity_convergence.jl
@@ -10,7 +10,7 @@ gamma = 2.0
 equations_euler = CompressibleEulerEquations2D(gamma)
 
 polydeg = 3
-solver_euler = DGSEM(polydeg, flux_hll)
+solver_euler = DGSEM(polydeg, FluxHLL(min_max_speed_naive))
 
 coordinates_min = (0.0, 0.0)
 coordinates_max = (2.0, 2.0)
diff --git a/examples/p4est_3d_dgsem/elixir_euler_source_terms_nonconforming_earth.jl b/examples/p4est_3d_dgsem/elixir_euler_source_terms_nonconforming_earth.jl
index 28a300cd681..28cdec12da5 100644
--- a/examples/p4est_3d_dgsem/elixir_euler_source_terms_nonconforming_earth.jl
+++ b/examples/p4est_3d_dgsem/elixir_euler_source_terms_nonconforming_earth.jl
@@ -68,7 +68,7 @@ boundary_condition = BoundaryConditionDirichlet(initial_condition)
 boundary_conditions = Dict(:inside => boundary_condition,
                            :outside => boundary_condition)
 
-surface_flux = flux_hll
+surface_flux = FluxHLL(min_max_speed_naive)
 # Note that a free stream is not preserved if N < 2 * N_geo, where N is the
 # polydeg of the solver and N_geo is the polydeg of the mesh.
 # However, the FSP error is negligible in this example.
diff --git a/examples/paper_self_gravitating_gas_dynamics/elixir_euler_convergence.jl b/examples/paper_self_gravitating_gas_dynamics/elixir_euler_convergence.jl
index aabfce0f66b..4f44d7b12ac 100644
--- a/examples/paper_self_gravitating_gas_dynamics/elixir_euler_convergence.jl
+++ b/examples/paper_self_gravitating_gas_dynamics/elixir_euler_convergence.jl
@@ -8,7 +8,7 @@ equations = CompressibleEulerEquations2D(2.0)
 
 initial_condition = initial_condition_eoc_test_coupled_euler_gravity
 
-solver = DGSEM(polydeg = 3, surface_flux = flux_hll)
+solver = DGSEM(polydeg = 3, surface_flux = FluxHLL(min_max_speed_naive))
 
 coordinates_min = (0.0, 0.0)
 coordinates_max = (2.0, 2.0)
diff --git a/examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_convergence.jl b/examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_convergence.jl
index ce1d2cd05bd..49b98803577 100644
--- a/examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_convergence.jl
+++ b/examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_convergence.jl
@@ -10,7 +10,7 @@ gamma = 2.0
 equations_euler = CompressibleEulerEquations2D(gamma)
 
 polydeg = 3
-solver_euler = DGSEM(polydeg, flux_hll)
+solver_euler = DGSEM(polydeg, FluxHLL(min_max_speed_naive))
 
 coordinates_min = (0.0, 0.0)
 coordinates_max = (2.0, 2.0)
diff --git a/examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_jeans_instability.jl b/examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_jeans_instability.jl
index f081f6bb91a..7461198fbb2 100644
--- a/examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_jeans_instability.jl
+++ b/examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_jeans_instability.jl
@@ -66,7 +66,7 @@ gamma = 5 / 3
 equations_euler = CompressibleEulerEquations2D(gamma)
 
 polydeg = 3
-solver_euler = DGSEM(polydeg, flux_hll)
+solver_euler = DGSEM(polydeg, FluxHLL(min_max_speed_naive))
 
 coordinates_min = (0.0, 0.0)
 coordinates_max = (1.0, 1.0)
diff --git a/examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_sedov_blast_wave.jl b/examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_sedov_blast_wave.jl
index b7be2320228..bc7ceb97c8b 100644
--- a/examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_sedov_blast_wave.jl
+++ b/examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_sedov_blast_wave.jl
@@ -85,7 +85,7 @@ function boundary_condition_sedov_self_gravity(u_inner, orientation, direction,
 end
 boundary_conditions = boundary_condition_sedov_self_gravity
 
-surface_flux = flux_hll
+surface_flux = FluxHLL(min_max_speed_naive)
 volume_flux = flux_chandrashekar
 polydeg = 3
 basis = LobattoLegendreBasis(polydeg)
diff --git a/examples/t8code_2d_dgsem/elixir_eulergravity_convergence.jl b/examples/t8code_2d_dgsem/elixir_eulergravity_convergence.jl
index 98a9a5521a9..cd10315945a 100644
--- a/examples/t8code_2d_dgsem/elixir_eulergravity_convergence.jl
+++ b/examples/t8code_2d_dgsem/elixir_eulergravity_convergence.jl
@@ -9,7 +9,7 @@ gamma = 2.0
 equations_euler = CompressibleEulerEquations2D(gamma)
 
 polydeg = 3
-solver_euler = DGSEM(polydeg, flux_hll)
+solver_euler = DGSEM(polydeg, FluxHLL(min_max_speed_naive))
 
 coordinates_min = (0.0, 0.0)
 coordinates_max = (2.0, 2.0)
diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_nonperiodic.jl b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_nonperiodic.jl
index e55fffc101e..9ed02c0e378 100644
--- a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_nonperiodic.jl
+++ b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_nonperiodic.jl
@@ -26,7 +26,9 @@ boundary_condition = BoundaryConditionDirichlet(initial_condition)
 # Get the DG approximation space
 
 volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
-solver = DGSEM(polydeg = 4, surface_flux = (flux_hll, flux_nonconservative_fjordholm_etal),
+solver = DGSEM(polydeg = 4,
+               surface_flux = (flux_hll,
+                               flux_nonconservative_fjordholm_etal),
                volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
 
 ###############################################################################
diff --git a/examples/tree_3d_dgsem/elixir_eulergravity_convergence.jl b/examples/tree_3d_dgsem/elixir_eulergravity_convergence.jl
index 21ef661d0b6..0a8c427bf8d 100644
--- a/examples/tree_3d_dgsem/elixir_eulergravity_convergence.jl
+++ b/examples/tree_3d_dgsem/elixir_eulergravity_convergence.jl
@@ -10,7 +10,7 @@ equations_euler = CompressibleEulerEquations3D(gamma)
 initial_condition = initial_condition_eoc_test_coupled_euler_gravity
 
 polydeg = 3
-solver_euler = DGSEM(polydeg, flux_hll)
+solver_euler = DGSEM(polydeg, FluxHLL(min_max_speed_naive))
 
 coordinates_min = (0.0, 0.0, 0.0)
 coordinates_max = (2.0, 2.0, 2.0)
diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_dirichlet.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_dirichlet.jl
index df1a69192ce..38e1279e220 100644
--- a/examples/unstructured_2d_dgsem/elixir_shallowwater_dirichlet.jl
+++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_dirichlet.jl
@@ -30,7 +30,9 @@ boundary_condition = Dict(:OuterCircle => boundary_condition_constant)
 # Get the DG approximation space
 
 volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
-solver = DGSEM(polydeg = 4, surface_flux = (flux_hll, flux_nonconservative_fjordholm_etal),
+solver = DGSEM(polydeg = 4,
+               surface_flux = (flux_hll,
+                               flux_nonconservative_fjordholm_etal),
                volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
 
 ###############################################################################
diff --git a/src/equations/numerical_fluxes.jl b/src/equations/numerical_fluxes.jl
index 87fcb412244..6794c71a32b 100644
--- a/src/equations/numerical_fluxes.jl
+++ b/src/equations/numerical_fluxes.jl
@@ -222,12 +222,12 @@ See [`FluxLaxFriedrichs`](@ref).
 const flux_lax_friedrichs = FluxLaxFriedrichs()
 
 """
-    FluxHLL(min_max_speed=min_max_speed_naive)
+    FluxHLL(min_max_speed=min_max_speed_davis)
 
 Create an HLL (Harten, Lax, van Leer) numerical flux where the minimum and maximum
 wave speeds are estimated as
 `λ_min, λ_max = min_max_speed(u_ll, u_rr, orientation_or_normal_direction, equations)`,
-defaulting to [`min_max_speed_naive`](@ref).
+defaulting to [`min_max_speed_davis`](@ref).
 Original paper:
 - Amiram Harten, Peter D. Lax, Bram van Leer (1983)
   On Upstream Differencing and Godunov-Type Schemes for Hyperbolic Conservation Laws
@@ -237,7 +237,7 @@ struct FluxHLL{MinMaxSpeed}
     min_max_speed::MinMaxSpeed
 end
 
-FluxHLL() = FluxHLL(min_max_speed_naive)
+FluxHLL() = FluxHLL(min_max_speed_davis)
 
 """
     min_max_speed_naive(u_ll, u_rr, orientation::Integer, equations)
@@ -246,10 +246,16 @@ FluxHLL() = FluxHLL(min_max_speed_naive)
 Simple and fast estimate(!) of the minimal and maximal wave speed of the Riemann problem with
 left and right states `u_ll, u_rr`, usually based only on the local wave speeds associated to
 `u_ll` and `u_rr`.
+Slightly more diffusive than [`min_max_speed_davis`](@ref).
 - Amiram Harten, Peter D. Lax, Bram van Leer (1983)
   On Upstream Differencing and Godunov-Type Schemes for Hyperbolic Conservation Laws
   [DOI: 10.1137/1025002](https://doi.org/10.1137/1025002)
 
+See eq. (10.37) from
+- Eleuterio F. Toro (2009)
+  Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction
+  [DOI: 10.1007/b79761](https://doi.org/10.1007/b79761)
+
 See also [`FluxHLL`](@ref), [`min_max_speed_davis`](@ref), [`min_max_speed_einfeldt`](@ref).
 """
 function min_max_speed_naive end
@@ -266,6 +272,10 @@ left and right states `u_ll, u_rr`, usually based only on the local wave speeds
   Simplified Second-Order Godunov-Type Methods
   [DOI: 10.1137/0909030](https://doi.org/10.1137/0909030)
 
+See eq. (10.38) from
+- Eleuterio F. Toro (2009)
+  Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction
+  [DOI: 10.1007/b79761](https://doi.org/10.1007/b79761)
 See also [`FluxHLL`](@ref), [`min_max_speed_naive`](@ref), [`min_max_speed_einfeldt`](@ref).
 """
 function min_max_speed_davis end
diff --git a/test/test_dgmulti_1d.jl b/test/test_dgmulti_1d.jl
index 0363086341f..e470de71efb 100644
--- a/test/test_dgmulti_1d.jl
+++ b/test/test_dgmulti_1d.jl
@@ -128,14 +128,12 @@ end
 @trixi_testset "elixir_euler_fdsbp_periodic.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_fdsbp_periodic.jl"),
                         l2=[
-                            9.146929180585711e-7,
-                            1.8997616878017292e-6,
-                            3.991417702211889e-6,
+                            9.146929178341782e-7, 1.8997616876521201e-6,
+                            3.991417701005622e-6,
                         ],
                         linf=[
-                            1.7321089884614338e-6,
-                            3.3252888855805907e-6,
-                            6.5252787737613005e-6,
+                            1.7321089882393892e-6, 3.3252888869128583e-6,
+                            6.525278767988141e-6,
                         ])
     show(stdout, semi.solver.basis)
     show(stdout, MIME"text/plain"(), semi.solver.basis)
diff --git a/test/test_dgmulti_2d.jl b/test/test_dgmulti_2d.jl
index 892c8ed37f0..ab6b505e208 100644
--- a/test/test_dgmulti_2d.jl
+++ b/test/test_dgmulti_2d.jl
@@ -17,6 +17,7 @@ isdir(outdir) && rm(outdir, recursive = true)
 @trixi_testset "elixir_euler_weakform.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"),
                         cells_per_dimension=(4, 4),
+                        surface_integral=SurfaceIntegralWeakForm(FluxHLL(min_max_speed_naive)),
                         # division by 2.0 corresponds to normalization by the square root of the size of the domain
                         l2=[
                             0.0013536930300254945,
@@ -44,6 +45,7 @@ end
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"),
                         cells_per_dimension=(4, 4),
                         approximation_type=SBP(),
+                        surface_integral=SurfaceIntegralWeakForm(FluxHLL(min_max_speed_naive)),
                         # division by 2.0 corresponds to normalization by the square root of the size of the domain
                         l2=[
                             0.0074706882014934735,
@@ -71,6 +73,7 @@ end
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"),
                         cells_per_dimension=(4, 4),
                         element_type=Quad(),
+                        surface_integral=SurfaceIntegralWeakForm(FluxHLL(min_max_speed_naive)),
                         # division by 2.0 corresponds to normalization by the square root of the size of the domain
                         l2=[
                             0.00031892254415307093,
@@ -184,16 +187,12 @@ end
 @trixi_testset "elixir_euler_bilinear.jl (Bilinear quadrilateral elements, SBP, flux differencing)" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_bilinear.jl"),
                         l2=[
-                            1.0259435706215337e-5,
-                            9.014090233720625e-6,
-                            9.014090233223014e-6,
-                            2.738953587401793e-5,
+                            1.0259432774540821e-5, 9.014087689495575e-6,
+                            9.01408768888544e-6, 2.738953324859446e-5,
                         ],
                         linf=[
-                            7.362609083649829e-5,
-                            6.874188055272512e-5,
-                            6.874188052830021e-5,
-                            0.0001912435192696904,
+                            7.362605996297233e-5, 6.874189724781488e-5,
+                            6.874189703509614e-5, 0.00019124355334110277,
                         ])
     # Ensure that we do not have excessive memory allocations 
     # (e.g., from type instabilities) 
@@ -208,16 +207,12 @@ end
 @trixi_testset "elixir_euler_curved.jl (Quadrilateral elements, SBP, flux differencing)" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_curved.jl"),
                         l2=[
-                            1.720476068165337e-5,
-                            1.592168205710526e-5,
-                            1.592168205812963e-5,
-                            4.894094865697305e-5,
+                            1.7204593127904542e-5, 1.5921547179522804e-5,
+                            1.5921547180107928e-5, 4.894071422525737e-5,
                         ],
                         linf=[
-                            0.00010525416930584619,
-                            0.00010003778091061122,
-                            0.00010003778085621029,
-                            0.00036426282101720275,
+                            0.00010525416937667842, 0.00010003778102718464,
+                            0.00010003778071832059, 0.0003642628211952825,
                         ])
     # Ensure that we do not have excessive memory allocations 
     # (e.g., from type instabilities) 
@@ -232,6 +227,7 @@ end
 @trixi_testset "elixir_euler_curved.jl (Quadrilateral elements, GaussSBP, flux differencing)" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_curved.jl"),
                         approximation_type=GaussSBP(),
+                        surface_integral=SurfaceIntegralWeakForm(FluxHLL(min_max_speed_naive)),
                         l2=[
                             3.4666312079259457e-6,
                             3.4392774480368986e-6,
@@ -259,6 +255,7 @@ end
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_curved.jl"),
                         element_type=Tri(), approximation_type=Polynomial(),
                         volume_integral=VolumeIntegralWeakForm(),
+                        surface_integral=SurfaceIntegralWeakForm(FluxHLL(min_max_speed_naive)),
                         l2=[
                             7.905498158659466e-6,
                             8.731690809663625e-6,
@@ -330,16 +327,12 @@ end
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform_periodic.jl"),
                         # division by 2.0 corresponds to normalization by the square root of the size of the domain
                         l2=[
-                            0.0014986508075708323,
-                            0.001528523420746786,
-                            0.0015285234207473158,
-                            0.004846505183839211,
-                        ] ./ 2.0,
+                            0.0007492755162295128, 0.0007641875305302599,
+                            0.0007641875305306243, 0.0024232389721009447,
+                        ],
                         linf=[
-                            0.0015062108658376872,
-                            0.0019373508504645365,
-                            0.0019373508504538783,
-                            0.004742686826709086,
+                            0.0015060064614331736, 0.0019371156800773726,
+                            0.0019371156800769285, 0.004742431684202408,
                         ])
     # Ensure that we do not have excessive memory allocations 
     # (e.g., from type instabilities) 
@@ -354,16 +347,12 @@ end
 @trixi_testset "elixir_euler_triangulate_pkg_mesh.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_triangulate_pkg_mesh.jl"),
                         l2=[
-                            2.344080455438114e-6,
-                            1.8610038753097983e-6,
-                            2.4095165666095305e-6,
-                            6.373308158814308e-6,
+                            2.344076909832665e-6, 1.8610002398709756e-6,
+                            2.4095132179484066e-6, 6.37330249340445e-6,
                         ],
                         linf=[
-                            2.5099852761334418e-5,
-                            2.2683684021362893e-5,
-                            2.6180448559287584e-5,
-                            5.5752932611508044e-5,
+                            2.509979394305084e-5, 2.2683711321080935e-5,
+                            2.6180377720841363e-5, 5.575278031910713e-5,
                         ])
     # Ensure that we do not have excessive memory allocations 
     # (e.g., from type instabilities) 
@@ -435,16 +424,12 @@ end
                                  "elixir_euler_rayleigh_taylor_instability.jl"),
                         cells_per_dimension=(8, 8), tspan=(0.0, 0.2),
                         l2=[
-                            0.0709665896982514,
-                            0.005182828752164663,
-                            0.013832655585206478,
-                            0.03247013800580221,
+                            0.07097806723891838, 0.005168550941966817,
+                            0.013820912272220933, 0.03243357220022434,
                         ],
                         linf=[
-                            0.4783963902824797,
-                            0.022527207050681054,
-                            0.040307056293369226,
-                            0.0852365428206836,
+                            0.4783395896753895, 0.02244629340135818,
+                            0.04023357731088538, 0.08515807256615027,
                         ])
     # Ensure that we do not have excessive memory allocations 
     # (e.g., from type instabilities) 
@@ -604,16 +589,12 @@ end
 @trixi_testset "elixir_euler_fdsbp_periodic.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_fdsbp_periodic.jl"),
                         l2=[
-                            1.3333320340010056e-6,
-                            2.044834627970641e-6,
-                            2.044834627855601e-6,
-                            5.282189803559564e-6,
+                            1.333332033888785e-6, 2.044834627786368e-6,
+                            2.0448346278315884e-6, 5.282189803437435e-6,
                         ],
                         linf=[
-                            2.7000151718858945e-6,
-                            3.988595028259212e-6,
-                            3.9885950273710336e-6,
-                            8.848583042286862e-6,
+                            2.7000151703315822e-6, 3.988595025372632e-6,
+                            3.9885950240403645e-6, 8.848583036513702e-6,
                         ])
     # Ensure that we do not have excessive memory allocations 
     # (e.g., from type instabilities) 
@@ -628,6 +609,7 @@ end
 @trixi_testset "elixir_euler_fdsbp_periodic.jl (arbitrary reference domain)" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_fdsbp_periodic.jl"),
                         xmin=-200.0, xmax=100.0, #= parameters for reference interval =#
+                        surface_flux=FluxHLL(min_max_speed_naive),
                         l2=[
                             1.333332034149886e-6,
                             2.0448346280892024e-6,
@@ -659,6 +641,7 @@ end
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_fdsbp_periodic.jl"),
                         approximation_type=D,
                         coordinates_min=(-3.0, -4.0), coordinates_max=(0.0, -1.0),
+                        surface_flux=FluxHLL(min_max_speed_naive),
                         l2=[
                             0.07318831033918516,
                             0.10039910610067465,
@@ -691,6 +674,7 @@ end
     global D = SummationByPartsOperators.couple_continuously(D_local, mesh_local)
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_fdsbp_periodic.jl"),
                         approximation_type=D,
+                        surface_flux=FluxHLL(min_max_speed_naive),
                         l2=[
                             1.5440402410017893e-5,
                             1.4913189903083485e-5,
diff --git a/test/test_dgmulti_3d.jl b/test/test_dgmulti_3d.jl
index 3a1db255484..fa70b11447c 100644
--- a/test/test_dgmulti_3d.jl
+++ b/test/test_dgmulti_3d.jl
@@ -17,20 +17,15 @@ isdir(outdir) && rm(outdir, recursive = true)
 # 3d tet/hex tests
 @trixi_testset "elixir_euler_weakform.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"),
-                        # division by sqrt(8.0) corresponds to normalization by the square root of the size of the domain
                         l2=[
-                            0.0010029534292051608,
-                            0.0011682205957721673,
-                            0.001072975385793516,
-                            0.000997247778892257,
-                            0.0039364354651358294,
-                        ] ./ sqrt(8),
+                            0.000354593110864001, 0.00041301573702385284,
+                            0.00037934556184883277, 0.0003525767114354012,
+                            0.0013917457634530887,
+                        ],
                         linf=[
-                            0.003660737033303718,
-                            0.005625620600749226,
-                            0.0030566354814669516,
-                            0.0041580358824311325,
-                            0.019326660236036464,
+                            0.0036608123230692513, 0.005625540942772123,
+                            0.0030565781898950206, 0.004158099048202857,
+                            0.01932716837214299,
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -74,6 +69,7 @@ end
 @trixi_testset "elixir_euler_weakform.jl (Hexahedral elements)" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"),
                         element_type=Hex(),
+                        surface_integral=SurfaceIntegralWeakForm(FluxHLL(min_max_speed_naive)),
                         # division by sqrt(8.0) corresponds to normalization by the square root of the size of the domain
                         l2=[
                             0.00030580190715769566,
@@ -102,18 +98,13 @@ end
 @trixi_testset "elixir_euler_curved.jl (Hex elements, SBP, flux differencing)" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_curved.jl"),
                         l2=[
-                            0.018354883045936066,
-                            0.024412704052042846,
-                            0.024408520416087945,
-                            0.01816314570880129,
-                            0.039342805507972006,
+                            0.01835488304593566, 0.024412704052042534,
+                            0.02440852041608929, 0.018163145708800853,
+                            0.03934280550797125,
                         ],
                         linf=[
-                            0.14862225990775757,
-                            0.28952368161864683,
-                            0.2912054484817035,
-                            0.1456603133854122,
-                            0.3315354586775472,
+                            0.14862225990793032, 0.2895236816183626, 0.291205448481636,
+                            0.14566031338563246, 0.33153545867790246,
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -129,18 +120,14 @@ end
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_curved.jl"),
                         approximation_type=GaussSBP(),
                         l2=[
-                            0.002631131519508634,
-                            0.0029144224044954105,
-                            0.002913889110662827,
-                            0.002615140832314194,
-                            0.006881528610614373,
+                            0.0026311315195097097, 0.002914422404496567,
+                            0.0029138891106640368, 0.002615140832315232,
+                            0.006881528610616624,
                         ],
                         linf=[
-                            0.020996114874140215,
-                            0.021314522450134543,
-                            0.021288322783006297,
-                            0.020273381695435244,
-                            0.052598740390024545,
+                            0.02099611487415931, 0.021314522450152307,
+                            0.021288322783027613, 0.020273381695449455,
+                            0.05259874039006007,
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -154,20 +141,15 @@ end
 
 @trixi_testset "elixir_euler_weakform_periodic.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform_periodic.jl"),
-                        # division by sqrt(8.0) corresponds to normalization by the square root of the size of the domain
                         l2=[
-                            0.0010317074322517949,
-                            0.0012277090547035293,
-                            0.0011273991123913515,
-                            0.0010418496196130177,
-                            0.004058878478404962,
-                        ] ./ sqrt(8),
+                            0.00036475807571383924, 0.00043404536371780537,
+                            0.0003985850214093045, 0.0003683451584072326,
+                            0.00143503620472638,
+                        ],
                         linf=[
-                            0.003227752881827861,
-                            0.005620317864620361,
-                            0.0030514833972379307,
-                            0.003987027618439498,
-                            0.019282224709831652,
+                            0.0032278615418719347, 0.005620238272054934,
+                            0.0030514261010661237, 0.0039871165455998,
+                            0.019282771780667396,
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -182,6 +164,7 @@ end
 @trixi_testset "elixir_euler_weakform_periodic.jl (Hexahedral elements)" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform_periodic.jl"),
                         element_type=Hex(),
+                        surface_integral=SurfaceIntegralWeakForm(FluxHLL(min_max_speed_naive)),
                         # division by sqrt(8.0) corresponds to normalization by the square root of the size of the domain
                         l2=[
                             0.00034230612468547436,
diff --git a/test/test_parabolic_3d.jl b/test/test_parabolic_3d.jl
index 1eaa9f51a56..863daeeaf35 100644
--- a/test/test_parabolic_3d.jl
+++ b/test/test_parabolic_3d.jl
@@ -400,6 +400,7 @@ end
     @test_trixi_include(joinpath(examples_dir(), "p4est_3d_dgsem",
                                  "elixir_navierstokes_taylor_green_vortex.jl"),
                         initial_refinement_level=2, tspan=(0.0, 0.25),
+                        surface_flux=FluxHLL(min_max_speed_naive),
                         l2=[
                             0.0001547509861140407,
                             0.015637861347119624,
diff --git a/test/test_special_elixirs.jl b/test/test_special_elixirs.jl
index ba670a6025e..277ade9bd5c 100644
--- a/test/test_special_elixirs.jl
+++ b/test/test_special_elixirs.jl
@@ -286,7 +286,8 @@ end
             equations = CompressibleEulerEquations1D(1.4)
             mesh = TreeMesh((-1.0,), (1.0,), initial_refinement_level = 3,
                             n_cells_max = 10^4)
-            solver = DGSEM(3, flux_hll, VolumeIntegralFluxDifferencing(flux_ranocha))
+            solver = DGSEM(3, FluxHLL(min_max_speed_naive),
+                           VolumeIntegralFluxDifferencing(flux_ranocha))
             initial_condition = (x, t, equations) -> begin
                 rho = 2 + sinpi(k * sum(x))
                 v1 = 0.1
diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl
index f5fb169033a..64a1faf05b8 100644
--- a/test/test_structured_2d.jl
+++ b/test/test_structured_2d.jl
@@ -606,16 +606,12 @@ end
     @test_trixi_include(joinpath(EXAMPLES_DIR,
                                  "elixir_euler_rayleigh_taylor_instability.jl"),
                         l2=[
-                            0.06365630381017849,
-                            0.007166887387738937,
-                            0.002878708825497772,
-                            0.010247678114070121,
+                            0.06365630515019809, 0.007166887172039836,
+                            0.0028787103533600804, 0.010247678008197966,
                         ],
                         linf=[
-                            0.4799214336153155,
-                            0.024595483032220266,
-                            0.02059808120543466,
-                            0.03190756362943725,
+                            0.47992143569849377, 0.02459548251933757,
+                            0.02059810091623976, 0.0319077000843877,
                         ],
                         cells_per_dimension=(8, 8),
                         tspan=(0.0, 0.3))
@@ -659,14 +655,12 @@ end
 @trixi_testset "elixir_eulerpolytropic_convergence.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulerpolytropic_convergence.jl"),
                         l2=[
-                            0.0016688820596537988,
-                            0.0025921681885685425,
-                            0.003280950351435014,
+                            0.00166898321776379, 0.00259202637930991,
+                            0.0032810744946276406,
                         ],
                         linf=[
-                            0.010994679664394269,
-                            0.01331197845637,
-                            0.020080117011346488,
+                            0.010994883201888683, 0.013309526619369905,
+                            0.020080326611175536,
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -678,18 +672,19 @@ end
     end
 end
 
-@trixi_testset "elixir_eulerpolytropic_convergence.jl: HLL(Davis)" begin
-    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulerpolytropic_convergence.jl"),
+@trixi_testset "elixir_eulerpolytropic_convergence.jl with FluxHLL(min_max_speed_naive)" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                 "elixir_eulerpolytropic_convergence.jl"),
                         solver=DGSEM(polydeg = 3,
-                                     surface_flux = FluxHLL(min_max_speed_davis),
+                                     surface_flux = FluxHLL(min_max_speed_naive),
                                      volume_integral = VolumeIntegralFluxDifferencing(volume_flux)),
                         l2=[
-                            0.0016689832177644243, 0.0025920263793104445,
-                            0.003281074494629298,
+                            0.001668882059653298, 0.002592168188567654,
+                            0.0032809503514328307,
                         ],
                         linf=[
-                            0.01099488320190023, 0.013309526619350365,
-                            0.02008032661117909,
+                            0.01099467966437917, 0.013311978456333584,
+                            0.020080117011337606,
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -727,14 +722,12 @@ end
     @test_trixi_include(joinpath(EXAMPLES_DIR,
                                  "elixir_eulerpolytropic_isothermal_wave.jl"),
                         l2=[
-                            0.004998778491726366,
-                            0.004998916000294425,
-                            9.259136963058664e-17,
+                            0.004998778512795407, 0.004998916021367992,
+                            8.991558055435833e-17,
                         ],
                         linf=[
-                            0.010001103673834888,
-                            0.010051165098399503,
-                            7.623942913643681e-16,
+                            0.010001103632831354, 0.010051165055185603,
+                            7.60697457718599e-16,
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -749,14 +742,12 @@ end
 @trixi_testset "elixir_eulerpolytropic_wave.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulerpolytropic_wave.jl"),
                         l2=[
-                            0.23642682112204072,
-                            0.20904264390331334,
-                            8.174982691297391e-17,
+                            0.23642871172548174, 0.2090519382039672,
+                            8.778842676292274e-17,
                         ],
                         linf=[
-                            0.4848250368349989,
-                            0.253350873815695,
-                            4.984552457753618e-16,
+                            0.4852276879687425, 0.25327870807625175,
+                            5.533921691832115e-16,
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
diff --git a/test/test_threaded.jl b/test/test_threaded.jl
index a8a1b1b425a..7365dcef21c 100644
--- a/test/test_threaded.jl
+++ b/test/test_threaded.jl
@@ -394,10 +394,10 @@ end
                                      "elixir_euler_curved.jl"),
                             alg=RDPK3SpFSAL49(thread = OrdinaryDiffEq.True()),
                             l2=[
-                                1.720476068165337e-5,
-                                1.592168205710526e-5,
-                                1.592168205812963e-5,
-                                4.894094865697305e-5,
+                                1.7204593127904542e-5,
+                                1.5921547179522804e-5,
+                                1.5921547180107928e-5,
+                                4.894071422525737e-5,
                             ],
                             linf=[
                                 0.00010525416930584619,
@@ -420,16 +420,16 @@ end
         @test_trixi_include(joinpath(examples_dir(), "dgmulti_2d",
                                      "elixir_euler_triangulate_pkg_mesh.jl"),
                             l2=[
-                                2.344080455438114e-6,
-                                1.8610038753097983e-6,
-                                2.4095165666095305e-6,
-                                6.373308158814308e-6,
+                                2.344076909832665e-6,
+                                1.8610002398709756e-6,
+                                2.4095132179484066e-6,
+                                6.37330249340445e-6,
                             ],
                             linf=[
-                                2.5099852761334418e-5,
-                                2.2683684021362893e-5,
-                                2.6180448559287584e-5,
-                                5.5752932611508044e-5,
+                                2.509979394305084e-5,
+                                2.2683711321080935e-5,
+                                2.6180377720841363e-5,
+                                5.575278031910713e-5,
                             ])
 
         # Ensure that we do not have excessive memory allocations
diff --git a/test/test_tree_1d_euler.jl b/test/test_tree_1d_euler.jl
index 39a1f6e30ba..f26500b411c 100644
--- a/test/test_tree_1d_euler.jl
+++ b/test/test_tree_1d_euler.jl
@@ -221,11 +221,11 @@ end
 
 @trixi_testset "elixir_euler_ec.jl with flux_hll" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
-                        l2=[0.07852272782240548, 0.10209790867523805, 0.293873048809011],
+                        l2=[0.07855251823583848, 0.10213903748267686, 0.293985892532479],
                         linf=[
-                            0.19244768908604093,
-                            0.2515941686151897,
-                            0.7258000837553769,
+                            0.192621556068018,
+                            0.25184744005299536,
+                            0.7264977555504792,
                         ],
                         maxiters=10,
                         surface_flux=flux_hll,
diff --git a/test/test_tree_1d_eulergravity.jl b/test/test_tree_1d_eulergravity.jl
index 9ab5b287d0b..17bc0c71a7a 100644
--- a/test/test_tree_1d_eulergravity.jl
+++ b/test/test_tree_1d_eulergravity.jl
@@ -13,14 +13,12 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem")
 @trixi_testset "elixir_eulergravity_convergence.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulergravity_convergence.jl"),
                         l2=[
-                            0.0002170799126638106,
-                            0.0002913792848717502,
-                            0.0006112320856262327,
+                            0.00021708496949694728, 0.0002913795242132917,
+                            0.0006112500956552259,
                         ],
                         linf=[
-                            0.0004977401033188222,
-                            0.0013594223337776157,
-                            0.002041891084400227,
+                            0.0004977733237385706, 0.0013594226727522418,
+                            0.0020418739554664,
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
diff --git a/test/test_tree_1d_shallowwater.jl b/test/test_tree_1d_shallowwater.jl
index f9be63b87fd..41ad5c32bbd 100644
--- a/test/test_tree_1d_shallowwater.jl
+++ b/test/test_tree_1d_shallowwater.jl
@@ -143,17 +143,18 @@ end
 @trixi_testset "elixir_shallowwater_source_terms.jl with flux_hll" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"),
                         l2=[
-                            0.0022758146627220154,
-                            0.015864082886204556,
+                            0.002275023323848826,
+                            0.015861093821754046,
                             4.436491725585346e-5,
                         ],
                         linf=[
-                            0.008457195427364006,
-                            0.057201667446161064,
+                            0.008461451098266792,
+                            0.05722331401673486,
                             9.098379777405796e-5,
                         ],
                         tspan=(0.0, 0.025),
-                        surface_flux=(flux_hll, flux_nonconservative_fjordholm_etal))
+                        surface_flux=(flux_hll,
+                                      flux_nonconservative_fjordholm_etal))
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     let
@@ -228,7 +229,7 @@ end
                             0.05720939349382359,
                             9.098379777405796e-5,
                         ],
-                        surface_flux=(FluxHydrostaticReconstruction(flux_hll,
+                        surface_flux=(FluxHydrostaticReconstruction(FluxHLL(min_max_speed_naive),
                                                                     hydrostatic_reconstruction_audusse_etal),
                                       flux_nonconservative_audusse_etal),
                         tspan=(0.0, 0.025))
@@ -255,7 +256,9 @@ end
                             3.469453422316143e-15,
                             3.844551077492042e-8,
                         ],
-                        tspan=(0.0, 0.25))
+                        tspan=(0.0, 0.25),
+                        surface_flux=(FluxHLL(min_max_speed_naive),
+                                      flux_nonconservative_fjordholm_etal),)
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     let
@@ -280,6 +283,8 @@ end
                             3.844551077492042e-8,
                         ],
                         tspan=(0.0, 0.25),
+                        surface_flux=(FluxHLL(min_max_speed_naive),
+                                      flux_nonconservative_fjordholm_etal),
                         boundary_condition=boundary_condition_slip_wall)
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
diff --git a/test/test_tree_2d_mhd.jl b/test/test_tree_2d_mhd.jl
index 1f8458075aa..66b47138a44 100644
--- a/test/test_tree_2d_mhd.jl
+++ b/test/test_tree_2d_mhd.jl
@@ -183,7 +183,7 @@ end
     end
 end
 
-@trixi_testset "elixir_mhd_orszag_tang.jl with flux_hll" begin
+@trixi_testset "elixir_mhd_orszag_tang.jl with flux_hlle" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_orszag_tang.jl"),
                         l2=[
                             0.10806619664693064,
diff --git a/test/test_tree_2d_shallowwater.jl b/test/test_tree_2d_shallowwater.jl
index 93a8cb63667..01742644736 100644
--- a/test/test_tree_2d_shallowwater.jl
+++ b/test/test_tree_2d_shallowwater.jl
@@ -195,6 +195,33 @@ end
 end
 
 @trixi_testset "elixir_shallowwater_source_terms.jl with flux_hll" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"),
+                        l2=[
+                            0.0018952610547425214,
+                            0.016943425162728183,
+                            0.017556784292859465,
+                            6.274146767717414e-5,
+                        ],
+                        linf=[
+                            0.0151635341334182,
+                            0.07967467926956129,
+                            0.08400050790965174,
+                            0.0001819675955490041,
+                        ],
+                        tspan=(0.0, 0.025),
+                        surface_flux=(flux_hll,
+                                      flux_nonconservative_fjordholm_etal))
+    # 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_shallowwater_source_terms.jl with FluxHLL(min_max_speed_naive)" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"),
                         l2=[
                             0.0018957692481057034,
@@ -209,7 +236,8 @@ end
                             0.0001819675955490041,
                         ],
                         tspan=(0.0, 0.025),
-                        surface_flux=(flux_hll, flux_nonconservative_fjordholm_etal))
+                        surface_flux=(FluxHLL(min_max_speed_naive),
+                                      flux_nonconservative_fjordholm_etal))
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     let
diff --git a/test/test_tree_3d_euler.jl b/test/test_tree_3d_euler.jl
index 02e657e001a..e9e2b82fec5 100644
--- a/test/test_tree_3d_euler.jl
+++ b/test/test_tree_3d_euler.jl
@@ -92,18 +92,14 @@ end
 @trixi_testset "elixir_euler_convergence.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_convergence.jl"),
                         l2=[
-                            0.0003637241020254405,
-                            0.0003955570866382718,
-                            0.0003955570866383613,
-                            0.00039555708663834417,
-                            0.0007811613481640202,
+                            0.0003637241020254673, 0.00039555708663848046,
+                            0.00039555708663832644, 0.0003955570866385083,
+                            0.0007811613481643962,
                         ],
                         linf=[
-                            0.0024000660244674066,
-                            0.0029635410025339315,
-                            0.0029635410025292686,
-                            0.002963541002525938,
-                            0.007191437359396424,
+                            0.0024000660244567484, 0.002963541002521053,
+                            0.0029635410025201647, 0.002963541002522385,
+                            0.007191437359379549,
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -462,6 +458,7 @@ end
                             2.397746252817731,
                         ],
                         maxiters=5, max_level=6,
+                        surface_flux=FluxHLL(min_max_speed_naive),
                         coverage_override=(maxiters = 2, initial_refinement_level = 1,
                                            base_level = 1, max_level = 3))
     # Ensure that we do not have excessive memory allocations
diff --git a/test/test_tree_3d_mhd.jl b/test/test_tree_3d_mhd.jl
index e75685f0b43..74107d462de 100644
--- a/test/test_tree_3d_mhd.jl
+++ b/test/test_tree_3d_mhd.jl
@@ -184,9 +184,9 @@ end
     end
 end
 
-@trixi_testset "elixir_mhd_alfven_wave.jl with Orszag-Tang setup + flux_hll" begin
+@trixi_testset "elixir_mhd_alfven_wave.jl with Orszag-Tang setup + flux_hlle" begin
     # OBS! This setup does not make much sense and is only used to exercise all components of the
-    # flux_hll implementation
+    # flux_hlle implementation
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
                         l2=[
                             0.004391143689111404,
diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl
index 04eb9f679aa..87d677e1623 100644
--- a/test/test_unstructured_2d.jl
+++ b/test/test_unstructured_2d.jl
@@ -17,16 +17,12 @@ isdir(outdir) && rm(outdir, recursive = true)
 @trixi_testset "elixir_euler_periodic.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_periodic.jl"),
                         l2=[
-                            0.00010978828464875207,
-                            0.00013010359527356914,
-                            0.00013010359527326057,
-                            0.0002987656724828824,
+                            0.0001099216141882387, 0.0001303795774982892,
+                            0.00013037957749794242, 0.0002993727892598759,
                         ],
                         linf=[
-                            0.00638626102818618,
-                            0.009804042508242183,
-                            0.009804042508253286,
-                            0.02183139311614468,
+                            0.006407280810928562, 0.009836067015418948,
+                            0.009836067015398076, 0.021903519038095176,
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -41,16 +37,12 @@ end
 @trixi_testset "elixir_euler_free_stream.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream.jl"),
                         l2=[
-                            3.3937971107485363e-14,
-                            2.447586447887882e-13,
-                            1.4585205789296455e-13,
-                            4.716993468962946e-13,
+                            3.3937365073416665e-14, 2.44759188939065e-13,
+                            1.4585198700082895e-13, 4.716940764877479e-13,
                         ],
                         linf=[
-                            8.804734719092266e-12,
-                            6.261270668606045e-11,
-                            2.93670088247211e-11,
-                            1.205400224080222e-10,
+                            8.804956763697191e-12, 6.261199891888225e-11,
+                            2.936639820205755e-11, 1.20543575121701e-10,
                         ],
                         tspan=(0.0, 0.1),
                         atol=3.0e-13)
@@ -78,7 +70,8 @@ end
                             0.29339040847600434,
                             0.5915610037764794,
                         ],
-                        tspan=(0.0, 0.25))
+                        tspan=(0.0, 0.25),
+                        surface_flux=FluxHLL(min_max_speed_naive))
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     let
@@ -407,15 +400,15 @@ end
 @trixi_testset "elixir_shallowwater_source_terms.jl with FluxHydrostaticReconstruction" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"),
                         l2=[
-                            0.0011197139793938152,
-                            0.015430259691310781,
-                            0.017081031802719724,
+                            0.001119678684752799,
+                            0.015429108794630785,
+                            0.01708275441241111,
                             5.089218476758271e-6,
                         ],
                         linf=[
-                            0.014300809338967824,
-                            0.12783372461225184,
-                            0.17625472321992852,
+                            0.014299564388827513,
+                            0.12785126473870534,
+                            0.17626788561725526,
                             2.6407324614341476e-5,
                         ],
                         surface_flux=(FluxHydrostaticReconstruction(flux_hll,
@@ -464,18 +457,19 @@ end
 @trixi_testset "elixir_shallowwater_source_terms.jl with flux_hll" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"),
                         l2=[
-                            0.0011197139793938727,
-                            0.015430259691311309,
-                            0.017081031802719554,
+                            0.0011196786847528799,
+                            0.015429108794631075,
+                            0.017082754412411742,
                             5.089218476759981e-6,
                         ],
                         linf=[
-                            0.014300809338967824,
-                            0.12783372461224918,
-                            0.17625472321993918,
+                            0.014299564388830177,
+                            0.12785126473870667,
+                            0.17626788561728546,
                             2.6407324614341476e-5,
                         ],
-                        surface_flux=(flux_hll, flux_nonconservative_fjordholm_etal),
+                        surface_flux=(flux_hll,
+                                      flux_nonconservative_fjordholm_etal),
                         tspan=(0.0, 0.025))
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -490,16 +484,12 @@ end
 @trixi_testset "elixir_shallowwater_dirichlet.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_dirichlet.jl"),
                         l2=[
-                            1.1577518608940115e-5,
-                            4.867189932537344e-13,
-                            4.647273240470541e-13,
-                            1.1577518608933468e-5,
+                            1.1577518608938916e-5, 4.859252379740366e-13,
+                            4.639600837197925e-13, 1.1577518608952174e-5,
                         ],
                         linf=[
-                            8.394063878602864e-5,
-                            1.1469760027632646e-10,
-                            1.1146619484429974e-10,
-                            8.394063879602065e-5,
+                            8.3940638787805e-5, 1.1446362498574484e-10,
+                            1.1124515748367981e-10, 8.39406387962427e-5,
                         ],
                         tspan=(0.0, 2.0))
     # Ensure that we do not have excessive memory allocations
@@ -516,16 +506,12 @@ end
     @test_trixi_include(joinpath(EXAMPLES_DIR,
                                  "elixir_shallowwater_wall_bc_shockcapturing.jl"),
                         l2=[
-                            0.04444388691670699,
-                            0.1527771788033111,
-                            0.1593763537203512,
-                            6.225080476986749e-8,
+                            0.0442113635677511, 0.1537465759364839, 0.16003586586203947,
+                            6.225080477067782e-8,
                         ],
                         linf=[
-                            0.6526506870169639,
-                            1.980765893182952,
-                            2.4807635459119757,
-                            3.982097158683473e-7,
+                            0.6347820607387928, 2.0078125433846736, 2.530726684667019,
+                            3.982097165344811e-7,
                         ],
                         tspan=(0.0, 0.05))
     # Ensure that we do not have excessive memory allocations
@@ -609,12 +595,12 @@ end
                             1.0066867437607972e-13,
                             6.889210012578449e-14,
                             1.568290814572709e-13],
-                        linf=[2.353373051988683e-10,
-                            2.801543719233024e-11,
-                            3.930469838486772e-11,
+                        linf=[5.6139981552405516e-11,
+                            2.842849566864203e-11,
+                            1.8290174930157832e-11,
                             4.61017890529547e-11],
                         tspan=(0.0, 0.1),
-                        atol=1.0e-11)
+                        atol=1.0e-10)
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     let