diff --git a/examples/postprocessing/restart_poiseuille_flow_2d.jl b/examples/postprocessing/restart_poiseuille_flow_2d.jl
new file mode 100644
index 0000000000..0332a1c31d
--- /dev/null
+++ b/examples/postprocessing/restart_poiseuille_flow_2d.jl
@@ -0,0 +1,31 @@
+# ==========================================================================================
+# Restart Example: Poiseuille Flow 2D
+#
+# This example demonstrates how to restart a simulation.
+# We first run a simulation of 2D Poiseuille flow up to t=0.3s, then restart from the
+# saved state and continue the simulation until t=0.6s.
+# ==========================================================================================
+using TrixiParticles
+
+trixi_include(@__MODULE__,
+              joinpath(examples_dir(), "fluid", "poiseuille_flow_2d.jl"),
+              tspan=(0.0, 0.3), sound_speed_factor=10, particle_spacing=4e-5)
+
+# Get latest iteration
+iter = saving_callback.condition.index[] - 1
+
+restart_file_fluid = joinpath("out", "fluid_1_$iter.vtu")
+restart_file_open_boundary = joinpath("out", "open_boundary_1_$iter.vtu")
+restart_file_boundary = joinpath("out", "boundary_1_$iter.vtu")
+
+ode_restart = semidiscretize(semi, (0.3, 0.6);
+                             restart_with=(restart_file_fluid,
+                                           restart_file_open_boundary,
+                                           restart_file_boundary))
+
+saving_callback = SolutionSavingCallback(dt=0.02, prefix="restart")
+
+callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback())
+
+sol_restart = solve(ode_restart, RDPK3SpFSAL35(), abstol=1e-5, reltol=1e-3, dtmax=1e-2,
+                    save_everystep=false, callback=callbacks)
diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl
index 215a2b3bf1..9cfd4647d7 100644
--- a/src/TrixiParticles.jl
+++ b/src/TrixiParticles.jl
@@ -61,6 +61,7 @@ include("general/semidiscretization.jl")
 include("general/gpu.jl")
 include("preprocessing/preprocessing.jl")
 include("io/io.jl")
+include("general/restart.jl")
 include("visualization/recipes_plots.jl")
 
 export Semidiscretization, semidiscretize, restart_with!
diff --git a/src/general/restart.jl b/src/general/restart.jl
new file mode 100644
index 0000000000..0f0f80dcc6
--- /dev/null
+++ b/src/general/restart.jl
@@ -0,0 +1,119 @@
+function set_initial_conditions!(v0_ode, u0_ode, semi, restart_with::Tuple{Vararg{String}})
+    # Check number of systems
+    if length(semi.systems) != length(restart_with)
+        throw(ArgumentError("Number of systems in `semi` does not match number of restart files provided " *
+                            "in `restart_with`"))
+    end
+
+    # Check that systems match
+    expected_system_names = system_names(semi.systems)
+    foreach_system(semi) do system
+        system_index = system_indices(system, semi)
+        filename = restart_with[system_index]
+        expected_system_name = expected_system_names[system_index]
+        if !occursin(expected_system_name, basename(splitext(filename)[1]))
+            throw(ArgumentError("Filename '$filename' for system $system_index does not contain expected name '$expected_system_name'. " *
+                                "Expected a VTK file for system of type $(nameof(typeof(system)))."))
+        end
+    end
+
+    # Set initial conditions
+    foreach_noalloc(semi.systems, restart_with) do (system, restart_file)
+        v0_system = wrap_v(v0_ode, system, semi)
+        u0_system = wrap_u(u0_ode, system, semi)
+
+        restart_data = vtk2trixi(restart_file; create_initial_condition=false)
+        v_restart = restart_v(system, restart_data)
+        u_restart = restart_u(system, restart_data)
+
+        v0_system .= Adapt.adapt(semi.parallelization_backend, v_restart)
+        u0_system .= Adapt.adapt(semi.parallelization_backend, u_restart)
+
+        restore_previous_state!(system, restart_file)
+    end
+end
+
+function time_span(tspan, restart_with::Tuple{Vararg{String}})
+    restart_data = vtk2trixi(first(restart_with))
+    t_restart = convert(eltype(tspan), restart_data.time)
+
+    if !isapprox(tspan[1], t_restart)
+        @info "Adjusting initial time from $(tspan[1]) to restart time $t_restart"
+    end
+
+    return (t_restart, tspan[2])
+end
+
+function write_density_and_pressure!(v_restart, system, density_calculator,
+                                     pressure, density)
+    return v_restart
+end
+
+function write_density_and_pressure!(v_restart, system,
+                                     density_calculator::ContinuityDensity,
+                                     pressure, density)
+    v_restart[size(v_restart, 1), :] = density
+
+    return v_restart
+end
+
+function write_density_and_pressure!(v_restart, system::EntropicallyDampedSPHSystem,
+                                     density_calculator::ContinuityDensity,
+                                     pressure, density)
+    v_restart[size(v_restart, 1), :] = density
+    v_restart[size(v_restart, 1) - 1, :] = pressure
+
+    return v_restart
+end
+
+restore_previous_state!(system, restart_file) = system
+
+function initialize_neighborhood_searches!(semi, u0_ode,
+                                           restart_with::Tuple{Vararg{String}})
+    foreach_system(semi) do system
+        foreach_system(semi) do neighbor
+            initialize_neighborhood_search!(semi, system, neighbor, u0_ode)
+        end
+    end
+
+    return semi
+end
+
+function initialize_neighborhood_search!(semi, system, neighbor, u0_ode)
+    # TODO Initialize after adapting to the GPU.
+    # Currently, this cannot use `semi.parallelization_backend`
+    # because data is still on the CPU.
+    PointNeighbors.initialize!(get_neighborhood_search(system, neighbor, semi),
+                               initial_restart_coordinates(system, u0_ode, semi),
+                               initial_restart_coordinates(neighbor, u0_ode, semi),
+                               eachindex_y=each_active_particle(neighbor),
+                               parallelization_backend=PolyesterBackend())
+    return semi
+end
+
+function initialize_neighborhood_search!(semi, system::TotalLagrangianSPHSystem,
+                                         neighbor::TotalLagrangianSPHSystem, u0_ode)
+    # For TLSPH, the self-interaction NHS is already initialized in the system constructor
+    return semi
+end
+
+function initial_restart_coordinates(system, u0_ode, semi)
+    # Transfer to CPU if data is on the GPU. Do nothing if already on CPU.
+    return transfer2cpu(wrap_u(u0_ode, system, semi))
+end
+
+function initial_restart_coordinates(system::Union{WallBoundarySystem,
+                                                   AbstractStructureSystem}, u0_ode, semi)
+    return initial_coordinates(system)
+end
+
+function initialize!(semi::Semidiscretization, restart_with::Tuple{Vararg{String}})
+    foreach_system(semi) do system
+        # Initialize this system
+        initialize_restart!(system, semi)
+    end
+
+    return semi
+end
+
+initialize_restart!(system, semi) = initialize!(system, semi)
diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl
index 7b91c0baf9..9d94107cff 100644
--- a/src/general/semidiscretization.jl
+++ b/src/general/semidiscretization.jl
@@ -184,7 +184,7 @@ end
 end
 
 """
-    semidiscretize(semi, tspan; reset_threads=true)
+    semidiscretize(semi, tspan; reset_threads=true, restart_with=nothing)
 
 Create an `ODEProblem` from the semidiscretization with the specified `tspan`.
 
@@ -193,6 +193,10 @@ Create an `ODEProblem` from the semidiscretization with the specified `tspan`.
 - `tspan`: The time span over which the simulation will be run.
 
 # Keywords
+- `restart_with`: Can be used to restart the simulation from VTK solution files (see [`SolutionSavingCallback`](@ref)).
+  This has to be a tuple of filenames, one for each system in the [`Semidiscretization`](@ref).
+  The order of the filenames has to match the order of the systems in the [`Semidiscretization`](@ref).
+  If no restart is desired, use `nothing` (default).
 - `reset_threads`: A boolean flag to reset Polyester.jl threads before the simulation (default: `true`).
   After an error within a threaded loop, threading might be disabled. Resetting the threads before the simulation
   ensures that threading is enabled again for the simulation.
@@ -219,7 +223,7 @@ timespan: (0.0, 1.0)
 u0: ([...], [...]) *this line is ignored by filter*
 ```
 """
-function semidiscretize(semi, tspan; reset_threads=true)
+function semidiscretize(semi, tspan; reset_threads=true, restart_with=nothing)
     (; systems) = semi
 
     # Check that all systems have the same eltype
@@ -267,14 +271,11 @@ function semidiscretize(semi, tspan; reset_threads=true)
     end
 
     # Set initial condition
-    foreach_system_wrapped(semi, v0_ode, u0_ode) do system, v0_system, u0_system
-        write_u0!(u0_system, system)
-        write_v0!(v0_system, system)
-    end
+    set_initial_conditions!(v0_ode, u0_ode, semi, restart_with)
 
     # TODO initialize after adapting to the GPU.
     # Requires https://github.com/trixi-framework/PointNeighbors.jl/pull/86.
-    initialize_neighborhood_searches!(semi)
+    initialize_neighborhood_searches!(semi, u0_ode, restart_with)
 
     if semi.parallelization_backend isa KernelAbstractions.GPU
         # Convert all arrays in the systems to the correct array type.
@@ -304,15 +305,31 @@ function semidiscretize(semi, tspan; reset_threads=true)
     end
 
     # Initialize all particle systems
-    foreach_system(semi_new) do system
-        # Initialize this system
-        initialize!(system, semi_new)
-    end
+    initialize!(semi_new, restart_with)
 
     # Reset callback flag that will be set by the `UpdateCallback`
     semi_new.update_callback_used[] = false
 
-    return DynamicalODEProblem(kick!, drift!, v0_ode, u0_ode, tspan, semi_new)
+    return DynamicalODEProblem(kick!, drift!, v0_ode, u0_ode,
+                               time_span(tspan, restart_with), semi_new)
+end
+
+function set_initial_conditions!(v0_ode, u0_ode, semi, restart_with::Nothing)
+    foreach_system_wrapped(semi, v0_ode, u0_ode) do system, v0_system, u0_system
+        write_u0!(u0_system, system)
+        write_v0!(v0_system, system)
+    end
+end
+
+time_span(tspan, restart_with::Nothing) = tspan
+
+function initialize!(semi::Semidiscretization, restart_with::Nothing)
+    foreach_system(semi) do system
+        # Initialize this system
+        initialize!(system, semi)
+    end
+
+    return semi
 end
 
 """
@@ -348,6 +365,39 @@ function restart_with!(semi, sol; reset_threads=true)
     return semi
 end
 
+function initialize_neighborhood_searches!(semi, u0_ode, restart_with::Nothing)
+    initialize_neighborhood_searches!(semi)
+end
+
+function initialize_neighborhood_searches!(semi)
+    foreach_system(semi) do system
+        foreach_system(semi) do neighbor
+            initialize_neighborhood_search!(semi, system, neighbor)
+        end
+    end
+
+    return semi
+end
+
+function initialize_neighborhood_search!(semi, system, neighbor)
+    # TODO Initialize after adapting to the GPU.
+    # Currently, this cannot use `semi.parallelization_backend`
+    # because data is still on the CPU.
+    PointNeighbors.initialize!(get_neighborhood_search(system, neighbor, semi),
+                               initial_coordinates(system),
+                               initial_coordinates(neighbor),
+                               eachindex_y=each_active_particle(neighbor),
+                               parallelization_backend=PolyesterBackend())
+
+    return semi
+end
+
+function initialize_neighborhood_search!(semi, system::TotalLagrangianSPHSystem,
+                                         neighbor::TotalLagrangianSPHSystem)
+    # For TLSPH, the self-interaction NHS is already initialized in the system constructor
+    return semi
+end
+
 # We have to pass `system` here for type stability,
 # since the type of `system` determines the return type.
 @inline function wrap_v(v_ode, system, semi)
diff --git a/src/io/write_vtk.jl b/src/io/write_vtk.jl
index 6747af2763..8e04e81b85 100644
--- a/src/io/write_vtk.jl
+++ b/src/io/write_vtk.jl
@@ -458,6 +458,16 @@ function write2vtk!(vtk, v, u, t, system::OpenBoundarySystem)
                       for particle in eachparticle(system)]
     vtk["pressure"] = [current_pressure(v, system, particle)
                        for particle in eachparticle(system)]
+    vtk["zone_id"] = [system.boundary_zone_indices[particle]
+                      for particle in eachparticle(system)]
+
+    if any(pm -> isa(pm, AbstractPressureModel), system.cache.pressure_reference_values)
+        for (i, pressure_model) in enumerate(system.cache.pressure_reference_values)
+            if pressure_model isa AbstractPressureModel
+                vtk["boundary_zone_pressure_$i"] = system.cache.pressure_reference_values[i].pressure[]
+            end
+        end
+    end
 
     if system.calculate_flow_rate
         Q_total = zero(eltype(system))
diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl
index 3e4754d4ad..89c78b5024 100644
--- a/src/schemes/boundary/open_boundary/system.jl
+++ b/src/schemes/boundary/open_boundary/system.jl
@@ -133,6 +133,9 @@ function initialize!(system::OpenBoundarySystem, semi)
     return system
 end
 
+# Skip during restart, as boundary zone indices are updated in `restore_previous_state!`
+initialize_restart!(system::OpenBoundarySystem, semi) = system
+
 function create_cache_open_boundary(boundary_model, fluid_system, initial_condition,
                                     density_diffusion, calculate_flow_rate, boundary_zones)
     reference_values = map(bz -> bz.reference_values, boundary_zones)
@@ -251,6 +254,8 @@ end
     return system.shifting_technique
 end
 
+@inline density_calculator(system::OpenBoundarySystem) = nothing
+
 system_sound_speed(system::OpenBoundarySystem) = system_sound_speed(system.fluid_system)
 
 @inline hydrodynamic_mass(system::OpenBoundarySystem, particle) = system.mass[particle]
@@ -680,6 +685,66 @@ function interpolate_velocity!(system::OpenBoundarySystem, boundary_zone,
     return system
 end
 
+function restart_u(system::OpenBoundarySystem, data)
+    coords_total = zeros(coordinates_eltype(system), u_nvariables(system),
+                         n_integrated_particles(system))
+    coords_total .= coordinates_eltype(system)(1e16)
+
+    coords_active = data.coordinates
+    for particle in axes(coords_active, 2)
+        for dim in 1:ndims(system)
+            coords_total[dim, particle] = coords_active[dim, particle]
+        end
+    end
+
+    system.buffer.active_particle .= false
+    system.buffer.active_particle[1:size(coords_active, 2)] .= true
+
+    update_system_buffer!(system.buffer)
+
+    return coords_total
+end
+
+function restart_v(system::OpenBoundarySystem, data)
+    v_total = zeros(eltype(system), v_nvariables(system),
+                    n_integrated_particles(system))
+
+    v_active = zeros(eltype(system), v_nvariables(system), size(data.velocity, 2))
+
+    v_active[1:ndims(system), :] = data.velocity
+    write_density_and_pressure!(v_active, system.fluid_system,
+                                density_calculator(system), data.pressure, data.density)
+
+    for particle in axes(v_active, 2)
+        for i in axes(v_active, 1)
+            v_total[i, particle] = v_active[i, particle]
+        end
+    end
+
+    return v_total
+end
+
+function restore_previous_state!(system::OpenBoundarySystem, file)
+    # We cannot simply use `update_boundary_zone_indices!` because rounding errors during file I/O
+    # may result in particles being located outside their intended boundary zone, even though they
+    # were written as active particles.
+    set_zero!(system.boundary_zone_indices)
+
+    values = vtk2trixi(file; create_initial_condition=false)
+    system.boundary_zone_indices[each_integrated_particle(system)] .= values.zone_id
+
+    if any(pm -> isa(pm, AbstractPressureModel), system.cache.pressure_reference_values)
+        for (i, pressure_model) in enumerate(system.cache.pressure_reference_values)
+            if pressure_model isa AbstractPressureModel
+                pressure_model.pressure[] = values[Symbol(:boundary_zone_pressure_, i)]
+                pressure_model.flow_rate[] = values[Symbol(:Q_, i)]
+            end
+        end
+    end
+
+    return system
+end
+
 # Open-boundary interpolation should reconstruct the surrounding fluid-like velocity field.
 # Therefore, only actual fluid systems and other open-boundary particles contribute;
 # rigid bodies, walls, and other non-fluid systems are intentionally excluded.
diff --git a/src/schemes/boundary/wall_boundary/dummy_particles.jl b/src/schemes/boundary/wall_boundary/dummy_particles.jl
index 6453fc0c2a..895eb2e0a7 100644
--- a/src/schemes/boundary/wall_boundary/dummy_particles.jl
+++ b/src/schemes/boundary/wall_boundary/dummy_particles.jl
@@ -328,6 +328,10 @@ function Base.show(io::IO, model::BoundaryModelDummyParticles)
     print(io, ")")
 end
 
+@inline function density_calculator(model::BoundaryModelDummyParticles)
+    return model.density_calculator
+end
+
 # Set the initial pressure to zero for visualization
 function initial_boundary_pressure(initial_density, density_calculator, _)
     return zero(initial_density)
diff --git a/src/schemes/boundary/wall_boundary/monaghan_kajtar.jl b/src/schemes/boundary/wall_boundary/monaghan_kajtar.jl
index 1df9a38cf9..f23f69a1fe 100644
--- a/src/schemes/boundary/wall_boundary/monaghan_kajtar.jl
+++ b/src/schemes/boundary/wall_boundary/monaghan_kajtar.jl
@@ -43,6 +43,10 @@ function Base.show(io::IO, model::BoundaryModelMonaghanKajtar)
     print(io, ")")
 end
 
+@inline function density_calculator(model::BoundaryModelMonaghanKajtar)
+    return nothing
+end
+
 @inline function pressure_acceleration(particle_system,
                                        neighbor_system::Union{WallBoundarySystem{<:BoundaryModelMonaghanKajtar},
                                                               TotalLagrangianSPHSystem{<:BoundaryModelMonaghanKajtar}},
diff --git a/src/schemes/boundary/wall_boundary/system.jl b/src/schemes/boundary/wall_boundary/system.jl
index 3a546f4013..115b6202d3 100644
--- a/src/schemes/boundary/wall_boundary/system.jl
+++ b/src/schemes/boundary/wall_boundary/system.jl
@@ -264,6 +264,23 @@ function restart_with!(system::WallBoundarySystem{<:BoundaryModelDummyParticles{
     return system
 end
 
+function restart_u(system::WallBoundarySystem, data)
+    if n_integrated_particles(system) > 0
+        throw(ArgumentError("`WallBoundarySystem` does not support integrated particle coordinates"))
+    end
+
+    return zeros(eltype(system), u_nvariables(system), n_integrated_particles(system))
+end
+
+function restart_v(system::WallBoundarySystem, data)
+    v_ode = zeros(eltype(system), v_nvariables(system), n_integrated_particles(system))
+
+    write_density_and_pressure!(v_ode, system, density_calculator(system),
+                                data.pressure, data.density)
+
+    return v_ode
+end
+
 # To incorporate the effect at boundaries in the viscosity term of the RHS, the neighbor
 # viscosity model has to be used.
 @inline function viscosity_model(system::WallBoundarySystem,
@@ -315,6 +332,10 @@ function system_correction(system::WallBoundarySystem{<:BoundaryModelDummyPartic
     return system.boundary_model.correction
 end
 
+@inline function density_calculator(system::WallBoundarySystem)
+    return density_calculator(system.boundary_model)
+end
+
 function system_data(system::WallBoundarySystem, dv_ode, du_ode, v_ode, u_ode, semi)
     dv = [current_acceleration(system, particle) for particle in eachparticle(system)]
     v = wrap_v(v_ode, system, semi)
diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl
index eb6338f1cc..be9c7ac6dd 100644
--- a/src/schemes/fluid/fluid.jl
+++ b/src/schemes/fluid/fluid.jl
@@ -253,6 +253,47 @@ end
     return nothing
 end
 
+function restart_u(system::AbstractFluidSystem, data)
+    coords_total = zeros(coordinates_eltype(system), u_nvariables(system),
+                         n_integrated_particles(system))
+    coords_total .= coordinates_eltype(system)(1e16)
+
+    coords_active = data.coordinates
+
+    for particle in axes(coords_active, 2)
+        for dim in 1:ndims(system)
+            coords_total[dim, particle] = coords_active[dim, particle]
+        end
+    end
+
+    if !isnothing(buffer(system))
+        system.buffer.active_particle .= false
+        system.buffer.active_particle[1:size(coords_active, 2)] .= true
+    end
+
+    update_system_buffer!(system.buffer)
+
+    return coords_total
+end
+
+function restart_v(system::AbstractFluidSystem, data)
+    velocity_total = zeros(eltype(system), v_nvariables(system),
+                           n_integrated_particles(system))
+    velocity_active = zeros(eltype(system), v_nvariables(system), size(data.velocity, 2))
+
+    velocity_active[1:ndims(system), :] = data.velocity
+    write_density_and_pressure!(velocity_active, system, density_calculator(system),
+                                data.pressure, data.density)
+
+    for particle in axes(velocity_active, 2)
+        for i in axes(velocity_active, 1)
+            velocity_total[i, particle] = velocity_active[i, particle]
+        end
+    end
+
+    return velocity_total
+end
+
 function check_configuration(fluid_system::AbstractFluidSystem, systems, nhs)
     if !(fluid_system isa ParticlePackingSystem) && !isnothing(fluid_system.surface_tension)
         foreach_system(systems) do neighbor
diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl
index f12004401a..cf5b1942c2 100644
--- a/src/schemes/structure/total_lagrangian_sph/system.jl
+++ b/src/schemes/structure/total_lagrangian_sph/system.jl
@@ -602,6 +602,16 @@ function restart_with!(system::TotalLagrangianSPHSystem, v, u)
     restart_with!(system, system.boundary_model, v, u)
 end
 
+function restart_u(system::AbstractStructureSystem, data)
+    # The integration array requires only the particles that need to be integrated
+    return data.coordinates[:, each_integrated_particle(system)]
+end
+
+function restart_v(system::AbstractStructureSystem, data)
+    # The integration array requires only the particles that need to be integrated
+    return data.velocity[:, each_integrated_particle(system)]
+end
+
 # An explanation of these equation can be found in
 # J. Lubliner, 2008. Plasticity theory.
 # See here below Equation 5.3.21 for the equation for the equivalent stress.
diff --git a/test/examples/gpu.jl b/test/examples/gpu.jl
index 8f84611647..12a192dc2c 100644
--- a/test/examples/gpu.jl
+++ b/test/examples/gpu.jl
@@ -732,4 +732,58 @@ end
             end
         end
     end
+
+    @testset verbose=true "Restart" begin
+        trixi_include_changeprecision(Float32, @__MODULE__,
+                                      joinpath(examples_dir(), "fluid",
+                                               "poiseuille_flow_2d.jl"),
+                                      tspan=(0.0f0, 0.6f0), sound_speed_factor=10,
+                                      particle_spacing=4.0f-5,
+                                      coordinates_eltype=Float32,
+                                      parallelization_backend=Main.parallelization_backend)
+
+        # Since this is an open boundary simulation, the number of active particles may
+        # differ. The results must be interpolated to enable comparison with the restart
+        # simulation. The fluid domain starts at `x = 10 * particle_spacing`.
+        n_interpolation_points = 10
+        start_point = [0.0f0 + 10 * particle_spacing, wall_distance / 2]
+        end_point = [flow_length - 10 * particle_spacing, wall_distance / 2]
+        result_full = interpolate_line(start_point, end_point, n_interpolation_points,
+                                       sol.prob.p, sol.prob.p.systems[1], sol,
+                                       cut_off_bnd=false)
+
+        # Run half simulation and safe checkpoint
+        trixi_include_changeprecision(Float32, @__MODULE__,
+                                      joinpath(examples_dir(), "fluid",
+                                               "poiseuille_flow_2d.jl"),
+                                      tspan=(0.0f0, 0.3f0), sound_speed_factor=10,
+                                      particle_spacing=4.0f-5,
+                                      coordinates_eltype=Float32,
+                                      parallelization_backend=Main.parallelization_backend)
+
+        # Transfer `semi` back to CPU
+        semi_new = TrixiParticles.Adapt.adapt(Array, sol.prob.p)
+
+        iter = round(Int, 0.3 / 0.02)
+        fluid_restart = joinpath("out", "fluid_1_$iter.vtu")
+        open_boundary_restart = joinpath("out", "open_boundary_1_$iter.vtu")
+        boundary_restart = joinpath("out", "boundary_1_$iter.vtu")
+
+        ode_restart = semidiscretize(semi_new, (0.3f0, 0.6f0);
+                                     restart_with=(fluid_restart, open_boundary_restart,
+                                                   boundary_restart))
+
+        sol_restart = solve(ode_restart, RDPK3SpFSAL35(), abstol=1.0f-5,
+                            reltol=1.0f-3, dtmax=1.0f-2, save_everystep=false,
+                            callback=UpdateCallback())
+
+        result_restart = interpolate_line(start_point, end_point,
+                                          n_interpolation_points, sol_restart.prob.p,
+                                          sol_restart.prob.p.systems[1],
+                                          sol_restart, cut_off_bnd=false)
+
+        @test isapprox(result_full.velocity, result_restart.velocity, rtol=8e-3)
+        @test isapprox(result_full.density, result_restart.density, rtol=8e-4)
+        @test isapprox(result_full.pressure, result_restart.pressure, rtol=8e-2)
+    end
 end
diff --git a/test/general/general.jl b/test/general/general.jl
index 2eb24d50a0..acb07de8b9 100644
--- a/test/general/general.jl
+++ b/test/general/general.jl
@@ -6,4 +6,5 @@ include("interpolation.jl")
 include("buffer.jl")
 include("util.jl")
 include("custom_quantities.jl")
+include("restart.jl")
 include("neighborhood_search.jl")
diff --git a/test/general/restart.jl b/test/general/restart.jl
new file mode 100644
index 0000000000..386038e6b1
--- /dev/null
+++ b/test/general/restart.jl
@@ -0,0 +1,73 @@
+@trixi_testset "Restart" begin
+    @trixi_testset "Half-Simulation Restart" begin
+        # Run full simulation
+        trixi_include(@__MODULE__,
+                      joinpath(examples_dir(), "fluid", "poiseuille_flow_2d.jl"),
+                      tspan=(0.0, 0.6), sound_speed_factor=10, particle_spacing=4e-5)
+
+        # Since this is an open boundary simulation, the number of active particles may
+        # differ. The results must be interpolated to enable comparison with the restart
+        # simulation. The fluid domain starts at `x = 10 * particle_spacing`.
+        n_interpolation_points = 10
+        start_point = [0.0 + 10 * particle_spacing, channel_height / 2]
+        end_point = [channel_length - 10 * particle_spacing, channel_height / 2]
+        result_full = interpolate_line(start_point, end_point, n_interpolation_points,
+                                       semi, fluid_system, sol, cut_off_bnd=false)
+
+        # Run half simulation
+        trixi_include(@__MODULE__,
+                      joinpath(examples_dir(), "fluid", "poiseuille_flow_2d.jl"),
+                      tspan=(0.0, 0.3), sound_speed_factor=10, particle_spacing=4e-5)
+
+        iter = round(Int, 0.3 / 0.02)
+        fluid_restart = joinpath("out", "fluid_1_$iter.vtu")
+        open_boundary_restart = joinpath("out", "open_boundary_1_$iter.vtu")
+        boundary_restart = joinpath("out", "boundary_1_$iter.vtu")
+
+        ode_restart = semidiscretize(semi, (0.3, 0.6);
+                                     restart_with=(fluid_restart, open_boundary_restart,
+                                                   boundary_restart))
+
+        sol_restart = solve(ode_restart, RDPK3SpFSAL35(), abstol=1e-5, reltol=1e-3,
+                            dtmax=1e-2, save_everystep=false, callback=UpdateCallback())
+
+        result_restart = interpolate_line(start_point, end_point,
+                                          n_interpolation_points, sol_restart.prob.p,
+                                          sol_restart.prob.p.systems[1],
+                                          sol_restart, cut_off_bnd=false)
+
+        @test isapprox(result_full.velocity, result_restart.velocity, rtol=1e-2)
+        @test isapprox(result_full.density, result_restart.density, rtol=8e-4)
+        @test isapprox(result_full.pressure, result_restart.pressure, rtol=8e-2)
+    end
+
+    @trixi_testset "Restore Previous State" begin
+        R1 = 1.7714
+        R2 = 106.66
+        C = 1.1808e-2
+        pressure_model = RCRWindkesselModel(; peripheral_resistance=R2,
+                                            compliance=C,
+                                            characteristic_resistance=R1)
+        # Run half simulation
+        trixi_include(@__MODULE__,
+                      joinpath(examples_dir(), "fluid", "poiseuille_flow_2d.jl"),
+                      outlet_reference_pressure=pressure_model,
+                      tspan=(0.0, 0.3), sound_speed_factor=10, particle_spacing=4e-5)
+
+        iter = round(Int, 0.3 / 0.02)
+        fluid_restart = joinpath("out", "fluid_1_$iter.vtu")
+        open_boundary_restart = joinpath("out", "open_boundary_1_$iter.vtu")
+        boundary_restart = joinpath("out", "boundary_1_$iter.vtu")
+
+        ode_restart = semidiscretize(semi, (0.3, 0.6);
+                                     restart_with=(fluid_restart, open_boundary_restart,
+                                                   boundary_restart))
+        restart_pressure = ode_restart.p.systems[2].cache.pressure_reference_values[2].pressure[]
+        restart_flow_rate = ode_restart.p.systems[2].cache.pressure_reference_values[2].flow_rate[]
+
+        @test isapprox(restart_pressure,
+                       open_boundary.cache.pressure_reference_values[2].pressure[])
+        @test isapprox(restart_flow_rate,
+                       open_boundary.cache.pressure_reference_values[2].flow_rate[])
+    end
+end