diff --git a/NEWS.md b/NEWS.md
index 857a55ef60..5ade32aba8 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -4,9 +4,16 @@ TrixiParticles.jl follows the interpretation of
 [semantic versioning (semver)](https://julialang.github.io/Pkg.jl/dev/compatibility/#Version-specifier-format-1)
 used in the Julia ecosystem. Notable changes will be documented in this file for human readability.
 
+## Version 0.5.2
+
+### Features
+- Added the ability to restart simulations from VTK solution files generated by `SolutionSavingCallback`.
+  Users can now pass the `restart_with` keyword argument to `semidiscretize` (#1190).
+
 ## Version 0.5.1
 
 ### Features
+
 - Implemented stage-level coupling for split integration (#1049).
 - Added `MechanicalWorkCalculatorCallback` (#940).
 
diff --git a/examples/postprocessing/restart_poiseuille_flow_2d.jl b/examples/postprocessing/restart_poiseuille_flow_2d.jl
new file mode 100644
index 0000000000..4564220805
--- /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.affect!.affect!.latest_saved_iter
+
+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 0bc56370f2..7101be8e49 100644
--- a/src/TrixiParticles.jl
+++ b/src/TrixiParticles.jl
@@ -62,6 +62,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/neighborhood_search.jl b/src/general/neighborhood_search.jl
index bc81241e13..a2b6c22fdf 100644
--- a/src/general/neighborhood_search.jl
+++ b/src/general/neighborhood_search.jl
@@ -261,6 +261,11 @@ end
 end
 
 # === Initialization ===
+
+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
diff --git a/src/general/restart.jl b/src/general/restart.jl
new file mode 100644
index 0000000000..99410c6a90
--- /dev/null
+++ b/src/general/restart.jl
@@ -0,0 +1,139 @@
+# Update `v0_ode` and `u0_ode` with the data from the restart files
+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
+
+# Compute a new `tspan` based on the restart files
+function time_span(tspan, restart_with::Tuple{Vararg{String}})
+    # Read restart times from all files
+    restart_times = [vtk2trixi(file).time for file in restart_with]
+    t_restart = convert(eltype(tspan), first(restart_times))
+
+    # Check if all restart files have the same time
+    if !all(isapprox(t, t_restart) for t in restart_times)
+        throw(ArgumentError("all restart files must start from the same time"))
+    end
+
+    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
+
+function write_density_and_pressure!(v_restart, system::EntropicallyDampedSPHSystem,
+                                     density_calculator::SummationDensity,
+                                     pressure, 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)
+
+
+# TODO
+# Store controller values for an adaptive time stepping scheme
diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl
index 463b3f0a9e..959ed70deb 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,14 @@ 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 can be either `nothing` (default, no restart) or a `Tuple` of filenames,
+  one for each system in the [`Semidiscretization`](@ref).
+  The order of the filenames must match the order of the systems in the [`Semidiscretization`](@ref).
+  Note that `semidiscretize` replaces the initial time (`tspan[1]`) with the timestamp read
+  from the VTK files. If the user-provided `tspan[1]` does not match the restart time,
+  it is adjusted and an info message is logged. If multiple files are provided, their
+  timestamps must match.
 - `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,9 +227,16 @@ 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
 
+    if restart_with isa String
+        restart_with = (restart_with,)
+    elseif !isnothing(restart_with) && !(restart_with isa NTuple{<:Any, String})
+        throw(ArgumentError("`restart_with` must be `nothing`, a string, or a tuple of strings, " *
+                            "got $(typeof(restart_with))"))
+    end
+
     # Check that all systems have the same eltype
     first_system = first(systems)
     if !all(system -> eltype(system) === eltype(first_system), systems)
@@ -267,14 +282,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,10 +316,7 @@ 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
@@ -318,7 +327,27 @@ function semidiscretize(semi, tspan; reset_threads=true)
     # here, since we cannot change `p` from within the callback (only its contents).
     p = @NamedTuple{semi::typeof(semi_new), split_integration_data::Any}((semi_new,
                                                                           nothing))
-    return DynamicalODEProblem(kick!, drift!, v0_ode, u0_ode, tspan, p)
+
+    return DynamicalODEProblem(kick!, drift!, v0_ode, u0_ode,
+                               time_span(tspan, restart_with), p)
+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
 
 """
diff --git a/src/io/read_vtk.jl b/src/io/read_vtk.jl
index 86bae1af77..4ed9dcb5f8 100644
--- a/src/io/read_vtk.jl
+++ b/src/io/read_vtk.jl
@@ -51,8 +51,19 @@ function vtk2trixi(file; element_type=nothing, coordinates_eltype=nothing,
     point_coords = ReadVTK.get_points(vtk_file)
 
     cELTYPE = isnothing(coordinates_eltype) ? eltype(point_coords) : coordinates_eltype
-    ELTYPE = isnothing(element_type) ?
-             eltype(first(ReadVTK.get_data(point_data["pressure"]))) : element_type
+
+    if !isnothing(element_type)
+        ELTYPE = element_type
+    else
+        # Try to get element type from pressure or density (whichever exists)
+        ELTYPE = if "pressure" in keys(point_data)
+            eltype(first(ReadVTK.get_data(point_data["pressure"])))
+        elseif "material_density" in keys(point_data)
+            eltype(first(ReadVTK.get_data(point_data["material_density"])))
+        else
+            error("neither 'pressure' nor 'material_density' field found in VTK file")
+        end
+    end
 
     results = Dict{Symbol, Any}()
 
@@ -87,7 +98,8 @@ function vtk2trixi(file; element_type=nothing, coordinates_eltype=nothing,
                                  results[:particle_spacing]
     results[:coordinates] = coordinates
     results[:time] = "time" in keys(field_data) ?
-                     first(ReadVTK.get_data(field_data["time"])) : zero(ELTYPE)
+                     convert(ELTYPE, first(ReadVTK.get_data(field_data["time"]))) :
+                     zero(ELTYPE)
 
     append!(used_keys, ["index", "ndims"])
     # Load any custom quantities
diff --git a/src/io/write_vtk.jl b/src/io/write_vtk.jl
index 0b19f32337..fa3fb80985 100644
--- a/src/io/write_vtk.jl
+++ b/src/io/write_vtk.jl
@@ -447,6 +447,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 7a56b2cfcc..4348f7df6f 100644
--- a/src/schemes/boundary/open_boundary/system.jl
+++ b/src/schemes/boundary/open_boundary/system.jl
@@ -143,6 +143,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)
@@ -261,6 +264,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]
@@ -705,6 +710,72 @@ 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)
+
+    # Since only active particles are written during saving, the loaded `data` contains
+    # only active particles. These are placed at the beginning of the array, leaving the
+    # inactive buffer particles at the end. Thus, we can safely activate the first N particles.
+    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))
+
+    # Since only active particles are written during saving, the loaded `data` contains
+    # only active particles. These are placed at the beginning of the array, leaving the
+    # inactive buffer particles at the end. Thus, we can safely activate the first N particles.
+    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 3f0addc408..c7dfdfe03c 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 7300cc03d9..2ea6e82cf6 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..750ad72ee4 100644
--- a/src/schemes/fluid/fluid.jl
+++ b/src/schemes/fluid/fluid.jl
@@ -253,6 +253,48 @@ end
     return nothing
 end
 
+function restart_u(system::AbstractFluidSystem, data)
+    inactive_coords = convert(coordinates_eltype(system), 1e16)
+    coords_total = fill(inactive_coords, u_nvariables(system),
+                        n_integrated_particles(system))
+
+    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
+
+    buf = buffer(system)
+    if !isnothing(buf)
+        buf.active_particle .= false
+        buf.active_particle[1:size(coords_active, 2)] .= true
+        update_system_buffer!(buf)
+    end
+
+    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/fluid/implicit_incompressible_sph/system.jl b/src/schemes/fluid/implicit_incompressible_sph/system.jl
index 6e54ce7b66..9b8fedc288 100644
--- a/src/schemes/fluid/implicit_incompressible_sph/system.jl
+++ b/src/schemes/fluid/implicit_incompressible_sph/system.jl
@@ -202,6 +202,8 @@ end
 # TODO: What do we do with the sound speed? This is needed for the viscosity.
 @inline system_sound_speed(system::ImplicitIncompressibleSPHSystem) = system.artificial_sound_speed
 
+@inline density_calculator(system::ImplicitIncompressibleSPHSystem) = nothing
+
 # Calculates the pressure values by solving a linear system with a relaxed Jacobi scheme
 function update_quantities!(system::ImplicitIncompressibleSPHSystem, v, u,
                             v_ode, u_ode, semi, t)
diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl
index fbd4dc0333..a1cd48fb54 100644
--- a/src/schemes/structure/total_lagrangian_sph/system.jl
+++ b/src/schemes/structure/total_lagrangian_sph/system.jl
@@ -634,6 +634,16 @@ function restart_with!(system::TotalLagrangianSPHSystem, v, u)
     restart_with!(system, system.boundary_model, v, u)
 end
 
+function restart_u(system::TotalLagrangianSPHSystem, 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::TotalLagrangianSPHSystem, 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 e3415c0e9d..ca8ab992de 100644
--- a/test/examples/gpu.jl
+++ b/test/examples/gpu.jl
@@ -896,4 +896,64 @@ 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, sol=nothing,
+                                      coordinates_eltype=Float32,
+                                      parallelization_backend=Main.parallelization_backend)
+
+        sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+                    dt=1.0f0, save_everystep=false,
+                    callback=CallbackSet(callbacks, StepsizeCallback(cfl=1.5f0)))
+
+        # 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, channel_height / 2]
+        end_point = [channel_length - 10 * particle_spacing, channel_height / 2]
+        result_full = interpolate_line(start_point, end_point, n_interpolation_points,
+                                       sol.prob.p.semi, sol.prob.p.semi.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, sol=nothing,
+                                      coordinates_eltype=Float32,
+                                      parallelization_backend=Main.parallelization_backend)
+
+        sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+                    dt=1.0f0, save_everystep=false,
+                    callback=CallbackSet(callbacks, StepsizeCallback(cfl=1.5f0)))
+
+        iter = saving_callback.affect!.affect!.latest_saved_iter
+        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.3f0, 0.6f0);
+                                     restart_with=(fluid_restart, open_boundary_restart,
+                                                   boundary_restart))
+
+        sol_restart = solve(ode_restart, CarpenterKennedy2N54(williamson_condition=false),
+                            dt=1.0f0, save_everystep=false,
+                            callback=CallbackSet(UpdateCallback(),
+                                                 StepsizeCallback(cfl=1.5f0)))
+
+        result_restart = interpolate_line(start_point, end_point,
+                                          n_interpolation_points, sol_restart.prob.p.semi,
+                                          sol_restart.prob.p.semi.systems[1],
+                                          sol_restart, cut_off_bnd=false)
+
+        @test isapprox(result_full.velocity, result_restart.velocity, rtol=7.0f-5)
+        @test isapprox(result_full.density, result_restart.density, rtol=5.0f-6)
+        @test isapprox(result_full.pressure, result_restart.pressure, rtol=5.0f-4)
+    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..d70f57b69e
--- /dev/null
+++ b/test/general/restart.jl
@@ -0,0 +1,113 @@
+@trixi_testset "Restart" begin
+    @trixi_testset "Poiseuille Flow 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,
+                      info_callback=nothing, sol=nothing)
+
+        sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+                    dt=1.0, save_everystep=false,
+                    callback=CallbackSet(callbacks, StepsizeCallback(cfl=1.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,
+                      info_callback=nothing, sol=nothing)
+
+        sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+                    dt=1.0, save_everystep=false,
+                    callback=CallbackSet(callbacks, StepsizeCallback(cfl=1.5)))
+
+        iter = saving_callback.affect!.affect!.latest_saved_iter
+        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, CarpenterKennedy2N54(williamson_condition=false),
+                            dt=1.0, save_everystep=false,
+                            callback=CallbackSet(UpdateCallback(),
+                                                 StepsizeCallback(cfl=1.5)))
+
+        result_restart = interpolate_line(start_point, end_point,
+                                          n_interpolation_points, sol_restart.prob.p.semi,
+                                          sol_restart.prob.p.semi.systems[1],
+                                          sol_restart, cut_off_bnd=false)
+
+        @test isapprox(result_full.velocity, result_restart.velocity, rtol=5e-5)
+        @test isapprox(result_full.density, result_restart.density, rtol=5e-6)
+        @test isapprox(result_full.pressure, result_restart.pressure, rtol=5e-4)
+    end
+
+    @trixi_testset "Poiseuille Flow 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, info_callback=nothing,
+                      tspan=(0.0, 0.3), sound_speed_factor=10, particle_spacing=4e-5)
+
+        iter = saving_callback.affect!.affect!.latest_saved_iter
+        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.semi.systems[2].cache.pressure_reference_values[2].pressure[]
+        restart_flow_rate = ode_restart.p.semi.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
+
+    @trixi_testset "Oscillating Beam Half-Simulation Restart" begin
+        # Run full simulation
+        trixi_include(@__MODULE__,
+                      joinpath(examples_dir(), "structure", "oscillating_beam_2d.jl"),
+                      tspan=(0.0, 1.0), n_particles_y=5, info_callback=nothing)
+
+        # Store the final solution for comparison
+        full_sol = sol
+        positions_full = full_sol.u[end].x[2]
+
+        # Run half simulation
+        trixi_include(@__MODULE__,
+                      joinpath(examples_dir(), "structure", "oscillating_beam_2d.jl"),
+                      tspan=(0.0, 0.5), n_particles_y=5)
+
+        # Get latest iteration and create restart
+        iter = saving_callback.affect!.affect!.latest_saved_iter
+        restart_file = joinpath("out", "structure_1_$iter.vtu")
+
+        ode_restart = semidiscretize(semi, (0.5, 1.0);
+                                     restart_with=restart_file)
+
+        sol_restart = solve(ode_restart, RDPK3SpFSAL49(), save_everystep=false, dt=1e-6)
+
+        # Compare the final particle velocities
+        @test isapprox(sol_restart.u[end].x[2], positions_full, rtol=8e-5)
+    end
+end