diff --git a/docs/literate/src/files/adding_nonconservative_equation.jl b/docs/literate/src/files/adding_nonconservative_equation.jl
index 1029c31408b..631dfa334a1 100644
--- a/docs/literate/src/files/adding_nonconservative_equation.jl
+++ b/docs/literate/src/files/adding_nonconservative_equation.jl
@@ -145,7 +145,7 @@ plot(sol)
 # above.
 
 error_1 = analysis_callback(sol).l2 |> first
-@test isapprox(error_1, 0.00029609575838969394) #src
+@test isapprox(error_1, 0.0002961027497380263) #src
 # Next, we increase the grid resolution by one refinement level and run the
 # simulation again.
 
diff --git a/examples/p4est_2d_dgsem/elixir_euler_riemannproblem_quadrants_amr.jl b/examples/p4est_2d_dgsem/elixir_euler_riemannproblem_quadrants_amr.jl
new file mode 100644
index 00000000000..138190819cd
--- /dev/null
+++ b/examples/p4est_2d_dgsem/elixir_euler_riemannproblem_quadrants_amr.jl
@@ -0,0 +1,128 @@
+using OrdinaryDiffEqSSPRK
+using Trixi
+
+###############################################################################
+## Semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+# Variant of the 4-quadrant Riemann problem considered in 
+# - Carsten W. Schulz-Rinne:
+#   Classification of the Riemann Problem for Two-Dimensional Gas Dynamics
+#   https://doi.org/10.1137/0524006
+# and 
+# - Carsten W. Schulz-Rinne, James P. Collins, and Harland M. Glaz
+#   Numerical Solution of the Riemann Problem for Two-Dimensional Gas Dynamics
+#   https://doi.org/10.1137/0914082
+
+# Specify the initial condition as a discontinuous initial condition (see docstring of 
+# `DiscontinuousFunction` for more information) which comes with a specialized 
+# initialization routine suited for Riemann problems.
+# In short, if a discontinuity is right at an interface, the boundary nodes (which are at the same location)
+# on that interface will be initialized with the left and right state of the discontinuity, i.e., 
+#                         { u_1, if element = left element and x_{element}^{(n)} = x_jump
+# u(x_jump, t, element) = {
+#                         { u_2, if element = right element and x_{element}^{(1)} = x_jump
+# This is realized by shifting the outer DG nodes inwards, i.e., on reference element
+# the outer nodes at `[-1, 1]` are shifted inwards to `[-1 + ε, 1 - ε]` with machine precision `ε`.
+struct InitialConditionRP <: DiscontinuousFunction end
+
+function (initial_condition_rp::InitialConditionRP)(x_, t,
+                                                    equations::CompressibleEulerEquations2D)
+    x, y = x_[1], x_[2]
+
+    if x >= 0.5 && y >= 0.5
+        rho, v1, v2, p = (0.5313, 0.0, 0.0, 0.4)
+    elseif x < 0.5 && y >= 0.5
+        rho, v1, v2, p = (1.0, 0.7276, 0.0, 1.0)
+    elseif x < 0.5 && y < 0.5
+        rho, v1, v2, p = (0.8, 0.0, 0.0, 1.0)
+    elseif x >= 0.5 && y < 0.5
+        rho, v1, v2, p = (1.0, 0.0, 0.7276, 1.0)
+    end
+
+    prim = SVector(rho, v1, v2, p)
+    return prim2cons(prim, equations)
+end
+# Note calling the constructor of the struct: `InitialConditionRP()` instead of `initial_condition_rp` !
+initial_condition = InitialConditionRP()
+
+# Extend domain by specifying free outflow
+@inline function boundary_condition_outflow(u_inner, normal_direction::AbstractVector, x, t,
+                                            surface_flux_function,
+                                            equations::CompressibleEulerEquations2D)
+    flux = Trixi.flux(u_inner, normal_direction, equations)
+
+    return flux
+end
+
+boundary_conditions = Dict(:x_neg => boundary_condition_outflow,
+                           :x_pos => boundary_condition_outflow,
+                           :y_neg => boundary_condition_outflow,
+                           :y_pos => boundary_condition_outflow)
+
+coordinates_min = (0.0, 0.0)
+coordinates_max = (1.0, 1.0)
+
+trees_per_dimension = (4, 4) # Set base discretization relatively coarse to allow large cells
+mesh = P4estMesh(trees_per_dimension, polydeg = 1,
+                 coordinates_min = coordinates_min, coordinates_max = coordinates_max,
+                 periodicity = (false, false),
+                 initial_refinement_level = 2)
+
+surface_flux = flux_hllc
+volume_flux = flux_ranocha
+
+polydeg = 3
+basis = LobattoLegendreBasis(polydeg)
+shock_indicator = IndicatorHennemannGassner(equations, basis,
+                                            alpha_max = 0.5,
+                                            alpha_min = 0.001,
+                                            alpha_smooth = true,
+                                            variable = Trixi.density)
+
+volume_integral = VolumeIntegralShockCapturingHG(shock_indicator;
+                                                 volume_flux_dg = volume_flux,
+                                                 volume_flux_fv = surface_flux)
+
+solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
+               volume_integral = volume_integral)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver;
+                                    boundary_conditions = boundary_conditions)
+
+###############################################################################
+## ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.25)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+stepsize_callback = StepsizeCallback(cfl = 1.8)
+
+amr_indicator = IndicatorLöhner(semi, variable = Trixi.density)
+amr_controller = ControllerThreeLevel(semi, amr_indicator,
+                                      base_level = 0,
+                                      med_level = 3, med_threshold = 0.02,
+                                      max_level = 6, max_threshold = 0.04)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval = 10,
+                           adapt_initial_condition = true,
+                           adapt_initial_condition_only_refine = true)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        amr_callback,
+                        stepsize_callback)
+
+###############################################################################
+## Run the simulation
+
+sol = solve(ode, SSPRK54();
+            dt = 1.0, save_everystep = false, callback = callbacks);
diff --git a/examples/p4est_3d_dgsem/elixir_euler_sedov.jl b/examples/p4est_3d_dgsem/elixir_euler_sedov.jl
index c7532de5a51..b1e3de4bdb0 100644
--- a/examples/p4est_3d_dgsem/elixir_euler_sedov.jl
+++ b/examples/p4est_3d_dgsem/elixir_euler_sedov.jl
@@ -6,15 +6,28 @@ using Trixi
 
 equations = CompressibleEulerEquations3D(1.4)
 
+# Specify the initial condition as a discontinuous initial condition (see docstring of 
+# `DiscontinuousFunction` for more information) which comes with a specialized 
+# initialization routine suited for the Riemann problems.
+# In short, if a discontinuity is right at an interface, the boundary nodes (which are at the same location)
+# on that interface will be initialized with the left and right state of the discontinuity, i.e., 
+#                         { u_1, if element = left element and x_{element}^{(n)} = x_jump
+# u(x_jump, t, element) = {
+#                         { u_2, if element = right element and x_{element}^{(1)} = x_jump
+# This is realized by shifting the outer DG nodes inwards, i.e., on reference element
+# the outer nodes are at `[-1, 1]` are shifted to `[-1 + ε, 1 - ε]` with machine precision `ε`.
+struct InitialConditionMediumSedovBlast <: DiscontinuousFunction end
+
 """
-    initial_condition_medium_sedov_blast_wave(x, t, equations::CompressibleEulerEquations3D)
+    (initial_condition_medium_sedov_blast_wave::InitialConditionMediumSedovBlast)(x, t, 
+                                                                                  equations::CompressibleEulerEquations3D)
 
 The Sedov blast wave setup based on Flash
 - https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000
 with smaller strength of the initial discontinuity.
 """
-function initial_condition_medium_sedov_blast_wave(x, t,
-                                                   equations::CompressibleEulerEquations3D)
+function (initial_condition_medium_sedov_blast_wave::InitialConditionMediumSedovBlast)(x, t,
+                                                                                       equations::CompressibleEulerEquations3D)
     # Set up polar coordinates
     inicenter = SVector(0.0, 0.0, 0.0)
     x_norm = x[1] - inicenter[1]
@@ -37,8 +50,9 @@ function initial_condition_medium_sedov_blast_wave(x, t,
 
     return prim2cons(SVector(rho, v1, v2, v3, p), equations)
 end
-
-initial_condition = initial_condition_medium_sedov_blast_wave
+# Note calling the constructor of the struct: `InitialConditionMediumSedovBlast()` instead of 
+# `initial_condition_medium_sedov_blast_wave` !
+initial_condition = InitialConditionMediumSedovBlast()
 
 surface_flux = flux_lax_friedrichs
 volume_flux = flux_ranocha
diff --git a/examples/structured_1d_dgsem/elixir_advection_shockcapturing.jl b/examples/structured_1d_dgsem/elixir_advection_shockcapturing.jl
index f69f5e4061c..c9e17a971f1 100644
--- a/examples/structured_1d_dgsem/elixir_advection_shockcapturing.jl
+++ b/examples/structured_1d_dgsem/elixir_advection_shockcapturing.jl
@@ -6,8 +6,21 @@ using Trixi
 
 advection_velocity = 1.0
 
+# Specify the initial condition as a discontinuous initial condition (see docstring of 
+# `DiscontinuousFunction` for more information) which comes with a specialized 
+# initialization routine suited for Riemann problems.
+# In short, if a discontinuity is right at an interface, the boundary nodes (which are at the same location)
+# on that interface will be initialized with the left and right state of the discontinuity, i.e., 
+#                         { u_1, if element = left element and x_{element}^{(n)} = x_jump
+# u(x_jump, t, element) = {
+#                         { u_2, if element = right element and x_{element}^{(1)} = x_jump
+# This is realized by shifting the outer DG nodes inwards, i.e., on reference element
+# the outer nodes at `[-1, 1]` are shifted inwards to `[-1 + ε, 1 - ε]` with machine precision `ε`.
+struct InitialConditionComposite <: DiscontinuousFunction end
+
 """
-    initial_condition_composite(x, t, equations::LinearScalarAdvectionEquation1D)
+    (initial_condition_composite::InitialConditionComposite)(x, t, 
+                                                             equations::LinearScalarAdvectionEquation1D)
 
 Wave form that is a combination of a Gaussian pulse, a square wave, a triangle wave,
 and half an ellipse with periodic boundary conditions.
@@ -16,7 +29,8 @@ Slight simplification from
   Efficient Implementation of Weighted ENO Schemes
   [DOI: 10.1006/jcph.1996.0130](https://doi.org/10.1006/jcph.1996.0130)
 """
-function initial_condition_composite(x, t, equations::LinearScalarAdvectionEquation1D)
+function (initial_condition_composite::InitialConditionComposite)(x, t,
+                                                                  equations::LinearScalarAdvectionEquation1D)
     xmin, xmax = -1.0, 1.0 # Only works if the domain is [-1.0,1.0]
     x_trans = x[1] - t
     L = xmax - xmin
@@ -42,8 +56,8 @@ function initial_condition_composite(x, t, equations::LinearScalarAdvectionEquat
 
     return SVector(value)
 end
-
-initial_condition = initial_condition_composite
+# Note calling the constructor of the struct: `InitialConditionComposite()` instead of `initial_condition_composite` !
+initial_condition = InitialConditionComposite()
 
 equations = LinearScalarAdvectionEquation1D(advection_velocity)
 
diff --git a/examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl b/examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl
index ff28db9c77e..45fa3924a46 100644
--- a/examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl
+++ b/examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl
@@ -14,15 +14,31 @@ cells_per_dimension = (64,)
 mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max,
                       periodicity = false)
 
+# Specify the initial condition as a discontinuous initial condition (see docstring of 
+# `DiscontinuousFunction` for more information) which comes with a specialized 
+# initialization routine suited for Riemann problems.
+# In short, if a discontinuity is right at an interface, the boundary nodes (which are at the same location)
+# on that interface will be initialized with the left and right state of the discontinuity, i.e., 
+#                         { u_1, if element = left element and x_{element}^{(n)} = x_jump
+# u(x_jump, t, element) = {
+#                         { u_2, if element = right element and x_{element}^{(1)} = x_jump
+# This is realized by shifting the outer DG nodes inwards, i.e., on reference element
+# the outer nodes at `[-1, 1]` are shifted inwards to `[-1 + ε, 1 - ε]` with machine precision `ε`.
+struct InitialConditionGreenlight <: DiscontinuousFunction end
+
 # Example inspired from http://www.clawpack.org/riemann_book/html/Traffic_flow.html#Example:-green-light
 # Green light that at x = 0 which switches at t = 0 from red to green.
 # To the left there are cars bumper to bumper, to the right there are no cars.
-function initial_condition_greenlight(x, t, equation::TrafficFlowLWREquations1D)
+function (initial_condition_greenlight::InitialConditionGreenlight)(x, t,
+                                                                    equation::TrafficFlowLWREquations1D)
     RealT = eltype(x)
     scalar = x[1] < 0 ? one(RealT) : zero(RealT)
 
     return SVector(scalar)
 end
+# Note calling the constructor of the struct: `InitialConditionGreenlight()` instead of
+# `initial_condition_greenlight` !
+initial_condition = InitialConditionGreenlight()
 
 ###############################################################################
 # Specify non-periodic boundary conditions
@@ -30,7 +46,7 @@ end
 # Assume that there are always cars waiting at the left
 function inflow(x, t, equations::TrafficFlowLWREquations1D)
     # -1.0 = coordinates_min
-    return initial_condition_greenlight(-1.0, t, equations)
+    return initial_condition(-1.0, t, equations)
 end
 boundary_condition_inflow = BoundaryConditionDirichlet(inflow)
 
@@ -47,8 +63,6 @@ end
 boundary_conditions = (x_neg = boundary_condition_inflow,
                        x_pos = boundary_condition_outflow)
 
-initial_condition = initial_condition_greenlight
-
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                     boundary_conditions = boundary_conditions)
 
diff --git a/examples/tree_1d_dgsem/elixir_burgers_rarefaction.jl b/examples/tree_1d_dgsem/elixir_burgers_rarefaction.jl
index a234521d1e6..5365a17615c 100644
--- a/examples/tree_1d_dgsem/elixir_burgers_rarefaction.jl
+++ b/examples/tree_1d_dgsem/elixir_burgers_rarefaction.jl
@@ -32,18 +32,34 @@ mesh = TreeMesh(coordinate_min, coordinate_max,
                 n_cells_max = 10_000,
                 periodicity = false)
 
+# Specify the initial condition as a discontinuous initial condition (see docstring of 
+# `DiscontinuousFunction` for more information) which comes with a specialized 
+# initialization routine suited for Riemann problems.
+# In short, if a discontinuity is right at an interface, the boundary nodes (which are at the same location)
+# on that interface will be initialized with the left and right state of the discontinuity, i.e., 
+#                         { u_1, if element = left element and x_{element}^{(n)} = x_jump
+# u(x_jump, t, element) = {
+#                         { u_2, if element = right element and x_{element}^{(1)} = x_jump
+# This is realized by shifting the outer DG nodes inwards, i.e., on reference element
+# the outer nodes at `[-1, 1]` are shifted inwards to `[-1 + ε, 1 - ε]` with machine precision `ε`.
+struct InitialConditionRarefaction <: DiscontinuousFunction end
+
 # Discontinuous initial condition (Riemann Problem) leading to a rarefaction fan.
-function initial_condition_rarefaction(x, t, equation::InviscidBurgersEquation1D)
+function (initial_condition_rarefaction::InitialConditionRarefaction)(x, t,
+                                                                      equation::InviscidBurgersEquation1D)
     RealT = eltype(x)
     scalar = x[1] < 0.5f0 ? convert(RealT, 0.5f0) : convert(RealT, 1.5f0)
 
     return SVector(scalar)
 end
+# Note calling the constructor of the struct: `InitialConditionRarefaction()` instead of 
+# `initial_condition_rarefaction` !
+initial_condition = InitialConditionRarefaction()
 
 ###############################################################################
 # Specify non-periodic boundary conditions
 
-boundary_condition_inflow = BoundaryConditionDirichlet(initial_condition_rarefaction)
+boundary_condition_inflow = BoundaryConditionDirichlet(initial_condition)
 
 function boundary_condition_outflow(u_inner, orientation, normal_direction, x, t,
                                     surface_flux_function,
@@ -57,8 +73,6 @@ end
 boundary_conditions = (x_neg = boundary_condition_inflow,
                        x_pos = boundary_condition_outflow)
 
-initial_condition = initial_condition_rarefaction
-
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                     boundary_conditions = boundary_conditions)
 
diff --git a/examples/tree_1d_dgsem/elixir_burgers_shock.jl b/examples/tree_1d_dgsem/elixir_burgers_shock.jl
index 4f48af081a7..b0d9f1f9470 100644
--- a/examples/tree_1d_dgsem/elixir_burgers_shock.jl
+++ b/examples/tree_1d_dgsem/elixir_burgers_shock.jl
@@ -32,18 +32,33 @@ mesh = TreeMesh(coordinate_min, coordinate_max,
                 n_cells_max = 10_000,
                 periodicity = false)
 
+# Specify the initial condition as a discontinuous initial condition (see docstring of 
+# `DiscontinuousFunction` for more information) which comes with a specialized 
+# initialization routine suited for Riemann problems.
+# In short, if a discontinuity is right at an interface, the boundary nodes (which are at the same location)
+# on that interface will be initialized with the left and right state of the discontinuity, i.e., 
+#                         { u_1, if element = left element and x_{element}^{(n)} = x_jump
+# u(x_jump, t, element) = {
+#                         { u_2, if element = right element and x_{element}^{(1)} = x_jump
+# This is realized by shifting the outer DG nodes inwards, i.e., on reference element
+# the outer nodes at `[-1, 1]` are shifted inwards to `[-1 + ε, 1 - ε]` with machine precision `ε`.
+struct InitialConditionShock <: DiscontinuousFunction end
+
 # Discontinuous initial condition (Riemann Problem) leading to a shock to test e.g. correct shock speed.
-function initial_condition_shock(x, t, equation::InviscidBurgersEquation1D)
+function (initial_condition_shock::InitialConditionShock)(x, t,
+                                                          equation::InviscidBurgersEquation1D)
     RealT = eltype(x)
     scalar = x[1] < 0.5f0 ? convert(RealT, 1.5f0) : convert(RealT, 0.5f0)
 
     return SVector(scalar)
 end
+# Note calling the constructor of the struct: `InitialConditionShock()` instead of `initial_condition_shock` !
+initial_condition = InitialConditionShock()
 
 ###############################################################################
 # Specify non-periodic boundary conditions
 
-boundary_condition_inflow = BoundaryConditionDirichlet(initial_condition_shock)
+boundary_condition_inflow = BoundaryConditionDirichlet(initial_condition)
 
 function boundary_condition_outflow(u_inner, orientation, normal_direction, x, t,
                                     surface_flux_function,
@@ -57,8 +72,6 @@ end
 boundary_conditions = (x_neg = boundary_condition_inflow,
                        x_pos = boundary_condition_outflow)
 
-initial_condition = initial_condition_shock
-
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                     boundary_conditions = boundary_conditions)
 
diff --git a/examples/tree_1d_dgsem/elixir_euler_blast_wave.jl b/examples/tree_1d_dgsem/elixir_euler_blast_wave.jl
index 6689ff6a8e6..068fd225f59 100644
--- a/examples/tree_1d_dgsem/elixir_euler_blast_wave.jl
+++ b/examples/tree_1d_dgsem/elixir_euler_blast_wave.jl
@@ -6,15 +6,29 @@ using Trixi
 
 equations = CompressibleEulerEquations1D(1.4)
 
+# Specify the initial condition as a discontinuous initial condition (see docstring of 
+# `DiscontinuousFunction` for more information) which comes with a specialized 
+# initialization routine suited for Riemann problems.
+# In short, if a discontinuity is right at an interface, the boundary nodes (which are at the same location)
+# on that interface will be initialized with the left and right state of the discontinuity, i.e., 
+#                         { u_1, if element = left element and x_{element}^{(n)} = x_jump
+# u(x_jump, t, element) = {
+#                         { u_2, if element = right element and x_{element}^{(1)} = x_jump
+# This is realized by shifting the outer DG nodes inwards, i.e., on reference element
+# the outer nodes at `[-1, 1]` are shifted inwards to `[-1 + ε, 1 - ε]` with machine precision `ε`.
+struct InitialConditionBlastWave <: DiscontinuousFunction end
+
 """
-    initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations1D)
+    (initial_condition_blast_wave::InitialConditionBlastWave)(x, t,
+                                                              equations::CompressibleEulerEquations1D)
 
 A medium blast wave taken from
 - Sebastian Hennemann, Gregor J. Gassner (2020)
   A provably entropy stable subcell shock capturing approach for high order split form DG
   [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044)
 """
-function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations1D)
+function (initial_condition_blast_wave::InitialConditionBlastWave)(x, t,
+                                                                   equations::CompressibleEulerEquations1D)
     # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave"
     # Set up polar coordinates
     RealT = eltype(x)
@@ -34,7 +48,9 @@ function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquation
 
     return prim2cons(SVector(rho, v1, p), equations)
 end
-initial_condition = initial_condition_blast_wave
+# Note calling the constructor of the struct: `InitialConditionBlastWave()` instead of
+# `initial_condition_blast_wave` !
+initial_condition = InitialConditionBlastWave()
 
 surface_flux = flux_lax_friedrichs
 volume_flux = flux_ranocha
diff --git a/examples/tree_1d_dgsem/elixir_euler_blast_wave_entropy_bounded.jl b/examples/tree_1d_dgsem/elixir_euler_blast_wave_entropy_bounded.jl
index 2d6e31d874c..f79090ced5b 100644
--- a/examples/tree_1d_dgsem/elixir_euler_blast_wave_entropy_bounded.jl
+++ b/examples/tree_1d_dgsem/elixir_euler_blast_wave_entropy_bounded.jl
@@ -6,15 +6,29 @@ using Trixi
 
 equations = CompressibleEulerEquations1D(1.4)
 
+# Specify the initial condition as a discontinuous initial condition (see docstring of 
+# `DiscontinuousFunction` for more information) which comes with a specialized 
+# initialization routine suited for Riemann problems.
+# In short, if a discontinuity is right at an interface, the boundary nodes (which are at the same location)
+# on that interface will be initialized with the left and right state of the discontinuity, i.e., 
+#                         { u_1, if element = left element and x_{element}^{(n)} = x_jump
+# u(x_jump, t, element) = {
+#                         { u_2, if element = right element and x_{element}^{(1)} = x_jump
+# This is realized by shifting the outer DG nodes inwards, i.e., on reference element
+# the outer nodes at `[-1, 1]` are shifted inwards to `[-1 + ε, 1 - ε]` with machine precision `ε`.
+struct InitialConditionBlastWave <: DiscontinuousFunction end
+
 """
-    initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations1D)
+    (initial_condition_blast_wave::InitialConditionBlastWave)(x, t,
+                                                              equations::CompressibleEulerEquations1D)
 
 A medium blast wave taken from
 - Sebastian Hennemann, Gregor J. Gassner (2020)
   A provably entropy stable subcell shock capturing approach for high order split form DG
   [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044)
 """
-function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations1D)
+function (initial_condition_blast_wave::InitialConditionBlastWave)(x, t,
+                                                                   equations::CompressibleEulerEquations1D)
     # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave"
     # Set up polar coordinates
     RealT = eltype(x)
@@ -34,7 +48,9 @@ function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquation
 
     return prim2cons(SVector(rho, v1, p), equations)
 end
-initial_condition = initial_condition_blast_wave
+# Note calling the constructor of the struct: `InitialConditionBlastWave()` instead of
+# `initial_condition_blast_wave` !
+initial_condition = InitialConditionBlastWave()
 
 # Note: We do not need to use the shock-capturing methodology here,
 # in contrast to the standard `euler_blast_wave.jl` example.
diff --git a/examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl
index be06997061b..98433fb1790 100644
--- a/examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl
+++ b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl
@@ -7,8 +7,21 @@ using Trixi
 
 equations = CompressibleEulerEquationsQuasi1D(1.4)
 
+# Specify the initial condition as a discontinuous initial condition (see docstring of 
+# `DiscontinuousFunction` for more information) which comes with a specialized 
+# initialization routine suited for Riemann problems.
+# In short, if a discontinuity is right at an interface, the boundary nodes (which are at the same location)
+# on that interface will be initialized with the left and right state of the discontinuity, i.e., 
+#                         { u_1, if element = left element and x_{element}^{(n)} = x_jump
+# u(x_jump, t, element) = {
+#                         { u_2, if element = right element and x_{element}^{(1)} = x_jump
+# This is realized by shifting the outer DG nodes inwards, i.e., on reference element
+# the outer nodes at `[-1, 1]` are shifted inwards to `[-1 + ε, 1 - ε]` with machine precision `ε`.
+struct InitialConditionDiscontinuity <: DiscontinuousFunction end
+
 """
-    initial_condition_discontinuity(x, t, equations::CompressibleEulerEquations1D)
+    (initial_condition_discontinuity::InitialConditionDiscontinuity)(x, t,
+                                                                     equations::CompressibleEulerEquations1D)
 
 A discontinuous initial condition taken from
 - Jesse Chan, Khemraj Shukla, Xinhui Wu, Ruofeng Liu, Prani Nalluri (2023)
@@ -16,8 +29,8 @@ A discontinuous initial condition taken from
     shallow water and compressible Euler equations
     [DOI: 10.48550/arXiv.2307.12089](https://doi.org/10.48550/arXiv.2307.12089)
 """
-function initial_condition_discontinuity(x, t,
-                                         equations::CompressibleEulerEquationsQuasi1D)
+function (initial_condition_discontinuity::InitialConditionDiscontinuity)(x, t,
+                                                                          equations::CompressibleEulerEquationsQuasi1D)
     RealT = eltype(x)
     rho = (x[1] < 0) ? RealT(3.4718) : RealT(2.0)
     v1 = (x[1] < 0) ? RealT(-2.5923) : RealT(-3.0)
@@ -26,8 +39,9 @@ function initial_condition_discontinuity(x, t,
 
     return prim2cons(SVector(rho, v1, p, a), equations)
 end
-
-initial_condition = initial_condition_discontinuity
+# Note calling the constructor of the struct: `InitialConditionDiscontinuity()` instead of
+# `initial_condition_discontinuity` !
+initial_condition = InitialConditionDiscontinuity()
 
 surface_flux = (flux_lax_friedrichs, flux_nonconservative_chan_etal)
 volume_flux = (flux_chan_etal, flux_nonconservative_chan_etal)
diff --git a/examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave_pure_fv.jl b/examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave_pure_fv.jl
index efa5e238d9b..badcc7555c2 100644
--- a/examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave_pure_fv.jl
+++ b/examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave_pure_fv.jl
@@ -6,13 +6,27 @@ using Trixi
 
 equations = CompressibleEulerEquations1D(1.4)
 
+# Specify the initial condition as a discontinuous initial condition (see docstring of 
+# `DiscontinuousFunction` for more information) which comes with a specialized 
+# initialization routine suited for Riemann problems.
+# In short, if a discontinuity is right at an interface, the boundary nodes (which are at the same location)
+# on that interface will be initialized with the left and right state of the discontinuity, i.e., 
+#                         { u_1, if element = left element and x_{element}^{(n)} = x_jump
+# u(x_jump, t, element) = {
+#                         { u_2, if element = right element and x_{element}^{(1)} = x_jump
+# This is realized by shifting the outer DG nodes inwards, i.e., on reference element
+# the outer nodes at `[-1, 1]` are shifted inwards to `[-1 + ε, 1 - ε]` with machine precision `ε`.
+struct InitialConditionSedovBlastWave <: DiscontinuousFunction end
+
 """
-    initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations1D)
+    (initial_condition_sedov_blast_wave::InitialConditionSedovBlastWave)(x, t, 
+                                                                         equations::CompressibleEulerEquations1D)
 
 The Sedov blast wave setup based on Flash
 - https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000
 """
-function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations1D)
+function (initial_condition_sedov_blast_wave::InitialConditionSedovBlastWave)(x, t,
+                                                                              equations::CompressibleEulerEquations1D)
     # Set up polar coordinates
     RealT = eltype(x)
     inicenter = SVector(0)
@@ -34,7 +48,9 @@ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEq
 
     return prim2cons(SVector(rho, v1, p), equations)
 end
-initial_condition = initial_condition_sedov_blast_wave
+# Note calling the constructor of the struct: `InitialConditionSedovBlastWave()` instead of
+# `initial_condition_sedov_blast_wave` !
+initial_condition = InitialConditionSedovBlastWave()
 
 surface_flux = flux_hllc
 basis = LobattoLegendreBasis(3)
diff --git a/examples/tree_1d_dgsem/elixir_eulermulti_two_interacting_blast_waves.jl b/examples/tree_1d_dgsem/elixir_eulermulti_two_interacting_blast_waves.jl
index ce7e7702d73..0fa941732ac 100644
--- a/examples/tree_1d_dgsem/elixir_eulermulti_two_interacting_blast_waves.jl
+++ b/examples/tree_1d_dgsem/elixir_eulermulti_two_interacting_blast_waves.jl
@@ -7,16 +7,28 @@ using Trixi
 equations = CompressibleEulerMulticomponentEquations1D(gammas = (1.4, 1.4, 1.4),
                                                        gas_constants = (0.4, 0.4, 0.4))
 
+# Specify the initial condition as a discontinuous initial condition (see docstring of 
+# `DiscontinuousFunction` for more information) which comes with a specialized 
+# initialization routine suited for Riemann problems.
+# In short, if a discontinuity is right at an interface, the boundary nodes (which are at the same location)
+# on that interface will be initialized with the left and right state of the discontinuity, i.e., 
+#                         { u_1, if element = left element and x_{element}^{(n)} = x_jump
+# u(x_jump, t, element) = {
+#                         { u_2, if element = right element and x_{element}^{(1)} = x_jump
+# This is realized by shifting the outer DG nodes inwards, i.e., on reference element
+# the outer nodes at `[-1, 1]` are shifted inwards to `[-1 + ε, 1 - ε]` with machine precision `ε`.
+struct InitialConditionTwoBlastWaves <: DiscontinuousFunction end
 """
-    initial_condition_two_interacting_blast_waves(x, t, equations::CompressibleEulerMulticomponentEquations1D)
+    (initial_condition_two_blast_waves::InitialConditionTwoBlastWaves)(x, t,
+                                                                       equations::CompressibleEulerMulticomponentEquations1D)
 
 A multicomponent two interacting blast wave test taken from
 - T. Plewa & E. Müller (1999)
   The consistent multi-fluid advection method
   [arXiv: 9807241](https://arxiv.org/pdf/astro-ph/9807241.pdf)
 """
-function initial_condition_two_interacting_blast_waves(x, t,
-                                                       equations::CompressibleEulerMulticomponentEquations1D)
+function (initial_condition_two_blast_waves::InitialConditionTwoBlastWaves)(x, t,
+                                                                            equations::CompressibleEulerMulticomponentEquations1D)
     RealT = eltype(x)
     rho1 = 0.5f0 * x[1]^2
     rho2 = 0.5f0 * (sin(20 * x[1]))^2
@@ -38,7 +50,9 @@ function initial_condition_two_interacting_blast_waves(x, t,
 
     return prim2cons(vcat(prim_other, prim_rho), equations)
 end
-initial_condition = initial_condition_two_interacting_blast_waves
+# Note calling the constructor of the struct: `InitialConditionTwoBlastWaves()` instead of 
+# `initial_condition_two_blast_waves` !
+initial_condition = InitialConditionTwoBlastWaves()
 
 function boundary_condition_two_interacting_blast_waves(u_inner, orientation, direction,
                                                         x, t, surface_flux_function,
diff --git a/examples/tree_1d_dgsem/elixir_mhd_briowu_shock_tube.jl b/examples/tree_1d_dgsem/elixir_mhd_briowu_shock_tube.jl
index c8626a8e6e9..07d5c0cddad 100644
--- a/examples/tree_1d_dgsem/elixir_mhd_briowu_shock_tube.jl
+++ b/examples/tree_1d_dgsem/elixir_mhd_briowu_shock_tube.jl
@@ -6,8 +6,21 @@ using Trixi
 gamma = 2
 equations = IdealGlmMhdEquations1D(gamma)
 
+# Specify the initial condition as a discontinuous initial condition (see docstring of 
+# `DiscontinuousFunction` for more information) which comes with a specialized 
+# initialization routine suited for Riemann problems.
+# In short, if a discontinuity is right at an interface, the boundary nodes (which are at the same location)
+# on that interface will be initialized with the left and right state of the discontinuity, i.e., 
+#                         { u_1, if element = left element and x_{element}^{(n)} = x_jump
+# u(x_jump, t, element) = {
+#                         { u_2, if element = right element and x_{element}^{(1)} = x_jump
+# This is realized by shifting the outer DG nodes inwards, i.e., on reference element
+# the outer nodes at `[-1, 1]` are shifted inwards to `[-1 + ε, 1 - ε]` with machine precision `ε`.
+struct InitialConditionBrioWuShockTube <: DiscontinuousFunction end
+
 """
-    initial_condition_briowu_shock_tube(x, t, equations::IdealGlmMhdEquations1D)
+    (initial_condition_briowu_shock_tube::InitialConditionBrioWuShockTube)(x, t,
+                                                                           equations::IdealGlmMhdEquations1D)
 
 Compound shock tube test case for one dimensional ideal MHD equations. It is basically an
 MHD extension of the Sod shock tube. Taken from Section V of the article
@@ -15,7 +28,8 @@ MHD extension of the Sod shock tube. Taken from Section V of the article
   An Upwind Differencing Scheme for the Equations of Ideal Magnetohydrodynamics
   [DOI: 10.1016/0021-9991(88)90120-9](https://doi.org/10.1016/0021-9991(88)90120-9)
 """
-function initial_condition_briowu_shock_tube(x, t, equations::IdealGlmMhdEquations1D)
+function (initial_condition_briowu_shock_tube::InitialConditionBrioWuShockTube)(x, t,
+                                                                                equations::IdealGlmMhdEquations1D)
     # domain must be set to [0, 1], γ = 2, final time = 0.12
     RealT = eltype(x)
     rho = x[1] < 0.5f0 ? 1.0f0 : 0.125f0
@@ -28,7 +42,9 @@ function initial_condition_briowu_shock_tube(x, t, equations::IdealGlmMhdEquatio
     B3 = 0
     return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3), equations)
 end
-initial_condition = initial_condition_briowu_shock_tube
+# Note calling the constructor of the struct: `InitialConditionBrioWuShockTube()`
+# instead of `initial_condition_briowu_shock_tube` !
+initial_condition = InitialConditionBrioWuShockTube()
 
 boundary_conditions = BoundaryConditionDirichlet(initial_condition)
 
diff --git a/examples/tree_1d_dgsem/elixir_mhd_ryujones_shock_tube.jl b/examples/tree_1d_dgsem/elixir_mhd_ryujones_shock_tube.jl
index ce55973fee5..a05208b8a61 100644
--- a/examples/tree_1d_dgsem/elixir_mhd_ryujones_shock_tube.jl
+++ b/examples/tree_1d_dgsem/elixir_mhd_ryujones_shock_tube.jl
@@ -6,8 +6,21 @@ using Trixi
 gamma = 5 / 3
 equations = IdealGlmMhdEquations1D(gamma)
 
+# Specify the initial condition as a discontinuous initial condition (see docstring of 
+# `DiscontinuousFunction` for more information) which comes with a specialized 
+# initialization routine suited for Riemann problems.
+# In short, if a discontinuity is right at an interface, the boundary nodes (which are at the same location)
+# on that interface will be initialized with the left and right state of the discontinuity, i.e., 
+#                         { u_1, if element = left element and x_{element}^{(n)} = x_jump
+# u(x_jump, t, element) = {
+#                         { u_2, if element = right element and x_{element}^{(1)} = x_jump
+# This is realized by shifting the outer DG nodes inwards, i.e., on reference element
+# the outer nodes at `[-1, 1]` are shifted inwards to `[-1 + ε, 1 - ε]` with machine precision `ε`.
+struct InitialConditionRyuJonesShockTube <: DiscontinuousFunction end
+
 """
-    initial_condition_ryujones_shock_tube(x, t, equations::IdealGlmMhdEquations1D)
+    (initial_condition_ryujones_shock_tube::InitialConditionRyuJonesShockTube)(x, t,
+                                                                               equations::IdealGlmMhdEquations1D)
 
 Ryu and Jones shock tube test case for one dimensional ideal MHD equations. Contains
 fast shocks, slow shocks, and rational discontinuities that propagate on either side
@@ -20,7 +33,8 @@ present in the one dimensional MHD equations. It is the second test from Section
 !!! note
     This paper has a typo in the initial conditions. Their variable `E` should be `p`.
 """
-function initial_condition_ryujones_shock_tube(x, t, equations::IdealGlmMhdEquations1D)
+function (initial_condition_ryujones_shock_tube::InitialConditionRyuJonesShockTube)(x, t,
+                                                                                    equations::IdealGlmMhdEquations1D)
     # domain must be set to [0, 1], γ = 5/3, final time = 0.2
     RealT = eltype(x)
     rho = x[1] <= 0.5f0 ? RealT(1.08) : one(RealT)
@@ -35,7 +49,9 @@ function initial_condition_ryujones_shock_tube(x, t, equations::IdealGlmMhdEquat
 
     return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3), equations)
 end
-initial_condition = initial_condition_ryujones_shock_tube
+# Note calling the constructor of the struct: `InitialConditionRyuJonesShockTube()` instead of
+# `initial_condition_ryujones_shock_tube` !
+initial_condition = InitialConditionRyuJonesShockTube()
 
 boundary_conditions = BoundaryConditionDirichlet(initial_condition)
 
diff --git a/examples/tree_1d_dgsem/elixir_mhdmulti_briowu_shock_tube.jl b/examples/tree_1d_dgsem/elixir_mhdmulti_briowu_shock_tube.jl
index 19b5274df9c..3fe3e4decd0 100644
--- a/examples/tree_1d_dgsem/elixir_mhdmulti_briowu_shock_tube.jl
+++ b/examples/tree_1d_dgsem/elixir_mhdmulti_briowu_shock_tube.jl
@@ -6,8 +6,20 @@ using Trixi
 equations = IdealGlmMhdMulticomponentEquations1D(gammas = (2.0, 2.0),
                                                  gas_constants = (2.0, 2.0))
 
+# Specify the initial condition as a discontinuous initial condition (see docstring of 
+# `DiscontinuousFunction` for more information) which comes with a specialized 
+# initialization routine suited for Riemann problems.
+# In short, if a discontinuity is right at an interface, the boundary nodes (which are at the same location)
+# on that interface will be initialized with the left and right state of the discontinuity, i.e., 
+#                         { u_1, if element = left element and x_{element}^{(n)} = x_jump
+# u(x_jump, t, element) = {
+#                         { u_2, if element = right element and x_{element}^{(1)} = x_jump
+# This is realized by shifting the outer DG nodes inwards, i.e., on reference element
+# the outer nodes at `[-1, 1]` are shifted inwards to `[-1 + ε, 1 - ε]` with machine precision `ε`.
+struct InitialConditionBrioWuShockTube <: DiscontinuousFunction end
 """
-    initial_condition_briowu_shock_tube(x, t, equations::IdealGlmMhdMulticomponentEquations1D)
+    (initial_condition_briowu_shock_tube::InitialConditionBrioWuShockTube)(x, t,
+                                                                           equations::IdealGlmMhdMulticomponentEquations1D)
 
 Compound shock tube test case for one dimensional ideal MHD equations. It is basically an
 MHD extension of the Sod shock tube. Taken from Section V of the article
@@ -15,8 +27,8 @@ MHD extension of the Sod shock tube. Taken from Section V of the article
   An Upwind Differencing Scheme for the Equations of Ideal Magnetohydrodynamics
   [DOI: 10.1016/0021-9991(88)90120-9](https://doi.org/10.1016/0021-9991(88)90120-9)
 """
-function initial_condition_briowu_shock_tube(x, t,
-                                             equations::IdealGlmMhdMulticomponentEquations1D)
+function (initial_condition_briowu_shock_tube::InitialConditionBrioWuShockTube)(x, t,
+                                                                                equations::IdealGlmMhdMulticomponentEquations1D)
     # domain must be set to [0, 1], γ = 2, final time = 0.12
     RealT = eltype(x)
     if x[1] < 0.5f0
@@ -43,7 +55,9 @@ function initial_condition_briowu_shock_tube(x, t,
     prim_other = SVector(v1, v2, v3, p, B1, B2, B3)
     return prim2cons(vcat(prim_other, prim_rho), equations)
 end
-initial_condition = initial_condition_briowu_shock_tube
+# Note calling the constructor of the struct: `InitialConditionBrioWuShockTube()`
+# instead of `initial_condition_briowu_shock_tube` !
+initial_condition = InitialConditionBrioWuShockTube()
 
 boundary_conditions = BoundaryConditionDirichlet(initial_condition)
 
diff --git a/examples/tree_2d_dgsem/elixir_euler_riemannproblem_quadrants_amr.jl b/examples/tree_2d_dgsem/elixir_euler_riemannproblem_quadrants_amr.jl
new file mode 100644
index 00000000000..48eff65c70e
--- /dev/null
+++ b/examples/tree_2d_dgsem/elixir_euler_riemannproblem_quadrants_amr.jl
@@ -0,0 +1,128 @@
+using OrdinaryDiffEqSSPRK
+using Trixi
+
+###############################################################################
+## Semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+# Variant of the 4-quadrant Riemann problem considered in 
+# - Carsten W. Schulz-Rinne:
+#   Classification of the Riemann Problem for Two-Dimensional Gas Dynamics
+#   https://doi.org/10.1137/0524006
+# and 
+# - Carsten W. Schulz-Rinne, James P. Collins, and Harland M. Glaz
+#   Numerical Solution of the Riemann Problem for Two-Dimensional Gas Dynamics
+#   https://doi.org/10.1137/0914082
+
+# Specify the initial condition as a discontinuous initial condition (see docstring of 
+# `DiscontinuousFunction` for more information) which comes with a specialized 
+# initialization routine suited for Riemann problems.
+# In short, if a discontinuity is right at an interface, the boundary nodes (which are at the same location)
+# on that interface will be initialized with the left and right state of the discontinuity, i.e., 
+#                         { u_1, if element = left element and x_{element}^{(n)} = x_jump
+# u(x_jump, t, element) = {
+#                         { u_2, if element = right element and x_{element}^{(1)} = x_jump
+# This is realized by shifting the outer DG nodes inwards, i.e., on reference element
+# the outer nodes at `[-1, 1]` are shifted inwards to `[-1 + ε, 1 - ε]` with machine precision `ε`.
+struct InitialConditionRP <: DiscontinuousFunction end
+
+function (initial_condition_rp::InitialConditionRP)(x_, t,
+                                                    equations::CompressibleEulerEquations2D)
+    x, y = x_[1], x_[2]
+
+    if x >= 0.5 && y >= 0.5
+        rho, v1, v2, p = (0.5313, 0.0, 0.0, 0.4)
+    elseif x < 0.5 && y >= 0.5
+        rho, v1, v2, p = (1.0, 0.7276, 0.0, 1.0)
+    elseif x < 0.5 && y < 0.5
+        rho, v1, v2, p = (0.8, 0.0, 0.0, 1.0)
+    elseif x >= 0.5 && y < 0.5
+        rho, v1, v2, p = (1.0, 0.0, 0.7276, 1.0)
+    end
+
+    prim = SVector(rho, v1, v2, p)
+    return prim2cons(prim, equations)
+end
+# Note calling the constructor of the struct: `InitialConditionRP()` instead of `initial_condition_rp` !
+initial_condition = InitialConditionRP()
+
+# Extend domain by specifying free outflow
+function boundary_condition_outflow(u_inner,
+                                    orientation::Integer, normal_direction,
+                                    x, t,
+                                    surface_flux_function,
+                                    equations::CompressibleEulerEquations2D)
+    flux = Trixi.flux(u_inner, orientation, equations)
+
+    return flux
+end
+
+boundary_conditions = (x_neg = boundary_condition_outflow,
+                       x_pos = boundary_condition_outflow,
+                       y_neg = boundary_condition_outflow,
+                       y_pos = boundary_condition_outflow)
+
+coordinates_min = (0.0, 0.0) # minimum coordinates (min(x), min(y))
+coordinates_max = (1.0, 1.0) # maximum coordinates (max(x), max(y))
+
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 4,
+                n_cells_max = 100_000,
+                periodicity = false)
+
+surface_flux = flux_hllc
+volume_flux = flux_ranocha
+
+polydeg = 3
+basis = LobattoLegendreBasis(polydeg)
+shock_indicator = IndicatorHennemannGassner(equations, basis,
+                                            alpha_max = 0.5,
+                                            alpha_min = 0.001,
+                                            alpha_smooth = true,
+                                            variable = Trixi.density)
+
+volume_integral = VolumeIntegralShockCapturingHG(shock_indicator;
+                                                 volume_flux_dg = volume_flux,
+                                                 volume_flux_fv = surface_flux)
+
+solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
+               volume_integral = volume_integral)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver;
+                                    boundary_conditions = boundary_conditions)
+
+###############################################################################
+## ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.25)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+stepsize_callback = StepsizeCallback(cfl = 1.8)
+
+amr_indicator = IndicatorLöhner(semi, variable = Trixi.density)
+amr_controller = ControllerThreeLevel(semi, amr_indicator,
+                                      base_level = 2,
+                                      med_level = 5, med_threshold = 0.02,
+                                      max_level = 8, max_threshold = 0.04)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval = 10,
+                           adapt_initial_condition = true)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        amr_callback,
+                        stepsize_callback)
+
+###############################################################################
+## Run the simulation
+
+sol = solve(ode, SSPRK54();
+            dt = 1.0, save_everystep = false, callback = callbacks);
diff --git a/examples/tree_3d_dgsem/elixir_euler_sedov_blast_wave.jl b/examples/tree_3d_dgsem/elixir_euler_sedov_blast_wave.jl
index 1413ba05966..3291ac93f3a 100644
--- a/examples/tree_3d_dgsem/elixir_euler_sedov_blast_wave.jl
+++ b/examples/tree_3d_dgsem/elixir_euler_sedov_blast_wave.jl
@@ -6,8 +6,20 @@ using Trixi
 
 equations = CompressibleEulerEquations3D(1.4)
 
+# Specify the initial condition as a discontinuous initial condition (see docstring of 
+# `DiscontinuousFunction` for more information) which comes with a specialized 
+# initialization routine suited for the Riemann problems.
+# In short, if a discontinuity is right at an interface, the boundary nodes (which are at the same location)
+# on that interface will be initialized with the left and right state of the discontinuity, i.e., 
+#                         { u_1, if element = left element and x_{element}^{(n)} = x_jump
+# u(x_jump, t, element) = {
+#                         { u_2, if element = right element and x_{element}^{(1)} = x_jump
+# This is realized by shifting the outer DG nodes inwards, i.e., on reference element
+# the outer nodes are at `[-1, 1]` are shifted to `[-1 + ε, 1 - ε]` with machine precision `ε`.
+struct InitialConditionSedovSelfGravity <: DiscontinuousFunction end
+
 """
-    initial_condition_sedov_self_gravity(x, t, equations::CompressibleEulerEquations3D)
+    (initial_condition_sedov_self_gravity::InitialConditionSedovSelfGravity)(x, t, equations::CompressibleEulerEquations3D)
 
 Adaptation of the Sedov blast wave with self-gravity taken from
 - Michael Schlottke-Lakemper, Andrew R. Winters, Hendrik Ranocha, Gregor J. Gassner (2020)
@@ -17,7 +29,8 @@ based on
 - https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000
 Should be used together with [`boundary_condition_sedov_self_gravity`](@ref).
 """
-function initial_condition_sedov_self_gravity(x, t, equations::CompressibleEulerEquations3D)
+function (initial_condition_sedov_self_gravity::InitialConditionSedovSelfGravity)(x, t,
+                                                                                  equations::CompressibleEulerEquations3D)
     # Calculate radius as distance from origin
     r = sqrt(x[1]^2 + x[2]^2 + x[3]^2)
 
@@ -47,7 +60,9 @@ function initial_condition_sedov_self_gravity(x, t, equations::CompressibleEuler
 
     return prim2cons(SVector(rho, v1, v2, v3, p), equations)
 end
-initial_condition = initial_condition_sedov_self_gravity
+# Note calling the constructor of the struct: `InitialConditionSedovSelfGravity()` instead of 
+# `initial_condition_sedov_self_gravity` !
+initial_condition = InitialConditionSedovSelfGravity()
 
 """
     boundary_condition_sedov_self_gravity(u_inner, orientation, direction, x, t,
diff --git a/src/Trixi.jl b/src/Trixi.jl
index 75d2d739f59..81a54dc4783 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -213,7 +213,8 @@ export splitting_steger_warming, splitting_vanleer_haenel,
 export initial_condition_constant,
        initial_condition_gauss,
        initial_condition_density_wave,
-       initial_condition_weak_blast_wave
+       initial_condition_weak_blast_wave,
+       DiscontinuousFunction
 
 export boundary_condition_do_nothing,
        boundary_condition_periodic,
diff --git a/src/basic_types.jl b/src/basic_types.jl
index 1aaf6d6d464..d0d63845ddd 100644
--- a/src/basic_types.jl
+++ b/src/basic_types.jl
@@ -121,4 +121,32 @@ const boundary_condition_do_nothing = BoundaryConditionDoNothing()
 function Base.show(io::IO, ::BoundaryConditionDoNothing)
     print(io, "boundary_condition_do_nothing")
 end
+
+@doc raw"""
+    DiscontinuousFunction
+
+An abstract type for convenient setting of discontinuous functions,
+mainly initial conditions for Riemann problems.
+For a struct that inherits from this type, i.e., 
+`struct initial_condition <: DiscontinuousFunction end`
+a specialized initialization is performed.
+
+The outer nodes (i.e., ±1 on the DG reference element) are shifted inwards
+by the smallest amount possible, i.e., `[-1 + ϵ, +1 - ϵ]`, when the initial condition
+is evaluated at the nodes.
+This avoids steep gradients in elements if a discontinuity is right at a cell boundary,
+i.e., if the jump location `x_jump` is at the position of an interface which is shared by
+the nodes ``x_{e-1}^{(n)} = x_{e}^{(1)}`` where ``e`` indexes the elements and ``n`` the DG nodes.
+
+In particular, this results in the typically desired behaviour for discontinuous initial conditions of the form
+```math
+    u(x, t) = 
+    \begin{cases}
+        u_1, & \text{if } x \leq x_{\text{jump}} \\
+        u_2, & \text{if } x > x_{\text{jump}}
+    \end{cases}
+```
+where ``x_\text{jump}`` is the location of the discontinuity.
+"""
+abstract type DiscontinuousFunction end
 end # @muladd
diff --git a/src/solvers/dg.jl b/src/solvers/dg.jl
index ad211b3c003..8adb5137c85 100644
--- a/src/solvers/dg.jl
+++ b/src/solvers/dg.jl
@@ -732,45 +732,276 @@ end
                  nelements(dg, cache)))
 end
 
-function compute_coefficients!(u, func, t, mesh::AbstractMesh{1}, equations, dg::DG,
-                               cache)
+# Standard function for computing coefficients of `func` (usually the initial condition)
+function compute_coefficients!(u, func, t,
+                               mesh::AbstractMesh{1}, equations, dg::DG, cache)
     @threaded for element in eachelement(dg, cache)
         for i in eachnode(dg)
-            x_node = get_node_coords(cache.elements.node_coordinates, equations, dg, i,
-                                     element)
-            # Changing the node positions passed to the initial condition by the minimum
-            # amount possible with the current type of floating point numbers allows setting
-            # discontinuous initial data in a simple way. In particular, a check like `if x < x_jump`
-            # works if the jump location `x_jump` is at the position of an interface.
+            x_node = get_node_coords(cache.elements.node_coordinates, equations, dg,
+                                     i, element)
+            u_node = func(x_node, t, equations)
+            set_node_vars!(u, u_node, equations, dg, i, element)
+        end
+    end
+end
+# Specialized function for computing coefficients of `func`, here a discontinuous initial condition
+# Shift the outer (i.e., ±1 on the reference element) nodes passed to the 
+# `func` (usually the initial condition) inwards by the smallest amount possible,
+# i.e., [-1 + ϵ, +1 - ϵ].
+# This avoids steep gradients in elements if a discontinuity is right at a cell boundary,
+# i.e., if the jump location `x_jump` is at the position of an interface which is shared by 
+# the nodes x_{e-1}^{(i)} = x_{e}^{(1)}.
+# In particular, this results in the typically desired behaviour for 
+# initial conditions of the form
+#           { u_1, if x <= x_jump
+# u(x, t) = {
+#           { u_2, if x > x_jump
+function compute_coefficients!(u, func::DiscontinuousFunction, t,
+                               mesh::AbstractMesh{1}, equations, dg::DG, cache)
+    @threaded for element in eachelement(dg, cache)
+        for i in eachnode(dg)
+            x_node = get_node_coords(cache.elements.node_coordinates, equations, dg,
+                                     i, element)
+            # For explanation of the implemented logic see comment above
             if i == 1
                 x_node = SVector(nextfloat(x_node[1]))
             elseif i == nnodes(dg)
                 x_node = SVector(prevfloat(x_node[1]))
             end
+
             u_node = func(x_node, t, equations)
             set_node_vars!(u, u_node, equations, dg, i, element)
         end
     end
 end
 
-function compute_coefficients!(u, func, t, mesh::AbstractMesh{2}, equations, dg::DG,
-                               cache)
+# Standard function for computing coefficients of `func` (usually the initial condition)
+function compute_coefficients!(u, func, t,
+                               mesh::AbstractMesh{2}, equations, dg::DG, cache)
+    @threaded for element in eachelement(dg, cache)
+        for j in eachnode(dg), i in eachnode(dg)
+            x_node = get_node_coords(cache.elements.node_coordinates, equations, dg,
+                                     i, j, element)
+            u_node = func(x_node, t, equations)
+            set_node_vars!(u, u_node, equations, dg, i, j, element)
+        end
+    end
+end
+# Specialized function for computing coefficients of `func`, here a discontinuous function
+# `TreeMesh` allows for an easy implementation
+function compute_coefficients!(u, func::DiscontinuousFunction, t,
+                               mesh::TreeMesh{2}, equations, dg::DG, cache)
+    @threaded for element in eachelement(dg, cache)
+        for j in eachnode(dg), i in eachnode(dg)
+            x_node = get_node_coords(cache.elements.node_coordinates, equations, dg,
+                                     i, j, element)
+            # For explanation of the implemented logic see comment above for the 1D case, i.e., 
+            # `compute_coefficients!` for `func::DiscontinuousFunction`, `mesh::AbstractMesh{1}`
+            # or the docstring of `DiscontinuousFunction`
+            if i == 1
+                x_node = SVector(nextfloat(x_node[1]), x_node[2])
+            elseif i == nnodes(dg)
+                x_node = SVector(prevfloat(x_node[1]), x_node[2])
+            end
+            if j == 1
+                x_node = SVector(x_node[1], nextfloat(x_node[2]))
+            elseif j == nnodes(dg)
+                x_node = SVector(x_node[1], prevfloat(x_node[2]))
+            end
+
+            u_node = func(x_node, t, equations)
+            set_node_vars!(u, u_node, equations, dg, i, j, element)
+        end
+    end
+end
+# Mesh-type general function for dragging the outer nodes inwards for 
+# discontinuous functions/initial conditions
+function compute_coefficients!(u, func::DiscontinuousFunction, t,
+                               mesh::AbstractMesh{2}, equations, dg::DG, cache)
     @threaded for element in eachelement(dg, cache)
+        # Helper variables to avoid querying `get_node_coords` too often
+        x1_min_at_node1 = true
+        x2_min_at_node1 = true
         for j in eachnode(dg), i in eachnode(dg)
-            x_node = get_node_coords(cache.elements.node_coordinates, equations, dg, i,
-                                     j, element)
+            x_node = get_node_coords(cache.elements.node_coordinates, equations, dg,
+                                     i, j, element)
+            # For explanation of the implemented logic see comment above for the 1D case, i.e., 
+            # `compute_coefficients!` for `func::DiscontinuousFunction`, `mesh::AbstractMesh{1}`
+            # or the docstring of `DiscontinuousFunction`
+            if i == 1
+                x_node_n = get_node_coords(cache.elements.node_coordinates, equations,
+                                           dg, nnodes(dg), j, element)
+
+                # For non-`TreeMesh`es the node orientation within an element is not necessarily
+                # consistent with the global coordinate system, thus we need to check what node
+                # lies on the "left" and which on the "right" side.
+                # NOTE: This still assumes that the first node index somewhat resembles the first coordinate direction!
+                if x_node[1] < x_node_n[1]
+                    x_node = SVector(nextfloat(x_node[1]), x_node[2])
+                else
+                    x1_min_at_node1 = false
+                    x_node = SVector(prevfloat(x_node[1]), x_node[2])
+                end
+            elseif i == nnodes(dg)
+                if x1_min_at_node1
+                    x_node = SVector(prevfloat(x_node[1]), x_node[2])
+                else
+                    x_node = SVector(nextfloat(x_node[1]), x_node[2])
+                end
+            end
+            if j == 1
+                x_node_n = get_node_coords(cache.elements.node_coordinates, equations,
+                                           dg, i, nnodes(dg), element)
+
+                # For non-`TreeMesh`es the node orientation within an element is not necessarily
+                # consistent with the global coordinate system, thus we need to check what node
+                # lies on the "bottom" and which on the "top" side.
+                # NOTE: This still assumes that the second node index somewhat resembles the second coordinate direction!
+                if x_node[2] < x_node_n[2]
+                    x_node = SVector(x_node[1], nextfloat(x_node[2]))
+                else
+                    x2_min_at_node1 = false
+                    x_node = SVector(x_node[1], prevfloat(x_node[2]))
+                end
+            elseif j == nnodes(dg)
+                if x2_min_at_node1
+                    x_node = SVector(x_node[1], prevfloat(x_node[2]))
+                else
+                    x_node = SVector(x_node[1], nextfloat(x_node[2]))
+                end
+            end
+
             u_node = func(x_node, t, equations)
             set_node_vars!(u, u_node, equations, dg, i, j, element)
         end
     end
 end
 
-function compute_coefficients!(u, func, t, mesh::AbstractMesh{3}, equations, dg::DG,
-                               cache)
+# Standard function for computing coefficients of `func` (usually the initial condition)
+function compute_coefficients!(u, func, t,
+                               mesh::AbstractMesh{3}, equations, dg::DG, cache)
+    @threaded for element in eachelement(dg, cache)
+        for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
+            x_node = get_node_coords(cache.elements.node_coordinates, equations, dg,
+                                     i, j, k, element)
+            u_node = func(x_node, t, equations)
+            set_node_vars!(u, u_node, equations, dg, i, j, k, element)
+        end
+    end
+end
+# Specialized function for computing coefficients of `func`, here a discontinuous function
+# `TreeMesh` allows for an easy implementation
+function compute_coefficients!(u, func::DiscontinuousFunction, t,
+                               mesh::TreeMesh{3}, equations, dg::DG, cache)
+    @threaded for element in eachelement(dg, cache)
+        for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
+            x_node = get_node_coords(cache.elements.node_coordinates, equations, dg,
+                                     i, j, k, element)
+            # For explanation of the implemented logic see comment above for the 1D case, i.e., 
+            # `compute_coefficients!` for `func::DiscontinuousFunction`, `mesh::AbstractMesh{1}`
+            # or the docstring of `DiscontinuousFunction`
+            if i == 1
+                x_node = SVector(nextfloat(x_node[1]), x_node[2], x_node[3])
+            elseif i == nnodes(dg)
+                x_node = SVector(prevfloat(x_node[1]), x_node[2], x_node[3])
+            end
+            if j == 1
+                x_node = SVector(x_node[1], nextfloat(x_node[2]), x_node[3])
+            elseif j == nnodes(dg)
+                x_node = SVector(x_node[1], prevfloat(x_node[2]), x_node[3])
+            end
+            if k == 1
+                x_node = SVector(x_node[1], x_node[2], nextfloat(x_node[3]))
+            elseif k == nnodes(dg)
+                x_node = SVector(x_node[1], x_node[2], prevfloat(x_node[3]))
+            end
+
+            u_node = func(x_node, t, equations)
+            set_node_vars!(u, u_node, equations, dg, i, j, k, element)
+        end
+    end
+end
+# Mesh-type general function for dragging the outer nodes inwards for 
+# discontinuous functions/initial conditions
+function compute_coefficients!(u, func::DiscontinuousFunction, t,
+                               mesh::AbstractMesh{3}, equations, dg::DG, cache)
     @threaded for element in eachelement(dg, cache)
+        # Helper variables to avoid querying `get_node_coords` too often
+        x1_min_at_node1 = true
+        x2_min_at_node1 = true
+        x3_min_at_node1 = true
         for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
-            x_node = get_node_coords(cache.elements.node_coordinates, equations, dg, i,
-                                     j, k, element)
+            x_node = get_node_coords(cache.elements.node_coordinates, equations, dg,
+                                     i, j, k, element)
+            # For explanation of the implemented logic see comment above for the 1D case, i.e., 
+            # `compute_coefficients!` for `func::DiscontinuousFunction`, `mesh::AbstractMesh{1}`
+            # or the docstring of `DiscontinuousFunction`
+            if i == 1
+                x_node_n = get_node_coords(cache.elements.node_coordinates, equations,
+                                           dg, nnodes(dg), j, k, element)
+
+                # For non-`TreeMesh`es the node orientation within an element is not necessarily
+                # consistent with the global coordinate system, thus we need to check what node
+                # lies on the "left" and which on the "right" side.
+                # NOTE: This still assumes that the first node index somewhat resembles the first coordinate direction!
+                if x_node[1] < x_node_n[1]
+                    x_node = SVector(nextfloat(x_node[1]), x_node[2], x_node[3])
+                else
+                    x1_min_at_node1 = false
+                    x_node = SVector(prevfloat(x_node[1]), x_node[2], x_node[3])
+                end
+            elseif i == nnodes(dg)
+                if x1_min_at_node1
+                    x_node = SVector(prevfloat(x_node[1]), x_node[2], x_node[3])
+                else
+                    x_node = SVector(nextfloat(x_node[1]), x_node[2], x_node[3])
+                end
+            end
+
+            if j == 1
+                x_node_n = get_node_coords(cache.elements.node_coordinates, equations,
+                                           dg, i, nnodes(dg), k, element)
+
+                # For non-`TreeMesh`es the node orientation within an element is not necessarily
+                # consistent with the global coordinate system, thus we need to check what node
+                # lies on the "bottom" and which on the "top" side.
+                # NOTE: This still assumes that the second node index somewhat resembles the second coordinate direction!
+                if x_node[2] < x_node_n[2]
+                    x_node = SVector(x_node[1], nextfloat(x_node[2]), x_node[3])
+                else
+                    x2_min_at_node1 = false
+                    x_node = SVector(x_node[1], prevfloat(x_node[2]), x_node[3])
+                end
+            elseif j == nnodes(dg)
+                if x2_min_at_node1
+                    x_node = SVector(x_node[1], prevfloat(x_node[2]), x_node[3])
+                else
+                    x_node = SVector(x_node[1], nextfloat(x_node[2]), x_node[3])
+                end
+            end
+
+            if k == 1
+                x_node_n = get_node_coords(cache.elements.node_coordinates, equations,
+                                           dg, i, j, nnodes(dg), element)
+
+                # For non-`TreeMesh`es the node orientation within an element is not necessarily
+                # consistent with the global coordinate system, thus we need to check what node
+                # lies on the "back" and which on the "front" side.
+                # NOTE: This still assumes that the third node index somewhat resembles the third coordinate direction!
+                if x_node[3] < x_node_n[3]
+                    x_node = SVector(x_node[1], x_node[2], nextfloat(x_node[3]))
+                else
+                    x3_min_at_node1 = false
+                    x_node = SVector(x_node[1], x_node[2], prevfloat(x_node[3]))
+                end
+            elseif k == nnodes(dg)
+                if x3_min_at_node1
+                    x_node = SVector(x_node[1], x_node[2], prevfloat(x_node[3]))
+                else
+                    x_node = SVector(x_node[1], x_node[2], nextfloat(x_node[3]))
+                end
+            end
+
             u_node = func(x_node, t, equations)
             set_node_vars!(u, u_node, equations, dg, i, j, k, element)
         end
diff --git a/src/solvers/dgsem_structured/dg_2d.jl b/src/solvers/dgsem_structured/dg_2d.jl
index 54349c6bf83..15bfbd612b6 100644
--- a/src/solvers/dgsem_structured/dg_2d.jl
+++ b/src/solvers/dgsem_structured/dg_2d.jl
@@ -437,7 +437,7 @@ end
     left_direction = right_direction - 1
 
     for i in eachnode(dg)
-        if orientation == 1
+        if orientation == 1 # (x-direction)
             u_ll = get_node_vars(u, equations, dg, nnodes(dg), i, left_element)
             u_rr = get_node_vars(u, equations, dg, 1, i, right_element)
 
@@ -451,7 +451,7 @@ end
             normal_direction = sign_jacobian *
                                get_contravariant_vector(1, contravariant_vectors,
                                                         1, i, right_element)
-        else # orientation == 2
+        else # orientation == 2 (y-direction)
             u_ll = get_node_vars(u, equations, dg, i, nnodes(dg), left_element)
             u_rr = get_node_vars(u, equations, dg, i, 1, right_element)
 
@@ -495,7 +495,7 @@ end
     left_direction = right_direction - 1
 
     for i in eachnode(dg)
-        if orientation == 1
+        if orientation == 1 # (x-direction)
             u_ll = get_node_vars(u, equations, dg, nnodes(dg), i, left_element)
             u_rr = get_node_vars(u, equations, dg, 1, i, right_element)
 
@@ -509,7 +509,7 @@ end
             normal_direction = sign_jacobian *
                                get_contravariant_vector(1, contravariant_vectors,
                                                         1, i, right_element)
-        else # orientation == 2
+        else # orientation == 2 (y-direction)
             u_ll = get_node_vars(u, equations, dg, i, nnodes(dg), left_element)
             u_rr = get_node_vars(u, equations, dg, i, 1, right_element)
 
diff --git a/src/solvers/dgsem_structured/dg_3d.jl b/src/solvers/dgsem_structured/dg_3d.jl
index 983a62d4144..994d78376ec 100644
--- a/src/solvers/dgsem_structured/dg_3d.jl
+++ b/src/solvers/dgsem_structured/dg_3d.jl
@@ -552,7 +552,7 @@ end
     left_direction = right_direction - 1
 
     for j in eachnode(dg), i in eachnode(dg)
-        if orientation == 1
+        if orientation == 1 # (x-direction)
             u_ll = get_node_vars(u, equations, dg, nnodes(dg), i, j, left_element)
             u_rr = get_node_vars(u, equations, dg, 1, i, j, right_element)
 
@@ -566,7 +566,7 @@ end
             normal_direction = sign_jacobian *
                                get_contravariant_vector(1, contravariant_vectors,
                                                         1, i, j, right_element)
-        elseif orientation == 2
+        elseif orientation == 2 # (y-direction)
             u_ll = get_node_vars(u, equations, dg, i, nnodes(dg), j, left_element)
             u_rr = get_node_vars(u, equations, dg, i, 1, j, right_element)
 
@@ -577,7 +577,7 @@ end
             normal_direction = sign_jacobian *
                                get_contravariant_vector(2, contravariant_vectors,
                                                         i, 1, j, right_element)
-        else # orientation == 3
+        else # orientation == 3 (z-direction)
             u_ll = get_node_vars(u, equations, dg, i, j, nnodes(dg), left_element)
             u_rr = get_node_vars(u, equations, dg, i, j, 1, right_element)
 
@@ -620,7 +620,7 @@ end
     left_direction = right_direction - 1
 
     for j in eachnode(dg), i in eachnode(dg)
-        if orientation == 1
+        if orientation == 1 # (x-direction)
             u_ll = get_node_vars(u, equations, dg, nnodes(dg), i, j, left_element)
             u_rr = get_node_vars(u, equations, dg, 1, i, j, right_element)
 
@@ -634,7 +634,7 @@ end
             normal_direction = sign_jacobian *
                                get_contravariant_vector(1, contravariant_vectors,
                                                         1, i, j, right_element)
-        elseif orientation == 2
+        elseif orientation == 2 # (y-direction)
             u_ll = get_node_vars(u, equations, dg, i, nnodes(dg), j, left_element)
             u_rr = get_node_vars(u, equations, dg, i, 1, j, right_element)
 
@@ -645,7 +645,7 @@ end
             normal_direction = sign_jacobian *
                                get_contravariant_vector(2, contravariant_vectors,
                                                         i, 1, j, right_element)
-        else # orientation == 3
+        else # orientation == 3 (z-direction)
             u_ll = get_node_vars(u, equations, dg, i, j, nnodes(dg), left_element)
             u_rr = get_node_vars(u, equations, dg, i, j, 1, right_element)
 
diff --git a/test/test_parabolic_1d.jl b/test/test_parabolic_1d.jl
index 8969b06eb13..e902f2558b8 100644
--- a/test/test_parabolic_1d.jl
+++ b/test/test_parabolic_1d.jl
@@ -17,7 +17,7 @@ isdir(outdir) && rm(outdir, recursive = true)
                                  "elixir_advection_diffusion.jl"),
                         initial_refinement_level=4, tspan=(0.0, 0.4), polydeg=3,
                         l2=[8.40483031802723e-6],
-                        linf=[2.8990878868540015e-5])
+                        linf=[2.899087786395471e-5])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     let
@@ -62,7 +62,7 @@ end
     @test_trixi_include(joinpath(examples_dir(), "tree_1d_dgsem",
                                  "elixir_advection_diffusion_restart.jl"),
                         l2=[1.0679933947301556e-5],
-                        linf=[3.910500545667439e-5])
+                        linf=[3.910500469289646e-5])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     let
@@ -165,9 +165,9 @@ end
                             0.001496626616191212
                         ],
                         linf=[
-                            0.0029963751636357117,
-                            0.0028639041695096433,
-                            0.012691132694550689
+                            0.0029963750445691772,
+                            0.002863904416491592,
+                            0.01269113341266781
                         ],
                         atol=1e-10)
     # Ensure that we do not have excessive memory allocations
@@ -216,13 +216,13 @@ end
                                                                                 Prandtl = prandtl_number()),
                         l2=[
                             2.5278845598681636e-5,
-                            2.5540145802666872e-5,
-                            0.0001211867535580826
+                            2.5538884774356816e-5,
+                            0.00012118293863523099
                         ],
                         linf=[
                             0.0001466387202588848,
-                            0.00019422419092429135,
-                            0.0009556449835592673
+                            0.0001942257659360519,
+                            0.0009556393581480194
                         ],
                         atol=1e-9)
     # Ensure that we do not have excessive memory allocations
@@ -250,7 +250,7 @@ end
                         linf=[
                             0.00011850494672183132,
                             0.00018987676556476442,
-                            0.0009597423024825247
+                            0.0009597509690966177
                         ],
                         atol=1e-9)
     # Ensure that we do not have excessive memory allocations
diff --git a/test/test_structured_1d.jl b/test/test_structured_1d.jl
index 40a238b262a..0619e9733d9 100644
--- a/test/test_structured_1d.jl
+++ b/test/test_structured_1d.jl
@@ -105,8 +105,8 @@ end
                         save_solution=SaveSolutionCallback(dt = 0.1 + 1.0e-8), # Adding a small epsilon to avoid floating-point precision issues
                         callbacks=CallbackSet(summary_callback, save_solution,
                                               analysis_callback, alive_callback),
-                        l2=[5.726144824784944e-7],
-                        linf=[3.43073006914274e-6])
+                        l2=[5.728387746154631e-7],
+                        linf=[3.434003242208661e-6])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     let
diff --git a/test/test_tree_1d.jl b/test/test_tree_1d.jl
index d8658062b90..a1d57f298d2 100644
--- a/test/test_tree_1d.jl
+++ b/test/test_tree_1d.jl
@@ -294,7 +294,7 @@ end
     sol = solve(ode, Tsit5(), abstol = 1.0e-6, reltol = 1.0e-6;
                 ode_default_options()..., callback = callbacks)
 
-    @test analysis_callback(sol).l2 ≈ [0.00029609575838969394, 5.5681704039507985e-6]
+    @test analysis_callback(sol).l2 ≈ [0.0002961027497380263, 5.573684084938363e-6]
 end
 
 # Clean up afterwards: delete Trixi.jl output directory
diff --git a/test/test_tree_1d_euler.jl b/test/test_tree_1d_euler.jl
index 56e76e80db3..20e39bbcdb6 100644
--- a/test/test_tree_1d_euler.jl
+++ b/test/test_tree_1d_euler.jl
@@ -76,7 +76,7 @@ end
                         linf=[
                             0.004090978306812376,
                             0.0004090978313582294,
-                            2.045489210189544e-5
+                            2.0454891760834926e-5
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
diff --git a/test/test_tree_1d_hypdiff.jl b/test/test_tree_1d_hypdiff.jl
index 562c940d44e..c35229b051c 100644
--- a/test/test_tree_1d_hypdiff.jl
+++ b/test/test_tree_1d_hypdiff.jl
@@ -13,7 +13,7 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem")
 @trixi_testset "elixir_hypdiff_nonperiodic.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_hypdiff_nonperiodic.jl"),
                         l2=[1.3655114953278825e-7, 1.0200345026471077e-6],
-                        linf=[7.173285075379177e-7, 4.507116828644797e-6])
+                        linf=[7.173285075379177e-7, 4.507116716734316e-6])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     let