diff --git a/Project.toml b/Project.toml
index da618dc78c1..17c02a2adac 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Trixi"
 uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
 authors = ["Michael Schlottke-Lakemper <m.schlottke-lakemper@hlrs.de>", "Gregor Gassner <ggassner@uni-koeln.de>", "Hendrik Ranocha <mail@ranocha.de>", "Andrew R. Winters <andrew.ross.winters@liu.se>", "Jesse Chan <jesse.chan@rice.edu>"]
-version = "0.6.1-pre"
+version = "0.6.2-pre"
 
 [deps]
 CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
diff --git a/examples/p4est_2d_dgsem/elixir_advection_restart.jl b/examples/p4est_2d_dgsem/elixir_advection_restart.jl
index 4f43e122ab3..0f573714c1f 100644
--- a/examples/p4est_2d_dgsem/elixir_advection_restart.jl
+++ b/examples/p4est_2d_dgsem/elixir_advection_restart.jl
@@ -31,7 +31,7 @@ save_solution.condition.save_initial_solution = false
 
 integrator = init(ode, CarpenterKennedy2N54(williamson_condition = false),
                   dt = dt, # solve needs some value here but it will be overwritten by the stepsize_callback
-                  save_everystep = false, callback = callbacks);
+                  save_everystep = false, callback = callbacks, maxiters = 100_000);
 
 # Get the last time index and work with that.
 load_timestep!(integrator, restart_filename)
diff --git a/examples/p4est_3d_dgsem/elixir_advection_restart.jl b/examples/p4est_3d_dgsem/elixir_advection_restart.jl
index cd97d69d692..b3dead42399 100644
--- a/examples/p4est_3d_dgsem/elixir_advection_restart.jl
+++ b/examples/p4est_3d_dgsem/elixir_advection_restart.jl
@@ -29,7 +29,7 @@ save_solution.condition.save_initial_solution = false
 
 integrator = init(ode, CarpenterKennedy2N54(williamson_condition = false),
                   dt = dt, # solve needs some value here but it will be overwritten by the stepsize_callback
-                  save_everystep = false, callback = callbacks);
+                  save_everystep = false, callback = callbacks, maxiters = 100_000);
 
 # Get the last time index and work with that.
 load_timestep!(integrator, restart_filename)
diff --git a/examples/structured_2d_dgsem/elixir_advection_restart.jl b/examples/structured_2d_dgsem/elixir_advection_restart.jl
index 19863faae8d..0accbdba702 100644
--- a/examples/structured_2d_dgsem/elixir_advection_restart.jl
+++ b/examples/structured_2d_dgsem/elixir_advection_restart.jl
@@ -30,7 +30,7 @@ save_solution.condition.save_initial_solution = false
 
 integrator = init(ode, CarpenterKennedy2N54(williamson_condition = false),
                   dt = dt, # solve needs some value here but it will be overwritten by the stepsize_callback
-                  save_everystep = false, callback = callbacks);
+                  save_everystep = false, callback = callbacks, maxiters = 100_000);
 
 # Get the last time index and work with that.
 load_timestep!(integrator, restart_filename)
diff --git a/examples/structured_3d_dgsem/elixir_advection_restart.jl b/examples/structured_3d_dgsem/elixir_advection_restart.jl
index e81ad5d6430..e516d794df8 100644
--- a/examples/structured_3d_dgsem/elixir_advection_restart.jl
+++ b/examples/structured_3d_dgsem/elixir_advection_restart.jl
@@ -29,7 +29,7 @@ save_solution.condition.save_initial_solution = false
 
 integrator = init(ode, CarpenterKennedy2N54(williamson_condition = false),
                   dt = dt, # solve needs some value here but it will be overwritten by the stepsize_callback
-                  save_everystep = false, callback = callbacks);
+                  save_everystep = false, callback = callbacks, maxiters = 100_000);
 
 # Get the last time index and work with that.
 load_timestep!(integrator, restart_filename)
diff --git a/examples/tree_2d_dgsem/elixir_advection_restart.jl b/examples/tree_2d_dgsem/elixir_advection_restart.jl
index 770629bb15e..e0d1003f524 100644
--- a/examples/tree_2d_dgsem/elixir_advection_restart.jl
+++ b/examples/tree_2d_dgsem/elixir_advection_restart.jl
@@ -29,7 +29,7 @@ save_solution.condition.save_initial_solution = false
 
 integrator = init(ode, alg,
                   dt = dt, # solve needs some value here but it will be overwritten by the stepsize_callback
-                  save_everystep = false, callback = callbacks; ode_default_options()...)
+                  callback = callbacks, maxiters = 100_000; ode_default_options()...)
 
 # Load saved context for adaptive time integrator
 if integrator.opts.adaptive
diff --git a/examples/tree_2d_dgsem/elixir_eulerpolytropic_convergence.jl b/examples/tree_2d_dgsem/elixir_eulerpolytropic_convergence.jl
new file mode 100644
index 00000000000..95ef38eb701
--- /dev/null
+++ b/examples/tree_2d_dgsem/elixir_eulerpolytropic_convergence.jl
@@ -0,0 +1,63 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the polytropic Euler equations
+
+gamma = 1.4
+kappa = 0.5     # Scaling factor for the pressure.
+equations = PolytropicEulerEquations2D(gamma, kappa)
+
+initial_condition = initial_condition_convergence_test
+
+volume_flux = flux_winters_etal
+solver = DGSEM(polydeg = 3, surface_flux = FluxHLL(min_max_speed_davis),
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+coordinates_min = (0.0, 0.0)
+coordinates_max = (1.0, 1.0)
+
+# Create a uniformly refined mesh with periodic boundaries
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 2,
+                periodicity = true,
+                n_cells_max = 30_000)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms = source_terms_convergence_test)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.1)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     extra_analysis_errors = (:l2_error_primitive,
+                                                              :linf_error_primitive))
+
+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.1)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_solution,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, # 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
diff --git a/examples/tree_3d_dgsem/elixir_advection_restart.jl b/examples/tree_3d_dgsem/elixir_advection_restart.jl
index f81e013fc0a..a53f4ebf5b1 100644
--- a/examples/tree_3d_dgsem/elixir_advection_restart.jl
+++ b/examples/tree_3d_dgsem/elixir_advection_restart.jl
@@ -27,7 +27,7 @@ save_solution.condition.save_initial_solution = false
 
 integrator = init(ode, CarpenterKennedy2N54(williamson_condition = false),
                   dt = dt, # solve needs some value here but it will be overwritten by the stepsize_callback
-                  save_everystep = false, callback = callbacks);
+                  save_everystep = false, callback = callbacks, maxiters = 100_000);
 
 # Get the last time index and work with that.
 load_timestep!(integrator, restart_filename)
diff --git a/examples/unstructured_2d_dgsem/elixir_euler_restart.jl b/examples/unstructured_2d_dgsem/elixir_euler_restart.jl
index 4d6af65b8a4..c1908691902 100644
--- a/examples/unstructured_2d_dgsem/elixir_euler_restart.jl
+++ b/examples/unstructured_2d_dgsem/elixir_euler_restart.jl
@@ -29,7 +29,7 @@ save_solution.condition.save_initial_solution = false
 
 integrator = init(ode, CarpenterKennedy2N54(williamson_condition = false),
                   dt = dt, # solve needs some value here but it will be overwritten by the stepsize_callback
-                  save_everystep = false, callback = callbacks);
+                  save_everystep = false, callback = callbacks, maxiters = 100_000);
 
 # Get the last time index and work with that.
 load_timestep!(integrator, restart_filename)
diff --git a/src/auxiliary/special_elixirs.jl b/src/auxiliary/special_elixirs.jl
index 919eddb18fa..5fdd9aea0c5 100644
--- a/src/auxiliary/special_elixirs.jl
+++ b/src/auxiliary/special_elixirs.jl
@@ -36,6 +36,16 @@ julia> redirect_stdout(devnull) do
 ```
 """
 function trixi_include(mod::Module, elixir::AbstractString; kwargs...)
+    # Check that all kwargs exist as assignments
+    code = read(elixir, String)
+    expr = Meta.parse("begin \n$code \nend")
+    expr = insert_maxiters(expr)
+
+    for (key, val) in kwargs
+        # This will throw an error when `key` is not found
+        find_assignment(expr, key)
+    end
+
     # Print information on potential wait time only in non-parallel case
     if !mpi_isparallel()
         @info "You just called `trixi_include`. Julia may now compile the code, please be patient."
@@ -243,6 +253,7 @@ end
 function find_assignment(expr, destination)
     # declare result to be able to assign to it in the closure
     local result
+    found = false
 
     # find explicit and keyword assignments
     walkexpr(expr) do x
@@ -250,12 +261,17 @@ function find_assignment(expr, destination)
             if (x.head === Symbol("=") || x.head === :kw) &&
                x.args[1] === Symbol(destination)
                 result = x.args[2]
+                found = true
                 # dump(x)
             end
         end
         return x
     end
 
+    if !found
+        throw(ArgumentError("assignment `$destination` not found in expression"))
+    end
+
     result
 end
 
@@ -274,17 +290,28 @@ function extract_initial_resolution(elixir, kwargs)
             return initial_refinement_level
         end
     catch e
-        if isa(e, UndefVarError)
-            # get cells_per_dimension from the elixir
-            cells_per_dimension = eval(find_assignment(expr, :cells_per_dimension))
-
-            if haskey(kwargs, :cells_per_dimension)
-                return kwargs[:cells_per_dimension]
-            else
-                return cells_per_dimension
+        # If `initial_refinement_level` is not found, we will get an `ArgumentError`
+        if isa(e, ArgumentError)
+            try
+                # get cells_per_dimension from the elixir
+                cells_per_dimension = eval(find_assignment(expr, :cells_per_dimension))
+
+                if haskey(kwargs, :cells_per_dimension)
+                    return kwargs[:cells_per_dimension]
+                else
+                    return cells_per_dimension
+                end
+            catch e2
+                # If `cells_per_dimension` is not found either
+                if isa(e2, ArgumentError)
+                    throw(ArgumentError("`convergence_test` requires the elixir to define " *
+                                        "`initial_refinement_level` or `cells_per_dimension`"))
+                else
+                    rethrow()
+                end
             end
         else
-            throw(e)
+            rethrow()
         end
     end
 end
diff --git a/src/equations/polytropic_euler_2d.jl b/src/equations/polytropic_euler_2d.jl
index d4902bbafb2..f5d2f7b0bad 100644
--- a/src/equations/polytropic_euler_2d.jl
+++ b/src/equations/polytropic_euler_2d.jl
@@ -120,7 +120,7 @@ function initial_condition_weak_blast_wave(x, t, equations::PolytropicEulerEquat
     return prim2cons(SVector(rho, v1, v2), equations)
 end
 
-# Calculate 1D flux for a single point in the normal direction
+# Calculate 2D flux for a single point in the normal direction
 # Note, this directional vector is not normalized
 @inline function flux(u, normal_direction::AbstractVector,
                       equations::PolytropicEulerEquations2D)
@@ -135,8 +135,28 @@ end
     return SVector(f1, f2, f3)
 end
 
+# Calculate 2D flux for a single point
+@inline function flux(u, orientation::Integer, equations::PolytropicEulerEquations2D)
+    _, v1, v2 = cons2prim(u, equations)
+    p = pressure(u, equations)
+
+    rho_v1 = u[2]
+    rho_v2 = u[3]
+
+    if orientation == 1
+        f1 = rho_v1
+        f2 = rho_v1 * v1 + p
+        f3 = rho_v1 * v2
+    else
+        f1 = rho_v2
+        f2 = rho_v2 * v1
+        f3 = rho_v2 * v2 + p
+    end
+    return SVector(f1, f2, f3)
+end
+
 """
-    flux_winters_etal(u_ll, u_rr, normal_direction,
+    flux_winters_etal(u_ll, u_rr, orientation_or_normal_direction,
                       equations::PolytropicEulerEquations2D)
 
 Entropy conserving two-point flux for isothermal or polytropic gases.
@@ -178,6 +198,37 @@ For details see Section 3.2 of the following reference
     return SVector(f1, f2, f3)
 end
 
+@inline function flux_winters_etal(u_ll, u_rr, orientation::Integer,
+                                   equations::PolytropicEulerEquations2D)
+    # Unpack left and right state
+    rho_ll, v1_ll, v2_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_rr = cons2prim(u_rr, equations)
+    p_ll = equations.kappa * rho_ll^equations.gamma
+    p_rr = equations.kappa * rho_rr^equations.gamma
+
+    # Compute the necessary mean values
+    if equations.gamma == 1.0 # isothermal gas
+        rho_mean = ln_mean(rho_ll, rho_rr)
+    else # equations.gamma > 1 # polytropic gas
+        rho_mean = stolarsky_mean(rho_ll, rho_rr, equations.gamma)
+    end
+    v1_avg = 0.5 * (v1_ll + v1_rr)
+    v2_avg = 0.5 * (v2_ll + v2_rr)
+    p_avg = 0.5 * (p_ll + p_rr)
+
+    if orientation == 1 # x-direction
+        f1 = rho_mean * 0.5 * (v1_ll + v1_rr)
+        f2 = f1 * v1_avg + p_avg
+        f3 = f1 * v2_avg
+    else # y-direction
+        f1 = rho_mean * 0.5 * (v2_ll + v2_rr)
+        f2 = f1 * v1_avg
+        f3 = f1 * v2_avg + p_avg
+    end
+
+    return SVector(f1, f2, f3)
+end
+
 @inline function min_max_speed_naive(u_ll, u_rr, normal_direction::AbstractVector,
                                      equations::PolytropicEulerEquations2D)
     rho_ll, v1_ll, v2_ll = cons2prim(u_ll, equations)
@@ -196,6 +247,53 @@ end
     return lambda_min, lambda_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::PolytropicEulerEquations2D)
+    rho_ll, v1_ll, v2_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_rr = cons2prim(u_rr, equations)
+    # Pressure for polytropic Euler
+    p_ll = equations.kappa * rho_ll^equations.gamma
+    p_rr = equations.kappa * rho_rr^equations.gamma
+
+    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::PolytropicEulerEquations2D)
+    rho_ll, v1_ll, v2_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_rr = cons2prim(u_rr, equations)
+    # Pressure for polytropic Euler
+    p_ll = equations.kappa * rho_ll^equations.gamma
+    p_rr = equations.kappa * rho_rr^equations.gamma
+
+    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
+
 @inline function max_abs_speeds(u, equations::PolytropicEulerEquations2D)
     rho, v1, v2 = cons2prim(u, equations)
     c = sqrt(equations.gamma * equations.kappa * rho^(equations.gamma - 1))
diff --git a/test/test_mpi_p4est_2d.jl b/test/test_mpi_p4est_2d.jl
index 1edbce8f6c8..da90537fcfd 100644
--- a/test/test_mpi_p4est_2d.jl
+++ b/test/test_mpi_p4est_2d.jl
@@ -69,7 +69,10 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "p4est_2d_dgsem")
     @trixi_testset "elixir_advection_restart.jl" begin
         @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart.jl"),
                             l2=[4.507575525876275e-6],
-                            linf=[6.21489667023134e-5])
+                            linf=[6.21489667023134e-5],
+                            # With the default `maxiters = 1` in coverage tests,
+                            # there would be no time steps after the restart.
+                            coverage_override=(maxiters = 100_000,))
     end
 
     @trixi_testset "elixir_euler_source_terms_nonconforming_unstructured_flag.jl" begin
diff --git a/test/test_mpi_p4est_3d.jl b/test/test_mpi_p4est_3d.jl
index 8082930b3b4..75f43650082 100644
--- a/test/test_mpi_p4est_3d.jl
+++ b/test/test_mpi_p4est_3d.jl
@@ -63,7 +63,10 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "p4est_3d_dgsem")
     @trixi_testset "elixir_advection_restart.jl" begin
         @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart.jl"),
                             l2=[0.002590388934758452],
-                            linf=[0.01840757696885409])
+                            linf=[0.01840757696885409],
+                            # With the default `maxiters = 1` in coverage tests,
+                            # there would be no time steps after the restart.
+                            coverage_override=(maxiters = 100_000,))
     end
 
     @trixi_testset "elixir_advection_cubed_sphere.jl" begin
diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl
index 07c6d02bbcd..db34aecc168 100644
--- a/test/test_p4est_2d.jl
+++ b/test/test_p4est_2d.jl
@@ -94,7 +94,10 @@ end
 @trixi_testset "elixir_advection_restart.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart.jl"),
                         l2=[4.507575525876275e-6],
-                        linf=[6.21489667023134e-5])
+                        linf=[6.21489667023134e-5],
+                        # With the default `maxiters = 1` in coverage tests,
+                        # there would be no time steps after the restart.
+                        coverage_override=(maxiters = 100_000,))
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     let
@@ -216,8 +219,8 @@ end
                         ],
                         surface_flux=flux_hlle,
                         tspan=(0.0, 0.3))
-    # Ensure that we do not have excessive memory allocations 
-    # (e.g., from type instabilities) 
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
     let
         t = sol.t[end]
         u_ode = sol.u[end]
diff --git a/test/test_p4est_3d.jl b/test/test_p4est_3d.jl
index 346c61c7448..f2467f30204 100644
--- a/test/test_p4est_3d.jl
+++ b/test/test_p4est_3d.jl
@@ -110,7 +110,10 @@ end
 @trixi_testset "elixir_advection_restart.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart.jl"),
                         l2=[0.002590388934758452],
-                        linf=[0.01840757696885409])
+                        linf=[0.01840757696885409],
+                        # With the default `maxiters = 1` in coverage tests,
+                        # there would be no time steps after the restart.
+                        coverage_override=(maxiters = 100_000,))
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     let
@@ -305,8 +308,8 @@ end
                         ],
                         tspan=(0.0, 0.3),
                         surface_flux=flux_hlle)
-    # Ensure that we do not have excessive memory allocations 
-    # (e.g., from type instabilities) 
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
     let
         t = sol.t[end]
         u_ode = sol.u[end]
diff --git a/test/test_special_elixirs.jl b/test/test_special_elixirs.jl
index 85671002ba6..ba670a6025e 100644
--- a/test/test_special_elixirs.jl
+++ b/test/test_special_elixirs.jl
@@ -92,8 +92,7 @@ coverage = occursin("--code-coverage", cmd) && !occursin("--code-coverage=none",
         # the convergence test logic
         @test_nowarn_mod convergence_test(@__MODULE__,
                                           joinpath(EXAMPLES_DIR, "tree_2d_dgsem",
-                                                   "elixir_advection_basic.jl"), 2,
-                                          tspan = (0.0, 0.01))
+                                                   "elixir_advection_basic.jl"), 2)
         @test_nowarn_mod convergence_test(@__MODULE__,
                                           joinpath(EXAMPLES_DIR, "tree_2d_dgsem",
                                                    "elixir_advection_extended.jl"), 2,
@@ -101,12 +100,10 @@ coverage = occursin("--code-coverage", cmd) && !occursin("--code-coverage=none",
                                           tspan = (0.0, 0.1))
         @test_nowarn_mod convergence_test(@__MODULE__,
                                           joinpath(EXAMPLES_DIR, "structured_2d_dgsem",
-                                                   "elixir_advection_basic.jl"), 2,
-                                          tspan = (0.0, 0.01))
+                                                   "elixir_advection_basic.jl"), 2)
         @test_nowarn_mod convergence_test(@__MODULE__,
                                           joinpath(EXAMPLES_DIR, "structured_2d_dgsem",
-                                                   "elixir_advection_coupled.jl"), 2,
-                                          tspan = (0.0, 0.01))
+                                                   "elixir_advection_coupled.jl"), 2)
         @test_nowarn_mod convergence_test(@__MODULE__,
                                           joinpath(EXAMPLES_DIR, "structured_2d_dgsem",
                                                    "elixir_advection_extended.jl"), 2,
@@ -345,8 +342,8 @@ end
     end
 
     @timed_testset "elixir_euler_ad.jl" begin
-        @test_trixi_include(joinpath(examples_dir(), "special_elixirs",
-                                     "elixir_euler_ad.jl"))
+        @test_nowarn_mod trixi_include(joinpath(examples_dir(), "special_elixirs",
+                                                "elixir_euler_ad.jl"))
     end
 end
 end
diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl
index 622d6daf2e2..e317783f95b 100644
--- a/test/test_structured_2d.jl
+++ b/test/test_structured_2d.jl
@@ -194,7 +194,10 @@ end
 @trixi_testset "elixir_advection_restart.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart.jl"),
                         l2=[4.219208035582454e-6],
-                        linf=[3.438434404412494e-5])
+                        linf=[3.438434404412494e-5],
+                        # With the default `maxiters = 1` in coverage tests,
+                        # there would be no time steps after the restart.
+                        coverage_override=(maxiters = 100_000,))
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     let
@@ -211,7 +214,10 @@ end
                         linf=[0.0015194252169410394],
                         rtol=5.0e-5, # Higher tolerance to make tests pass in CI (in particular with macOS)
                         elixir_file="elixir_advection_waving_flag.jl",
-                        restart_file="restart_000021.h5")
+                        restart_file="restart_000021.h5",
+                        # With the default `maxiters = 1` in coverage tests,
+                        # there would be no time steps after the restart.
+                        coverage_override=(maxiters = 100_000,))
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     let
@@ -227,7 +233,10 @@ end
                         l2=[7.841217436552029e-15],
                         linf=[1.0857981180834031e-13],
                         elixir_file="elixir_advection_free_stream.jl",
-                        restart_file="restart_000036.h5")
+                        restart_file="restart_000036.h5",
+                        # With the default `maxiters = 1` in coverage tests,
+                        # there would be no time steps after the restart.
+                        coverage_override=(maxiters = 100_000,))
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     let
@@ -825,6 +834,29 @@ end
     end
 end
 
+@trixi_testset "elixir_eulerpolytropic_convergence.jl: HLL(Davis)" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulerpolytropic_convergence.jl"),
+                        solver=DGSEM(polydeg = 3,
+                                     surface_flux = FluxHLL(min_max_speed_davis),
+                                     volume_integral = VolumeIntegralFluxDifferencing(volume_flux)),
+                        l2=[
+                            0.0016689832177644243, 0.0025920263793104445,
+                            0.003281074494629298,
+                        ],
+                        linf=[
+                            0.01099488320190023, 0.013309526619350365,
+                            0.02008032661117909,
+                        ])
+    # 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_eulerpolytropic_ec.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulerpolytropic_ec.jl"),
                         l2=[
diff --git a/test/test_structured_3d.jl b/test/test_structured_3d.jl
index 0213e1a9813..a52c459d6be 100644
--- a/test/test_structured_3d.jl
+++ b/test/test_structured_3d.jl
@@ -62,7 +62,10 @@ end
 @trixi_testset "elixir_advection_restart.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart.jl"),
                         l2=[0.0025903889347585777],
-                        linf=[0.018407576968841655])
+                        linf=[0.018407576968841655],
+                        # With the default `maxiters = 1` in coverage tests,
+                        # there would be no time steps after the restart.
+                        coverage_override=(maxiters = 100_000,))
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     let
diff --git a/test/test_threaded.jl b/test/test_threaded.jl
index 478c90b476a..dbbcbf4c7ce 100644
--- a/test/test_threaded.jl
+++ b/test/test_threaded.jl
@@ -208,7 +208,10 @@ end
                             linf=[0.0015194252169410394],
                             rtol=5.0e-5, # Higher tolerance to make tests pass in CI (in particular with macOS)
                             elixir_file="elixir_advection_waving_flag.jl",
-                            restart_file="restart_000021.h5")
+                            restart_file="restart_000021.h5",
+                            # With the default `maxiters = 1` in coverage tests,
+                            # there would be no time steps after the restart.
+                            coverage_override=(maxiters = 100_000,))
 
         # Ensure that we do not have excessive memory allocations
         # (e.g., from type instabilities)
diff --git a/test/test_tree_2d_eulerpolytropic.jl b/test/test_tree_2d_eulerpolytropic.jl
new file mode 100644
index 00000000000..545cf7274ff
--- /dev/null
+++ b/test/test_tree_2d_eulerpolytropic.jl
@@ -0,0 +1,35 @@
+module TestExamples2DEulerMulticomponent
+
+using Test
+using Trixi
+
+include("test_trixi.jl")
+
+EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem")
+
+@testset "Polytropic Euler" begin
+#! format: noindent
+
+@trixi_testset "elixir_eulerpolytropic_convergence.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                 "elixir_eulerpolytropic_convergence.jl"),
+                        l2=[
+                            0.0016689832177626373, 0.0025920263793094526,
+                            0.003281074494626679,
+                        ],
+                        linf=[
+                            0.010994883201896677, 0.013309526619350365,
+                            0.02008032661117376,
+                        ])
+    # 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
+end
+
+end # module
diff --git a/test/test_tree_2d_part2.jl b/test/test_tree_2d_part2.jl
index d8e86d14f18..622f12109ff 100644
--- a/test/test_tree_2d_part2.jl
+++ b/test/test_tree_2d_part2.jl
@@ -26,6 +26,9 @@ isdir(outdir) && rm(outdir, recursive = true)
     # Compressible Euler Multicomponent
     include("test_tree_2d_eulermulti.jl")
 
+    # Compressible Polytropic Euler
+    include("test_tree_2d_eulerpolytropic.jl")
+
     # Compressible Euler coupled with acoustic perturbation equations
     include("test_tree_2d_euleracoustics.jl")
 
diff --git a/test/test_tree_3d_advection.jl b/test/test_tree_3d_advection.jl
index 56278629417..ae53a2df52f 100644
--- a/test/test_tree_3d_advection.jl
+++ b/test/test_tree_3d_advection.jl
@@ -27,7 +27,10 @@ end
 @trixi_testset "elixir_advection_restart.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart.jl"),
                         l2=[0.00016017848135651983],
-                        linf=[0.0014175368788298393])
+                        linf=[0.0014175368788298393],
+                        # With the default `maxiters = 1` in coverage tests,
+                        # there would be no time steps after the restart.
+                        coverage_override=(maxiters = 100_000,))
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     let
diff --git a/test/test_trixi.jl b/test/test_trixi.jl
index 245efbc0175..cebe2164ae6 100644
--- a/test/test_trixi.jl
+++ b/test/test_trixi.jl
@@ -33,7 +33,7 @@ macro test_trixi_include(elixir, args...)
     local kwargs = Pair{Symbol, Any}[]
     for arg in args
         if (arg.head == :(=) &&
-            !(arg.args[1] in (:l2, :linf, :atol, :rtol, :coverage_override))
+            !(arg.args[1] in (:l2, :linf, :atol, :rtol, :coverage_override, :skip_coverage))
             && !(coverage && arg.args[1] in keys(coverage_override)))
             push!(kwargs, Pair(arg.args...))
         end
diff --git a/test/test_unit.jl b/test/test_unit.jl
index d4c11da635c..af91a365458 100644
--- a/test/test_unit.jl
+++ b/test/test_unit.jl
@@ -788,6 +788,54 @@ end
     end
 end
 
+@timed_testset "Consistency check for HLL flux with Davis wave speed estimates: Polytropic CEE" begin
+    flux_hll = FluxHLL(min_max_speed_davis)
+
+    gamma = 1.4
+    kappa = 0.5     # Scaling factor for the pressure.
+    equations = PolytropicEulerEquations2D(gamma, kappa)
+    u = SVector(1.1, -0.5, 2.34)
+
+    orientations = [1, 2]
+    for orientation in orientations
+        @test flux_hll(u, u, orientation, equations) ≈ flux(u, orientation, equations)
+    end
+
+    normal_directions = [SVector(1.0, 0.0),
+        SVector(0.0, 1.0),
+        SVector(0.5, -0.5),
+        SVector(-1.2, 0.3)]
+
+    for normal_direction in normal_directions
+        @test flux_hll(u, u, normal_direction, equations) ≈
+              flux(u, normal_direction, equations)
+    end
+end
+
+@timed_testset "Consistency check for Winters flux: Polytropic CEE" begin
+    for gamma in [1.4, 1.0, 5 / 3]
+        kappa = 0.5     # Scaling factor for the pressure.
+        equations = PolytropicEulerEquations2D(gamma, kappa)
+        u = SVector(1.1, -0.5, 2.34)
+
+        orientations = [1, 2]
+        for orientation in orientations
+            @test flux_winters_etal(u, u, orientation, equations) ≈
+                  flux(u, orientation, equations)
+        end
+
+        normal_directions = [SVector(1.0, 0.0),
+            SVector(0.0, 1.0),
+            SVector(0.5, -0.5),
+            SVector(-1.2, 0.3)]
+
+        for normal_direction in normal_directions
+            @test flux_winters_etal(u, u, normal_direction, equations) ≈
+                  flux(u, normal_direction, equations)
+        end
+    end
+end
+
 @timed_testset "Consistency check for HLL flux with Davis wave speed estimates: LEE" begin
     flux_hll = FluxHLL(min_max_speed_davis)
 
@@ -1332,6 +1380,103 @@ end
 
     @test mesh.boundary_faces[:entire_boundary] == [1, 2]
 end
+
+@testset "trixi_include" begin
+    @trixi_testset "Basic" begin
+        example = """
+            x = 4
+            """
+
+        filename = tempname()
+        try
+            open(filename, "w") do file
+                write(file, example)
+            end
+
+            # Use `@trixi_testset`, which wraps code in a temporary module, and call
+            # `trixi_include` with `@__MODULE__` in order to isolate this test.
+            @test_warn "You just called" trixi_include(@__MODULE__, filename)
+            @test @isdefined x
+            @test x == 4
+
+            @test_warn "You just called" trixi_include(@__MODULE__, filename, x = 7)
+            @test x == 7
+
+            @test_throws "assignment `y` not found in expression" trixi_include(@__MODULE__,
+                                                                                filename,
+                                                                                y = 3)
+        finally
+            rm(filename, force = true)
+        end
+    end
+
+    @trixi_testset "With `solve` Without `maxiters`" begin
+        # `trixi_include` assumes this to be the `solve` function of OrdinaryDiffEq,
+        # and therefore tries to insert the kwarg `maxiters`, which will fail here.
+        example = """
+            solve() = 0
+            x = solve()
+            """
+
+        filename = tempname()
+        try
+            open(filename, "w") do file
+                write(file, example)
+            end
+
+            # Use `@trixi_testset`, which wraps code in a temporary module, and call
+            # `trixi_include` with `@__MODULE__` in order to isolate this test.
+            @test_throws "no method matching solve(; maxiters::Int64)" trixi_include(@__MODULE__,
+                                                                                     filename)
+
+            @test_throws "no method matching solve(; maxiters::Int64)" trixi_include(@__MODULE__,
+                                                                                     filename,
+                                                                                     maxiters = 3)
+        finally
+            rm(filename, force = true)
+        end
+    end
+
+    @trixi_testset "With `solve` with `maxiters`" begin
+        # We need another example file that we include with `Base.include` first, in order to
+        # define the `solve` method without `trixi_include` trying to insert `maxiters` kwargs.
+        # Then, we can test that `trixi_include` inserts the kwarg in the `solve()` call.
+        example1 = """
+            solve(; maxiters=0) = maxiters
+            """
+
+        example2 = """
+            x = solve()
+            """
+
+        filename1 = tempname()
+        filename2 = tempname()
+        try
+            open(filename1, "w") do file
+                write(file, example1)
+            end
+            open(filename2, "w") do file
+                write(file, example2)
+            end
+
+            # Use `@trixi_testset`, which wraps code in a temporary module, and call
+            # `Base.include` and `trixi_include` with `@__MODULE__` in order to isolate this test.
+            Base.include(@__MODULE__, filename1)
+            @test_warn "You just called" trixi_include(@__MODULE__, filename2)
+            @test @isdefined x
+            # This is the default `maxiters` inserted by `trixi_include`
+            @test x == 10^5
+
+            @test_warn "You just called" trixi_include(@__MODULE__, filename2,
+                                                       maxiters = 7)
+            # Test that `maxiters` got overwritten
+            @test x == 7
+        finally
+            rm(filename1, force = true)
+            rm(filename2, force = true)
+        end
+    end
+end
 end
 
 end #module
diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl
index d4416ac5b6a..5341d86a7d1 100644
--- a/test/test_unstructured_2d.jl
+++ b/test/test_unstructured_2d.jl
@@ -129,7 +129,10 @@ end
                             0.005243995459478956,
                             0.004685630332338153,
                             0.01750217718347713,
-                        ])
+                        ],
+                        # With the default `maxiters = 1` in coverage tests,
+                        # there would be no time steps after the restart.
+                        coverage_override=(maxiters = 100_000,))
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     let
diff --git a/test/test_visualization.jl b/test/test_visualization.jl
index 48164a70fb3..6444dc91d5d 100644
--- a/test/test_visualization.jl
+++ b/test/test_visualization.jl
@@ -243,7 +243,6 @@ end
     @test_nowarn_mod trixi_include(@__MODULE__,
                                    joinpath(examples_dir(), "tree_2d_dgsem",
                                             "elixir_advection_basic.jl"),
-                                   tspan = (0, 0.1),
                                    analysis_callback = Trixi.TrivialCallback())
     @test adapt_to_mesh_level(sol, 5) isa Tuple
 
@@ -259,7 +258,6 @@ end
     @test_nowarn_mod trixi_include(@__MODULE__,
                                    joinpath(examples_dir(), "tree_3d_dgsem",
                                             "elixir_advection_basic.jl"),
-                                   tspan = (0, 0.1),
                                    analysis_callback = Trixi.TrivialCallback(),
                                    initial_refinement_level = 1)
     @test PlotData2D(sol) isa Trixi.PlotData2DCartesian
@@ -288,8 +286,7 @@ end
 
     @test_nowarn_mod trixi_include(@__MODULE__,
                                    joinpath(examples_dir(), "structured_3d_dgsem",
-                                            "elixir_advection_basic.jl"),
-                                   tspan = (0, 0.1))
+                                            "elixir_advection_basic.jl"))
 
     @testset "1D plot from 3D solution and general mesh" begin
         @testset "Create 1D plot as slice" begin