diff --git a/AUTHORS.md b/AUTHORS.md
index 973e311920..abaa3e7e03 100644
--- a/AUTHORS.md
+++ b/AUTHORS.md
@@ -24,6 +24,7 @@ are listed in alphabetical order:
 
 * Maximilian D. Bertrand
 * Benjamin Bolm
+* Simon Candelaresi
 * Jesse Chan
 * Lars Christmann
 * Christof Czernik
diff --git a/docs/make.jl b/docs/make.jl
index 5069e4dc49..57629577dd 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -92,6 +92,7 @@ makedocs(
         "Getting started" => [
             "Overview" => "overview.md",
             "Visualization" => "visualization.md",
+            "Restart simulation" => "restart.md",
         ],
         "Tutorials" => tutorials,
         "Basic building blocks" => [
diff --git a/docs/src/restart.md b/docs/src/restart.md
new file mode 100644
index 0000000000..d24d93cb29
--- /dev/null
+++ b/docs/src/restart.md
@@ -0,0 +1,89 @@
+# [Restart simulation](@id restart)
+
+You can continue running an already finished simulation by first
+preparing the simulation for the restart and then performing the restart.
+Here we suppose that in the first run your simulation stops at time 1.0
+and then you want it to run further to time 2.0.
+
+## [Prepare the simulation for a restart](@id restart_preparation)
+In you original elixir you need to specify to write out restart files.
+Those will later be read for the restart of your simulation.
+This is done almost the same way as writing the snapshots using the
+[`SaveSolutionCallback`](@ref) callback.
+For the restart files it is called [`SaveRestartCallback`](@ref):
+```julia
+save_restart = SaveRestartCallback(interval=100,
+                                   save_final_restart=true)
+```
+Make this part of your `CallbackSet`.
+
+An example is
+[```examples/examples/structured_2d_dgsem/elixir_advection_extended.jl```](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/structured_2d_dgsem/elixir_advection_extended.jl).
+
+
+## [Perform the simulation restart](@id restart_perform)
+Since all of the information about the simulation can be obtained from the
+last snapshot, the restart can be done with relatively few lines
+in an extra elixir file.
+However, some might prefer to keep everything in one elixir and
+conditionals like ```if restart``` with a boolean variable ```restart``` that is user defined.
+
+First we need to define from which file we want to restart, e.g.
+```julia
+restart_file = "restart_000021.h5"
+restart_filename = joinpath("out", restart_file)
+```
+
+Then we load the mesh file:
+```julia
+mesh = load_mesh(restart_filename)
+```
+
+This is then needed for the semidiscretization:
+```julia
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+```
+
+We then define a new time span for the simulation that takes as starting
+time the one form the snapshot:
+```julia
+tspan = (load_time(restart_filename), 2.0)
+```
+
+We now also take the last ```dt```, so that our solver does not need to first find
+one to fulfill the CFL condition:
+```julia
+dt = load_dt(restart_filename)
+```
+
+The ODE that we will pass to the solver is now:
+```julia
+ode = semidiscretize(semi, tspan, restart_filename)
+```
+
+You should now define a [`SaveSolutionCallback`](@ref) similar to the
+[original simulation](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/structured_2d_dgsem/elixir_advection_extended.jl),
+but with ```save_initial_solution=false```, otherwise our initial snapshot will be overwritten.
+If you are using one file for the original simulation and the restart
+you can reuse your [`SaveSolutionCallback`](@ref), but need to set
+```julia
+save_solution.condition.save_initial_solution = false
+```
+
+Before we compute the solution using 
+[OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl)
+we need to set the integrator
+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
+```
+
+Now we can compute the solution:
+```julia
+sol = solve!(integrator)
+```
+
+An example is in `[``examples/structured_2d_dgsem/elixir_advection_restart.jl```](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/structured_2d_dgsem/elixir_advection_restart.jl).
diff --git a/examples/p4est_2d_dgsem/elixir_advection_restart.jl b/examples/p4est_2d_dgsem/elixir_advection_restart.jl
index 1906fb2896..79a35199b8 100644
--- a/examples/p4est_2d_dgsem/elixir_advection_restart.jl
+++ b/examples/p4est_2d_dgsem/elixir_advection_restart.jl
@@ -24,13 +24,23 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                     boundary_conditions=boundary_conditions)
 
 tspan = (load_time(restart_filename), 2.0)
+dt = load_dt(restart_filename)
 ode = semidiscretize(semi, tspan, restart_filename);
 
+# Do not overwrite the initial snapshot written by elixir_advection_extended.jl.
+save_solution.condition.save_initial_solution = false
+
+integrator = init(ode, CarpenterKennedy2N54(williamson_condition=false),
+                  dt=dt, # solve needs some value here but it will be overwritten by the stepsize_callback
+                  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
+
 
 ###############################################################################
 # run the simulation
 
-sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
-            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
-            save_everystep=false, callback=callbacks);
+sol = solve!(integrator)
 summary_callback() # print the timer summary
diff --git a/examples/p4est_3d_dgsem/elixir_advection_restart.jl b/examples/p4est_3d_dgsem/elixir_advection_restart.jl
index 71b37e9f39..b27eaab62e 100644
--- a/examples/p4est_3d_dgsem/elixir_advection_restart.jl
+++ b/examples/p4est_3d_dgsem/elixir_advection_restart.jl
@@ -21,13 +21,23 @@ mesh = load_mesh(restart_filename)
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, solver)
 
 tspan = (load_time(restart_filename), 2.0)
+dt = load_dt(restart_filename)
 ode = semidiscretize(semi, tspan, restart_filename);
 
+# Do not overwrite the initial snapshot written by elixir_advection_extended.jl.
+save_solution.condition.save_initial_solution = false
+
+integrator = init(ode, CarpenterKennedy2N54(williamson_condition=false),
+                  dt=dt, # solve needs some value here but it will be overwritten by the stepsize_callback
+                  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
+
 
 ###############################################################################
 # run the simulation
 
-sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
-            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
-            save_everystep=false, callback=callbacks);
+sol = solve!(integrator)
 summary_callback() # print the timer summary
diff --git a/examples/structured_2d_dgsem/elixir_advection_restart.jl b/examples/structured_2d_dgsem/elixir_advection_restart.jl
index 2c2a0ef8f5..98c44fac71 100644
--- a/examples/structured_2d_dgsem/elixir_advection_restart.jl
+++ b/examples/structured_2d_dgsem/elixir_advection_restart.jl
@@ -23,13 +23,22 @@ mesh = load_mesh(restart_filename)
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
 
 tspan = (load_time(restart_filename), 2.0)
+dt = load_dt(restart_filename)
 ode = semidiscretize(semi, tspan, restart_filename);
 
+# Do not overwrite the initial snapshot written by elixir_advection_extended.jl.
+save_solution.condition.save_initial_solution = false
+
+integrator = init(ode, CarpenterKennedy2N54(williamson_condition=false),
+                 dt=dt, # solve needs some value here but it will be overwritten by the stepsize_callback
+                 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
 
 ###############################################################################
 # run the simulation
 
-sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
-            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
-            save_everystep=false, callback=callbacks);
+sol = solve!(integrator)
 summary_callback() # print the timer summary
diff --git a/examples/structured_3d_dgsem/elixir_advection_restart.jl b/examples/structured_3d_dgsem/elixir_advection_restart.jl
index 39e1a67516..39d28848c7 100644
--- a/examples/structured_3d_dgsem/elixir_advection_restart.jl
+++ b/examples/structured_3d_dgsem/elixir_advection_restart.jl
@@ -21,13 +21,23 @@ mesh = load_mesh(restart_filename)
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, solver)
 
 tspan = (load_time(restart_filename), 2.0)
+dt = load_dt(restart_filename)
 ode = semidiscretize(semi, tspan, restart_filename);
 
+# Do not overwrite the initial snapshot written by elixir_advection_extended.jl.
+save_solution.condition.save_initial_solution = false
+
+integrator = init(ode, CarpenterKennedy2N54(williamson_condition=false),
+                  dt=dt, # solve needs some value here but it will be overwritten by the stepsize_callback
+                  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
+
 
 ###############################################################################
 # run the simulation
 
-sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
-            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
-            save_everystep=false, callback=callbacks);
+sol = solve!(integrator)
 summary_callback() # print the timer summary
diff --git a/examples/tree_2d_dgsem/elixir_advection_restart.jl b/examples/tree_2d_dgsem/elixir_advection_restart.jl
index 2cb45c0b47..4ceb593257 100644
--- a/examples/tree_2d_dgsem/elixir_advection_restart.jl
+++ b/examples/tree_2d_dgsem/elixir_advection_restart.jl
@@ -20,13 +20,23 @@ mesh = load_mesh(restart_filename)
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
 
 tspan = (load_time(restart_filename), 2.0)
+dt = load_dt(restart_filename)
 ode = semidiscretize(semi, tspan, restart_filename);
 
+# Do not overwrite the initial snapshot written by elixir_advection_extended.jl.
+save_solution.condition.save_initial_solution = false
+
+integrator = init(ode, CarpenterKennedy2N54(williamson_condition=false),
+                  dt=dt, # solve needs some value here but it will be overwritten by the stepsize_callback
+                  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
 
 ###############################################################################
 # run the simulation
 
-sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
-            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
-            save_everystep=false, callback=callbacks);
+sol = solve!(integrator)
+
 summary_callback() # print the timer summary
diff --git a/examples/tree_3d_dgsem/elixir_advection_restart.jl b/examples/tree_3d_dgsem/elixir_advection_restart.jl
index 83bf4418b9..3061f16587 100644
--- a/examples/tree_3d_dgsem/elixir_advection_restart.jl
+++ b/examples/tree_3d_dgsem/elixir_advection_restart.jl
@@ -20,13 +20,23 @@ mesh = load_mesh(restart_filename)
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
 
 tspan = (load_time(restart_filename), 2.0)
+dt = load_dt(restart_filename)
 ode = semidiscretize(semi, tspan, restart_filename);
 
+# Do not overwrite the initial snapshot written by elixir_advection_extended.jl.
+save_solution.condition.save_initial_solution = false
+
+integrator = init(ode, CarpenterKennedy2N54(williamson_condition=false),
+                  dt=dt, # solve needs some value here but it will be overwritten by the stepsize_callback
+                  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
+
 
 ###############################################################################
 # run the simulation
 
-sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
-            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
-            save_everystep=false, callback=callbacks);
+sol = solve!(integrator)
 summary_callback() # print the timer summary
diff --git a/examples/unstructured_2d_dgsem/elixir_euler_restart.jl b/examples/unstructured_2d_dgsem/elixir_euler_restart.jl
index 2ac6765202..b85cc2c6d7 100644
--- a/examples/unstructured_2d_dgsem/elixir_euler_restart.jl
+++ b/examples/unstructured_2d_dgsem/elixir_euler_restart.jl
@@ -22,14 +22,24 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                     boundary_conditions=boundary_conditions)
 
 tspan = (load_time(restart_filename), 1.0)
+dt = load_dt(restart_filename)
 ode = semidiscretize(semi, tspan, restart_filename);
 
+# Do not overwrite the initial snapshot written by elixir_advection_extended.jl.
+save_solution.condition.save_initial_solution = false
+
+integrator = init(ode, CarpenterKennedy2N54(williamson_condition=false),
+                  dt=dt, # solve needs some value here but it will be overwritten by the stepsize_callback
+                  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
+
 
 ###############################################################################
 # run the simulation
 
-sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
-            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
-            save_everystep=false, callback=callbacks);
+sol = solve!(integrator)
 summary_callback() # print the timer summary
 
diff --git a/src/Trixi.jl b/src/Trixi.jl
index 66878f4b45..6fc62f5052 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -241,7 +241,7 @@ export SummaryCallback, SteadyStateCallback, AnalysisCallback, AliveCallback,
        GlmSpeedCallback, LBMCollisionCallback, EulerAcousticsCouplingCallback,
        TrivialCallback, AnalysisCallbackCoupled
 
-export load_mesh, load_time
+export load_mesh, load_time, 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 e23f58f26e..f567a5c7fd 100644
--- a/src/callbacks_step/save_restart.jl
+++ b/src/callbacks_step/save_restart.jl
@@ -130,6 +130,28 @@ function load_time(restart_file::AbstractString)
     end
 end
 
+"""
+    load_timestep(restart_file::AbstractString)
+
+Load the time step number (`iter` in OrdinaryDiffEq.jl) saved in a `restart_file`.
+"""
+function load_timestep(restart_file::AbstractString)
+    h5open(restart_file, "r") do file
+        read(attributes(file)["timestep"])
+    end
+end
+
+"""
+    load_dt(restart_file::AbstractString)
+
+Load the time step size (`dt` in OrdinaryDiffEq.jl) saved in a `restart_file`.
+"""
+function load_dt(restart_file::AbstractString)
+    h5open(restart_file, "r") do file
+        read(attributes(file)["dt"])
+    end
+end
+
 function load_restart_file(semi::AbstractSemidiscretization, restart_file)
     load_restart_file(mesh_equations_solver_cache(semi)..., restart_file)
 end
diff --git a/src/meshes/mesh_io.jl b/src/meshes/mesh_io.jl
index ede85d8010..da67fe23e0 100644
--- a/src/meshes/mesh_io.jl
+++ b/src/meshes/mesh_io.jl
@@ -286,6 +286,7 @@ function load_mesh_serial(mesh_file::AbstractString; n_cells_max, RealT)
 
         mesh = StructuredMesh(size, mapping; RealT = RealT, unsaved_changes = false,
                               mapping_as_string = mapping_as_string)
+        mesh.current_filename = mesh_file
     elseif mesh_type == "UnstructuredMesh2D"
         mesh_filename, periodicity_ = h5open(mesh_file, "r") do file
             return read(attributes(file)["mesh_filename"]),