diff --git a/src/semidiscretization/semidiscretization_hyperbolic.jl b/src/semidiscretization/semidiscretization_hyperbolic.jl
index 2a563c0222..b51b10feab 100644
--- a/src/semidiscretization/semidiscretization_hyperbolic.jl
+++ b/src/semidiscretization/semidiscretization_hyperbolic.jl
@@ -409,4 +409,25 @@ function rhs!(du_ode, u_ode, semi::SemidiscretizationHyperbolic, t)
 
     return nothing
 end
+
+# `rhs!` for partitioned Runge-Kutta methods, such as the Paired Explicit Runge-Kutta (PERK) methods
+function rhs!(du_ode, u_ode, semi::SemidiscretizationHyperbolic, t,
+              element_indices, interface_indices, boundary_indices, mortar_indices)
+    @unpack mesh, equations, boundary_conditions, source_terms, solver, cache = semi
+
+    u = wrap_array(u_ode, mesh, equations, solver, cache)
+    du = wrap_array(du_ode, mesh, equations, solver, cache)
+
+    # TODO: Taal decide, do we need to pass the mesh?
+    time_start = time_ns()
+    @trixi_timeit timer() "rhs! (part.)" rhs!(du, u, t, mesh, equations,
+                                              boundary_conditions, source_terms,
+                                              solver, cache,
+                                              element_indices, interface_indices,
+                                              boundary_indices, mortar_indices)
+    runtime = time_ns() - time_start
+    put!(semi.performance_counter, runtime)
+
+    return nothing
+end
 end # @muladd
diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl
index c24b8f2a29..c5c6720e14 100644
--- a/src/solvers/dgsem_tree/dg_1d.jl
+++ b/src/solvers/dgsem_tree/dg_1d.jl
@@ -70,33 +70,37 @@ end
 function rhs!(du, u, t,
               mesh::TreeMesh{1}, equations,
               boundary_conditions, source_terms::Source,
-              dg::DG, cache) where {Source}
+              dg::DG, cache,
+              element_indices = eachelement(dg, cache),
+              interface_indices = eachinterface(dg, cache),
+              boundary_indices = eachboundary(dg, cache),
+              mortar_indices = nothing) where {Source}
     # Reset du
-    @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache)
+    @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache, element_indices)
 
     # Calculate volume integral
     @trixi_timeit timer() "volume integral" begin
         calc_volume_integral!(du, u, mesh,
                               have_nonconservative_terms(equations), equations,
-                              dg.volume_integral, dg, cache)
+                              dg.volume_integral, dg, cache, element_indices)
     end
 
     # Prolong solution to interfaces
     @trixi_timeit timer() "prolong2interfaces" begin
-        prolong2interfaces!(cache, u, mesh, equations, dg)
+        prolong2interfaces!(cache, u, mesh, equations, dg, interface_indices)
     end
 
     # Calculate interface fluxes
     @trixi_timeit timer() "interface flux" begin
         calc_interface_flux!(cache.elements.surface_flux_values, mesh,
                              have_nonconservative_terms(equations), equations,
-                             dg.surface_integral, dg, cache)
+                             dg.surface_integral, dg, cache, interface_indices)
     end
 
     # Prolong solution to boundaries
     @trixi_timeit timer() "prolong2boundaries" begin
         prolong2boundaries!(cache, u, mesh, equations,
-                            dg.surface_integral, dg)
+                            dg.surface_integral, dg, boundary_indices)
     end
 
     # Calculate boundary fluxes
@@ -108,15 +112,16 @@ function rhs!(du, u, t,
     # Calculate surface integrals
     @trixi_timeit timer() "surface integral" begin
         calc_surface_integral!(du, u, mesh, equations,
-                               dg.surface_integral, dg, cache)
+                               dg.surface_integral, dg, cache, element_indices)
     end
 
     # Apply Jacobian from mapping to reference element
-    @trixi_timeit timer() "Jacobian" apply_jacobian!(du, mesh, equations, dg, cache)
+    @trixi_timeit timer() "Jacobian" apply_jacobian!(du, mesh, equations, dg, cache,
+                                                     element_indices)
 
     # Calculate source terms
     @trixi_timeit timer() "source terms" begin
-        calc_sources!(du, u, t, source_terms, equations, dg, cache)
+        calc_sources!(du, u, t, source_terms, equations, dg, cache, element_indices)
     end
 
     return nothing
@@ -245,11 +250,12 @@ end
 function calc_volume_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1}},
                                have_nonconservative_terms, equations,
                                volume_integral::VolumeIntegralPureLGLFiniteVolumeO2,
-                               dg::DGSEM, cache)
+                               dg::DGSEM, cache,
+                               element_indices = eachelement(dg, cache))
     @unpack x_interfaces, volume_flux_fv, reconstruction_mode, slope_limiter = volume_integral
 
     # Calculate LGL second-order FV volume integral
-    @threaded for element in eachelement(dg, cache)
+    @threaded for element in element_indices
         fvO2_kernel!(du, u, mesh,
                      have_nonconservative_terms, equations,
                      volume_flux_fv, dg, cache, element,
@@ -408,12 +414,13 @@ end
     return nothing
 end
 
-function prolong2interfaces!(cache, u, mesh::TreeMesh{1}, equations, dg::DG)
+function prolong2interfaces!(cache, u, mesh::TreeMesh{1}, equations, dg::DG,
+                             interface_indices = eachinterface(dg, cache))
     @unpack interfaces = cache
     @unpack neighbor_ids = interfaces
     interfaces_u = interfaces.u
 
-    @threaded for interface in eachinterface(dg, cache)
+    @threaded for interface in interface_indices
         left_element = neighbor_ids[1, interface]
         right_element = neighbor_ids[2, interface]
 
@@ -430,11 +437,12 @@ end
 function calc_interface_flux!(surface_flux_values,
                               mesh::TreeMesh{1},
                               have_nonconservative_terms::False, equations,
-                              surface_integral, dg::DG, cache)
+                              surface_integral, dg::DG, cache,
+                              interface_indices = eachinterface(dg, cache))
     @unpack surface_flux = surface_integral
     @unpack u, neighbor_ids, orientations = cache.interfaces
 
-    @threaded for interface in eachinterface(dg, cache)
+    @threaded for interface in interface_indices
         # Get neighboring elements
         left_id = neighbor_ids[1, interface]
         right_id = neighbor_ids[2, interface]
@@ -450,6 +458,8 @@ function calc_interface_flux!(surface_flux_values,
 
         # Copy flux to left and right element storage
         for v in eachvariable(equations)
+            # TODO: Accumulate fluxes, i.e., `+=`? 
+            # Reset would then be handled by `reset_du!`
             surface_flux_values[v, left_direction, left_id] = flux[v]
             surface_flux_values[v, right_direction, right_id] = flux[v]
         end
@@ -461,11 +471,12 @@ end
 function calc_interface_flux!(surface_flux_values,
                               mesh::TreeMesh{1},
                               have_nonconservative_terms::True, equations,
-                              surface_integral, dg::DG, cache)
+                              surface_integral, dg::DG, cache,
+                              interface_indices = eachinterface(dg, cache))
     surface_flux, nonconservative_flux = surface_integral.surface_flux
     @unpack u, neighbor_ids, orientations = cache.interfaces
 
-    @threaded for interface in eachinterface(dg, cache)
+    @threaded for interface in interface_indices
         # Get neighboring elements
         left_id = neighbor_ids[1, interface]
         right_id = neighbor_ids[2, interface]
@@ -501,11 +512,12 @@ function calc_interface_flux!(surface_flux_values,
 end
 
 function prolong2boundaries!(cache, u,
-                             mesh::TreeMesh{1}, equations, surface_integral, dg::DG)
+                             mesh::TreeMesh{1}, equations, surface_integral, dg::DG,
+                             boundary_indices = eachboundary(dg, cache))
     @unpack boundaries = cache
     @unpack neighbor_sides = boundaries
 
-    @threaded for boundary in eachboundary(dg, cache)
+    @threaded for boundary in boundary_indices
         element = boundaries.neighbor_ids[boundary]
 
         # boundary in x-direction
@@ -616,7 +628,8 @@ function calc_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:A
 end
 
 function calc_surface_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1}},
-                                equations, surface_integral, dg::DGSEM, cache)
+                                equations, surface_integral, dg::DGSEM, cache,
+                                element_indices = eachelement(dg, cache))
     @unpack boundary_interpolation = dg.basis
     @unpack surface_flux_values = cache.elements
 
@@ -626,7 +639,7 @@ function calc_surface_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1
     # into FMAs (see comment at the top of the file).
     factor_1 = boundary_interpolation[1, 1]
     factor_2 = boundary_interpolation[nnodes(dg), 2]
-    @threaded for element in eachelement(dg, cache)
+    @threaded for element in element_indices
         for v in eachvariable(equations)
             # surface at -x
             du[v, 1, element] = (du[v, 1, element] -
@@ -642,10 +655,11 @@ function calc_surface_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1
 end
 
 function apply_jacobian!(du, mesh::TreeMesh{1},
-                         equations, dg::DG, cache)
+                         equations, dg::DG, cache,
+                         element_indices = eachelement(dg, cache))
     @unpack inverse_jacobian = cache.elements
 
-    @threaded for element in eachelement(dg, cache)
+    @threaded for element in element_indices
         factor = -inverse_jacobian[element]
 
         for i in eachnode(dg)
@@ -660,15 +674,17 @@ end
 
 # Need dimension specific version to avoid error at dispatching
 function calc_sources!(du, u, t, source_terms::Nothing,
-                       equations::AbstractEquations{1}, dg::DG, cache)
+                       equations::AbstractEquations{1}, dg::DG, cache,
+                       element_indices = eachelement(dg, cache))
     return nothing
 end
 
 function calc_sources!(du, u, t, source_terms,
-                       equations::AbstractEquations{1}, dg::DG, cache)
+                       equations::AbstractEquations{1}, dg::DG, cache,
+                       element_indices = eachelement(dg, cache))
     @unpack node_coordinates = cache.elements
 
-    @threaded for element in eachelement(dg, cache)
+    @threaded for element in element_indices
         for i in eachnode(dg)
             u_local = get_node_vars(u, equations, dg, i, element)
             x_local = get_node_coords(node_coordinates, equations, dg,
diff --git a/src/solvers/dgsem_tree/dg_1d_parabolic.jl b/src/solvers/dgsem_tree/dg_1d_parabolic.jl
index d6e1147eb6..32e8e29320 100644
--- a/src/solvers/dgsem_tree/dg_1d_parabolic.jl
+++ b/src/solvers/dgsem_tree/dg_1d_parabolic.jl
@@ -16,27 +16,34 @@
 function rhs_parabolic!(du, u, t, mesh::TreeMesh{1},
                         equations_parabolic::AbstractEquationsParabolic,
                         boundary_conditions_parabolic, source_terms,
-                        dg::DG, parabolic_scheme, cache, cache_parabolic)
+                        dg::DG, parabolic_scheme, cache, cache_parabolic,
+                        element_indices = eachelement(dg, cache),
+                        interface_indices = eachinterface(dg, cache),
+                        boundary_indices = eachboundary(dg, cache),
+                        mortar_indices = nothing)
     @unpack viscous_container = cache_parabolic
     @unpack u_transformed, gradients, flux_viscous = viscous_container
 
     # Convert conservative variables to a form more suitable for viscous flux calculations
     @trixi_timeit timer() "transform variables" begin
         transform_variables!(u_transformed, u, mesh, equations_parabolic,
-                             dg, parabolic_scheme, cache, cache_parabolic)
+                             dg, parabolic_scheme, cache, cache_parabolic,
+                             element_indices)
     end
 
     # Compute the gradients of the transformed variables
     @trixi_timeit timer() "calculate gradient" begin
         calc_gradient!(gradients, u_transformed, t, mesh, equations_parabolic,
                        boundary_conditions_parabolic, dg, parabolic_scheme,
-                       cache, cache_parabolic)
+                       cache, cache_parabolic,
+                       element_indices, interface_indices, boundary_indices)
     end
 
     # Compute and store the viscous fluxes
     @trixi_timeit timer() "calculate viscous fluxes" begin
         calc_viscous_fluxes!(flux_viscous, gradients, u_transformed, mesh,
-                             equations_parabolic, dg, cache, cache_parabolic)
+                             equations_parabolic, dg, cache, cache_parabolic,
+                             element_indices)
     end
 
     # The remainder of this function is essentially a regular rhs! for
@@ -54,30 +61,31 @@ function rhs_parabolic!(du, u, t, mesh::TreeMesh{1},
     # TODO: parabolic; reconsider current data structure reuse strategy
 
     # Reset du
-    @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache)
+    @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache, element_indices)
 
     # Calculate volume integral
     @trixi_timeit timer() "volume integral" begin
-        calc_volume_integral!(du, flux_viscous, mesh, equations_parabolic, dg, cache)
+        calc_volume_integral!(du, flux_viscous, mesh, equations_parabolic, dg, cache,
+                              element_indices)
     end
 
     # Prolong solution to interfaces
     @trixi_timeit timer() "prolong2interfaces" begin
         prolong2interfaces!(cache_parabolic, flux_viscous, mesh, equations_parabolic,
-                            dg, cache)
+                            dg, cache, interface_indices)
     end
 
     # Calculate interface fluxes
     @trixi_timeit timer() "interface flux" begin
         calc_interface_flux!(cache_parabolic.elements.surface_flux_values,
                              mesh, equations_parabolic, dg, parabolic_scheme,
-                             cache_parabolic)
+                             cache_parabolic, interface_indices)
     end
 
     # Prolong solution to boundaries
     @trixi_timeit timer() "prolong2boundaries" begin
         prolong2boundaries!(cache_parabolic, flux_viscous, mesh, equations_parabolic,
-                            dg.surface_integral, dg, cache)
+                            dg.surface_integral, dg, cache, boundary_indices)
     end
 
     # Calculate boundary fluxes
@@ -91,12 +99,14 @@ function rhs_parabolic!(du, u, t, mesh::TreeMesh{1},
     # Calculate surface integrals
     @trixi_timeit timer() "surface integral" begin
         calc_surface_integral!(du, u, mesh, equations_parabolic,
-                               dg.surface_integral, dg, cache_parabolic)
+                               dg.surface_integral, dg, cache_parabolic,
+                               element_indices)
     end
 
     # Apply Jacobian from mapping to reference element
     @trixi_timeit timer() "Jacobian" begin
-        apply_jacobian_parabolic!(du, mesh, equations_parabolic, dg, cache_parabolic)
+        apply_jacobian_parabolic!(du, mesh, equations_parabolic, dg, cache_parabolic,
+                                  element_indices)
     end
 
     return nothing
@@ -107,10 +117,11 @@ end
 # TODO: can we avoid copying data?
 function transform_variables!(u_transformed, u, mesh::TreeMesh{1},
                               equations_parabolic::AbstractEquationsParabolic,
-                              dg::DG, parabolic_scheme, cache, cache_parabolic)
+                              dg::DG, parabolic_scheme, cache, cache_parabolic,
+                              element_indices = eachelement(dg, cache))
     transformation = gradient_variable_transformation(equations_parabolic)
 
-    @threaded for element in eachelement(dg, cache)
+    @threaded for element in element_indices
         # Calculate volume terms in one element
         for i in eachnode(dg)
             u_node = get_node_vars(u, equations_parabolic, dg, i, element)
@@ -127,10 +138,11 @@ end
 function calc_volume_integral!(du, flux_viscous,
                                mesh::TreeMesh{1},
                                equations_parabolic::AbstractEquationsParabolic,
-                               dg::DGSEM, cache)
+                               dg::DGSEM, cache,
+                               element_indices = eachelement(dg, cache))
     @unpack derivative_dhat = dg.basis
 
-    @threaded for element in eachelement(dg, cache)
+    @threaded for element in element_indices
         # Calculate volume terms in one element
         for i in eachnode(dg)
             flux_1_node = get_node_vars(flux_viscous, equations_parabolic, dg, i,
@@ -150,12 +162,13 @@ end
 function prolong2interfaces!(cache_parabolic, flux_viscous,
                              mesh::TreeMesh{1},
                              equations_parabolic::AbstractEquationsParabolic,
-                             dg::DG, cache)
+                             dg::DG, cache::NamedTuple,
+                             interface_indices = eachinterface(dg, cache))
     @unpack interfaces = cache_parabolic
     @unpack neighbor_ids = interfaces
     interfaces_u = interfaces.u
 
-    @threaded for interface in eachinterface(dg, cache)
+    @threaded for interface in interface_indices
         left_element = neighbor_ids[1, interface]
         right_element = neighbor_ids[2, interface]
 
@@ -174,10 +187,10 @@ end
 function calc_interface_flux!(surface_flux_values,
                               mesh::TreeMesh{1}, equations_parabolic,
                               dg::DG, parabolic_scheme,
-                              cache_parabolic)
+                              cache_parabolic, interface_indices)
     @unpack neighbor_ids, orientations = cache_parabolic.interfaces
 
-    @threaded for interface in eachinterface(dg, cache_parabolic)
+    @threaded for interface in interface_indices
         # Get neighboring elements
         left_id = neighbor_ids[1, interface]
         right_id = neighbor_ids[2, interface]
@@ -210,12 +223,13 @@ end
 function prolong2boundaries!(cache_parabolic, flux_viscous,
                              mesh::TreeMesh{1},
                              equations_parabolic::AbstractEquationsParabolic,
-                             surface_integral, dg::DG, cache)
+                             surface_integral, dg::DG, cache::NamedTuple,
+                             boundary_indices = eachboundary(dg, cache))
     @unpack boundaries = cache_parabolic
     @unpack neighbor_sides, neighbor_ids = boundaries
     boundaries_u = boundaries.u
 
-    @threaded for boundary in eachboundary(dg, cache_parabolic)
+    @threaded for boundary in boundary_indices
         element = neighbor_ids[boundary]
 
         if neighbor_sides[boundary] == 1
@@ -237,8 +251,9 @@ end
 
 function calc_viscous_fluxes!(flux_viscous, gradients, u_transformed, mesh::TreeMesh{1},
                               equations_parabolic::AbstractEquationsParabolic,
-                              dg::DG, cache, cache_parabolic)
-    @threaded for element in eachelement(dg, cache)
+                              dg::DG, cache, cache_parabolic,
+                              element_indices = eachelement(dg, cache))
+    @threaded for element in element_indices
         for i in eachnode(dg)
             # Get solution and gradients
             u_node = get_node_vars(u_transformed, equations_parabolic, dg, i, element)
@@ -417,10 +432,11 @@ end
 function calc_gradient_volume_integral!(gradients, u_transformed,
                                         mesh::TreeMesh{1}, # for dispatch only
                                         equations_parabolic::AbstractEquationsParabolic,
-                                        dg::DGSEM, cache)
+                                        dg::DGSEM, cache,
+                                        element_indices = eachelement(dg, cache))
     @unpack derivative_dhat = dg.basis
 
-    @threaded for element in eachelement(dg, cache)
+    @threaded for element in element_indices
         # Calculate volume terms in one element,
         # corresponds to `kernel` functions for the hyperbolic part of the flux
         for i in eachnode(dg)
@@ -441,10 +457,11 @@ end
 function calc_gradient_interface_flux!(surface_flux_values,
                                        mesh::TreeMesh{1}, equations_parabolic, dg::DG,
                                        parabolic_scheme,
-                                       cache, cache_parabolic)
+                                       cache, cache_parabolic,
+                                       interface_indices = eachinterface(dg, cache))
     @unpack neighbor_ids, orientations = cache_parabolic.interfaces
 
-    @threaded for interface in eachinterface(dg, cache_parabolic)
+    @threaded for interface in interface_indices
         # Get neighboring elements
         left_id = neighbor_ids[1, interface]
         right_id = neighbor_ids[2, interface]
@@ -474,7 +491,8 @@ end
 function calc_gradient_surface_integral!(gradients,
                                          mesh::TreeMesh{1}, # for dispatch only
                                          equations_parabolic::AbstractEquationsParabolic,
-                                         dg::DGSEM, cache, cache_parabolic)
+                                         dg::DGSEM, cache, cache_parabolic,
+                                         element_indices = eachelement(dg, cache))
     @unpack boundary_interpolation = dg.basis
     @unpack surface_flux_values = cache_parabolic.elements
 
@@ -484,7 +502,7 @@ function calc_gradient_surface_integral!(gradients,
     # into FMAs (see comment at the top of the file).
     factor_1 = boundary_interpolation[1, 1]
     factor_2 = boundary_interpolation[nnodes(dg), 2]
-    @threaded for element in eachelement(dg, cache)
+    @threaded for element in element_indices
         for v in eachvariable(equations_parabolic)
             # surface at -x
             gradients[v, 1, element] = (gradients[v, 1, element] -
@@ -504,11 +522,14 @@ end
 # Calculate the gradient of the transformed variables
 function calc_gradient!(gradients, u_transformed, t, mesh::TreeMesh{1},
                         equations_parabolic, boundary_conditions_parabolic,
-                        dg::DG, parabolic_scheme, cache, cache_parabolic)
+                        dg::DG, parabolic_scheme, cache, cache_parabolic,
+                        element_indices = eachelement(dg, cache),
+                        interface_indices = eachinterface(dg, cache),
+                        boundary_indices = eachboundary(dg, cache))
 
     # Reset du
     @trixi_timeit timer() "reset gradients" begin
-        reset_du!(gradients, dg, cache)
+        reset_du!(gradients, dg, cache, element_indices)
     end
 
     # Calculate volume integral
@@ -521,14 +542,14 @@ function calc_gradient!(gradients, u_transformed, t, mesh::TreeMesh{1},
     @trixi_timeit timer() "prolong2interfaces" prolong2interfaces!(cache_parabolic,
                                                                    u_transformed, mesh,
                                                                    equations_parabolic,
-                                                                   dg)
-
+                                                                   dg, cache,
+                                                                   interface_indices)
     # Calculate interface fluxes
     @trixi_timeit timer() "interface flux" begin
         @unpack surface_flux_values = cache_parabolic.elements
         calc_gradient_interface_flux!(surface_flux_values,
                                       mesh, equations_parabolic, dg, parabolic_scheme,
-                                      cache, cache_parabolic)
+                                      cache, cache_parabolic, interface_indices)
     end
 
     # Prolong solution to boundaries
@@ -536,7 +557,8 @@ function calc_gradient!(gradients, u_transformed, t, mesh::TreeMesh{1},
                                                                    u_transformed, mesh,
                                                                    equations_parabolic,
                                                                    dg.surface_integral,
-                                                                   dg)
+                                                                   dg, cache,
+                                                                   boundary_indices)
 
     # Calculate boundary fluxes
     @trixi_timeit timer() "boundary flux" calc_gradient_boundary_flux!(cache_parabolic,
@@ -550,13 +572,13 @@ function calc_gradient!(gradients, u_transformed, t, mesh::TreeMesh{1},
     # Calculate surface integrals
     @trixi_timeit timer() "surface integral" begin
         calc_gradient_surface_integral!(gradients, mesh, equations_parabolic,
-                                        dg, cache, cache_parabolic)
+                                        dg, cache, cache_parabolic, element_indices)
     end
 
     # Apply Jacobian from mapping to reference element
     @trixi_timeit timer() "Jacobian" begin
         apply_jacobian_parabolic!(gradients, mesh, equations_parabolic, dg,
-                                  cache_parabolic)
+                                  cache_parabolic, element_indices)
     end
 
     return nothing
@@ -593,10 +615,11 @@ end
 #   `du/dt + df/dx = dg/dx + source(x,t)`,
 # where f(u) is the inviscid flux and g(u) is the viscous flux.
 function apply_jacobian_parabolic!(du, mesh::TreeMesh{1},
-                                   equations::AbstractEquationsParabolic, dg::DG, cache)
+                                   equations::AbstractEquationsParabolic, dg::DG, cache,
+                                   element_indices = eachelement(dg, cache))
     @unpack inverse_jacobian = cache.elements
 
-    @threaded for element in eachelement(dg, cache)
+    @threaded for element in element_indices
         factor = inverse_jacobian[element]
 
         for i in eachnode(dg)
diff --git a/src/solvers/solvers.jl b/src/solvers/solvers.jl
index 9eb7bdfa0b..8f26814446 100644
--- a/src/solvers/solvers.jl
+++ b/src/solvers/solvers.jl
@@ -6,10 +6,10 @@
 #! format: noindent
 
 # Used by both `dg::DGSEM` and `dg::FDSBP`
-function reset_du!(du, dg, cache)
+function reset_du!(du, dg, cache, element_indices = eachelement(dg, cache))
     # du .= zero(eltype(du)) doesn't scale when using multiple threads.
     # See https://github.com/trixi-framework/Trixi.jl/pull/924 for a performance comparison.
-    @threaded for element in eachelement(dg, cache)
+    @threaded for element in element_indices
         du[.., element] .= zero(eltype(du))
     end