diff --git a/docs/src/restart.md b/docs/src/restart.md
index 767269ff27d..c7cbcd11852 100644
--- a/docs/src/restart.md
+++ b/docs/src/restart.md
@@ -77,8 +77,7 @@ and its time step number, e.g.:
 ```julia
 integrator = init(ode, CarpenterKennedy2N54(williamson_condition=false),
                   dt=dt, save_everystep=false, callback=callbacks);
-integrator.iter = load_timestep(restart_filename)
-integrator.stats.naccept = integrator.iter
+load_timestep!(integrator, restart_filename)
 ```
 
 Now we can compute the solution:
diff --git a/examples/p4est_2d_dgsem/elixir_advection_restart.jl b/examples/p4est_2d_dgsem/elixir_advection_restart.jl
index 79a35199b83..52917616a6a 100644
--- a/examples/p4est_2d_dgsem/elixir_advection_restart.jl
+++ b/examples/p4est_2d_dgsem/elixir_advection_restart.jl
@@ -35,8 +35,7 @@ integrator = init(ode, CarpenterKennedy2N54(williamson_condition=false),
                   save_everystep=false, callback=callbacks);
 
 # Get the last time index and work with that.
-integrator.iter = load_timestep(restart_filename)
-integrator.stats.naccept = integrator.iter
+load_timestep!(integrator, restart_filename)
 
 
 ###############################################################################
diff --git a/examples/p4est_3d_dgsem/elixir_advection_restart.jl b/examples/p4est_3d_dgsem/elixir_advection_restart.jl
index b27eaab62e2..26d10cf8826 100644
--- a/examples/p4est_3d_dgsem/elixir_advection_restart.jl
+++ b/examples/p4est_3d_dgsem/elixir_advection_restart.jl
@@ -32,8 +32,7 @@ integrator = init(ode, CarpenterKennedy2N54(williamson_condition=false),
                   save_everystep=false, callback=callbacks);
 
 # Get the last time index and work with that.
-integrator.iter = load_timestep(restart_filename)
-integrator.stats.naccept = integrator.iter
+load_timestep!(integrator, restart_filename)
 
 
 ###############################################################################
diff --git a/examples/structured_2d_dgsem/elixir_advection_restart.jl b/examples/structured_2d_dgsem/elixir_advection_restart.jl
index 98c44fac71a..82eaa21333a 100644
--- a/examples/structured_2d_dgsem/elixir_advection_restart.jl
+++ b/examples/structured_2d_dgsem/elixir_advection_restart.jl
@@ -34,8 +34,7 @@ integrator = init(ode, CarpenterKennedy2N54(williamson_condition=false),
                  save_everystep=false, callback=callbacks);
 
 # Get the last time index and work with that.
-integrator.iter = load_timestep(restart_filename)
-integrator.stats.naccept = integrator.iter
+load_timestep!(integrator, restart_filename)
 
 ###############################################################################
 # run the simulation
diff --git a/examples/structured_3d_dgsem/elixir_advection_restart.jl b/examples/structured_3d_dgsem/elixir_advection_restart.jl
index 39d28848c77..921c5310340 100644
--- a/examples/structured_3d_dgsem/elixir_advection_restart.jl
+++ b/examples/structured_3d_dgsem/elixir_advection_restart.jl
@@ -32,8 +32,7 @@ integrator = init(ode, CarpenterKennedy2N54(williamson_condition=false),
                   save_everystep=false, callback=callbacks);
 
 # Get the last time index and work with that.
-integrator.iter = load_timestep(restart_filename)
-integrator.stats.naccept = integrator.iter
+load_timestep!(integrator, restart_filename)
 
 
 ###############################################################################
diff --git a/examples/tree_2d_dgsem/elixir_advection_restart.jl b/examples/tree_2d_dgsem/elixir_advection_restart.jl
index 72efb7d0c84..771ec5aefe7 100644
--- a/examples/tree_2d_dgsem/elixir_advection_restart.jl
+++ b/examples/tree_2d_dgsem/elixir_advection_restart.jl
@@ -32,8 +32,7 @@ integrator = init(ode, alg,
                   save_everystep=false, callback=callbacks)
 
 # Get the last time index and work with that.
-integrator.iter = load_timestep(restart_filename)
-integrator.stats.naccept = integrator.iter
+load_timestep!(integrator, restart_filename)
 
 ###############################################################################
 # run the simulation
diff --git a/examples/tree_3d_dgsem/elixir_advection_restart.jl b/examples/tree_3d_dgsem/elixir_advection_restart.jl
index 3061f165874..b7835ed061f 100644
--- a/examples/tree_3d_dgsem/elixir_advection_restart.jl
+++ b/examples/tree_3d_dgsem/elixir_advection_restart.jl
@@ -31,8 +31,7 @@ integrator = init(ode, CarpenterKennedy2N54(williamson_condition=false),
                   save_everystep=false, callback=callbacks);
 
 # Get the last time index and work with that.
-integrator.iter = load_timestep(restart_filename)
-integrator.stats.naccept = integrator.iter
+load_timestep!(integrator, restart_filename)
 
 
 ###############################################################################
diff --git a/examples/unstructured_2d_dgsem/elixir_euler_restart.jl b/examples/unstructured_2d_dgsem/elixir_euler_restart.jl
index b85cc2c6d70..6653f8662d9 100644
--- a/examples/unstructured_2d_dgsem/elixir_euler_restart.jl
+++ b/examples/unstructured_2d_dgsem/elixir_euler_restart.jl
@@ -33,8 +33,7 @@ integrator = init(ode, CarpenterKennedy2N54(williamson_condition=false),
                   save_everystep=false, callback=callbacks);
 
 # Get the last time index and work with that.
-integrator.iter = load_timestep(restart_filename)
-integrator.stats.naccept = integrator.iter
+load_timestep!(integrator, restart_filename)
 
 
 ###############################################################################
diff --git a/src/Trixi.jl b/src/Trixi.jl
index ec4d20558e5..be43c45b93d 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -253,7 +253,7 @@ export SummaryCallback, SteadyStateCallback, AnalysisCallback, AliveCallback,
        GlmSpeedCallback, LBMCollisionCallback, EulerAcousticsCouplingCallback,
        TrivialCallback, AnalysisCallbackCoupled
 
-export load_mesh, load_time, load_timestep, load_dt
+export load_mesh, load_time, load_timestep, load_timestep!, load_dt
 
 export ControllerThreeLevel, ControllerThreeLevelCombined,
        IndicatorLöhner, IndicatorLoehner, IndicatorMax,
diff --git a/src/callbacks_step/save_restart.jl b/src/callbacks_step/save_restart.jl
index f567a5c7fda..06817a9b730 100644
--- a/src/callbacks_step/save_restart.jl
+++ b/src/callbacks_step/save_restart.jl
@@ -141,6 +141,18 @@ function load_timestep(restart_file::AbstractString)
     end
 end
 
+"""
+    load_timestep!(integrator, restart_file::AbstractString)
+
+Load the time step number saved in a `restart_file` and assign it to both the time step
+number and and the number of accepted steps
+(`iter` and `stats.naccept` in OrdinaryDiffEq.jl, respectively) in `integrator`.
+"""
+function load_timestep!(integrator, restart_file::AbstractString)
+    integrator.iter = load_timestep(restart_file)
+    integrator.stats.naccept = integrator.iter
+end
+
 """
     load_dt(restart_file::AbstractString)