diff --git a/examples/1d/elixir_euler_density_wave.jl b/examples/1d/elixir_euler_density_wave.jl
new file mode 100644
index 0000000000..91fd897099
--- /dev/null
+++ b/examples/1d/elixir_euler_density_wave.jl
@@ -0,0 +1,53 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+equations = CompressibleEulerEquations1D(1.4)
+
+initial_condition = initial_condition_density_wave
+
+surface_flux = flux_lax_friedrichs
+solver = DGSEM(5, surface_flux)
+
+coordinates_min = -1
+coordinates_max =  1
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=2,
+                n_cells_max=30_000)
+
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 2.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+stepsize_callback = StepsizeCallback(cfl=0.8)
+
+save_solution = SaveSolutionCallback(interval=100,
+                                     save_initial_solution=true,
+                                     save_final_solution=true,
+                                     solution_variables=:primitive)
+# TODO: Taal, restart
+# restart_interval = 10
+
+analysis_interval = 2000
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
+
+callbacks = CallbackSet(summary_callback, stepsize_callback, save_solution, analysis_callback, alive_callback)
+
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), dt=stepsize_callback(ode),
+            save_everystep=false, callback=callbacks);
+summary_callback() # print the timer summary
diff --git a/examples/1d/elixir_euler_sedov_blast_wave_shockcapturing_amr.jl b/examples/1d/elixir_euler_sedov_blast_wave_shockcapturing_amr.jl
new file mode 100644
index 0000000000..5f1f089002
--- /dev/null
+++ b/examples/1d/elixir_euler_sedov_blast_wave_shockcapturing_amr.jl
@@ -0,0 +1,75 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations1D(1.4)
+
+initial_condition = initial_condition_sedov_blast_wave
+
+surface_flux = flux_lax_friedrichs
+volume_flux  = flux_chandrashekar
+basis = LobattoLegendreBasis(3)
+indicator_sc = IndicatorHennemannGassner(equations, basis,
+                                         alpha_max=0.5,
+                                         alpha_min=0.001,
+                                         alpha_smooth=true,
+                                         variable=density_pressure)
+volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+                                                 volume_flux_dg=volume_flux,
+                                                 volume_flux_fv=surface_flux)
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+coordinates_min = (-2,)
+coordinates_max = ( 2,)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=6,
+                n_cells_max=10_000)
+
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 12.5)
+ode = semidiscretize(semi, tspan)
+
+amr_indicator = IndicatorHennemannGassner(semi,
+                                          alpha_max=0.5,
+                                          alpha_min=0.001,
+                                          alpha_smooth=true,
+                                          variable=density_pressure)
+amr_controller = ControllerThreeLevel(semi, amr_indicator,
+                                      base_level=4,
+                                      max_level=6, max_threshold=0.01)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval=5,
+                           adapt_initial_condition=true,
+                           adapt_initial_condition_only_refine=true)
+
+summary_callback = SummaryCallback()
+
+stepsize_callback = StepsizeCallback(cfl=0.8)
+
+save_solution = SaveSolutionCallback(interval=100,
+                                     save_initial_solution=true,
+                                     save_final_solution=true,
+                                     solution_variables=:primitive)
+
+analysis_interval = 100
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
+
+callbacks = CallbackSet(summary_callback, stepsize_callback, save_solution, analysis_callback, alive_callback)
+
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), dt=stepsize_callback(ode),
+            save_everystep=false, callback=callbacks);
+summary_callback() # print the timer summary
diff --git a/examples/2d/elixir_advection_amr_nonperiodic.jl b/examples/2d/elixir_advection_amr_nonperiodic.jl
new file mode 100644
index 0000000000..f42942b201
--- /dev/null
+++ b/examples/2d/elixir_advection_amr_nonperiodic.jl
@@ -0,0 +1,74 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advectionvelocity = (1.0, 1.0)
+# advectionvelocity = (0.2, -0.3)
+equations = LinearScalarAdvectionEquation2D(advectionvelocity)
+
+initial_condition = initial_condition_gauss
+
+# you can either use a single function to impose the BCs weakly in all
+# 2*ndims == 4 directions or you can pass a tuple containing BCs for each direction
+# boundary_conditions = boundary_conditions_convergence_test
+boundary_conditions = boundary_conditions_gauss
+
+surface_flux = flux_lax_friedrichs
+solver = DGSEM(3, surface_flux)
+
+coordinates_min = (-5, -5)
+coordinates_max = ( 5,  5)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=4,
+                n_cells_max=30_000,
+                periodicity=false)
+
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions=boundary_conditions)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 10.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable=first),
+                                      base_level=4,
+                                      med_level=5, med_threshold=0.1,
+                                      max_level=6, max_threshold=0.6)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval=5,
+                           adapt_initial_condition=true,
+                           adapt_initial_condition_only_refine=true)
+
+stepsize_callback = StepsizeCallback(cfl=1.6)
+
+save_solution = SaveSolutionCallback(interval=100,
+                                     save_initial_solution=true,
+                                     save_final_solution=true,
+                                     solution_variables=:primitive)
+# TODO: Taal, IO
+# restart_interval = 10
+
+analysis_interval = 100
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval,
+                                     extra_analysis_integrals=(entropy,))
+
+# TODO: Taal decide, first AMR or save solution etc.
+callbacks = CallbackSet(summary_callback, amr_callback, stepsize_callback, save_solution, analysis_callback, alive_callback);
+
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), dt=stepsize_callback(ode),
+            save_everystep=false, callback=callbacks);
+summary_callback() # print the timer summary
diff --git a/examples/2d/parameters_advection_amr_nonperiodic.toml b/examples/2d/parameters_advection_amr_nonperiodic.toml
index 298bb81b07..892d59af2d 100644
--- a/examples/2d/parameters_advection_amr_nonperiodic.toml
+++ b/examples/2d/parameters_advection_amr_nonperiodic.toml
@@ -11,7 +11,7 @@ advectionvelocity = [1.0, 1.0] # Need only for linarscalaradvection
 initial_condition = "initial_condition_gauss"
 boundary_conditions = "boundary_condition_gauss"
 amr_indicator = "gauss"
-# surface_flux = "lax_riedrichs_flux"
+# surface_flux = "flux_lax_friedrichs"
 # source_terms =
 t_start = 0.0
 t_end   = 5.0
diff --git a/src/callbacks/stepsize_dg1d.jl b/src/callbacks/stepsize_dg1d.jl
index 0dd3f58062..c6bfb81ab6 100644
--- a/src/callbacks/stepsize_dg1d.jl
+++ b/src/callbacks/stepsize_dg1d.jl
@@ -1,7 +1,9 @@
 
 function max_dt(u::AbstractArray{<:Any,3}, t, mesh::TreeMesh{1},
                 constant_speed::Val{false}, equations, dg::DG, cache)
-  max_λ1 = max_λ2 = nextfloat(zero(t))
+  # Use `nextfloat` to avoid division by zero if lambda is ~zero (e.g., linear scalar advection with
+  # vanishing velocity)
+  max_λ1 = nextfloat(zero(t))
 
   for element in eachelement(dg, cache)
     inv_jacobian = cache.elements.inverse_jacobian[element]
@@ -18,6 +20,8 @@ end
 
 function max_dt(u::AbstractArray{<:Any,3}, t, mesh::TreeMesh{1},
                 constant_speed::Val{true}, equations, dg::DG, cache)
+  # Use `nextfloat` to avoid division by zero if lambda is ~zero (e.g., linear scalar advection with
+  # vanishing velocity)
   max_λ1 = nextfloat(zero(t))
 
   for element in eachelement(dg, cache)
diff --git a/src/callbacks/stepsize_dg2d.jl b/src/callbacks/stepsize_dg2d.jl
index 1edf02e651..931679cdcc 100644
--- a/src/callbacks/stepsize_dg2d.jl
+++ b/src/callbacks/stepsize_dg2d.jl
@@ -1,6 +1,8 @@
 
 function max_dt(u::AbstractArray{<:Any,4}, t, mesh::TreeMesh{2},
                 constant_speed::Val{false}, equations, dg::DG, cache)
+  # Use `nextfloat` to avoid division by zero if lambda is ~zero (e.g., linear scalar advection with
+  # vanishing velocity)
   max_λ1 = max_λ2 = nextfloat(zero(t))
 
   for element in eachelement(dg, cache)
@@ -19,6 +21,8 @@ end
 
 function max_dt(u::AbstractArray{<:Any,4}, t, mesh::TreeMesh{2},
                 constant_speed::Val{true}, equations, dg::DG, cache)
+  # Use `nextfloat` to avoid division by zero if lambda is ~zero (e.g., linear scalar advection with
+  # vanishing velocity)
   max_λ1 = max_λ2 = nextfloat(zero(t))
 
   for element in eachelement(dg, cache)
diff --git a/src/run.jl b/src/run.jl
index b65e56a4b3..21bfdc6eae 100644
--- a/src/run.jl
+++ b/src/run.jl
@@ -26,10 +26,10 @@ end
 # of `Trixi`. However, users will want to evaluate in the global scope of `Main` or something
 # similar to manage dependencies on their own.
 """
-    trixi_include([mod::Module=Main,] parameters_file; kwargs...)
+    trixi_include([mod::Module=Main,] elixir_file; kwargs...)
 
-`include` the file `parameters_file` and evaluate its content in the global scope of module `mod`.
-You can override specific assignments in `parameters_file` by supplying keyword arguments.
+`include` the file `elixir_file` and evaluate its content in the global scope of module `mod`.
+You can override specific assignments in `elixir_file` by supplying keyword arguments.
 It's basic purpose is to make it easier to modify some parameters while running Trixi from the
 REPL. Additionally, this is used in tests to reduce the computational burden for CI while still
 providing examples with sensible default values for users.
@@ -43,11 +43,11 @@ julia> sol.t[end]
 0.1
 ```
 """
-function trixi_include(mod::Module, parameters_file::AbstractString; kwargs...)
-  Base.include(ex -> replace_assignments(ex; kwargs...), mod, parameters_file)
+function trixi_include(mod::Module, elixir_file::AbstractString; kwargs...)
+  Base.include(ex -> replace_assignments(ex; kwargs...), mod, elixir_file)
 end
 
-trixi_include(parameters_file::AbstractString; kwargs...) = trixi_include(Main, parameters_file; kwargs...)
+trixi_include(elixir_file::AbstractString; kwargs...) = trixi_include(Main, elixir_file; kwargs...)
 
 
 
diff --git a/test/test_trixi.jl b/test/test_trixi.jl
index ca8e7e7d61..59baade3e8 100644
--- a/test/test_trixi.jl
+++ b/test/test_trixi.jl
@@ -30,23 +30,23 @@ end
 
 
 """
-    test_trixi_include(parameters_file; l2=nothing, linf=nothing, atol=10*eps(), rtol=0.001, parameters...)
+    test_trixi_include(elixir_file; l2=nothing, linf=nothing, atol=10*eps(), rtol=0.001, parameters...)
 
-Test Trixi by calling `trixi_include(parameters_file; parameters...)`.
+Test Trixi by calling `trixi_include(elixir_file; parameters...)`.
 By default, only the absence of error output is checked.
 If `l2` or `linf` are specified, in addition the resulting L2/Linf errors
 are compared approximately against these reference values, using `atol, rtol`
 as absolute/relative tolerance.
 """
-function test_trixi_include(parameters_file; l2=nothing, linf=nothing,
-                                             atol=200*eps(), rtol=0.001,
-                                             kwargs...)
+function test_trixi_include(elixir_file; l2=nothing, linf=nothing,
+                                         atol=200*eps(), rtol=0.001,
+                                         kwargs...)
 
   println("#"^80)
-  println(parameters_file)
+  println(elixir_file)
 
   # evaluate examples in the scope of the module they're called from
-  @test_nowarn trixi_include(@__MODULE__, parameters_file; kwargs...)
+  @test_nowarn trixi_include(@__MODULE__, elixir_file; kwargs...)
 
   # If present, compare L2 and Linf errors against reference values
   if !isnothing(l2) || !isnothing(linf)