diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 0c636ee8b0b..725cd348de9 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -81,11 +81,7 @@ jobs:
           - mpi
           - threaded
         include:
-          - version: '1.8'
-            os: ubuntu-latest
-            arch: x64
-            trixi_test: threaded_legacy
-          - version: '1.9'
+          - version: '1.10'
             os: ubuntu-latest
             arch: x64
             trixi_test: threaded_legacy
diff --git a/Project.toml b/Project.toml
index e0ce86c23ba..28b4562391c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -79,7 +79,7 @@ FillArrays = "0.13.2, 1"
 ForwardDiff = "0.10.24"
 HDF5 = "0.16.10, 0.17"
 IfElse = "0.1"
-LinearAlgebra = "1"
+LinearAlgebra = "1.8.5, 1.10"
 LinearMaps = "2.7, 3.0"
 LoopVectorization = "0.12.151"
 MPI = "0.20.6"
@@ -94,7 +94,7 @@ PrecompileTools = "1.1"
 Preferences = "1.3"
 Printf = "1"
 RecipesBase = "1.1"
-RecursiveArrayTools = "2.38.10"
+RecursiveArrayTools = "3.27.1"
 Reexport = "1.0"
 Requires = "1.1"
 SciMLBase = "1.90, 2"
@@ -115,7 +115,7 @@ TriplotBase = "0.1"
 TriplotRecipes = "0.1"
 TrixiBase = "0.1.3"
 UUIDs = "1.6"
-julia = "1.8"
+julia = "1.10"
 
 [extras]
 Convex = "f65535da-76fb-5f13-bab9-19810c17039a"
diff --git a/examples/dgmulti_2d/elixir_navierstokes_lid_driven_cavity.jl b/examples/dgmulti_2d/elixir_navierstokes_lid_driven_cavity.jl
index 6be384138e2..6ee0c3b76e1 100644
--- a/examples/dgmulti_2d/elixir_navierstokes_lid_driven_cavity.jl
+++ b/examples/dgmulti_2d/elixir_navierstokes_lid_driven_cavity.jl
@@ -5,10 +5,10 @@ using Trixi
 # semidiscretization of the ideal compressible Navier-Stokes equations
 
 prandtl_number() = 0.72
-mu = 0.001
+mu() = 0.001
 
 equations = CompressibleEulerEquations2D(1.4)
-equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu,
+equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu(),
                                                           Prandtl = prandtl_number())
 
 # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
diff --git a/examples/dgmulti_3d/elixir_navierstokes_taylor_green_vortex.jl b/examples/dgmulti_3d/elixir_navierstokes_taylor_green_vortex.jl
index 5b844fbe3d2..7b568a84a96 100644
--- a/examples/dgmulti_3d/elixir_navierstokes_taylor_green_vortex.jl
+++ b/examples/dgmulti_3d/elixir_navierstokes_taylor_green_vortex.jl
@@ -6,10 +6,10 @@ using Trixi
 # semidiscretization of the compressible Navier-Stokes equations
 
 prandtl_number() = 0.72
-mu = 6.25e-4 # equivalent to Re = 1600
+mu() = 6.25e-4 # equivalent to Re = 1600
 
 equations = CompressibleEulerEquations3D(1.4)
-equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu,
+equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu(),
                                                           Prandtl = prandtl_number())
 
 """
diff --git a/src/Trixi.jl b/src/Trixi.jl
index 7d557ddde38..1dcf4f7d11e 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -61,6 +61,7 @@ using OffsetArrays: OffsetArray, OffsetVector
 using P4est
 using T8code
 using RecipesBase: RecipesBase
+using RecursiveArrayTools: VectorOfArray
 using Requires: @require
 using Static: Static, One, True, False
 @reexport using StaticArrays: SVector
diff --git a/src/auxiliary/precompile.jl b/src/auxiliary/precompile.jl
index dfda8ece687..3fa5f9575e7 100644
--- a/src/auxiliary/precompile.jl
+++ b/src/auxiliary/precompile.jl
@@ -523,9 +523,9 @@ function _precompile_manual_()
                                                  typeof(Trixi.initialize_summary_callback),
                                                  typeof(SciMLBase.FINALIZE_DEFAULT)}
         @assert Base.precompile(Tuple{typeof(show), Base.TTY, summary_callback_type})
-        @assert Base.precompile(Tuple{typeof(show), IOContext{Base.TTY}, MIME"text/plain",
-                                      summary_callback_type})
-        @assert Base.precompile(Tuple{summary_callback_type, Base.TTY})
+        # @assert Base.precompile(Tuple{typeof(show), IOContext{Base.TTY}, MIME"text/plain",
+        #                               summary_callback_type})
+        # @assert Base.precompile(Tuple{summary_callback_type, Base.TTY})
 
         # TODO: SteadyStateCallback, AnalysisCallback
 
diff --git a/src/semidiscretization/semidiscretization.jl b/src/semidiscretization/semidiscretization.jl
index dd5c3c4791d..f41c7ea4a7f 100644
--- a/src/semidiscretization/semidiscretization.jl
+++ b/src/semidiscretization/semidiscretization.jl
@@ -277,6 +277,11 @@ function _jacobian_ad_forward(semi, t0, u0_ode, du_ode, config)
     return J
 end
 
+# unpack u if it is wrapped in VectorOfArray (mainly for DGMulti solvers)
+jacobian_ad_forward(semi::AbstractSemidiscretization, t0, u0_ode::VectorOfArray) = jacobian_ad_forward(semi,
+                                                                                                       t0,
+                                                                                                       parent(u0_ode))
+
 # This version is specialized to `StructArray`s used by some `DGMulti` solvers.
 # We need to convert the numerical solution vectors since ForwardDiff cannot
 # handle arrays of `SVector`s.
diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl
index 2d588c5c79d..a48a6288eba 100644
--- a/src/solvers/dgmulti/dg.jl
+++ b/src/solvers/dgmulti/dg.jl
@@ -128,6 +128,14 @@ end
 # interface with semidiscretization_hyperbolic
 wrap_array(u_ode, mesh::DGMultiMesh, equations, dg::DGMulti, cache) = u_ode
 wrap_array_native(u_ode, mesh::DGMultiMesh, equations, dg::DGMulti, cache) = u_ode
+
+# used to initialize `u_ode` in `semidiscretize`
+function allocate_coefficients(mesh::DGMultiMesh, equations, dg::DGMulti, cache)
+    return VectorOfArray(allocate_nested_array(real(dg), nvariables(equations),
+                                               size(mesh.md.x), dg))
+end
+wrap_array(u_ode::VectorOfArray, mesh::DGMultiMesh, equations, dg::DGMulti, cache) = parent(u_ode)
+
 function digest_boundary_conditions(boundary_conditions::NamedTuple{Keys, ValueTypes},
                                     mesh::DGMultiMesh,
                                     dg::DGMulti,
@@ -199,10 +207,6 @@ function create_cache(mesh::DGMultiMesh{NDIMS}, equations, dg::DGMultiWeakForm,
             local_values_threaded, flux_threaded, rotated_flux_threaded)
 end
 
-function allocate_coefficients(mesh::DGMultiMesh, equations, dg::DGMulti, cache)
-    return allocate_nested_array(real(dg), nvariables(equations), size(mesh.md.x), dg)
-end
-
 function compute_coefficients!(u, initial_condition, t,
                                mesh::DGMultiMesh, equations, dg::DGMulti, cache)
     md = mesh.md
diff --git a/src/visualization/types.jl b/src/visualization/types.jl
index b294ce25607..cf7c9857730 100644
--- a/src/visualization/types.jl
+++ b/src/visualization/types.jl
@@ -614,6 +614,20 @@ function PlotData1D(u, mesh::TreeMesh, equations, solver, cache;
                       orientation_x)
 end
 
+# unwrap u if it is VectorOfArray 
+PlotData1D(u::VectorOfArray, mesh, equations, dg::DGMulti{1}, cache; kwargs...) = PlotData1D(parent(u),
+                                                                                             mesh,
+                                                                                             equations,
+                                                                                             dg,
+                                                                                             cache;
+                                                                                             kwargs...)
+PlotData2D(u::VectorOfArray, mesh, equations, dg::DGMulti{2}, cache; kwargs...) = PlotData2D(parent(u),
+                                                                                             mesh,
+                                                                                             equations,
+                                                                                             dg,
+                                                                                             cache;
+                                                                                             kwargs...)
+
 function PlotData1D(u, mesh, equations, solver, cache;
                     solution_variables = nothing, nvisnodes = nothing,
                     slice = :x, point = (0.0, 0.0, 0.0), curve = nothing)
diff --git a/test/test_dgmulti_1d.jl b/test/test_dgmulti_1d.jl
index 1c3cd604df1..72d343d565b 100644
--- a/test/test_dgmulti_1d.jl
+++ b/test/test_dgmulti_1d.jl
@@ -25,7 +25,8 @@ isdir(outdir) && rm(outdir, recursive = true)
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -41,7 +42,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -65,7 +67,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -89,7 +92,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -106,7 +110,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -130,7 +135,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -161,7 +167,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -183,7 +190,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -207,14 +215,14 @@ end
                         cells_per_dimension=(8,),
                         approximation_type=SBP(),
                         l2=[
-                            3.03001101100507e-6,
-                            1.692177335948727e-5,
-                            3.002634351734614e-16,
-                            1.1636653574178203e-15
+                            3.0300196635805022e-6,
+                            1.6921833812545857e-5,
+                            2.9844594164368975e-16,
+                            1.1012004949980629e-15
                         ],
                         linf=[
-                            1.2043401988570679e-5,
-                            5.346847010329059e-5,
+                            1.2043309307818717e-5,
+                            5.346754311919e-5,
                             9.43689570931383e-16,
                             2.220446049250313e-15
                         ])
@@ -224,7 +232,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -250,7 +259,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 end
diff --git a/test/test_dgmulti_2d.jl b/test/test_dgmulti_2d.jl
index 2e41591d52c..b1c14837fef 100644
--- a/test/test_dgmulti_2d.jl
+++ b/test/test_dgmulti_2d.jl
@@ -37,7 +37,8 @@ isdir(outdir) && rm(outdir, recursive = true)
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -65,7 +66,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -93,7 +95,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -121,7 +124,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -150,7 +154,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -180,19 +185,24 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
 @trixi_testset "elixir_euler_bilinear.jl (Bilinear quadrilateral elements, SBP, flux differencing)" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_bilinear.jl"),
                         l2=[
-                            1.0259432774540821e-5, 9.014087689495575e-6,
-                            9.01408768888544e-6, 2.738953324859446e-5
+                            1.0267413589968656e-5,
+                            9.03069720963081e-6,
+                            9.030697209721065e-6,
+                            2.7436672091049314e-5
                         ],
                         linf=[
-                            7.362605996297233e-5, 6.874189724781488e-5,
-                            6.874189703509614e-5, 0.00019124355334110277
+                            7.36251369879426e-5,
+                            6.874041557969335e-5,
+                            6.874041552329402e-5,
+                            0.00019123932693609902
                         ])
     # Ensure that we do not have excessive memory allocations 
     # (e.g., from type instabilities) 
@@ -200,19 +210,24 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
 @trixi_testset "elixir_euler_curved.jl (Quadrilateral elements, SBP, flux differencing)" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_curved.jl"),
                         l2=[
-                            1.7204593127904542e-5, 1.5921547179522804e-5,
-                            1.5921547180107928e-5, 4.894071422525737e-5
+                            1.7209164346836478e-5,
+                            1.5928649356474767e-5,
+                            1.5928649356802847e-5,
+                            4.8963394546089164e-5
                         ],
                         linf=[
-                            0.00010525416937667842, 0.00010003778102718464,
-                            0.00010003778071832059, 0.0003642628211952825
+                            0.00010525404319428056,
+                            0.00010003768703326088,
+                            0.00010003768694910598,
+                            0.0003642622844113319
                         ])
     # Ensure that we do not have excessive memory allocations 
     # (e.g., from type instabilities) 
@@ -220,7 +235,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -229,16 +245,16 @@ end
                         approximation_type=GaussSBP(),
                         surface_integral=SurfaceIntegralWeakForm(FluxHLL(min_max_speed_naive)),
                         l2=[
-                            3.4666312079259457e-6,
-                            3.4392774480368986e-6,
-                            3.439277447953705e-6,
-                            1.0965598424665836e-5
+                            3.4664508443541302e-6,
+                            3.4389354928807557e-6,
+                            3.438935492692069e-6,
+                            1.0965259031107001e-5
                         ],
                         linf=[
-                            1.1327280377004811e-5,
-                            1.1343911926253725e-5,
-                            1.1343911906935844e-5,
-                            3.679582619220412e-5
+                            1.1326776948594741e-5,
+                            1.1343379410666543e-5,
+                            1.1343379308081936e-5,
+                            3.679395547040443e-5
                         ],
                         rtol=2 * sqrt(eps()))
     # Ensure that we do not have excessive memory allocations 
@@ -247,7 +263,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -257,16 +274,16 @@ end
                         volume_integral=VolumeIntegralWeakForm(),
                         surface_integral=SurfaceIntegralWeakForm(FluxHLL(min_max_speed_naive)),
                         l2=[
-                            7.905498158659466e-6,
-                            8.731690809663625e-6,
-                            8.731690811576996e-6,
-                            2.9113296018693953e-5
+                            7.906577233455084e-6,
+                            8.733496764560955e-6,
+                            8.733496764569368e-6,
+                            2.9118523221923932e-5
                         ],
                         linf=[
-                            3.298811230090237e-5,
-                            4.032272476939269e-5,
-                            4.032272526011127e-5,
-                            0.00012013725458537294
+                            3.29875525286738e-5,
+                            4.032296588629691e-5,
+                            4.032296651601541e-5,
+                            0.00012013778890995752
                         ])
     # Ensure that we do not have excessive memory allocations 
     # (e.g., from type instabilities) 
@@ -274,23 +291,24 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
 @trixi_testset "elixir_euler_hohqmesh.jl (Quadrilateral elements, SBP, flux differencing)" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_hohqmesh.jl"),
                         l2=[
-                            0.0008153911341517156,
-                            0.0007768159701964676,
-                            0.00047902606811690694,
-                            0.0015551846076348535
+                            0.0008153911341539523,
+                            0.0007768159702011952,
+                            0.0004790260681142826,
+                            0.0015551846076274918
                         ],
                         linf=[
-                            0.0029301131365355726,
-                            0.0034427051471457304,
-                            0.0028721569841545502,
-                            0.011125365074589944
+                            0.002930113136531798,
+                            0.003442705146861069,
+                            0.002872156984277563,
+                            0.011125365075300486
                         ])
     # Ensure that we do not have excessive memory allocations 
     # (e.g., from type instabilities) 
@@ -298,7 +316,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -319,7 +338,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -340,7 +360,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -360,7 +381,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -387,7 +409,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -415,7 +438,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -424,12 +448,16 @@ end
                                  "elixir_euler_rayleigh_taylor_instability.jl"),
                         cells_per_dimension=(8, 8), tspan=(0.0, 0.2),
                         l2=[
-                            0.07097806723891838, 0.005168550941966817,
-                            0.013820912272220933, 0.03243357220022434
+                            0.07097806924106471,
+                            0.005168545523460976,
+                            0.013820905434253445,
+                            0.03243358478653133
                         ],
                         linf=[
-                            0.4783395896753895, 0.02244629340135818,
-                            0.04023357731088538, 0.08515807256615027
+                            0.4783395366569936,
+                            0.022446258588973853,
+                            0.04023354591166624,
+                            0.08515791118082117
                         ])
     # Ensure that we do not have excessive memory allocations 
     # (e.g., from type instabilities) 
@@ -437,7 +465,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -462,7 +491,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -487,7 +517,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -512,7 +543,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -546,7 +578,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -582,7 +615,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -602,7 +636,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -628,7 +663,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -660,7 +696,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -693,7 +730,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -716,7 +754,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -739,7 +778,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -815,7 +855,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -824,15 +865,15 @@ end
                         cells_per_dimension=8, element_type=Quad(),
                         approximation_type=SBP(),
                         l2=[
-                            0.0020316462913319046,
-                            0.023669019044882247,
-                            0.03446194752754684,
-                            1.9333465252381796e-15
+                            0.0020316463892983217,
+                            0.02366902012965938,
+                            0.03446194535725363,
+                            1.921676942941478e-15
                         ],
                         linf=[
-                            0.010385010095182778,
-                            0.08750628939565086,
-                            0.12088392994348407,
+                            0.010384996665098178,
+                            0.08750632767286826,
+                            0.12088391569555768,
                             9.325873406851315e-15
                         ])
     # Ensure that we do not have excessive memory allocations 
@@ -841,7 +882,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -850,16 +892,16 @@ end
                         cells_per_dimension=8, element_type=Tri(),
                         approximation_type=SBP(),
                         l2=[
-                            0.004180680322490383,
-                            0.07026192411558974,
-                            0.11815151697006446,
-                            2.329788936151192e-15
+                            0.004180679992535108,
+                            0.07026193567927695,
+                            0.11815151184746633,
+                            2.3786840926019625e-15
                         ],
                         linf=[
-                            0.02076003852980346,
-                            0.29169601664914424,
-                            0.5674183379872275,
-                            1.1546319456101628e-14
+                            0.020760033097378283,
+                            0.29169608872805686,
+                            0.567418412384793,
+                            1.1102230246251565e-14
                         ])
     # Ensure that we do not have excessive memory allocations 
     # (e.g., from type instabilities) 
@@ -867,7 +909,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -878,16 +921,16 @@ end
                         # The last l2, linf error are the L2 projection error in approximating `b`, so they are not
                         # zero for general non-collocated quadrature rules (e.g., for `element_type=Tri()`, `polydeg > 2`).
                         l2=[
-                            0.0008309356912456799,
-                            0.01522451288799231,
-                            0.016033969387208476,
-                            1.2820247308150876e-5
+                            0.0008309358577296097,
+                            0.015224511207450263,
+                            0.016033971785878454,
+                            1.282024730815488e-5
                         ],
                         linf=[
-                            0.001888045014140971,
-                            0.05466838692127718,
-                            0.06345885709961152,
-                            3.3989933098554914e-5
+                            0.0018880416154898327,
+                            0.05466845626696504,
+                            0.06345896594568323,
+                            3.398993309877696e-5
                         ])
     # Ensure that we do not have excessive memory allocations 
     # (e.g., from type instabilities) 
@@ -895,7 +938,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -908,16 +952,16 @@ end
                         # a `(polydeg + 1)`-point Gauss quadrature rule in each coordinate (in general, StartUpDG.jl defaults
                         # to the quadrature rule with the fewest number of points which exactly integrates the mass matrix).
                         l2=[
-                            7.460461950323111e-5,
-                            0.003685589808444905,
-                            0.0039101604749887785,
-                            2.0636891126652983e-15
+                            7.460473151203597e-5,
+                            0.0036855901000765463,
+                            0.003910160802530521,
+                            6.743418333559633e-15
                         ],
                         linf=[
-                            0.000259995400729629,
-                            0.0072236204211630906,
-                            0.010364675200833062,
-                            1.021405182655144e-14
+                            0.0002599957400737374,
+                            0.007223608258381642,
+                            0.010364657535841815,
+                            2.042810365310288e-14
                         ])
     # Ensure that we do not have excessive memory allocations 
     # (e.g., from type instabilities) 
@@ -925,7 +969,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 end
diff --git a/test/test_dgmulti_3d.jl b/test/test_dgmulti_3d.jl
index 5af8c9ee911..33df79a4b38 100644
--- a/test/test_dgmulti_3d.jl
+++ b/test/test_dgmulti_3d.jl
@@ -33,7 +33,8 @@ isdir(outdir) && rm(outdir, recursive = true)
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -62,7 +63,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -91,7 +93,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -117,7 +120,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -140,7 +144,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -162,7 +167,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -191,7 +197,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -222,7 +229,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -230,18 +238,18 @@ end
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_taylor_green_vortex.jl"),
                         polydeg=3, tspan=(0.0, 1.0), cells_per_dimension=(2, 2, 2),
                         l2=[
-                            0.0003612827827560599,
-                            0.06219350883951729,
-                            0.062193508839503864,
-                            0.08121963221634831,
-                            0.07082703570808184
+                            0.00036128264902931644,
+                            0.06219350570157671,
+                            0.062193505701565316,
+                            0.08121963725209637,
+                            0.0708269605813566
                         ],
                         linf=[
-                            0.0007893509649821162,
-                            0.1481953939988877,
-                            0.14819539399791176,
-                            0.14847291108358926,
-                            0.21313533492212855
+                            0.0007893500666786846,
+                            0.14819541663164099,
+                            0.14819541663231595,
+                            0.148472950090691,
+                            0.2131352319423172
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -249,7 +257,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -258,18 +267,18 @@ end
                         polydeg=3, approximation_type=GaussSBP(), tspan=(0.0, 1.0),
                         cells_per_dimension=(2, 2, 2),
                         l2=[
-                            0.00036128278275524326,
-                            0.062193508839511434,
-                            0.06219350883949677,
-                            0.08121963221635205,
-                            0.07082703570765223
+                            0.0003612826490291416,
+                            0.06219350570157282,
+                            0.06219350570156088,
+                            0.08121963725209767,
+                            0.07082696058040763
                         ],
                         linf=[
-                            0.000789350964946367,
-                            0.14819539399525805,
-                            0.14819539399590542,
-                            0.14847291107658706,
-                            0.21313533492059378
+                            0.0007893500666356079,
+                            0.14819541663140268,
+                            0.14819541663222624,
+                            0.14847295009030398,
+                            0.21313523192395678
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -277,7 +286,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -311,7 +321,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -347,7 +358,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -373,7 +385,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 
@@ -387,7 +400,8 @@ end
         t = sol.t[end]
         u_ode = sol.u[end]
         du_ode = similar(u_ode)
-        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
     end
 end
 end
diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl
index ceefb65e99b..6271534a5ac 100644
--- a/test/test_parabolic_2d.jl
+++ b/test/test_parabolic_2d.jl
@@ -124,16 +124,16 @@ end
                                  "elixir_navierstokes_convergence.jl"),
                         cells_per_dimension=(4, 4), tspan=(0.0, 0.1),
                         l2=[
-                            0.0015355076812510957,
-                            0.0033843168272696756,
-                            0.0036531858107443434,
-                            0.009948436427519214
+                            0.0015355076237431118,
+                            0.003384316785885901,
+                            0.0036531858026850757,
+                            0.009948436101649498
                         ],
                         linf=[
-                            0.005522560467190019,
-                            0.013425258500730508,
-                            0.013962115643482154,
-                            0.027483102120502423
+                            0.005522560543588462,
+                            0.013425258431728926,
+                            0.013962115936715924,
+                            0.027483099961148838
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -150,16 +150,16 @@ end
                                  "elixir_navierstokes_convergence_curved.jl"),
                         cells_per_dimension=(4, 4), tspan=(0.0, 0.1),
                         l2=[
-                            0.004255101916146187,
-                            0.011118488923215765,
-                            0.011281831283462686,
-                            0.03573656447388509
+                            0.0042551020940351444,
+                            0.011118489080358264,
+                            0.011281831362358863,
+                            0.035736565778376306
                         ],
                         linf=[
-                            0.015071710669706473,
-                            0.04103132025858458,
-                            0.03990424085750277,
-                            0.1309401718598764
+                            0.015071709836357083,
+                            0.04103131887989486,
+                            0.03990424032494211,
+                            0.13094018584692968
                         ],)
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -176,16 +176,16 @@ end
                                  "elixir_navierstokes_lid_driven_cavity.jl"),
                         cells_per_dimension=(4, 4), tspan=(0.0, 0.5),
                         l2=[
-                            0.00022156125227115747,
-                            0.028318325921401,
-                            0.009509168701070296,
-                            0.028267900513550506
+                            0.0002215612357465129,
+                            0.028318325887331217,
+                            0.009509168805093485,
+                            0.028267893004691534
                         ],
                         linf=[
-                            0.001562278941298234,
-                            0.14886653390744856,
-                            0.0716323565533752,
-                            0.19472785105241996
+                            0.0015622793960574644,
+                            0.1488665309341318,
+                            0.07163235778907852,
+                            0.19472797949052278
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
diff --git a/test/test_parabolic_3d.jl b/test/test_parabolic_3d.jl
index 2690a08cbb9..f8f3e9bbc15 100644
--- a/test/test_parabolic_3d.jl
+++ b/test/test_parabolic_3d.jl
@@ -17,18 +17,18 @@ isdir(outdir) && rm(outdir, recursive = true)
                                  "elixir_navierstokes_convergence.jl"),
                         cells_per_dimension=(4, 4, 4), tspan=(0.0, 0.1),
                         l2=[
-                            0.0005532847115849239,
-                            0.000659263490965341,
-                            0.0007776436127362806,
-                            0.0006592634909662951,
-                            0.0038073628897809185
+                            0.0005532846479614563,
+                            0.000659263463988067,
+                            0.0007776436003494915,
+                            0.000659263463988129,
+                            0.0038073624941206956
                         ],
                         linf=[
-                            0.0017039861523615585,
-                            0.002628561703560073,
-                            0.003531057425112172,
-                            0.0026285617036090336,
-                            0.015587829540351095
+                            0.001703986341275776,
+                            0.0026285618026252733,
+                            0.00353105737957371,
+                            0.002628561802588858,
+                            0.015587831432887
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -45,18 +45,18 @@ end
                                  "elixir_navierstokes_convergence_curved.jl"),
                         cells_per_dimension=(4, 4, 4), tspan=(0.0, 0.1),
                         l2=[
-                            0.0014027227251207474,
-                            0.0021322235533273513,
-                            0.0027873741447455194,
-                            0.0024587473070627423,
-                            0.00997836818019202
+                            0.0014027227340680359,
+                            0.0021322235583299425,
+                            0.002787374145873934,
+                            0.002458747307842109,
+                            0.009978368214450204
                         ],
                         linf=[
-                            0.006341750402837576,
-                            0.010306014252246865,
-                            0.01520740250924979,
-                            0.010968264045485565,
-                            0.047454389831591115
+                            0.006341750448945582,
+                            0.010306014425485621,
+                            0.015207402553448324,
+                            0.010968264060799426,
+                            0.04745438898236998
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -73,18 +73,18 @@ end
                                  "elixir_navierstokes_taylor_green_vortex.jl"),
                         cells_per_dimension=(4, 4, 4), tspan=(0.0, 0.25),
                         l2=[
-                            0.0001825713444029892,
-                            0.015589736382772248,
-                            0.015589736382771884,
-                            0.021943924667273653,
-                            0.01927370280244222
+                            0.00018257125088549987,
+                            0.015589736346235174,
+                            0.015589736346235415,
+                            0.021943924698669025,
+                            0.019273688367502154
                         ],
                         linf=[
-                            0.0006268463584697681,
-                            0.03218881662749007,
-                            0.03218881662697948,
-                            0.053872495395614256,
-                            0.05183822000984151
+                            0.0006268461326666142,
+                            0.03218881686243058,
+                            0.03218881686357877,
+                            0.053872494644958,
+                            0.05183811394229565
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
diff --git a/test/test_threaded.jl b/test/test_threaded.jl
index 83544aa0072..d85e0332ef5 100644
--- a/test/test_threaded.jl
+++ b/test/test_threaded.jl
@@ -39,6 +39,7 @@ Trixi.MPI.Barrier(Trixi.mpi_comm())
             t = sol.t[end]
             u_ode = sol.u[end]
             du_ode = similar(u_ode)
+            Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
             @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
         end
     end
@@ -65,6 +66,7 @@ Trixi.MPI.Barrier(Trixi.mpi_comm())
             t = sol.t[end]
             u_ode = sol.u[end]
             du_ode = similar(u_ode)
+            Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
             @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
         end
     end
@@ -81,6 +83,7 @@ Trixi.MPI.Barrier(Trixi.mpi_comm())
             t = sol.t[end]
             u_ode = sol.u[end]
             du_ode = similar(u_ode)
+            Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
             @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
         end
     end
@@ -108,6 +111,7 @@ Trixi.MPI.Barrier(Trixi.mpi_comm())
             t = sol.t[end]
             u_ode = sol.u[end]
             du_ode = similar(u_ode)
+            Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
             @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
         end
     end
@@ -134,6 +138,7 @@ Trixi.MPI.Barrier(Trixi.mpi_comm())
             t = sol.t[end]
             u_ode = sol.u[end]
             du_ode = similar(u_ode)
+            Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
             @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
         end
     end
@@ -152,6 +157,7 @@ Trixi.MPI.Barrier(Trixi.mpi_comm())
             t = sol.t[end]
             u_ode = sol.u[end]
             du_ode = similar(u_ode)
+            Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
             @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
         end
     end
@@ -169,6 +175,7 @@ Trixi.MPI.Barrier(Trixi.mpi_comm())
             t = sol.t[end]
             u_ode = sol.u[end]
             du_ode = similar(u_ode)
+            Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
             @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
         end
     end
@@ -196,6 +203,7 @@ Trixi.MPI.Barrier(Trixi.mpi_comm())
             t = sol.t[end]
             u_ode = sol.u[end]
             du_ode = similar(u_ode)
+            Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
             @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
         end
     end
@@ -220,6 +228,7 @@ end
             t = sol.t[end]
             u_ode = sol.u[end]
             du_ode = similar(u_ode)
+            Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
             @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
         end
     end
@@ -245,6 +254,7 @@ end
             t = sol.t[end]
             u_ode = sol.u[end]
             du_ode = similar(u_ode)
+            Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
             @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
         end
     end
@@ -270,6 +280,7 @@ end
             t = sol.t[end]
             u_ode = sol.u[end]
             du_ode = similar(u_ode)
+            Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
             @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
         end
     end
@@ -298,6 +309,7 @@ end
             t = sol.t[end]
             u_ode = sol.u[end]
             du_ode = similar(u_ode)
+            Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
             @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
         end
     end
@@ -385,6 +397,7 @@ end
             t = sol.t[end]
             u_ode = sol.u[end]
             du_ode = similar(u_ode)
+            Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
             @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
         end
     end
@@ -394,16 +407,16 @@ end
                                      "elixir_euler_curved.jl"),
                             alg=RDPK3SpFSAL49(thread = OrdinaryDiffEq.True()),
                             l2=[
-                                1.7204593127904542e-5,
-                                1.5921547179522804e-5,
-                                1.5921547180107928e-5,
-                                4.894071422525737e-5
+                                1.720916434676505e-5,
+                                1.5928649356300228e-5,
+                                1.5928649356913923e-5,
+                                4.896339454587786e-5
                             ],
                             linf=[
-                                0.00010525416930584619,
-                                0.00010003778091061122,
-                                0.00010003778085621029,
-                                0.00036426282101720275
+                                0.00010525404319960963,
+                                0.00010003768703459315,
+                                0.00010003768696797977,
+                                0.0003642622844073351
                             ])
 
         # Ensure that we do not have excessive memory allocations
@@ -412,6 +425,7 @@ end
             t = sol.t[end]
             u_ode = sol.u[end]
             du_ode = similar(u_ode)
+            Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
             @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
         end
     end
@@ -438,6 +452,7 @@ end
             t = sol.t[end]
             u_ode = sol.u[end]
             du_ode = similar(u_ode)
+            Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
             @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
         end
     end
@@ -464,6 +479,7 @@ end
             t = sol.t[end]
             u_ode = sol.u[end]
             du_ode = similar(u_ode)
+            Trixi.rhs!(du_ode, u_ode, semi, t) # run once first to deal with spurious allocations
             @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
         end
     end
diff --git a/test/test_visualization.jl b/test/test_visualization.jl
index 5c7e5dbbd1f..dd7021e24ed 100644
--- a/test/test_visualization.jl
+++ b/test/test_visualization.jl
@@ -98,7 +98,12 @@ test_examples_2d = Dict("TreeMesh" => ("tree_2d_dgsem",
 
         semi = sol.prob.p
         if mesh == "DGMulti"
-            scalar_data = StructArrays.component(sol.u[end], 1)
+            if sol.u[end] isa Trixi.VectorOfArray
+                u = parent(sol.u[end])
+            else
+                u = sol.u[end]
+            end
+            scalar_data = StructArrays.component(u, 1)
             @test_nowarn_mod Plots.plot(ScalarPlotData2D(scalar_data, semi))
         else
             cache = semi.cache