From f1400fdcb534ccaec50d5b7e9c41b4fd30d283ea Mon Sep 17 00:00:00 2001
From: Marcel Schurer <marcel.schurer@outlook.de>
Date: Fri, 17 Oct 2025 16:03:25 +0200
Subject: [PATCH 01/62] Update vtk2trixi function to return `NamedTuple`
 instead of `InitialCondition`

---
 src/io/read_vtk.jl | 27 ++++++++++++---------------
 1 file changed, 12 insertions(+), 15 deletions(-)

diff --git a/src/io/read_vtk.jl b/src/io/read_vtk.jl
index 90daa721a3..99d049fe8d 100644
--- a/src/io/read_vtk.jl
+++ b/src/io/read_vtk.jl
@@ -1,7 +1,7 @@
 """
 	vtk2trixi(file::String)
 
-Load VTK file and convert data to an [`InitialCondition`](@ref).
+Load VTK file and convert data to an NamedTuple.
 
 # Arguments
 - `file`: Name of the VTK file to be loaded.
@@ -18,17 +18,11 @@ rectangular = RectangularShape(0.1, (10, 10), (0, 0), density=1.5, velocity=(1.0
 # Write the `InitialCondition` to a vtk file
 trixi2vtk(rectangular; filename="rectangular", output_directory="out")
 
-# Read the vtk file and convert it to `InitialCondition`
-ic = vtk2trixi(joinpath("out", "rectangular.vtu"))
+# Read the vtk file and convert the data to an `NamedTuple`
+data = vtk2trixi(joinpath("out", "rectangular.vtu"))
 
 # output
-┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
-│ InitialCondition{Float64}                                                                        │
-│ ═════════════════════════                                                                        │
-│ #dimensions: ……………………………………………… 2                                                                │
-│ #particles: ………………………………………………… 100                                                              │
-│ particle spacing: ………………………………… 0.1                                                              │
-└──────────────────────────────────────────────────────────────────────────────────────────────────┘
+NamedTuple{data...}
 """
 function vtk2trixi(file)
     vtk_file = ReadVTK.VTKFile(file)
@@ -39,6 +33,8 @@ function vtk2trixi(file)
 
     # Retrieve fields
     ndims = first(ReadVTK.get_data(field_data["ndims"]))
+    time = "time" in keys(field_data) ? first(ReadVTK.get_data(field_data["time"])) : 0.0
+
     coordinates = ReadVTK.get_points(vtk_file)[1:ndims, :]
 
     fields = ["velocity", "density", "pressure", "mass", "particle_spacing"]
@@ -62,9 +58,10 @@ function vtk2trixi(file)
                        first(results["particle_spacing"]) :
                        results["particle_spacing"]
 
-    return InitialCondition(; coordinates, particle_spacing=particle_spacing,
-                            velocity=results["velocity"],
-                            mass=results["mass"],
-                            density=results["density"],
-                            pressure=results["pressure"])
+
+    return (time, coordinates, particle_spacing,
+            velocity=results["velocity"],
+            mass=results["mass"],
+            density=results["density"],
+            pressure=results["pressure"])
 end

From 217f71103261ea5a83a8e6984a021ffecf526b61 Mon Sep 17 00:00:00 2001
From: Marcel Schurer <marcel.schurer@outlook.de>
Date: Fri, 17 Oct 2025 18:08:43 +0200
Subject: [PATCH 02/62] Refactor vtk2trixi function to accept custom fields

---
 src/io/read_vtk.jl | 34 +++++++++++++++++-----------------
 1 file changed, 17 insertions(+), 17 deletions(-)

diff --git a/src/io/read_vtk.jl b/src/io/read_vtk.jl
index 99d049fe8d..66efb764ac 100644
--- a/src/io/read_vtk.jl
+++ b/src/io/read_vtk.jl
@@ -24,44 +24,44 @@ data = vtk2trixi(joinpath("out", "rectangular.vtu"))
 # output
 NamedTuple{data...}
 """
-function vtk2trixi(file)
+function vtk2trixi(file; custom_fields...)
     vtk_file = ReadVTK.VTKFile(file)
 
     # Retrieve data fields (e.g., pressure, velocity, ...)
     point_data = ReadVTK.get_point_data(vtk_file)
     field_data = ReadVTK.get_field_data(vtk_file)
 
+    results = Dict{Symbol, Any}()
+
     # Retrieve fields
     ndims = first(ReadVTK.get_data(field_data["ndims"]))
-    time = "time" in keys(field_data) ? first(ReadVTK.get_data(field_data["time"])) : 0.0
-
     coordinates = ReadVTK.get_points(vtk_file)[1:ndims, :]
 
-    fields = ["velocity", "density", "pressure", "mass", "particle_spacing"]
-    results = Dict{String, Array{Float64}}()
-
+    fields = [:velocity, :density, :pressure, :mass, :particle_spacing]
     for field in fields
         # Look for any key that contains the field name
         all_keys = keys(point_data)
-        idx = findfirst(k -> occursin(field, k), all_keys)
+        idx = findfirst(k -> occursin(string(field), k), all_keys)
         if idx !== nothing
             results[field] = ReadVTK.get_data(point_data[all_keys[idx]])
         else
             # Use zeros as default values when a field is missing
-            results[field] = field in ["mass"] ?
-                             zeros(size(coordinates, 2)) : zero(coordinates)
+            results[field] = string(field) in ["mass"] ?
+                                    zeros(size(coordinates, 2)) : zero(coordinates)
             @info "No '$field' field found in VTK file. Will be set to zero."
         end
     end
 
-    particle_spacing = allequal(results["particle_spacing"]) ?
-                       first(results["particle_spacing"]) :
-                       results["particle_spacing"]
+    results[:particle_spacing] = allequal(results[:particle_spacing]) ?
+                                          first(results[:particle_spacing]) :
+                                          results[:particle_spacing]
+    results[:coordinates] = coordinates
+    results[:time] = "time" in keys(field_data) ?
+                                    first(ReadVTK.get_data(field_data["time"])) : 0.0
 
+    for (key, field) in custom_fields
+        results[key] = first(ReadVTK.get_data(field_data[field]))
+    end
 
-    return (time, coordinates, particle_spacing,
-            velocity=results["velocity"],
-            mass=results["mass"],
-            density=results["density"],
-            pressure=results["pressure"])
+    return NamedTuple(results)
 end

From 0fa590038776711b3023bb5b82758c587e2a765d Mon Sep 17 00:00:00 2001
From: Marcel Schurer <marcel.schurer@outlook.de>
Date: Mon, 20 Oct 2025 17:32:41 +0200
Subject: [PATCH 03/62] Enhance vtk2trixi function to support scalar and vector
 custom quantities and update tests

---
 src/io/read_vtk.jl  |  36 ++++++----
 test/io/read_vtk.jl | 164 +++++++++++++++++++++++++++++++++-----------
 2 files changed, 149 insertions(+), 51 deletions(-)

diff --git a/src/io/read_vtk.jl b/src/io/read_vtk.jl
index 66efb764ac..abd610a18b 100644
--- a/src/io/read_vtk.jl
+++ b/src/io/read_vtk.jl
@@ -1,10 +1,15 @@
 """
-	vtk2trixi(file::String)
+	vtk2trixi(file::String; custom_quantities...)
 
 Load VTK file and convert data to an NamedTuple.
 
 # Arguments
-- `file`: Name of the VTK file to be loaded.
+- `file`:                 Name of the VTK file to be loaded.
+- `custom_quantities...`: Additional custom quantities to be loaded from the VTK file.
+                          Each custom quantity must be explicitly listed in the
+                          `custom_quantities` during the simulation to ensure it is
+                          included in the VTK output and can be successfully loaded.
+                          See [Custom Quantities](@ref custom_quantities) for details.
 
 !!! warning "Experimental Implementation"
     This is an experimental feature and may change in any future releases.
@@ -15,16 +20,18 @@ Load VTK file and convert data to an NamedTuple.
 rectangular = RectangularShape(0.1, (10, 10), (0, 0), density=1.5, velocity=(1.0, -2.0),
                                pressure=1000.0)
 
-# Write the `InitialCondition` to a vtk file
-trixi2vtk(rectangular; filename="rectangular", output_directory="out")
+# Write the `InitialCondition` with custom quantity to a vtk file
+trixi2vtk(rectangular; filename="rectangular", output_directory="out",
+          my_custom_quantity=3.0)
 
 # Read the vtk file and convert the data to an `NamedTuple`
-data = vtk2trixi(joinpath("out", "rectangular.vtu"))
+data = vtk2trixi(joinpath("out", "rectangular.vtu");
+                 my_custom_quantity="my_custom_quantity")
 
 # output
 NamedTuple{data...}
 """
-function vtk2trixi(file; custom_fields...)
+function vtk2trixi(file; custom_quantities...)
     vtk_file = ReadVTK.VTKFile(file)
 
     # Retrieve data fields (e.g., pressure, velocity, ...)
@@ -47,20 +54,25 @@ function vtk2trixi(file; custom_fields...)
         else
             # Use zeros as default values when a field is missing
             results[field] = string(field) in ["mass"] ?
-                                    zeros(size(coordinates, 2)) : zero(coordinates)
+                             zeros(size(coordinates, 2)) : zero(coordinates)
             @info "No '$field' field found in VTK file. Will be set to zero."
         end
     end
 
     results[:particle_spacing] = allequal(results[:particle_spacing]) ?
-                                          first(results[:particle_spacing]) :
-                                          results[:particle_spacing]
+                                 first(results[:particle_spacing]) :
+                                 results[:particle_spacing]
     results[:coordinates] = coordinates
     results[:time] = "time" in keys(field_data) ?
-                                    first(ReadVTK.get_data(field_data["time"])) : 0.0
+                     first(ReadVTK.get_data(field_data["time"])) : 0.0
 
-    for (key, field) in custom_fields
-        results[key] = first(ReadVTK.get_data(field_data[field]))
+    for (key, quantity) in custom_quantities
+        if quantity in keys(point_data)
+            results[key] = ReadVTK.get_data(point_data[quantity])
+        end
+        if quantity in keys(field_data)
+            results[key] = first(ReadVTK.get_data(field_data[quantity]))
+        end
     end
 
     return NamedTuple(results)
diff --git a/test/io/read_vtk.jl b/test/io/read_vtk.jl
index d82e8c50c3..c603c597e0 100644
--- a/test/io/read_vtk.jl
+++ b/test/io/read_vtk.jl
@@ -3,19 +3,50 @@
         coordinates = fill(1.0, 2, 12)
         velocity = fill(2.0, 2, 12)
 
-        expected_ic = InitialCondition(; coordinates=coordinates, velocity=velocity,
+        expected_ic = InitialCondition(coordinates=coordinates, velocity=velocity,
                                        density=1000.0, pressure=900.0, mass=50.0)
 
-        @testset verbose=true "`InitialCondition`" begin
-            trixi2vtk(expected_ic; filename="tmp_initial_condition",
-                      output_directory=tmp_dir)
-
-            test_ic = vtk2trixi(joinpath(tmp_dir, "tmp_initial_condition.vtu"))
+        expected_cq_scalar = 3.0
+        expected_cq_vector = fill(expected_cq_scalar,
+                                  size(expected_ic.coordinates, 2))
+        scalar_cq_function(system, data, t) = expected_cq_scalar
+        vector_cq_function(system, data,
+                           t) = fill(expected_cq_scalar, nparticles(system))
 
-            @test isapprox(expected_ic.coordinates, test_ic.coordinates, rtol=1e-5)
-            @test isapprox(expected_ic.velocity, test_ic.velocity, rtol=1e-5)
-            @test isapprox(expected_ic.density, test_ic.density, rtol=1e-5)
-            @test isapprox(expected_ic.pressure, test_ic.pressure, rtol=1e-5)
+        @testset verbose=true "`InitialCondition`" begin
+            @testset verbose=true "Scalar custom quantity" begin
+                trixi2vtk(expected_ic; filename="tmp_initial_condition_scalar",
+                          output_directory=tmp_dir,
+                          cq_scalar=expected_cq_scalar)
+
+                # Load file containing test data
+                test_data = vtk2trixi(joinpath(tmp_dir, "tmp_initial_condition_scalar.vtu");
+                                      cq_scalar="cq_scalar")
+
+                @test isapprox(expected_ic.coordinates, test_data.coordinates, rtol=1e-5)
+                @test isapprox(expected_ic.velocity, test_data.velocity, rtol=1e-5)
+                @test isapprox(expected_ic.density, test_data.density, rtol=1e-5)
+                @test isapprox(expected_ic.pressure, test_data.pressure, rtol=1e-5)
+                @test isapprox(expected_cq_scalar, test_data.cq_scalar,
+                               rtol=1e-5)
+            end
+
+            @testset verbose=true "Vector custom quantity" begin
+                trixi2vtk(expected_ic; filename="tmp_initial_condition_vector",
+                          output_directory=tmp_dir,
+                          cq_vector=expected_cq_vector)
+
+                # Load file containing test data
+                test_data = vtk2trixi(joinpath(tmp_dir, "tmp_initial_condition_vector.vtu");
+                                      cq_vector="cq_vector")
+
+                @test isapprox(expected_ic.coordinates, test_data.coordinates, rtol=1e-5)
+                @test isapprox(expected_ic.velocity, test_data.velocity, rtol=1e-5)
+                @test isapprox(expected_ic.density, test_data.density, rtol=1e-5)
+                @test isapprox(expected_ic.pressure, test_data.pressure, rtol=1e-5)
+                @test isapprox(expected_cq_vector, test_data.cq_vector,
+                               rtol=1e-5)
+            end
         end
 
         @testset verbose=true "`AbstractFluidSystem`" begin
@@ -29,27 +60,52 @@
             semi = Semidiscretization(fluid_system)
 
             # Create random ODE solutions
-            dvdu_ode = nothing
             v = fill(2.0, ndims(fluid_system), nparticles(fluid_system))
             pressure = fill(3.0, nparticles(fluid_system))
             v_ode = vec([v; pressure'])
             u = fill(1.0, ndims(fluid_system), nparticles(fluid_system))
             u_ode = vec(u)
+            dv_ode = zero(v_ode)
+            du_ode = zero(u_ode)
             x = (; v_ode, u_ode)
             vu_ode = (; x)
-
-            # Write out `AbstractFluidSystem` Simulation-File
-            trixi2vtk(fluid_system, dvdu_ode, vu_ode, semi, 0.0,
-                      nothing; system_name="tmp_file_fluid", output_directory=tmp_dir,
-                      iter=1)
-
-            # Load `AbstractFluidSystem` Simulation-File
-            test = vtk2trixi(joinpath(tmp_dir, "tmp_file_fluid_1.vtu"))
-
-            @test isapprox(u, test.coordinates, rtol=1e-5)
-            @test isapprox(v, test.velocity, rtol=1e-5)
-            @test isapprox(pressure, test.pressure, rtol=1e-5)
-            @test isapprox(fluid_system.cache.density, test.density, rtol=1e-5)
+            dvdu_ode = (; x=(; v_ode=dv_ode, u_ode=du_ode))
+
+            @testset verbose=true "Scalar custom quantity" begin
+                trixi2vtk(fluid_system, dvdu_ode, vu_ode, semi, 0.0,
+                          nothing; system_name="tmp_file_fluid_scalar",
+                          output_directory=tmp_dir,
+                          iter=1, cq_scalar=scalar_cq_function)
+
+                # Load file containing test data
+                test_data = vtk2trixi(joinpath(tmp_dir, "tmp_file_fluid_scalar_1.vtu");
+                                      cq_scalar="cq_scalar")
+
+                @test isapprox(u, test_data.coordinates, rtol=1e-5)
+                @test isapprox(v, test_data.velocity, rtol=1e-5)
+                @test isapprox(pressure, test_data.pressure, rtol=1e-5)
+                @test isapprox(fluid_system.cache.density, test_data.density, rtol=1e-5)
+                @test isapprox(expected_cq_scalar, test_data.cq_scalar,
+                               rtol=1e-5)
+            end
+
+            @testset verbose=true "Vector custom quantity" begin
+                trixi2vtk(fluid_system, dvdu_ode, vu_ode, semi, 0.0,
+                          nothing; system_name="tmp_file_fluid_vector",
+                          output_directory=tmp_dir,
+                          iter=1, cq_vector=vector_cq_function)
+
+                # Load file containing test data
+                test_data = vtk2trixi(joinpath(tmp_dir, "tmp_file_fluid_vector_1.vtu");
+                                      cq_vector="cq_vector")
+
+                @test isapprox(u, test_data.coordinates, rtol=1e-5)
+                @test isapprox(v, test_data.velocity, rtol=1e-5)
+                @test isapprox(pressure, test_data.pressure, rtol=1e-5)
+                @test isapprox(fluid_system.cache.density, test_data.density, rtol=1e-5)
+                @test isapprox(fill(expected_cq_scalar, nparticles(fluid_system)),
+                               test_data.cq_vector, rtol=1e-5)
+            end
         end
 
         @testset verbose=true "`WallBoundarySystem`" begin
@@ -67,25 +123,55 @@
             semi = Semidiscretization(boundary_system)
 
             # Create dummy ODE solutions
-            dvdu_ode = nothing
             v_ode = zeros(ndims(boundary_system) * nparticles(boundary_system))
             u_ode = zeros(ndims(boundary_system) * nparticles(boundary_system))
+            dv_ode = zero(v_ode)
+            du_ode = zero(u_ode)
             x = (; v_ode, u_ode)
             vu_ode = (; x)
-
-            # Write out `WallBoundarySystem` Simulation-File
-            trixi2vtk(boundary_system, dvdu_ode, vu_ode, semi, 0.0,
-                      nothing; system_name="tmp_file_boundary", output_directory=tmp_dir,
-                      iter=1)
-
-            # Load `WallBoundarySystem` Simulation-File
-            test = vtk2trixi(joinpath(tmp_dir, "tmp_file_boundary_1.vtu"))
-
-            @test isapprox(boundary_system.coordinates, test.coordinates, rtol=1e-5)
-            # The velocity is always zero for `WallBoundarySystem`
-            @test isapprox(zeros(size(test.velocity)), test.velocity, rtol=1e-5)
-            @test isapprox(boundary_model.pressure, test.pressure, rtol=1e-5)
-            @test isapprox(boundary_model.cache.density, test.density, rtol=1e-5)
+            dvdu_ode = (; x=(; v_ode=dv_ode, u_ode=du_ode))
+
+            @testset verbose=true "Scalar custom quantity" begin
+                trixi2vtk(boundary_system, dvdu_ode, vu_ode, semi, 0.0,
+                          nothing; system_name="tmp_file_boundary_scalar",
+                          output_directory=tmp_dir,
+                          iter=1, cq_scalar=scalar_cq_function)
+
+                # Load file containing test data
+                test_data = vtk2trixi(joinpath(tmp_dir, "tmp_file_boundary_scalar_1.vtu");
+                                      cq_scalar="cq_scalar")
+
+                @test isapprox(boundary_system.coordinates, test_data.coordinates,
+                               rtol=1e-5)
+                # The velocity is always zero for `WallBoundarySystem`
+                @test isapprox(zeros(size(test_data.velocity)), test_data.velocity,
+                               rtol=1e-5)
+                @test isapprox(boundary_model.pressure, test_data.pressure, rtol=1e-5)
+                @test isapprox(boundary_model.cache.density, test_data.density, rtol=1e-5)
+                @test isapprox(expected_cq_scalar, test_data.cq_scalar,
+                               rtol=1e-5)
+            end
+
+            @testset verbose=true "Vector custom quantity" begin
+                trixi2vtk(boundary_system, dvdu_ode, vu_ode, semi, 0.0,
+                          nothing; system_name="tmp_file_boundary_vector",
+                          output_directory=tmp_dir,
+                          iter=1, cq_vector=vector_cq_function)
+
+                # Load file containing test data
+                test_data = vtk2trixi(joinpath(tmp_dir, "tmp_file_boundary_vector_1.vtu");
+                                      cq_vector="cq_vector")
+
+                @test isapprox(boundary_system.coordinates, test_data.coordinates,
+                               rtol=1e-5)
+                # The velocity is always zero for `WallBoundarySystem`
+                @test isapprox(zeros(size(test_data.velocity)), test_data.velocity,
+                               rtol=1e-5)
+                @test isapprox(boundary_model.pressure, test_data.pressure, rtol=1e-5)
+                @test isapprox(boundary_model.cache.density, test_data.density, rtol=1e-5)
+                @test isapprox(fill(expected_cq_scalar, nparticles(boundary_system)),
+                               test_data.cq_vector, rtol=1e-5)
+            end
         end
     end
 end

From 38a9b2aebfdf845a357fa903ad02d0380f5b7fe3 Mon Sep 17 00:00:00 2001
From: Marcel Schurer <marcel.schurer@outlook.de>
Date: Tue, 21 Oct 2025 14:31:30 +0200
Subject: [PATCH 04/62] Fix typo in documentation for vtk2trixi function

---
 src/io/read_vtk.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/io/read_vtk.jl b/src/io/read_vtk.jl
index abd610a18b..e189d81bfc 100644
--- a/src/io/read_vtk.jl
+++ b/src/io/read_vtk.jl
@@ -1,7 +1,7 @@
 """
 	vtk2trixi(file::String; custom_quantities...)
 
-Load VTK file and convert data to an NamedTuple.
+Load VTK file and convert data to a NamedTuple.
 
 # Arguments
 - `file`:                 Name of the VTK file to be loaded.

From 5d1192aaa014c3c47e7d318b5fce22b9b607e15a Mon Sep 17 00:00:00 2001
From: Marcel Schurer <marcel.schurer@outlook.de>
Date: Thu, 23 Oct 2025 11:30:38 +0200
Subject: [PATCH 05/62] Fix capitalization in comments

---
 src/io/read_vtk.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/io/read_vtk.jl b/src/io/read_vtk.jl
index e189d81bfc..02e341f1ae 100644
--- a/src/io/read_vtk.jl
+++ b/src/io/read_vtk.jl
@@ -20,11 +20,11 @@ Load VTK file and convert data to a NamedTuple.
 rectangular = RectangularShape(0.1, (10, 10), (0, 0), density=1.5, velocity=(1.0, -2.0),
                                pressure=1000.0)
 
-# Write the `InitialCondition` with custom quantity to a vtk file
+# Write the `InitialCondition` with custom quantity to a VTK file
 trixi2vtk(rectangular; filename="rectangular", output_directory="out",
           my_custom_quantity=3.0)
 
-# Read the vtk file and convert the data to an `NamedTuple`
+# Read the VTK file and convert the data to a `NamedTuple`
 data = vtk2trixi(joinpath("out", "rectangular.vtu");
                  my_custom_quantity="my_custom_quantity")
 

From 7e51a814bfc6e2d7ea693f96d0724bfa8e0d73aa Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Sat, 22 Nov 2025 08:54:34 +0100
Subject: [PATCH 06/62] add ELTYPE

---
 src/io/read_vtk.jl | 13 ++++++++-----
 1 file changed, 8 insertions(+), 5 deletions(-)

diff --git a/src/io/read_vtk.jl b/src/io/read_vtk.jl
index 90daa721a3..d87a586bf2 100644
--- a/src/io/read_vtk.jl
+++ b/src/io/read_vtk.jl
@@ -29,8 +29,10 @@ ic = vtk2trixi(joinpath("out", "rectangular.vtu"))
 │ #particles: ………………………………………………… 100                                                              │
 │ particle spacing: ………………………………… 0.1                                                              │
 └──────────────────────────────────────────────────────────────────────────────────────────────────┘
+```
 """
-function vtk2trixi(file)
+function vtk2trixi(file; element_type=Float64)
+    ELTYPE = element_type
     vtk_file = ReadVTK.VTKFile(file)
 
     # Retrieve data fields (e.g., pressure, velocity, ...)
@@ -39,7 +41,7 @@ function vtk2trixi(file)
 
     # Retrieve fields
     ndims = first(ReadVTK.get_data(field_data["ndims"]))
-    coordinates = ReadVTK.get_points(vtk_file)[1:ndims, :]
+    coordinates = convert.(ELTYPE, ReadVTK.get_points(vtk_file)[1:ndims, :])
 
     fields = ["velocity", "density", "pressure", "mass", "particle_spacing"]
     results = Dict{String, Array{Float64}}()
@@ -49,11 +51,11 @@ function vtk2trixi(file)
         all_keys = keys(point_data)
         idx = findfirst(k -> occursin(field, k), all_keys)
         if idx !== nothing
-            results[field] = ReadVTK.get_data(point_data[all_keys[idx]])
+            results[field] = convert.(ELTYPE, ReadVTK.get_data(point_data[all_keys[idx]]))
         else
             # Use zeros as default values when a field is missing
             results[field] = field in ["mass"] ?
-                             zeros(size(coordinates, 2)) : zero(coordinates)
+                             zeros(ELTYPE, size(coordinates, 2)) : zero(coordinates)
             @info "No '$field' field found in VTK file. Will be set to zero."
         end
     end
@@ -62,7 +64,8 @@ function vtk2trixi(file)
                        first(results["particle_spacing"]) :
                        results["particle_spacing"]
 
-    return InitialCondition(; coordinates, particle_spacing=particle_spacing,
+    return InitialCondition(; coordinates,
+                            particle_spacing=convert(ELTYPE, particle_spacing),
                             velocity=results["velocity"],
                             mass=results["mass"],
                             density=results["density"],

From 2b2fb72f724f4437fcdc35b15c0354f60ba9d0d5 Mon Sep 17 00:00:00 2001
From: Marcel Schurer <marcel.schurer@outlook.de>
Date: Wed, 3 Dec 2025 08:40:05 +0100
Subject: [PATCH 07/62] Fix formatting in documentation

---
 src/io/read_vtk.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/io/read_vtk.jl b/src/io/read_vtk.jl
index 02e341f1ae..bdf1b7033e 100644
--- a/src/io/read_vtk.jl
+++ b/src/io/read_vtk.jl
@@ -30,6 +30,7 @@ data = vtk2trixi(joinpath("out", "rectangular.vtu");
 
 # output
 NamedTuple{data...}
+```
 """
 function vtk2trixi(file; custom_quantities...)
     vtk_file = ReadVTK.VTKFile(file)

From b043bd11e4cc43cebb23bdba1388f5174f59a653 Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Mon, 8 Dec 2025 17:52:26 +0100
Subject: [PATCH 08/62] add `coordinates_eltype`

---
 src/io/read_vtk.jl | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/io/read_vtk.jl b/src/io/read_vtk.jl
index 326773d4b4..884588033b 100644
--- a/src/io/read_vtk.jl
+++ b/src/io/read_vtk.jl
@@ -33,8 +33,9 @@ ic = vtk2trixi(joinpath("out", "rectangular.vtu"))
 └──────────────────────────────────────────────────────────────────────────────────────────────────┘
 ```
 """
-function vtk2trixi(file; element_type=Float64)
+function vtk2trixi(file; element_type=Float64, coordinates_eltype=Float64)
     ELTYPE = element_type
+    cELTYPE = coordinates_eltype
     vtk_file = ReadVTK.VTKFile(file)
 
     # Retrieve data fields (e.g., pressure, velocity, ...)
@@ -43,7 +44,7 @@ function vtk2trixi(file; element_type=Float64)
 
     # Retrieve fields
     ndims = first(ReadVTK.get_data(field_data["ndims"]))
-    coordinates = convert.(ELTYPE, ReadVTK.get_points(vtk_file)[1:ndims, :])
+    coordinates = convert.(cELTYPE, ReadVTK.get_points(vtk_file)[1:ndims, :])
 
     fields = ["velocity", "density", "pressure", "mass", "particle_spacing"]
     results = Dict{String, Array{Float64}}()

From a9b0f2f9df1d2b9f247ff6053ab27a6d082766ed Mon Sep 17 00:00:00 2001
From: Marcel Schurer <marcel.schurer@outlook.de>
Date: Tue, 9 Dec 2025 14:03:28 +0100
Subject: [PATCH 09/62] Corrected failing doctest output

---
 src/io/read_vtk.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/io/read_vtk.jl b/src/io/read_vtk.jl
index bdf1b7033e..4189ea474a 100644
--- a/src/io/read_vtk.jl
+++ b/src/io/read_vtk.jl
@@ -29,7 +29,7 @@ data = vtk2trixi(joinpath("out", "rectangular.vtu");
                  my_custom_quantity="my_custom_quantity")
 
 # output
-NamedTuple{data...}
+
 ```
 """
 function vtk2trixi(file; custom_quantities...)

From 0359425131784b5ebc3064a15fdc03b2613d615f Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Wed, 10 Dec 2025 09:30:19 +0100
Subject: [PATCH 10/62] first prototype

---
 .../boundary/open_boundary/boundary_zones.jl  | 23 ++++-
 src/schemes/boundary/open_boundary/system.jl  | 99 ++++++++++++++++---
 2 files changed, 102 insertions(+), 20 deletions(-)

diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl
index 7088fea58e..be5979ba6e 100644
--- a/src/schemes/boundary/open_boundary/boundary_zones.jl
+++ b/src/schemes/boundary/open_boundary/boundary_zones.jl
@@ -147,7 +147,7 @@ bidirectional_flow = BoundaryZone(; boundary_face=face_vertices, face_normal,
 !!! warning "Experimental Implementation"
     This is an experimental feature and may change in any future releases.
 """
-struct BoundaryZone{IC, S, ZO, ZW, FD, FN, ELTYPE, R}
+struct BoundaryZone{IC, S, ZO, ZW, FD, FN, ELTYPE, R, C}
     initial_condition :: IC
     spanning_set      :: S
     zone_origin       :: ZO
@@ -156,6 +156,7 @@ struct BoundaryZone{IC, S, ZO, ZW, FD, FN, ELTYPE, R}
     face_normal       :: FN
     rest_pressure     :: ELTYPE # Only required for `BoundaryModelDynamicalPressureZhang`
     reference_values  :: R
+    cache             :: C
     # Note that the following can't be static type parameters, as all boundary zones in a system
     # must have the same type, so that we can loop over them in a type-stable way.
     average_inflow_velocity :: Bool
@@ -167,7 +168,7 @@ end
 function BoundaryZone(; boundary_face, face_normal, density, particle_spacing,
                       initial_condition=nothing, extrude_geometry=nothing,
                       open_boundary_layers::Integer, average_inflow_velocity=true,
-                      boundary_type=BidirectionalFlow(),
+                      boundary_type=BidirectionalFlow(), sample_points=nothing,
                       rest_pressure=zero(eltype(density)),
                       reference_density=nothing, reference_pressure=nothing,
                       reference_velocity=nothing)
@@ -262,10 +263,12 @@ function BoundaryZone(; boundary_face, face_normal, density, particle_spacing,
         ic.velocity .= stack(velocity_ref.(coordinates_svector, 0))
     end
 
+    cache = (; create_cache_boundary_zone(ic, sample_points)...)
+
     return BoundaryZone(ic, spanning_set_, zone_origin, zone_width,
                         flow_direction, face_normal_, rest_pressure, reference_values,
-                        average_inflow_velocity, prescribed_density, prescribed_pressure,
-                        prescribed_velocity)
+                        cache, average_inflow_velocity, prescribed_density,
+                        prescribed_pressure, prescribed_velocity)
 end
 
 function boundary_type_name(boundary_zone::BoundaryZone)
@@ -301,6 +304,16 @@ function Base.show(io::IO, ::MIME"text/plain", boundary_zone::BoundaryZone)
     end
 end
 
+create_cache_boundary_zone(initial_condition, sample_points::Nothing) = (;)
+
+function create_cache_boundary_zone(initial_condition, sample_points::Matrix)
+    # TODO: Check matrix (ndims etc.)
+    shepard_coefficient = zeros(eltype(initial_condition), axes(sample_points, 2))
+    dA = initial_condition.particle_spacing
+    return (; sample_points=sample_points, sample_velocity=copy(sample_points),
+            shepard_coefficient, dA)
+end
+
 function set_up_boundary_zone(boundary_face, face_normal, density, particle_spacing,
                               initial_condition, extrude_geometry, open_boundary_layers,
                               boundary_type)
@@ -453,7 +466,7 @@ function update_boundary_zone_indices!(system, u, boundary_zones, semi)
         # - Floating-point rounding when a particle lies almost exactly on the `boundary_face`
         #   during transition, causing a reset just outside the zone
         #   (fixed in https://github.com/trixi-framework/TrixiParticles.jl/pull/997).
-        @assert system.boundary_zone_indices[particle] != 0 "No boundary zone found for active buffer particle"
+        @assert system.boundary_zone_indices[particle]!=0 "No boundary zone found for active buffer particle"
     end
 
     return system
diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl
index 9396c3f1fb..2c4c5c6c0c 100644
--- a/src/schemes/boundary/open_boundary/system.jl
+++ b/src/schemes/boundary/open_boundary/system.jl
@@ -63,7 +63,7 @@ end
 
 function OpenBoundarySystem(boundary_zones::Union{BoundaryZone, Nothing}...;
                             fluid_system::AbstractFluidSystem, buffer_size::Integer,
-                            boundary_model,
+                            boundary_model, calculate_flow_rate=false,
                             pressure_acceleration=boundary_model isa
                                                   BoundaryModelDynamicalPressureZhang ?
                                                   fluid_system.pressure_acceleration_formulation :
@@ -85,7 +85,7 @@ function OpenBoundarySystem(boundary_zones::Union{BoundaryZone, Nothing}...;
     cache = (;
              create_cache_shifting(initial_conditions, shifting_technique)...,
              create_cache_open_boundary(boundary_model, fluid_system, initial_conditions,
-                                        boundary_zones_)...)
+                                        calculate_flow_rate, boundary_zones_)...)
 
     fluid_system_index = Ref(0)
 
@@ -110,6 +110,7 @@ function OpenBoundarySystem(boundary_zones::Union{BoundaryZone, Nothing}...;
                                                   zone.face_normal,
                                                   zone.rest_pressure,
                                                   nothing,
+                                                  zone.cache,
                                                   zone.average_inflow_velocity,
                                                   zone.prescribed_density,
                                                   zone.prescribed_pressure,
@@ -132,7 +133,7 @@ function initialize!(system::OpenBoundarySystem, semi)
 end
 
 function create_cache_open_boundary(boundary_model, fluid_system, initial_condition,
-                                    boundary_zones)
+                                    calculate_flow_rate, boundary_zones)
     reference_values = map(bz -> bz.reference_values, boundary_zones)
     ELTYPE = eltype(initial_condition)
 
@@ -141,6 +142,15 @@ function create_cache_open_boundary(boundary_model, fluid_system, initial_condit
     density_reference_values = map(ref -> ref.reference_density, reference_values)
     velocity_reference_values = map(ref -> ref.reference_velocity, reference_values)
 
+    cache = (; pressure_reference_values=pressure_reference_values,
+             density_reference_values=density_reference_values,
+             velocity_reference_values=velocity_reference_values, calculate_flow_rate)
+
+    if calculate_flow_rate
+        boundary_zones_flow_rate = zeros(ELTYPE, length(boundary_zones))
+        cache = (; boundary_zones_flow_rate, cache...)
+    end
+
     if boundary_model isa BoundaryModelCharacteristicsLastiwka
         characteristics = zeros(ELTYPE, 3, nparticles(initial_condition))
         previous_characteristics = zeros(ELTYPE, 3, nparticles(initial_condition))
@@ -148,10 +158,8 @@ function create_cache_open_boundary(boundary_model, fluid_system, initial_condit
         return (; characteristics=characteristics,
                 previous_characteristics=previous_characteristics,
                 pressure=copy(initial_condition.pressure),
-                density=copy(initial_condition.density),
-                pressure_reference_values=pressure_reference_values,
-                density_reference_values=density_reference_values,
-                velocity_reference_values=velocity_reference_values)
+                density=copy(initial_condition.density), cache...)
+
     elseif boundary_model isa BoundaryModelDynamicalPressureZhang
         # A separate array for the boundary pressure is required,
         # since it is specified independently from the computed pressure for the momentum equation.
@@ -171,10 +179,7 @@ function create_cache_open_boundary(boundary_model, fluid_system, initial_condit
         cache = (; density_calculator=ContinuityDensity(),
                  density_diffusion=density_diffusion_,
                  pressure_boundary=pressure_boundary,
-                 density_rest=density_rest,
-                 pressure_reference_values=pressure_reference_values,
-                 density_reference_values=density_reference_values,
-                 velocity_reference_values=velocity_reference_values)
+                 density_rest=density_rest, cache...)
 
         if fluid_system isa EntropicallyDampedSPHSystem
             # Density and pressure is stored in `v`
@@ -186,10 +191,7 @@ function create_cache_open_boundary(boundary_model, fluid_system, initial_condit
     else
         return (;
                 pressure=copy(initial_condition.pressure),
-                density=copy(initial_condition.density),
-                pressure_reference_values=pressure_reference_values,
-                density_reference_values=density_reference_values,
-                velocity_reference_values=velocity_reference_values)
+                density=copy(initial_condition.density), cache...)
     end
 end
 
@@ -305,6 +307,10 @@ function update_open_boundary_eachstep!(system::OpenBoundarySystem, v_ode, u_ode
         u = wrap_u(u_ode, system, semi)
         v = wrap_v(v_ode, system, semi)
 
+        # This must be called before `update_pressure_model!` and `check_domain!`
+        # to ensure quantities remain consistent with the current simulation state.
+        calculate_flow_rate!(system, v, u, v_ode, u_ode, semi)
+
         @trixi_timeit timer() "check domain" check_domain!(system, v, u, v_ode, u_ode, semi)
 
         update_pressure_model!(system, v, u, semi, integrator.dt)
@@ -323,6 +329,28 @@ end
 
 update_open_boundary_eachstep!(system, v_ode, u_ode, semi, t, integrator) = system
 
+function calculate_flow_rate!(system, v, u, v_ode, u_ode, semi)
+    system.cache.calculate_flow_rate || return system
+
+    for boundary_zone in system.boundary_zones
+        interpolate_velocity!(system, boundary_zone, v, u, v_ode, u_ode, semi)
+    end
+
+    foreach_enumerate(system.boundary_zones) do (zone_id, boundary_zone)
+        (; face_normal) = boundary_zone
+        (; sample_velocity, dA) = boundary_zone.cache
+        # Compute volumetric flow rate: Q = ∫ v ⋅ n dA
+        current_flow_rate = sum(axes(sample_velocity, 2)) do point
+            vn = dot(current_velocity(sample_velocity, system, point), -face_normal)
+            return vn * dA
+        end
+
+        system.cache.boundary_zones_flow_rate[zone_id] = current_flow_rate
+    end
+
+    return system
+end
+
 function check_domain!(system, v, u, v_ode, u_ode, semi)
     (; boundary_zones, boundary_candidates, fluid_candidates, fluid_system) = system
 
@@ -570,3 +598,44 @@ end
 
     return system
 end
+
+function interpolate_velocity!(system::OpenBoundarySystem, boundary_zone, v, u,
+                               v_ode, u_ode, semi)
+    (; sample_points, sample_velocity, shepard_coefficient) = boundary_zone.cache
+    smoothing_length = initial_smoothing_length(system)
+    smoothing_kernel = system_smoothing_kernel(system)
+
+    set_zero!(shepard_coefficient)
+    set_zero!(sample_velocity)
+
+    foreach_system(semi) do neighbor_system
+        v_neighbor = wrap_v(v_ode, neighbor_system, semi)
+        u_neighbor = wrap_u(u_ode, neighbor_system, semi)
+        neighbor_coords = current_coordinates(u_neighbor, neighbor_system)
+
+        foreach_point_neighbor(system, neighbor_system, sample_points, neighbor_coords,
+                               semi,
+                               points=axes(sample_points, 2)) do point, neighbor,
+                                                                 pos_diff, distance
+            m_b = hydrodynamic_mass(neighbor_system, neighbor)
+            volume_b = m_b / current_density(v_neighbor, neighbor_system, neighbor)
+            W_ab = kernel(smoothing_kernel, distance, smoothing_length)
+            shepard_coefficient[point] += volume_b * W_ab
+
+            velocity_neighbor = viscous_velocity(v_neighbor, neighbor_system, neighbor)
+            for i in axes(velocity_neighbor, 1)
+                sample_velocity[i, point] += velocity_neighbor[i] * volume_b * W_ab
+            end
+        end
+    end
+
+    @threaded semi for point in axes(sample_points, 2)
+        if shepard_coefficient[point] > eps()
+            for i in axes(sample_velocity, 1)
+                sample_velocity[i, point] /= shepard_coefficient[point]
+            end
+        end
+    end
+
+    return system
+end

From fd15671fcb56a9a37a276d75a396457b10b94cd5 Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Wed, 10 Dec 2025 11:03:16 +0100
Subject: [PATCH 11/62] write Q

---
 src/io/write_vtk.jl                                  | 8 ++++++++
 src/schemes/boundary/open_boundary/boundary_zones.jl | 2 +-
 2 files changed, 9 insertions(+), 1 deletion(-)

diff --git a/src/io/write_vtk.jl b/src/io/write_vtk.jl
index 07fa559d2e..d1146e67b1 100644
--- a/src/io/write_vtk.jl
+++ b/src/io/write_vtk.jl
@@ -409,6 +409,14 @@ function write2vtk!(vtk, v, u, t, system::OpenBoundarySystem)
     vtk["pressure"] = [current_pressure(v, system, particle)
                        for particle in eachparticle(system)]
 
+    if system.cache.calculate_flow_rate
+        for i in eachindex(system.cache.boundary_zones_flow_rate)
+            vtk["Q_$i"] = system.cache.boundary_zones_flow_rate[i]
+        end
+
+        vtk["Q_total"] = sum(system.cache.boundary_zones_flow_rate)
+    end
+
     return vtk
 end
 
diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl
index be5979ba6e..2bde137093 100644
--- a/src/schemes/boundary/open_boundary/boundary_zones.jl
+++ b/src/schemes/boundary/open_boundary/boundary_zones.jl
@@ -466,7 +466,7 @@ function update_boundary_zone_indices!(system, u, boundary_zones, semi)
         # - Floating-point rounding when a particle lies almost exactly on the `boundary_face`
         #   during transition, causing a reset just outside the zone
         #   (fixed in https://github.com/trixi-framework/TrixiParticles.jl/pull/997).
-        @assert system.boundary_zone_indices[particle]!=0 "No boundary zone found for active buffer particle"
+        @assert system.boundary_zone_indices[particle] != 0 "No boundary zone found for active buffer particle"
     end
 
     return system

From 364e24ea6894836aee569d2e11d7a934c0a7642b Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Wed, 10 Dec 2025 11:42:06 +0100
Subject: [PATCH 12/62] gpu fix

---
 src/io/write_vtk.jl                          |  6 ++++--
 src/schemes/boundary/open_boundary/system.jl | 20 +++++++++++++-------
 2 files changed, 17 insertions(+), 9 deletions(-)

diff --git a/src/io/write_vtk.jl b/src/io/write_vtk.jl
index d1146e67b1..ff2f2f85eb 100644
--- a/src/io/write_vtk.jl
+++ b/src/io/write_vtk.jl
@@ -410,11 +410,13 @@ function write2vtk!(vtk, v, u, t, system::OpenBoundarySystem)
                        for particle in eachparticle(system)]
 
     if system.cache.calculate_flow_rate
+        Q_total = zero(eltype(system))
         for i in eachindex(system.cache.boundary_zones_flow_rate)
-            vtk["Q_$i"] = system.cache.boundary_zones_flow_rate[i]
+            vtk["Q_$i"] = system.cache.boundary_zones_flow_rate[i][]
+            Q_total += system.cache.boundary_zones_flow_rate[i][]
         end
 
-        vtk["Q_total"] = sum(system.cache.boundary_zones_flow_rate)
+        vtk["Q_total"] = Q_total
     end
 
     return vtk
diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl
index 2c4c5c6c0c..12cde7de5e 100644
--- a/src/schemes/boundary/open_boundary/system.jl
+++ b/src/schemes/boundary/open_boundary/system.jl
@@ -147,7 +147,8 @@ function create_cache_open_boundary(boundary_model, fluid_system, initial_condit
              velocity_reference_values=velocity_reference_values, calculate_flow_rate)
 
     if calculate_flow_rate
-        boundary_zones_flow_rate = zeros(ELTYPE, length(boundary_zones))
+        boundary_zones_flow_rate = ntuple(i -> Ref(zero(ELTYPE)),
+                                          Val(length(boundary_zones)))
         cache = (; boundary_zones_flow_rate, cache...)
     end
 
@@ -329,23 +330,28 @@ end
 
 update_open_boundary_eachstep!(system, v_ode, u_ode, semi, t, integrator) = system
 
-function calculate_flow_rate!(system, v, u, v_ode, u_ode, semi)
-    system.cache.calculate_flow_rate || return system
+function calculate_flow_rate!(system::OpenBoundarySystem{<:Any, ELTYPE, NDIMS}, v, u, v_ode,
+                              u_ode, semi) where {ELTYPE, NDIMS}
+    (; calculate_flow_rate, boundary_zones_flow_rate) = system.cache
+    calculate_flow_rate || return system
 
     for boundary_zone in system.boundary_zones
         interpolate_velocity!(system, boundary_zone, v, u, v_ode, u_ode, semi)
     end
 
-    foreach_enumerate(system.boundary_zones) do (zone_id, boundary_zone)
+    foreach_enumerate(boundary_zones_flow_rate) do (zone_id, boundary_zone_flow_rate)
+        boundary_zone = system.boundary_zones[zone_id]
         (; face_normal) = boundary_zone
         (; sample_velocity, dA) = boundary_zone.cache
+
         # Compute volumetric flow rate: Q = ∫ v ⋅ n dA
-        current_flow_rate = sum(axes(sample_velocity, 2)) do point
-            vn = dot(current_velocity(sample_velocity, system, point), -face_normal)
+        velocities = reinterpret(reshape, SVector{NDIMS, ELTYPE}, sample_velocity)
+        current_flow_rate = sum(velocities) do velocity
+            vn = dot(velocity, -face_normal)
             return vn * dA
         end
 
-        system.cache.boundary_zones_flow_rate[zone_id] = current_flow_rate
+        boundary_zone_flow_rate[] = current_flow_rate
     end
 
     return system

From c1b0abf6536801fa6669908a5f32030b31d693f5 Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Wed, 10 Dec 2025 14:46:11 +0100
Subject: [PATCH 13/62] add checks

---
 .../boundary/open_boundary/boundary_zones.jl  | 54 ++++++++++++++++---
 src/schemes/boundary/open_boundary/system.jl  | 14 +++--
 2 files changed, 56 insertions(+), 12 deletions(-)

diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl
index 2bde137093..bfefefca4c 100644
--- a/src/schemes/boundary/open_boundary/boundary_zones.jl
+++ b/src/schemes/boundary/open_boundary/boundary_zones.jl
@@ -263,7 +263,8 @@ function BoundaryZone(; boundary_face, face_normal, density, particle_spacing,
         ic.velocity .= stack(velocity_ref.(coordinates_svector, 0))
     end
 
-    cache = (; create_cache_boundary_zone(ic, sample_points)...)
+    cache = (;
+             create_cache_boundary_zone(ic, boundary_face, face_normal_, sample_points)...)
 
     return BoundaryZone(ic, spanning_set_, zone_origin, zone_width,
                         flow_direction, face_normal_, rest_pressure, reference_values,
@@ -300,18 +301,55 @@ function Base.show(io::IO, ::MIME"text/plain", boundary_zone::BoundaryZone)
         summary_line(io, "boundary type", boundary_type_name(boundary_zone))
         summary_line(io, "#particles", nparticles(boundary_zone.initial_condition))
         summary_line(io, "width", round(boundary_zone.zone_width, digits=6))
+        if hasproperty(boundary_zone.cache, :cross_sectional_area)
+            summary_line(io, "cross sectional area",
+                         boundary_zone.cache.cross_sectional_area)
+        end
         summary_footer(io)
     end
 end
 
-create_cache_boundary_zone(initial_condition, sample_points::Nothing) = (;)
+function create_cache_boundary_zone(initial_condition, boundary_face, face_normal,
+                                    sample_points::Nothing)
+    return (; sample_points)
+end
+
+function create_cache_boundary_zone(initial_condition, boundary_face, face_normal,
+                                    sample_points)
+    (; particle_spacing) = initial_condition
+    area_increment = particle_spacing^(ndims(initial_condition) - 1)
+    if sample_points === :default
+        sample_points_ = extrude_geometry(boundary_face; particle_spacing, density=Inf,
+                                          direction=(-face_normal), n_extrude=1).coordinates
+    else
+        if !(sample_points isa Matrix && size(sample_points, 1) == ndims(initial_condition))
+            throw(ArgumentError("`sample_points` must be a matrix with " *
+                                "`ndims(initial_condition)` rows"))
+        end
+
+        # TODO: Check if regular grid?
+
+        sample_points_ = sample_points
+    end
+
+    discrete_face_area = area_increment * size(sample_points_, 2)
+
+    if ndims(initial_condition) == 3
+        v1, v2, v3 = boundary_face
+        face_area = norm(cross(v2 - v1, v3 - v1))
+    elseif ndims(initial_condition) == 2
+        v1, v2 = boundary_face
+        face_area = norm(v2 - v1)
+    end
+
+    if discrete_face_area > face_area
+        @warn "The sampled area of the boundary face " *
+              "($(discrete_face_area)) is larger than the actual face area ($(face_area)). "
+    end
 
-function create_cache_boundary_zone(initial_condition, sample_points::Matrix)
-    # TODO: Check matrix (ndims etc.)
-    shepard_coefficient = zeros(eltype(initial_condition), axes(sample_points, 2))
-    dA = initial_condition.particle_spacing
-    return (; sample_points=sample_points, sample_velocity=copy(sample_points),
-            shepard_coefficient, dA)
+    return (; sample_points=sample_points_, sample_velocity=copy(sample_points_),
+            shepard_coefficient=zeros(eltype(initial_condition), size(sample_points_, 2)),
+            area_increment, cross_sectional_area=discrete_face_area)
 end
 
 function set_up_boundary_zone(boundary_face, face_normal, density, particle_spacing,
diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl
index 12cde7de5e..4bed038600 100644
--- a/src/schemes/boundary/open_boundary/system.jl
+++ b/src/schemes/boundary/open_boundary/system.jl
@@ -147,6 +147,12 @@ function create_cache_open_boundary(boundary_model, fluid_system, initial_condit
              velocity_reference_values=velocity_reference_values, calculate_flow_rate)
 
     if calculate_flow_rate
+        if any(zone -> isnothing(zone.cache.sample_points), boundary_zones)
+            throw(ArgumentError("`sample_points` must be specified for all boundary zones when " *
+                                "`calculate_flow_rate` is true.\n" *
+                                "Use `sample_points=:default` to automatically generate sample points."))
+        end
+
         boundary_zones_flow_rate = ntuple(i -> Ref(zero(ELTYPE)),
                                           Val(length(boundary_zones)))
         cache = (; boundary_zones_flow_rate, cache...)
@@ -332,8 +338,8 @@ update_open_boundary_eachstep!(system, v_ode, u_ode, semi, t, integrator) = syst
 
 function calculate_flow_rate!(system::OpenBoundarySystem{<:Any, ELTYPE, NDIMS}, v, u, v_ode,
                               u_ode, semi) where {ELTYPE, NDIMS}
-    (; calculate_flow_rate, boundary_zones_flow_rate) = system.cache
-    calculate_flow_rate || return system
+    system.cache.calculate_flow_rate || return system
+    (; boundary_zones_flow_rate) = system.cache
 
     for boundary_zone in system.boundary_zones
         interpolate_velocity!(system, boundary_zone, v, u, v_ode, u_ode, semi)
@@ -342,13 +348,13 @@ function calculate_flow_rate!(system::OpenBoundarySystem{<:Any, ELTYPE, NDIMS},
     foreach_enumerate(boundary_zones_flow_rate) do (zone_id, boundary_zone_flow_rate)
         boundary_zone = system.boundary_zones[zone_id]
         (; face_normal) = boundary_zone
-        (; sample_velocity, dA) = boundary_zone.cache
+        (; sample_velocity, area_increment) = boundary_zone.cache
 
         # Compute volumetric flow rate: Q = ∫ v ⋅ n dA
         velocities = reinterpret(reshape, SVector{NDIMS, ELTYPE}, sample_velocity)
         current_flow_rate = sum(velocities) do velocity
             vn = dot(velocity, -face_normal)
-            return vn * dA
+            return vn * area_increment
         end
 
         boundary_zone_flow_rate[] = current_flow_rate

From eeebcc22a371b01d4d51d50ddf90ad31c3ae7a40 Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Wed, 10 Dec 2025 15:14:33 +0100
Subject: [PATCH 14/62] add docs

---
 src/schemes/boundary/open_boundary/boundary_zones.jl | 6 +++++-
 src/schemes/boundary/open_boundary/system.jl         | 4 +++-
 2 files changed, 8 insertions(+), 2 deletions(-)

diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl
index bfefefca4c..043beb8196 100644
--- a/src/schemes/boundary/open_boundary/boundary_zones.jl
+++ b/src/schemes/boundary/open_boundary/boundary_zones.jl
@@ -81,6 +81,10 @@ There are three ways to specify the actual shape of the boundary zone:
                        Per default it is set to zero (assuming a gauge pressure system).
                        - For `EntropicallyDampedSPHSystem`: Use the initial pressure from the `InitialCondition`
                        - For `WeaklyCompressibleSPHSystem`: Use the background pressure from the equation of state
+- `sample_points=:default`: Either `:default` to automatically generate sample points on the boundary face (default),
+                            or a matrix of dimensions `(ndims, n_points)` containing sample points
+                            on the boundary face used to compute the volumetric flow rate.
+                            Set to `nothing` to skip sampling.
 
 !!! note "Note"
     The reference values (`reference_velocity`, `reference_pressure`, `reference_density`)
@@ -168,7 +172,7 @@ end
 function BoundaryZone(; boundary_face, face_normal, density, particle_spacing,
                       initial_condition=nothing, extrude_geometry=nothing,
                       open_boundary_layers::Integer, average_inflow_velocity=true,
-                      boundary_type=BidirectionalFlow(), sample_points=nothing,
+                      boundary_type=BidirectionalFlow(), sample_points=:default,
                       rest_pressure=zero(eltype(density)),
                       reference_density=nothing, reference_pressure=nothing,
                       reference_velocity=nothing)
diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl
index 4bed038600..36889bc0e2 100644
--- a/src/schemes/boundary/open_boundary/system.jl
+++ b/src/schemes/boundary/open_boundary/system.jl
@@ -1,7 +1,7 @@
 @doc raw"""
     OpenBoundarySystem(boundary_zone::BoundaryZone;
                        fluid_system::AbstractFluidSystem, buffer_size::Integer,
-                       boundary_model)
+                       boundary_model, calculate_flow_rate=false)
 
 Open boundary system for in- and outflow particles.
 
@@ -11,6 +11,8 @@ Open boundary system for in- and outflow particles.
 # Keywords
 - `fluid_system`: The corresponding fluid system
 - `boundary_model`: Boundary model (see [Open Boundary Models](@ref open_boundary_models))
+- `calculate_flow_rate=false`: Set to `true` to calculate the volumetric flow rate through each boundary zone.
+                               Default is `false`.
 - `buffer_size`: Number of buffer particles.
 - `pressure_acceleration`: Pressure acceleration formulation for the system. Required only
                            when using [`BoundaryModelDynamicalPressureZhang`](@ref).

From d8b5b0091aba86fabd02d5b014d697674cbdcc7b Mon Sep 17 00:00:00 2001
From: Marcel Schurer <marcel.schurer@outlook.de>
Date: Wed, 10 Dec 2025 15:50:22 +0100
Subject: [PATCH 15/62] Corrected failing doctest output

---
 src/io/read_vtk.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/io/read_vtk.jl b/src/io/read_vtk.jl
index 4189ea474a..f98d35d10c 100644
--- a/src/io/read_vtk.jl
+++ b/src/io/read_vtk.jl
@@ -15,7 +15,7 @@ Load VTK file and convert data to a NamedTuple.
     This is an experimental feature and may change in any future releases.
 
 # Example
-```jldoctest; output = false
+```jldoctest; output = false, filter = r"density = \\[.*\\]|pressure = \\[.*\\]|mass = \\[.*\\]|velocity = \\[.*\\]|coordinates = \\[.*\\]"
 # Create a rectangular shape
 rectangular = RectangularShape(0.1, (10, 10), (0, 0), density=1.5, velocity=(1.0, -2.0),
                                pressure=1000.0)
@@ -29,7 +29,7 @@ data = vtk2trixi(joinpath("out", "rectangular.vtu");
                  my_custom_quantity="my_custom_quantity")
 
 # output
-
+(particle_spacing = 0.1, density = [...], time = 0.0, pressure = [...], mass = [...], my_custom_quantity = 3.0, velocity = [...], coordinates = [...])
 ```
 """
 function vtk2trixi(file; custom_quantities...)

From d21a5492311db0c13e75aa7fa77866cef49814d5 Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Wed, 10 Dec 2025 17:29:05 +0100
Subject: [PATCH 16/62] adapt pressure model

---
 .../boundary/open_boundary/pressure_model.jl  | 39 +++++--------------
 src/schemes/boundary/open_boundary/system.jl  |  4 ++
 src/util.jl                                   | 14 +++++++
 .../boundary/open_boundary/pressure_model.jl  |  6 +--
 4 files changed, 30 insertions(+), 33 deletions(-)

diff --git a/src/schemes/boundary/open_boundary/pressure_model.jl b/src/schemes/boundary/open_boundary/pressure_model.jl
index 81ce8f9ff4..7e180396a1 100644
--- a/src/schemes/boundary/open_boundary/pressure_model.jl
+++ b/src/schemes/boundary/open_boundary/pressure_model.jl
@@ -67,51 +67,32 @@ function update_pressure_model!(system, v, u, semi, dt)
 
     if any(pm -> isa(pm, AbstractPressureModel), system.cache.pressure_reference_values)
         @trixi_timeit timer() "update pressure model" begin
-            calculate_flow_rate_and_pressure!(system, v, u, dt)
+            calculate_pressure!(system, dt)
         end
     end
 
     return system
 end
 
-function calculate_flow_rate_and_pressure!(system, v, u, dt)
-    (; pressure_reference_values) = system.cache
-    foreach_enumerate(pressure_reference_values) do (zone_id, pressure_model)
-        boundary_zone = system.boundary_zones[zone_id]
-        calculate_flow_rate_and_pressure!(pressure_model, system, boundary_zone, v, u, dt)
+function calculate_pressure!(system, dt)
+    (; pressure_reference_values, boundary_zones_flow_rate) = system.cache
+
+    foreach_noalloc(pressure_reference_values,
+                    boundary_zones_flow_rate) do (pressure_model, flow_rate)
+        calculate_pressure!(pressure_model, system, flow_rate[], dt)
     end
 
     return system
 end
 
-function calculate_flow_rate_and_pressure!(pressure_model, system, boundary_zone, v, u, dt)
+function calculate_pressure!(pressure_model, system, current_flow_rate, dt)
     return pressure_model
 end
 
-function calculate_flow_rate_and_pressure!(pressure_model::RCRWindkesselModel, system,
-                                           boundary_zone, v, u, dt)
-    (; particle_spacing) = system.initial_condition
+function calculate_pressure!(pressure_model::RCRWindkesselModel, system,
+                             current_flow_rate, dt)
     (; characteristic_resistance, peripheral_resistance, compliance,
      flow_rate, pressure) = pressure_model
-    (; face_normal) = boundary_zone
-
-    # Find particles within the current boundary zone
-    candidates = findall(particle -> boundary_zone ==
-                                     current_boundary_zone(system, particle),
-                         each_integrated_particle(system))
-
-    # Assuming negligible transverse velocity gradients within the boundary zone,
-    # the full area of the zone is taken as the representative cross-sectional
-    # area for volumetric flow-rate estimation.
-    cross_sectional_area = length(candidates) * particle_spacing^(ndims(system) - 1)
-
-    # Division inside the `sum` closure to maintain GPU compatibility
-    velocity_avg = sum(candidates) do particle
-        return dot(current_velocity(v, system, particle), -face_normal) / length(candidates)
-    end
-
-    # Compute volumetric flow rate: Q = v * A
-    current_flow_rate = velocity_avg * cross_sectional_area
 
     previous_pressure = pressure[]
     previous_flow_rate = flow_rate[]
diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl
index 36889bc0e2..7038f48534 100644
--- a/src/schemes/boundary/open_boundary/system.jl
+++ b/src/schemes/boundary/open_boundary/system.jl
@@ -144,6 +144,10 @@ function create_cache_open_boundary(boundary_model, fluid_system, initial_condit
     density_reference_values = map(ref -> ref.reference_density, reference_values)
     velocity_reference_values = map(ref -> ref.reference_velocity, reference_values)
 
+    if any(pr -> isa(pr, AbstractPressureModel), pressure_reference_values)
+        calculate_flow_rate = true
+    end
+
     cache = (; pressure_reference_values=pressure_reference_values,
              density_reference_values=density_reference_values,
              velocity_reference_values=velocity_reference_values, calculate_flow_rate)
diff --git a/src/util.jl b/src/util.jl
index 18389cd44b..be032a8b95 100644
--- a/src/util.jl
+++ b/src/util.jl
@@ -11,6 +11,20 @@ end
 
 @inline foreach_noalloc(func, collection::Tuple{}) = nothing
 
+@inline function foreach_noalloc(func, collection1, collection2)
+    element1 = first(collection1)
+    remaining_collection1 = Base.tail(collection1)
+    element2 = first(collection2)
+    remaining_collection2 = Base.tail(collection2)
+
+    func((element1, element2))
+
+    # Process remaining collection
+    foreach_noalloc(func, remaining_collection1, remaining_collection2)
+end
+
+@inline foreach_noalloc(func, collection1::Tuple{}, collection2::Tuple{}) = nothing
+
 # Same as `foreach(enumerate(something))`, but without allocations.
 # Note that compile times may increase if this is used with big tuples.
 @inline foreach_enumerate(func, collection) = foreach_enumerate(func, collection, 1)
diff --git a/test/schemes/boundary/open_boundary/pressure_model.jl b/test/schemes/boundary/open_boundary/pressure_model.jl
index f35c2e867f..c0e1baffe7 100644
--- a/test/schemes/boundary/open_boundary/pressure_model.jl
+++ b/test/schemes/boundary/open_boundary/pressure_model.jl
@@ -81,15 +81,13 @@
                                         boundary_model=nothing,
                                         fluid_system=FluidSystemMockRCR(nothing, nothing))
             system.boundary_zone_indices .= 1
-
-            u = system.initial_condition.coordinates
             v = system.initial_condition.velocity
 
             times = collect(tspan[1]:dt:tspan[2])
             p_calculated = empty(times)
             for t in times
-                v[1, :] .= func(t)
-                TrixiParticles.calculate_flow_rate_and_pressure!(system, v, u, dt)
+                system.cache.boundary_zones_flow_rate[1][] = func(t)
+                TrixiParticles.calculate_pressure!(system, dt)
 
                 # Store only values after the seventh cycle
                 if t >= 7T

From b4020c8df060d51bb3f012ef4c977aabe14c7834 Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Thu, 11 Dec 2025 08:07:32 +0100
Subject: [PATCH 17/62] use `foreach_no_alloc`

---
 examples/fluid/poiseuille_flow_2d.jl         |  3 +-
 src/schemes/boundary/open_boundary/system.jl | 40 ++++++++++----------
 src/util.jl                                  | 14 +++++++
 3 files changed, 37 insertions(+), 20 deletions(-)

diff --git a/examples/fluid/poiseuille_flow_2d.jl b/examples/fluid/poiseuille_flow_2d.jl
index dd324a5bca..1afacd81a8 100644
--- a/examples/fluid/poiseuille_flow_2d.jl
+++ b/examples/fluid/poiseuille_flow_2d.jl
@@ -29,7 +29,7 @@ open_boundary_layers = 10
 
 # ==========================================================================================
 # ==== Experiment Setup
-tspan = (0.0, 2.0)
+tspan = (0.0, 1.0)
 wcsph = true
 
 domain_size = (flow_length, wall_distance)
@@ -136,6 +136,7 @@ outflow = BoundaryZone(; boundary_face=face_out, face_normal=(.-(flow_direction)
 
 open_boundary = OpenBoundarySystem(inflow, outflow; fluid_system,
                                    boundary_model=open_boundary_model,
+                                   calculate_flow_rate=true,
                                    buffer_size=n_buffer_particles)
 
 # ==========================================================================================
diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl
index 36889bc0e2..648a97cb03 100644
--- a/src/schemes/boundary/open_boundary/system.jl
+++ b/src/schemes/boundary/open_boundary/system.jl
@@ -318,7 +318,7 @@ function update_open_boundary_eachstep!(system::OpenBoundarySystem, v_ode, u_ode
 
         # This must be called before `update_pressure_model!` and `check_domain!`
         # to ensure quantities remain consistent with the current simulation state.
-        calculate_flow_rate!(system, v, u, v_ode, u_ode, semi)
+        calculate_flow_rate!(system, v_ode, u_ode, semi)
 
         @trixi_timeit timer() "check domain" check_domain!(system, v, u, v_ode, u_ode, semi)
 
@@ -338,28 +338,30 @@ end
 
 update_open_boundary_eachstep!(system, v_ode, u_ode, semi, t, integrator) = system
 
-function calculate_flow_rate!(system::OpenBoundarySystem{<:Any, ELTYPE, NDIMS}, v, u, v_ode,
-                              u_ode, semi) where {ELTYPE, NDIMS}
+function calculate_flow_rate!(system::OpenBoundarySystem{<:Any, ELTYPE, NDIMS},
+                              v_ode, u_ode, semi) where {ELTYPE, NDIMS}
     system.cache.calculate_flow_rate || return system
-    (; boundary_zones_flow_rate) = system.cache
 
-    for boundary_zone in system.boundary_zones
-        interpolate_velocity!(system, boundary_zone, v, u, v_ode, u_ode, semi)
-    end
+    @trixi_timeit timer() "flow rate calculation" begin
+        (; boundary_zones) = system
+        (; boundary_zones_flow_rate) = system.cache
 
-    foreach_enumerate(boundary_zones_flow_rate) do (zone_id, boundary_zone_flow_rate)
-        boundary_zone = system.boundary_zones[zone_id]
-        (; face_normal) = boundary_zone
-        (; sample_velocity, area_increment) = boundary_zone.cache
+        foreach_noalloc(boundary_zones,
+                        boundary_zones_flow_rate) do (boundary_zone, flow_rate)
+            (; face_normal) = boundary_zone
+            (; sample_velocity, area_increment) = boundary_zone.cache
 
-        # Compute volumetric flow rate: Q = ∫ v ⋅ n dA
-        velocities = reinterpret(reshape, SVector{NDIMS, ELTYPE}, sample_velocity)
-        current_flow_rate = sum(velocities) do velocity
-            vn = dot(velocity, -face_normal)
-            return vn * area_increment
-        end
+            interpolate_velocity!(system, boundary_zone, v_ode, u_ode, semi)
 
-        boundary_zone_flow_rate[] = current_flow_rate
+            # Compute volumetric flow rate: Q = ∫ v ⋅ n dA
+            velocities = reinterpret(reshape, SVector{NDIMS, ELTYPE}, sample_velocity)
+            current_flow_rate = sum(velocities) do velocity
+                vn = dot(velocity, -face_normal)
+                return vn * area_increment
+            end
+
+            flow_rate[] = current_flow_rate
+        end
     end
 
     return system
@@ -613,7 +615,7 @@ end
     return system
 end
 
-function interpolate_velocity!(system::OpenBoundarySystem, boundary_zone, v, u,
+function interpolate_velocity!(system::OpenBoundarySystem, boundary_zone,
                                v_ode, u_ode, semi)
     (; sample_points, sample_velocity, shepard_coefficient) = boundary_zone.cache
     smoothing_length = initial_smoothing_length(system)
diff --git a/src/util.jl b/src/util.jl
index 18389cd44b..be032a8b95 100644
--- a/src/util.jl
+++ b/src/util.jl
@@ -11,6 +11,20 @@ end
 
 @inline foreach_noalloc(func, collection::Tuple{}) = nothing
 
+@inline function foreach_noalloc(func, collection1, collection2)
+    element1 = first(collection1)
+    remaining_collection1 = Base.tail(collection1)
+    element2 = first(collection2)
+    remaining_collection2 = Base.tail(collection2)
+
+    func((element1, element2))
+
+    # Process remaining collection
+    foreach_noalloc(func, remaining_collection1, remaining_collection2)
+end
+
+@inline foreach_noalloc(func, collection1::Tuple{}, collection2::Tuple{}) = nothing
+
 # Same as `foreach(enumerate(something))`, but without allocations.
 # Note that compile times may increase if this is used with big tuples.
 @inline foreach_enumerate(func, collection) = foreach_enumerate(func, collection, 1)

From 55be6a601973cc7a9eff94cafe9d6b0a5114feaf Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Thu, 11 Dec 2025 08:28:07 +0100
Subject: [PATCH 18/62] fix unit tests

---
 .../boundary/open_boundary/boundary_zones.jl  |  2 +-
 .../boundary/open_boundary/boundary_zone.jl   | 62 ++++++++++++++-----
 2 files changed, 49 insertions(+), 15 deletions(-)

diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl
index 043beb8196..55b86be737 100644
--- a/src/schemes/boundary/open_boundary/boundary_zones.jl
+++ b/src/schemes/boundary/open_boundary/boundary_zones.jl
@@ -346,7 +346,7 @@ function create_cache_boundary_zone(initial_condition, boundary_face, face_norma
         face_area = norm(v2 - v1)
     end
 
-    if discrete_face_area > face_area
+    if discrete_face_area > (face_area + eps(face_area))
         @warn "The sampled area of the boundary face " *
               "($(discrete_face_area)) is larger than the actual face area ($(face_area)). "
     end
diff --git a/test/schemes/boundary/open_boundary/boundary_zone.jl b/test/schemes/boundary/open_boundary/boundary_zone.jl
index 619c048bbe..3e1fa13118 100644
--- a/test/schemes/boundary/open_boundary/boundary_zone.jl
+++ b/test/schemes/boundary/open_boundary/boundary_zone.jl
@@ -1,7 +1,7 @@
 @testset verbose=true "Boundary Zone" begin
     @testset "`show`" begin
         inflow = BoundaryZone(; boundary_face=([0.0, 0.0], [0.0, 1.0]),
-                              particle_spacing=0.05,
+                              particle_spacing=0.05, sample_points=nothing,
                               face_normal=(1.0, 0.0), density=1.0,
                               reference_density=0.0,
                               reference_pressure=0.0,
@@ -22,7 +22,7 @@
         @test repr("text/plain", inflow) == show_box
 
         outflow = BoundaryZone(; boundary_face=([0.0, 0.0], [0.0, 1.0]),
-                               particle_spacing=0.05,
+                               particle_spacing=0.05, sample_points=nothing,
                                reference_density=0.0,
                                reference_pressure=0.0,
                                reference_velocity=[0.0, 0.0],
@@ -41,6 +41,48 @@
         └──────────────────────────────────────────────────────────────────────────────────────────────────┘"""
 
         @test repr("text/plain", outflow) == show_box
+
+        bidirectional = BoundaryZone(; boundary_face=([0.0, 0.0], [0.0, 1.0]),
+                                     particle_spacing=0.05, sample_points=nothing,
+                                     reference_density=0.0,
+                                     reference_pressure=0.0,
+                                     reference_velocity=[0.0, 0.0],
+                                     face_normal=(1.0, 0.0), density=1.0,
+                                     open_boundary_layers=4)
+
+        show_compact = "BoundaryZone() with 80 particles"
+        @test repr(bidirectional) == show_compact
+        show_box = """
+        ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
+        │ BoundaryZone                                                                                     │
+        │ ════════════                                                                                     │
+        │ boundary type: ………………………………………… bidirectional_flow                                               │
+        │ #particles: ………………………………………………… 80                                                               │
+        │ width: ……………………………………………………………… 0.2                                                              │
+        └──────────────────────────────────────────────────────────────────────────────────────────────────┘"""
+
+        @test repr("text/plain", bidirectional) == show_box
+
+        zone = BoundaryZone(; boundary_face=([0.0, 0.0], [1 / sqrt(2), 1 / sqrt(2)]),
+                            particle_spacing=0.05, reference_density=0.0,
+                            reference_pressure=0.0,
+                            reference_velocity=[0.0, 0.0],
+                            face_normal=normalize([-1.0, 1.0]), density=1.0,
+                            open_boundary_layers=4)
+
+        show_compact = "BoundaryZone() with 80 particles"
+        @test repr(zone) == show_compact
+        show_box = """
+        ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
+        │ BoundaryZone                                                                                     │
+        │ ════════════                                                                                     │
+        │ boundary type: ………………………………………… bidirectional_flow                                               │
+        │ #particles: ………………………………………………… 80                                                               │
+        │ width: ……………………………………………………………… 0.2                                                              │
+        │ cross sectional area: ……………………… 1.0                                                              │
+        └──────────────────────────────────────────────────────────────────────────────────────────────────┘"""
+
+        @test repr("text/plain", zone) == show_box
     end
 
     @testset verbose=true "Illegal Inputs" begin
@@ -127,9 +169,7 @@
                     outflow
                 ]
 
-                @testset verbose=true "$(TrixiParticles.boundary_type_name(boundary_zone))" for boundary_zone in
-                                                                                                boundary_zones
-
+                @testset verbose=true "$(TrixiParticles.boundary_type_name(boundary_zone))" for boundary_zone in boundary_zones
                     zone_width = open_boundary_layers *
                                  boundary_zone.initial_condition.particle_spacing
                     sign_ = (TrixiParticles.boundary_type_name(boundary_zone) == "inflow") ?
@@ -192,9 +232,7 @@
                     outflow
                 ]
 
-                @testset verbose=true "$(TrixiParticles.boundary_type_name(boundary_zone))" for boundary_zone in
-                                                                                                boundary_zones
-
+                @testset verbose=true "$(TrixiParticles.boundary_type_name(boundary_zone))" for boundary_zone in boundary_zones
                     zone_width = open_boundary_layers *
                                  boundary_zone.initial_condition.particle_spacing
                     sign_ = (TrixiParticles.boundary_type_name(boundary_zone) == "inflow") ?
@@ -232,9 +270,7 @@
             outflow
         ]
 
-        @testset verbose=true "$(TrixiParticles.boundary_type_name(boundary_zone))" for boundary_zone in
-                                                                                        boundary_zones
-
+        @testset verbose=true "$(TrixiParticles.boundary_type_name(boundary_zone))" for boundary_zone in boundary_zones
             perturb_ = TrixiParticles.boundary_type_name(boundary_zone) == "inflow" ?
                        sqrt(eps()) :
                        -sqrt(eps())
@@ -280,9 +316,7 @@
             outflow
         ]
 
-        @testset verbose=true "$(TrixiParticles.boundary_type_name(boundary_zone))" for boundary_zone in
-                                                                                        boundary_zones
-
+        @testset verbose=true "$(TrixiParticles.boundary_type_name(boundary_zone))" for boundary_zone in boundary_zones
             perturb_ = TrixiParticles.boundary_type_name(boundary_zone) == "inflow" ?
                        eps() : -eps()
             point4 = boundary_zone.spanning_set[1] + boundary_zone.zone_origin

From a431873732519aaae7b4d748d612e3cea2b5e7be Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Thu, 11 Dec 2025 08:28:07 +0100
Subject: [PATCH 19/62] fix unit tests

---
 examples/fluid/poiseuille_flow_2d.jl          |  2 +-
 .../boundary/open_boundary/boundary_zones.jl  |  2 +-
 .../boundary/open_boundary/boundary_zone.jl   | 62 ++++++++++++++-----
 3 files changed, 50 insertions(+), 16 deletions(-)

diff --git a/examples/fluid/poiseuille_flow_2d.jl b/examples/fluid/poiseuille_flow_2d.jl
index 1afacd81a8..5e57c7c0bd 100644
--- a/examples/fluid/poiseuille_flow_2d.jl
+++ b/examples/fluid/poiseuille_flow_2d.jl
@@ -29,7 +29,7 @@ open_boundary_layers = 10
 
 # ==========================================================================================
 # ==== Experiment Setup
-tspan = (0.0, 1.0)
+tspan = (0.0, 2.0)
 wcsph = true
 
 domain_size = (flow_length, wall_distance)
diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl
index 043beb8196..55b86be737 100644
--- a/src/schemes/boundary/open_boundary/boundary_zones.jl
+++ b/src/schemes/boundary/open_boundary/boundary_zones.jl
@@ -346,7 +346,7 @@ function create_cache_boundary_zone(initial_condition, boundary_face, face_norma
         face_area = norm(v2 - v1)
     end
 
-    if discrete_face_area > face_area
+    if discrete_face_area > (face_area + eps(face_area))
         @warn "The sampled area of the boundary face " *
               "($(discrete_face_area)) is larger than the actual face area ($(face_area)). "
     end
diff --git a/test/schemes/boundary/open_boundary/boundary_zone.jl b/test/schemes/boundary/open_boundary/boundary_zone.jl
index 619c048bbe..3e1fa13118 100644
--- a/test/schemes/boundary/open_boundary/boundary_zone.jl
+++ b/test/schemes/boundary/open_boundary/boundary_zone.jl
@@ -1,7 +1,7 @@
 @testset verbose=true "Boundary Zone" begin
     @testset "`show`" begin
         inflow = BoundaryZone(; boundary_face=([0.0, 0.0], [0.0, 1.0]),
-                              particle_spacing=0.05,
+                              particle_spacing=0.05, sample_points=nothing,
                               face_normal=(1.0, 0.0), density=1.0,
                               reference_density=0.0,
                               reference_pressure=0.0,
@@ -22,7 +22,7 @@
         @test repr("text/plain", inflow) == show_box
 
         outflow = BoundaryZone(; boundary_face=([0.0, 0.0], [0.0, 1.0]),
-                               particle_spacing=0.05,
+                               particle_spacing=0.05, sample_points=nothing,
                                reference_density=0.0,
                                reference_pressure=0.0,
                                reference_velocity=[0.0, 0.0],
@@ -41,6 +41,48 @@
         └──────────────────────────────────────────────────────────────────────────────────────────────────┘"""
 
         @test repr("text/plain", outflow) == show_box
+
+        bidirectional = BoundaryZone(; boundary_face=([0.0, 0.0], [0.0, 1.0]),
+                                     particle_spacing=0.05, sample_points=nothing,
+                                     reference_density=0.0,
+                                     reference_pressure=0.0,
+                                     reference_velocity=[0.0, 0.0],
+                                     face_normal=(1.0, 0.0), density=1.0,
+                                     open_boundary_layers=4)
+
+        show_compact = "BoundaryZone() with 80 particles"
+        @test repr(bidirectional) == show_compact
+        show_box = """
+        ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
+        │ BoundaryZone                                                                                     │
+        │ ════════════                                                                                     │
+        │ boundary type: ………………………………………… bidirectional_flow                                               │
+        │ #particles: ………………………………………………… 80                                                               │
+        │ width: ……………………………………………………………… 0.2                                                              │
+        └──────────────────────────────────────────────────────────────────────────────────────────────────┘"""
+
+        @test repr("text/plain", bidirectional) == show_box
+
+        zone = BoundaryZone(; boundary_face=([0.0, 0.0], [1 / sqrt(2), 1 / sqrt(2)]),
+                            particle_spacing=0.05, reference_density=0.0,
+                            reference_pressure=0.0,
+                            reference_velocity=[0.0, 0.0],
+                            face_normal=normalize([-1.0, 1.0]), density=1.0,
+                            open_boundary_layers=4)
+
+        show_compact = "BoundaryZone() with 80 particles"
+        @test repr(zone) == show_compact
+        show_box = """
+        ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
+        │ BoundaryZone                                                                                     │
+        │ ════════════                                                                                     │
+        │ boundary type: ………………………………………… bidirectional_flow                                               │
+        │ #particles: ………………………………………………… 80                                                               │
+        │ width: ……………………………………………………………… 0.2                                                              │
+        │ cross sectional area: ……………………… 1.0                                                              │
+        └──────────────────────────────────────────────────────────────────────────────────────────────────┘"""
+
+        @test repr("text/plain", zone) == show_box
     end
 
     @testset verbose=true "Illegal Inputs" begin
@@ -127,9 +169,7 @@
                     outflow
                 ]
 
-                @testset verbose=true "$(TrixiParticles.boundary_type_name(boundary_zone))" for boundary_zone in
-                                                                                                boundary_zones
-
+                @testset verbose=true "$(TrixiParticles.boundary_type_name(boundary_zone))" for boundary_zone in boundary_zones
                     zone_width = open_boundary_layers *
                                  boundary_zone.initial_condition.particle_spacing
                     sign_ = (TrixiParticles.boundary_type_name(boundary_zone) == "inflow") ?
@@ -192,9 +232,7 @@
                     outflow
                 ]
 
-                @testset verbose=true "$(TrixiParticles.boundary_type_name(boundary_zone))" for boundary_zone in
-                                                                                                boundary_zones
-
+                @testset verbose=true "$(TrixiParticles.boundary_type_name(boundary_zone))" for boundary_zone in boundary_zones
                     zone_width = open_boundary_layers *
                                  boundary_zone.initial_condition.particle_spacing
                     sign_ = (TrixiParticles.boundary_type_name(boundary_zone) == "inflow") ?
@@ -232,9 +270,7 @@
             outflow
         ]
 
-        @testset verbose=true "$(TrixiParticles.boundary_type_name(boundary_zone))" for boundary_zone in
-                                                                                        boundary_zones
-
+        @testset verbose=true "$(TrixiParticles.boundary_type_name(boundary_zone))" for boundary_zone in boundary_zones
             perturb_ = TrixiParticles.boundary_type_name(boundary_zone) == "inflow" ?
                        sqrt(eps()) :
                        -sqrt(eps())
@@ -280,9 +316,7 @@
             outflow
         ]
 
-        @testset verbose=true "$(TrixiParticles.boundary_type_name(boundary_zone))" for boundary_zone in
-                                                                                        boundary_zones
-
+        @testset verbose=true "$(TrixiParticles.boundary_type_name(boundary_zone))" for boundary_zone in boundary_zones
             perturb_ = TrixiParticles.boundary_type_name(boundary_zone) == "inflow" ?
                        eps() : -eps()
             point4 = boundary_zone.spanning_set[1] + boundary_zone.zone_origin

From 9339bcbe7a217f84e59a819ca91b17e8fbefb647 Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Thu, 11 Dec 2025 08:28:07 +0100
Subject: [PATCH 20/62] fix unit tests

---
 examples/fluid/poiseuille_flow_2d.jl          |  2 +-
 .../boundary/open_boundary/boundary_zones.jl  |  2 +-
 .../boundary/open_boundary/boundary_zone.jl   | 46 ++++++++++++++++++-
 3 files changed, 46 insertions(+), 4 deletions(-)

diff --git a/examples/fluid/poiseuille_flow_2d.jl b/examples/fluid/poiseuille_flow_2d.jl
index 1afacd81a8..5e57c7c0bd 100644
--- a/examples/fluid/poiseuille_flow_2d.jl
+++ b/examples/fluid/poiseuille_flow_2d.jl
@@ -29,7 +29,7 @@ open_boundary_layers = 10
 
 # ==========================================================================================
 # ==== Experiment Setup
-tspan = (0.0, 1.0)
+tspan = (0.0, 2.0)
 wcsph = true
 
 domain_size = (flow_length, wall_distance)
diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl
index 043beb8196..55b86be737 100644
--- a/src/schemes/boundary/open_boundary/boundary_zones.jl
+++ b/src/schemes/boundary/open_boundary/boundary_zones.jl
@@ -346,7 +346,7 @@ function create_cache_boundary_zone(initial_condition, boundary_face, face_norma
         face_area = norm(v2 - v1)
     end
 
-    if discrete_face_area > face_area
+    if discrete_face_area > (face_area + eps(face_area))
         @warn "The sampled area of the boundary face " *
               "($(discrete_face_area)) is larger than the actual face area ($(face_area)). "
     end
diff --git a/test/schemes/boundary/open_boundary/boundary_zone.jl b/test/schemes/boundary/open_boundary/boundary_zone.jl
index 619c048bbe..8542fd554d 100644
--- a/test/schemes/boundary/open_boundary/boundary_zone.jl
+++ b/test/schemes/boundary/open_boundary/boundary_zone.jl
@@ -1,7 +1,7 @@
 @testset verbose=true "Boundary Zone" begin
     @testset "`show`" begin
         inflow = BoundaryZone(; boundary_face=([0.0, 0.0], [0.0, 1.0]),
-                              particle_spacing=0.05,
+                              particle_spacing=0.05, sample_points=nothing,
                               face_normal=(1.0, 0.0), density=1.0,
                               reference_density=0.0,
                               reference_pressure=0.0,
@@ -22,7 +22,7 @@
         @test repr("text/plain", inflow) == show_box
 
         outflow = BoundaryZone(; boundary_face=([0.0, 0.0], [0.0, 1.0]),
-                               particle_spacing=0.05,
+                               particle_spacing=0.05, sample_points=nothing,
                                reference_density=0.0,
                                reference_pressure=0.0,
                                reference_velocity=[0.0, 0.0],
@@ -41,6 +41,48 @@
         └──────────────────────────────────────────────────────────────────────────────────────────────────┘"""
 
         @test repr("text/plain", outflow) == show_box
+
+        bidirectional = BoundaryZone(; boundary_face=([0.0, 0.0], [0.0, 1.0]),
+                                     particle_spacing=0.05, sample_points=nothing,
+                                     reference_density=0.0,
+                                     reference_pressure=0.0,
+                                     reference_velocity=[0.0, 0.0],
+                                     face_normal=(1.0, 0.0), density=1.0,
+                                     open_boundary_layers=4)
+
+        show_compact = "BoundaryZone() with 80 particles"
+        @test repr(bidirectional) == show_compact
+        show_box = """
+        ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
+        │ BoundaryZone                                                                                     │
+        │ ════════════                                                                                     │
+        │ boundary type: ………………………………………… bidirectional_flow                                               │
+        │ #particles: ………………………………………………… 80                                                               │
+        │ width: ……………………………………………………………… 0.2                                                              │
+        └──────────────────────────────────────────────────────────────────────────────────────────────────┘"""
+
+        @test repr("text/plain", bidirectional) == show_box
+
+        zone = BoundaryZone(; boundary_face=([0.0, 0.0], [1 / sqrt(2), 1 / sqrt(2)]),
+                            particle_spacing=0.05, reference_density=0.0,
+                            reference_pressure=0.0,
+                            reference_velocity=[0.0, 0.0],
+                            face_normal=normalize([-1.0, 1.0]), density=1.0,
+                            open_boundary_layers=4)
+
+        show_compact = "BoundaryZone() with 80 particles"
+        @test repr(zone) == show_compact
+        show_box = """
+        ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
+        │ BoundaryZone                                                                                     │
+        │ ════════════                                                                                     │
+        │ boundary type: ………………………………………… bidirectional_flow                                               │
+        │ #particles: ………………………………………………… 80                                                               │
+        │ width: ……………………………………………………………… 0.2                                                              │
+        │ cross sectional area: ……………………… 1.0                                                              │
+        └──────────────────────────────────────────────────────────────────────────────────────────────────┘"""
+
+        @test repr("text/plain", zone) == show_box
     end
 
     @testset verbose=true "Illegal Inputs" begin

From d8cda93e7d425300e555074c107aa4f3214ccc57 Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Thu, 11 Dec 2025 09:39:52 +0100
Subject: [PATCH 21/62] fix doc tests

---
 src/schemes/boundary/open_boundary/boundary_zones.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl
index 55b86be737..6095ec71d1 100644
--- a/src/schemes/boundary/open_boundary/boundary_zones.jl
+++ b/src/schemes/boundary/open_boundary/boundary_zones.jl
@@ -145,6 +145,7 @@ bidirectional_flow = BoundaryZone(; boundary_face=face_vertices, face_normal,
 │ boundary type: ………………………………………… bidirectional_flow                                               │
 │ #particles: ………………………………………………… 234                                                              │
 │ width: ……………………………………………………………… 0.4                                                              │
+│ cross sectional area: ……………………… 1.0000000000000002                                               │
 └──────────────────────────────────────────────────────────────────────────────────────────────────┘
 ```
 

From be8a23ed944d8ed6f1960e01102074c67f7b34ba Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Thu, 11 Dec 2025 09:58:50 +0100
Subject: [PATCH 22/62] better docs

---
 src/schemes/boundary/open_boundary/boundary_zones.jl | 11 +++++++----
 src/schemes/boundary/open_boundary/system.jl         |  4 ++--
 2 files changed, 9 insertions(+), 6 deletions(-)

diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl
index 6095ec71d1..0a9577c16b 100644
--- a/src/schemes/boundary/open_boundary/boundary_zones.jl
+++ b/src/schemes/boundary/open_boundary/boundary_zones.jl
@@ -82,8 +82,10 @@ There are three ways to specify the actual shape of the boundary zone:
                        - For `EntropicallyDampedSPHSystem`: Use the initial pressure from the `InitialCondition`
                        - For `WeaklyCompressibleSPHSystem`: Use the background pressure from the equation of state
 - `sample_points=:default`: Either `:default` to automatically generate sample points on the boundary face (default),
-                            or a matrix of dimensions `(ndims, n_points)` containing sample points
+                            or a matrix of dimensions `(ndims, npoints)` containing sample points
                             on the boundary face used to compute the volumetric flow rate.
+                            Each sample point represents a discrete area of `particle_spacing^(ndims-1)`.
+                            Therefore, `sample_points` must form a regular grid.
                             Set to `nothing` to skip sampling.
 
 !!! note "Note"
@@ -332,8 +334,6 @@ function create_cache_boundary_zone(initial_condition, boundary_face, face_norma
                                 "`ndims(initial_condition)` rows"))
         end
 
-        # TODO: Check if regular grid?
-
         sample_points_ = sample_points
     end
 
@@ -347,6 +347,9 @@ function create_cache_boundary_zone(initial_condition, boundary_face, face_norma
         face_area = norm(v2 - v1)
     end
 
+    # We only check if the discretized area exceeds the boundary face area.
+    # For 3D boundary zones with complex or non-rectangular flow profiles
+    # (e.g., pipe flow), the cross-sectional area can legitimately be smaller then the boundary face area.
     if discrete_face_area > (face_area + eps(face_area))
         @warn "The sampled area of the boundary face " *
               "($(discrete_face_area)) is larger than the actual face area ($(face_area)). "
@@ -509,7 +512,7 @@ function update_boundary_zone_indices!(system, u, boundary_zones, semi)
         # - Floating-point rounding when a particle lies almost exactly on the `boundary_face`
         #   during transition, causing a reset just outside the zone
         #   (fixed in https://github.com/trixi-framework/TrixiParticles.jl/pull/997).
-        @assert system.boundary_zone_indices[particle] != 0 "No boundary zone found for active buffer particle"
+        @assert system.boundary_zone_indices[particle]!=0 "No boundary zone found for active buffer particle"
     end
 
     return system
diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl
index 648a97cb03..0e09c669c0 100644
--- a/src/schemes/boundary/open_boundary/system.jl
+++ b/src/schemes/boundary/open_boundary/system.jl
@@ -12,8 +12,8 @@ Open boundary system for in- and outflow particles.
 - `fluid_system`: The corresponding fluid system
 - `boundary_model`: Boundary model (see [Open Boundary Models](@ref open_boundary_models))
 - `calculate_flow_rate=false`: Set to `true` to calculate the volumetric flow rate through each boundary zone.
-                               Default is `false`.
-- `buffer_size`: Number of buffer particles.
+                               Default is `false`. To ensure accurate flow rate computation, the cross-sectional area
+                               of the `boundary_zone` must be properly sampled. See `sample_points` in [`BoundaryZone`](@ref).
 - `pressure_acceleration`: Pressure acceleration formulation for the system. Required only
                            when using [`BoundaryModelDynamicalPressureZhang`](@ref).
                            Defaults to the formulation from `fluid_system` if applicable; otherwise, `nothing`.

From 8aff117fcb73723e8b58ebae76f8dc2f3931bdd7 Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Fri, 12 Dec 2025 15:35:33 +0100
Subject: [PATCH 23/62] use stored type

---
 src/io/read_vtk.jl  | 10 ++++----
 test/io/read_vtk.jl | 56 +++++++++++++++++++++++++++++++++++++--------
 2 files changed, 53 insertions(+), 13 deletions(-)

diff --git a/src/io/read_vtk.jl b/src/io/read_vtk.jl
index 884588033b..d7a578c25e 100644
--- a/src/io/read_vtk.jl
+++ b/src/io/read_vtk.jl
@@ -33,18 +33,20 @@ ic = vtk2trixi(joinpath("out", "rectangular.vtu"))
 └──────────────────────────────────────────────────────────────────────────────────────────────────┘
 ```
 """
-function vtk2trixi(file; element_type=Float64, coordinates_eltype=Float64)
-    ELTYPE = element_type
-    cELTYPE = coordinates_eltype
+function vtk2trixi(file; element_type=:default, coordinates_eltype=Float64)
     vtk_file = ReadVTK.VTKFile(file)
 
     # Retrieve data fields (e.g., pressure, velocity, ...)
     point_data = ReadVTK.get_point_data(vtk_file)
     field_data = ReadVTK.get_field_data(vtk_file)
+    point_coords = ReadVTK.get_points(vtk_file)
+
+    cELTYPE = coordinates_eltype
+    ELTYPE = element_type === :default ? eltype(point_coords) : element_type
 
     # Retrieve fields
     ndims = first(ReadVTK.get_data(field_data["ndims"]))
-    coordinates = convert.(cELTYPE, ReadVTK.get_points(vtk_file)[1:ndims, :])
+    coordinates = convert.(cELTYPE, point_coords[1:ndims, :])
 
     fields = ["velocity", "density", "pressure", "mass", "particle_spacing"]
     results = Dict{String, Array{Float64}}()
diff --git a/test/io/read_vtk.jl b/test/io/read_vtk.jl
index d82e8c50c3..1331fdc9b9 100644
--- a/test/io/read_vtk.jl
+++ b/test/io/read_vtk.jl
@@ -7,15 +7,53 @@
                                        density=1000.0, pressure=900.0, mass=50.0)
 
         @testset verbose=true "`InitialCondition`" begin
-            trixi2vtk(expected_ic; filename="tmp_initial_condition",
-                      output_directory=tmp_dir)
-
-            test_ic = vtk2trixi(joinpath(tmp_dir, "tmp_initial_condition.vtu"))
-
-            @test isapprox(expected_ic.coordinates, test_ic.coordinates, rtol=1e-5)
-            @test isapprox(expected_ic.velocity, test_ic.velocity, rtol=1e-5)
-            @test isapprox(expected_ic.density, test_ic.density, rtol=1e-5)
-            @test isapprox(expected_ic.pressure, test_ic.pressure, rtol=1e-5)
+            @testset verbose=true "`Float64`" begin
+                trixi2vtk(expected_ic; filename="tmp_initial_condition_64",
+                          output_directory=tmp_dir)
+                file = joinpath(tmp_dir, "tmp_initial_condition_64.vtu")
+                test_ic = vtk2trixi(file)
+
+                @test isapprox(expected_ic.coordinates, test_ic.coordinates, rtol=1e-5)
+                @test isapprox(expected_ic.velocity, test_ic.velocity, rtol=1e-5)
+                @test isapprox(expected_ic.density, test_ic.density, rtol=1e-5)
+                @test isapprox(expected_ic.pressure, test_ic.pressure, rtol=1e-5)
+                @test eltype(test_ic) === Float64
+                @test eltype(test_ic.coordinates) === Float64
+            end
+
+            @testset verbose=true "`Float32`" begin
+                expected_ic_32 = InitialCondition(;
+                                                  coordinates=convert.(Float32,
+                                                                       coordinates),
+                                                  velocity=convert.(Float32, velocity),
+                                                  density=1000.0f0, pressure=900.0f0,
+                                                  mass=50.0f0, particle_spacing=0.1f0)
+                trixi2vtk(expected_ic_32; filename="tmp_initial_condition_32",
+                          output_directory=tmp_dir)
+                file = joinpath(tmp_dir, "tmp_initial_condition_32.vtu")
+                test_ic = vtk2trixi(file)
+
+                @test isapprox(expected_ic_32.coordinates, test_ic.coordinates, rtol=1e-5)
+                @test isapprox(expected_ic_32.velocity, test_ic.velocity, rtol=1e-5)
+                @test isapprox(expected_ic_32.density, test_ic.density, rtol=1e-5)
+                @test isapprox(expected_ic_32.pressure, test_ic.pressure, rtol=1e-5)
+                @test eltype(test_ic) === Float32
+                @test eltype(test_ic.coordinates) === Float64
+            end
+
+            @testset verbose=true "Custom Element Type" begin
+                trixi2vtk(expected_ic; filename="tmp_initial_condition_64",
+                          output_directory=tmp_dir)
+                file = joinpath(tmp_dir, "tmp_initial_condition_64.vtu")
+                test_ic = vtk2trixi(file, element_type=Float32, coordinates_eltype=Float32)
+
+                @test isapprox(expected_ic.coordinates, test_ic.coordinates, rtol=1e-5)
+                @test isapprox(expected_ic.velocity, test_ic.velocity, rtol=1e-5)
+                @test isapprox(expected_ic.density, test_ic.density, rtol=1e-5)
+                @test isapprox(expected_ic.pressure, test_ic.pressure, rtol=1e-5)
+                @test eltype(test_ic) === Float32
+                @test eltype(test_ic.coordinates) === Float32
+            end
         end
 
         @testset verbose=true "`AbstractFluidSystem`" begin

From 868b44922c8300cde607af1b8fb14361ec0d996b Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Fri, 12 Dec 2025 15:39:53 +0100
Subject: [PATCH 24/62] add docs

---
 src/io/read_vtk.jl | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/src/io/read_vtk.jl b/src/io/read_vtk.jl
index d7a578c25e..0638c5c21a 100644
--- a/src/io/read_vtk.jl
+++ b/src/io/read_vtk.jl
@@ -1,11 +1,16 @@
 """
-	vtk2trixi(file::String)
+	vtk2trixi(file::String; element_type=:default, coordinates_eltype=Float64)
 
 Load VTK file and convert data to an [`InitialCondition`](@ref).
 
 # Arguments
 - `file`: Name of the VTK file to be loaded.
 
+# Keywords
+- `element_type`: Element type for particle fields (`:default` keeps the type
+    stored in the VTK file, otherwise converted to the given type).
+- `coordinates_eltype`: Element type for particle coordinates (defaults to `Float64`).
+
 !!! warning "Experimental Implementation"
     This is an experimental feature and may change in any future releases.
 

From 6e390c2524dd8d358c8dc360e7b9b685505dcf1d Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Fri, 12 Dec 2025 17:15:04 +0100
Subject: [PATCH 25/62] add tests

---
 .../boundary/open_boundary/flow_rate.jl       | 64 +++++++++++++++++++
 .../boundary/open_boundary/open_boundary.jl   |  1 +
 2 files changed, 65 insertions(+)
 create mode 100644 test/schemes/boundary/open_boundary/flow_rate.jl

diff --git a/test/schemes/boundary/open_boundary/flow_rate.jl b/test/schemes/boundary/open_boundary/flow_rate.jl
new file mode 100644
index 0000000000..9bbaf0e277
--- /dev/null
+++ b/test/schemes/boundary/open_boundary/flow_rate.jl
@@ -0,0 +1,64 @@
+@testset verbose=true "Calculate Flow Rate" begin
+    particle_spacing = 0.01
+
+    # Define a parabolic velocity profile
+    velocity_function(pos) = [4 * (pos[2] - 0.5) * (1.5 - pos[2]), 0.0]
+
+    # Create fluid domain with the specified velocity profile
+    n_particles_y = round(Int, 2 / particle_spacing)
+    fluid = RectangularShape(particle_spacing, (4, n_particles_y), (0.0, 0.0),
+                             density=1000.0, velocity=velocity_function)
+
+    smoothing_length = 1.3 * particle_spacing
+    smoothing_kernel = WendlandC2Kernel{ndims(fluid)}()
+    fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length,
+                                               1.0)
+    fluid_system.cache.density .= fluid.density
+
+    # Use a smaller cross-sectional area to test user-defined area functionality
+    # and to perform interpolation within an embedded domain.
+    # (In a simulation with solid boundaries, wall velocities would also be included.)
+    sample_points = RectangularShape(particle_spacing, (1, round(Int, n_particles_y / 2)),
+                                     (0.0, 0.5), density=1.0).coordinates
+
+    zone = BoundaryZone(; boundary_face=([0.0, 0.0], [0.0, 2.0]),
+                        face_normal=(-1.0, 0.0), open_boundary_layers=10, density=1000.0,
+                        particle_spacing, sample_points=sample_points,
+                        reference_velocity=(pos, t) -> velocity_function(pos))
+
+    open_boundary = OpenBoundarySystem(zone; fluid_system,
+                                       boundary_model=BoundaryModelMirroringTafuni(),
+                                       calculate_flow_rate=true, buffer_size=0)
+
+    semi = Semidiscretization(fluid_system, open_boundary)
+    TrixiParticles.initialize_neighborhood_searches!(semi)
+    TrixiParticles.initialize!(open_boundary, semi)
+
+    # Set up ODE for initial conditions
+    ode = semidiscretize(semi, (0, 1))
+    v_ode, u_ode = ode.u0.x
+
+    boundary_zone = first(open_boundary.boundary_zones)
+
+    @testset verbose=true "Velocity Interpolation" begin
+        TrixiParticles.interpolate_velocity!(open_boundary, boundary_zone,
+                                             v_ode, u_ode, semi)
+
+        coords = reinterpret(reshape, SVector{2, Float64},
+                             boundary_zone.cache.sample_points)
+        vels_analytic = first.(velocity_function.(coords))
+        vels_interpolated = first.(reinterpret(reshape, SVector{2, Float64},
+                                               boundary_zone.cache.sample_velocity))
+
+        @test isapprox(vels_interpolated, vels_analytic, rtol=1e-3)
+    end
+
+    @testset verbose=true "Flow Rate Calculation" begin
+        TrixiParticles.calculate_flow_rate!(open_boundary, v_ode, u_ode, semi)
+
+        Q_analytic = zone.cache.cross_sectional_area * (2 / 3)
+        Q_calculated = first(open_boundary.cache.boundary_zones_flow_rate)[]
+
+        @test isapprox(Q_analytic, Q_calculated; rtol=1e-3)
+    end
+end
diff --git a/test/schemes/boundary/open_boundary/open_boundary.jl b/test/schemes/boundary/open_boundary/open_boundary.jl
index 78c73345e4..74fd8098df 100644
--- a/test/schemes/boundary/open_boundary/open_boundary.jl
+++ b/test/schemes/boundary/open_boundary/open_boundary.jl
@@ -2,3 +2,4 @@ include("characteristic_variables.jl")
 include("mirroring.jl")
 include("boundary_zone.jl")
 include("pressure_model.jl")
+include("flow_rate.jl")

From 6121e7d8b3c071027f46949a578b340e40c895bc Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Fri, 12 Dec 2025 17:43:17 +0100
Subject: [PATCH 26/62] fix formatting

---
 src/schemes/boundary/open_boundary/boundary_zones.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl
index 0a9577c16b..7f15f92def 100644
--- a/src/schemes/boundary/open_boundary/boundary_zones.jl
+++ b/src/schemes/boundary/open_boundary/boundary_zones.jl
@@ -349,7 +349,7 @@ function create_cache_boundary_zone(initial_condition, boundary_face, face_norma
 
     # We only check if the discretized area exceeds the boundary face area.
     # For 3D boundary zones with complex or non-rectangular flow profiles
-    # (e.g., pipe flow), the cross-sectional area can legitimately be smaller then the boundary face area.
+    # (e.g., pipe flow), the cross-sectional area can legitimately be smaller than the boundary face area.
     if discrete_face_area > (face_area + eps(face_area))
         @warn "The sampled area of the boundary face " *
               "($(discrete_face_area)) is larger than the actual face area ($(face_area)). "
@@ -512,7 +512,7 @@ function update_boundary_zone_indices!(system, u, boundary_zones, semi)
         # - Floating-point rounding when a particle lies almost exactly on the `boundary_face`
         #   during transition, causing a reset just outside the zone
         #   (fixed in https://github.com/trixi-framework/TrixiParticles.jl/pull/997).
-        @assert system.boundary_zone_indices[particle]!=0 "No boundary zone found for active buffer particle"
+        @assert system.boundary_zone_indices[particle] != 0 "No boundary zone found for active buffer particle"
     end
 
     return system

From 7369f71c127077a8d059ca29fb3a4cd12eccd0e5 Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Sat, 13 Dec 2025 09:54:42 +0100
Subject: [PATCH 27/62] fix eltype

---
 src/schemes/boundary/open_boundary/boundary_zones.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl
index 7f15f92def..8e42b7088a 100644
--- a/src/schemes/boundary/open_boundary/boundary_zones.jl
+++ b/src/schemes/boundary/open_boundary/boundary_zones.jl
@@ -326,15 +326,15 @@ function create_cache_boundary_zone(initial_condition, boundary_face, face_norma
     (; particle_spacing) = initial_condition
     area_increment = particle_spacing^(ndims(initial_condition) - 1)
     if sample_points === :default
-        sample_points_ = extrude_geometry(boundary_face; particle_spacing, density=Inf,
-                                          direction=(-face_normal), n_extrude=1).coordinates
+        points = extrude_geometry(boundary_face; particle_spacing, density=Inf,
+                                  direction=(-face_normal), n_extrude=1).coordinates
+        sample_points_ = convert.(eltype(initial_condition), points)
     else
         if !(sample_points isa Matrix && size(sample_points, 1) == ndims(initial_condition))
             throw(ArgumentError("`sample_points` must be a matrix with " *
                                 "`ndims(initial_condition)` rows"))
         end
-
-        sample_points_ = sample_points
+        sample_points_ = convert.(eltype(initial_condition), sample_points)
     end
 
     discrete_face_area = area_increment * size(sample_points_, 2)

From a9123267463a4a0a38f2355c6e226b34aa740e3c Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Tue, 16 Dec 2025 09:09:06 +0100
Subject: [PATCH 28/62] revise #959

---
 src/io/read_vtk.jl  | 13 ++++++----
 test/io/read_vtk.jl | 61 ++++++++++++++++++++-------------------------
 2 files changed, 35 insertions(+), 39 deletions(-)

diff --git a/src/io/read_vtk.jl b/src/io/read_vtk.jl
index acce879a04..e30a28d76c 100644
--- a/src/io/read_vtk.jl
+++ b/src/io/read_vtk.jl
@@ -2,7 +2,7 @@
 	vtk2trixi(file::String; element_type=:default, coordinates_eltype=Float64,
               custom_quantities...)
 
-Load VTK file and convert data to a NamedTuple.
+Load VTK file and convert data to a `NamedTuple`.
 
 # Arguments
 - `file`:                 Name of the VTK file to be loaded.
@@ -76,14 +76,17 @@ function vtk2trixi(file; element_type=:default, coordinates_eltype=Float64,
                                  results[:particle_spacing]
     results[:coordinates] = coordinates
     results[:time] = "time" in keys(field_data) ?
-                     first(ReadVTK.get_data(field_data["time"])) : 0.0
+                     first(ReadVTK.get_data(field_data["time"])) : zero(ELTYPE)
 
-    for (key, quantity) in custom_quantities
+    for (key, quantity_) in custom_quantities
+        quantity = string(quantity_)
         if quantity in keys(point_data)
             results[key] = ReadVTK.get_data(point_data[quantity])
-        end
-        if quantity in keys(field_data)
+        elseif quantity in keys(field_data)
             results[key] = first(ReadVTK.get_data(field_data[quantity]))
+        else
+            throw(ArgumentError("Custom quantity '$quantity' not found in VTK file. " *
+                                "Make sure it was included during the simulation."))
         end
     end
 
diff --git a/test/io/read_vtk.jl b/test/io/read_vtk.jl
index a4ac246dde..43bac75794 100644
--- a/test/io/read_vtk.jl
+++ b/test/io/read_vtk.jl
@@ -7,11 +7,9 @@
                                          density=1000.0, pressure=900.0, mass=50.0)
 
         expected_cq_scalar = 3.0
-        expected_cq_vector = fill(expected_cq_scalar,
-                                  size(expected_data.coordinates, 2))
+        expected_cq_vector = fill(expected_cq_scalar, nparticles(expected_data))
         scalar_cq_function(system, data, t) = expected_cq_scalar
-        vector_cq_function(system, data,
-        t) = fill(expected_cq_scalar, nparticles(system))
+        vector_cq_function(system, data, t) = fill(expected_cq_scalar, nparticles(system))
 
         @testset verbose=true "`InitialCondition`" begin
             @testset verbose=true "`Float64`" begin
@@ -24,6 +22,7 @@
                 @test isapprox(expected_data.velocity, data.velocity, rtol=1e-5)
                 @test isapprox(expected_data.density, data.density, rtol=1e-5)
                 @test isapprox(expected_data.pressure, data.pressure, rtol=1e-5)
+                @test all(key -> eltype(data[key]) === Float64, keys(data))
                 @test eltype(data.coordinates) === Float64
             end
 
@@ -43,6 +42,8 @@
                 @test isapprox(expected_ic_32.velocity, data.velocity, rtol=1e-5)
                 @test isapprox(expected_ic_32.density, data.density, rtol=1e-5)
                 @test isapprox(expected_ic_32.pressure, data.pressure, rtol=1e-5)
+                @test all(key -> eltype(data[key]) === Float32,
+                          setdiff(keys(data), (:coordinates,)))
                 @test eltype(data.coordinates) === Float64
             end
 
@@ -56,41 +57,39 @@
                 @test isapprox(expected_data.velocity, data.velocity, rtol=1e-5)
                 @test isapprox(expected_data.density, data.density, rtol=1e-5)
                 @test isapprox(expected_data.pressure, data.pressure, rtol=1e-5)
+                @test all(key -> eltype(data[key]) === Float32, keys(data))
                 @test eltype(data.coordinates) === Float32
             end
 
-            @testset verbose=true "Scalar custom quantity" begin
+            @testset verbose=true "Scalar Custom Quantity" begin
                 trixi2vtk(expected_data; filename="tmp_initial_condition_scalar",
-                          output_directory=tmp_dir,
-                          cq_scalar=expected_cq_scalar)
+                          output_directory=tmp_dir, cq_scalar=expected_cq_scalar)
 
                 # Load file containing test data
                 test_data = vtk2trixi(joinpath(tmp_dir, "tmp_initial_condition_scalar.vtu");
-                                      cq_scalar="cq_scalar")
+                                      cq_scalar=:cq_scalar)
 
                 @test isapprox(expected_data.coordinates, test_data.coordinates, rtol=1e-5)
                 @test isapprox(expected_data.velocity, test_data.velocity, rtol=1e-5)
                 @test isapprox(expected_data.density, test_data.density, rtol=1e-5)
                 @test isapprox(expected_data.pressure, test_data.pressure, rtol=1e-5)
-                @test isapprox(expected_cq_scalar, test_data.cq_scalar,
-                               rtol=1e-5)
+                @test isapprox(expected_cq_scalar, test_data.cq_scalar, rtol=1e-5)
             end
 
-            @testset verbose=true "Vector custom quantity" begin
+            @testset verbose=true "Vector Custom Quantity" begin
                 trixi2vtk(expected_data; filename="tmp_initial_condition_vector",
                           output_directory=tmp_dir,
                           cq_vector=expected_cq_vector)
 
                 # Load file containing test data
                 test_data = vtk2trixi(joinpath(tmp_dir, "tmp_initial_condition_vector.vtu");
-                                      cq_vector="cq_vector")
+                                      cq_vector=:cq_vector)
 
                 @test isapprox(expected_data.coordinates, test_data.coordinates, rtol=1e-5)
                 @test isapprox(expected_data.velocity, test_data.velocity, rtol=1e-5)
                 @test isapprox(expected_data.density, test_data.density, rtol=1e-5)
                 @test isapprox(expected_data.pressure, test_data.pressure, rtol=1e-5)
-                @test isapprox(expected_cq_vector, test_data.cq_vector,
-                               rtol=1e-5)
+                @test isapprox(expected_cq_vector, test_data.cq_vector, rtol=1e-5)
             end
         end
 
@@ -116,33 +115,30 @@
             vu_ode = (; x)
             dvdu_ode = (; x=(; v_ode=dv_ode, u_ode=du_ode))
 
-            @testset verbose=true "Scalar custom quantity" begin
+            @testset verbose=true "Scalar Custom Quantity" begin
                 trixi2vtk(fluid_system, dvdu_ode, vu_ode, semi, 0.0,
                           nothing; system_name="tmp_file_fluid_scalar",
-                          output_directory=tmp_dir,
-                          iter=1, cq_scalar=scalar_cq_function)
+                          output_directory=tmp_dir, iter=1, cq_scalar=scalar_cq_function)
 
                 # Load file containing test data
                 test_data = vtk2trixi(joinpath(tmp_dir, "tmp_file_fluid_scalar_1.vtu");
-                                      cq_scalar="cq_scalar")
+                                      cq_scalar=:cq_scalar)
 
                 @test isapprox(u, test_data.coordinates, rtol=1e-5)
                 @test isapprox(v, test_data.velocity, rtol=1e-5)
                 @test isapprox(pressure, test_data.pressure, rtol=1e-5)
                 @test isapprox(fluid_system.cache.density, test_data.density, rtol=1e-5)
-                @test isapprox(expected_cq_scalar, test_data.cq_scalar,
-                               rtol=1e-5)
+                @test isapprox(expected_cq_scalar, test_data.cq_scalar, rtol=1e-5)
             end
 
-            @testset verbose=true "Vector custom quantity" begin
+            @testset verbose=true "Vector Custom Quantity" begin
                 trixi2vtk(fluid_system, dvdu_ode, vu_ode, semi, 0.0,
                           nothing; system_name="tmp_file_fluid_vector",
-                          output_directory=tmp_dir,
-                          iter=1, cq_vector=vector_cq_function)
+                          output_directory=tmp_dir, iter=1, cq_vector=vector_cq_function)
 
                 # Load file containing test data
                 test_data = vtk2trixi(joinpath(tmp_dir, "tmp_file_fluid_vector_1.vtu");
-                                      cq_vector="cq_vector")
+                                      cq_vector=:cq_vector)
 
                 @test isapprox(u, test_data.coordinates, rtol=1e-5)
                 @test isapprox(v, test_data.velocity, rtol=1e-5)
@@ -176,15 +172,14 @@
             vu_ode = (; x)
             dvdu_ode = (; x=(; v_ode=dv_ode, u_ode=du_ode))
 
-            @testset verbose=true "Scalar custom quantity" begin
+            @testset verbose=true "Scalar Custom Quantity" begin
                 trixi2vtk(boundary_system, dvdu_ode, vu_ode, semi, 0.0,
                           nothing; system_name="tmp_file_boundary_scalar",
-                          output_directory=tmp_dir,
-                          iter=1, cq_scalar=scalar_cq_function)
+                          output_directory=tmp_dir, iter=1, cq_scalar=scalar_cq_function)
 
                 # Load file containing test data
                 test_data = vtk2trixi(joinpath(tmp_dir, "tmp_file_boundary_scalar_1.vtu");
-                                      cq_scalar="cq_scalar")
+                                      cq_scalar=:cq_scalar)
 
                 @test isapprox(boundary_system.coordinates, test_data.coordinates,
                                rtol=1e-5)
@@ -193,19 +188,17 @@
                                rtol=1e-5)
                 @test isapprox(boundary_model.pressure, test_data.pressure, rtol=1e-5)
                 @test isapprox(boundary_model.cache.density, test_data.density, rtol=1e-5)
-                @test isapprox(expected_cq_scalar, test_data.cq_scalar,
-                               rtol=1e-5)
+                @test isapprox(expected_cq_scalar, test_data.cq_scalar, rtol=1e-5)
             end
 
-            @testset verbose=true "Vector custom quantity" begin
+            @testset verbose=true "Vector Custom Quantity" begin
                 trixi2vtk(boundary_system, dvdu_ode, vu_ode, semi, 0.0,
                           nothing; system_name="tmp_file_boundary_vector",
-                          output_directory=tmp_dir,
-                          iter=1, cq_vector=vector_cq_function)
+                          output_directory=tmp_dir, iter=1, cq_vector=vector_cq_function)
 
                 # Load file containing test data
                 test_data = vtk2trixi(joinpath(tmp_dir, "tmp_file_boundary_vector_1.vtu");
-                                      cq_vector="cq_vector")
+                                      cq_vector=:cq_vector)
 
                 @test isapprox(boundary_system.coordinates, test_data.coordinates,
                                rtol=1e-5)

From 29b5b1a283d2d7f7582a8a92dca93f44e9e9ff73 Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Tue, 6 Jan 2026 11:51:14 +0100
Subject: [PATCH 29/62] fix buffer

---
 src/general/buffer.jl                        | 4 ++--
 src/schemes/boundary/open_boundary/system.jl | 8 ++++----
 test/general/buffer.jl                       | 9 +++------
 3 files changed, 9 insertions(+), 12 deletions(-)

diff --git a/src/general/buffer.jl b/src/general/buffer.jl
index 1e2c196799..48ab0d92be 100644
--- a/src/general/buffer.jl
+++ b/src/general/buffer.jl
@@ -40,10 +40,10 @@ end
 # Dispatch by system type to handle systems that provide a buffer.
 @inline buffer(system) = nothing
 
-@inline update_system_buffer!(buffer::Nothing, semi) = buffer
+@inline update_system_buffer!(buffer::Nothing) = buffer
 
 # TODO `resize` allocates. Find a non-allocating version
-@inline function update_system_buffer!(buffer::SystemBuffer, semi)
+@inline function update_system_buffer!(buffer::SystemBuffer)
     (; active_particle) = buffer
 
     # TODO: Parallelize (see https://github.com/trixi-framework/TrixiParticles.jl/issues/810)
diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl
index a1e7140c52..b74da7153d 100644
--- a/src/schemes/boundary/open_boundary/system.jl
+++ b/src/schemes/boundary/open_boundary/system.jl
@@ -405,8 +405,8 @@ function check_domain!(system, v, u, v_ode, u_ode, semi)
                           v, u, v_fluid, u_fluid)
     end
 
-    update_system_buffer!(system.buffer, semi)
-    update_system_buffer!(fluid_system.buffer, semi)
+    update_system_buffer!(system.buffer)
+    update_system_buffer!(fluid_system.buffer)
 
     fluid_candidates .= false
 
@@ -436,8 +436,8 @@ function check_domain!(system, v, u, v_ode, u_ode, semi)
                           v, u, v_fluid, u_fluid)
     end
 
-    update_system_buffer!(system.buffer, semi)
-    update_system_buffer!(fluid_system.buffer, semi)
+    update_system_buffer!(system.buffer)
+    update_system_buffer!(fluid_system.buffer)
 
     # Since particles have been transferred, the neighborhood searches must be updated
     update_nhs!(semi, u_ode)
diff --git a/test/general/buffer.jl b/test/general/buffer.jl
index f9cdf25cb2..9eb6d87fac 100644
--- a/test/general/buffer.jl
+++ b/test/general/buffer.jl
@@ -28,16 +28,14 @@
         particle_id = findfirst(==(false), system_buffer.buffer.active_particle)
         system_buffer.buffer.active_particle[particle_id] = true
 
-        TrixiParticles.update_system_buffer!(system_buffer.buffer,
-                                             DummySemidiscretization())
+        TrixiParticles.update_system_buffer!(system_buffer.buffer)
 
         @test TrixiParticles.each_integrated_particle(system_buffer) == 1:(n_particles + 1)
 
         TrixiParticles.deactivate_particle!(system_buffer, particle_id,
                                             ones(2, particle_id))
 
-        TrixiParticles.update_system_buffer!(system_buffer.buffer,
-                                             DummySemidiscretization())
+        TrixiParticles.update_system_buffer!(system_buffer.buffer)
 
         @test TrixiParticles.each_integrated_particle(system_buffer) == 1:n_particles
 
@@ -45,8 +43,7 @@
         TrixiParticles.deactivate_particle!(system_buffer, particle_id,
                                             ones(2, particle_id))
 
-        TrixiParticles.update_system_buffer!(system_buffer.buffer,
-                                             DummySemidiscretization())
+        TrixiParticles.update_system_buffer!(system_buffer.buffer)
 
         @test TrixiParticles.each_integrated_particle(system_buffer) ==
               setdiff(1:n_particles, particle_id)

From daceea26d9c78bad972761a46c5f9c5a6eb9c18f Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Tue, 6 Jan 2026 12:15:49 +0100
Subject: [PATCH 30/62] first poc

---
 src/general/restart.jl                        | 94 +++++++++++++++++++
 src/general/semidiscretization.jl             | 25 +++--
 src/schemes/boundary/open_boundary/system.jl  | 45 +++++++++
 .../boundary/wall_boundary/dummy_particles.jl |  4 +
 .../boundary/wall_boundary/monaghan_kajtar.jl |  4 +
 src/schemes/boundary/wall_boundary/system.jl  | 17 ++++
 src/schemes/fluid/fluid.jl                    | 45 +++++++++
 .../structure/total_lagrangian_sph/system.jl  | 11 +++
 8 files changed, 236 insertions(+), 9 deletions(-)
 create mode 100644 src/general/restart.jl

diff --git a/src/general/restart.jl b/src/general/restart.jl
new file mode 100644
index 0000000000..df3eaef261
--- /dev/null
+++ b/src/general/restart.jl
@@ -0,0 +1,94 @@
+struct RestartCondition{V, U}
+    v_restart::V
+    u_restart::U
+    t_restart::Real
+end
+
+function RestartCondition(system::AbstractSystem, restart_file; precondition_values=nothing)
+    restart_data = vtk2trixi(restart_file)
+    v_restart = restart_v(system, restart_data)
+    u_restart = restart_u(system, restart_data)
+
+    # TODO
+    t_restart = 0.0
+
+    precondition_system!(system, precondition_values)
+
+    return RestartCondition(v_restart, u_restart, t_restart)
+end
+
+function set_intial_conditions!(v0_ode, u0_ode, semi, restart_conditions)
+    if length(semi.systems) != length(restart_conditions)
+        throw(ArgumentError("Number of systems in `semi` does not match number of `restart_conditions`"))
+    end
+
+    foreach_noalloc(semi.systems, restart_conditions) do (system, restart_condition)
+        v0_system = wrap_v(v0_ode, system, semi)
+        u0_system = wrap_u(u0_ode, system, semi)
+
+        u0_system .= restart_condition.u_restart
+        v0_system .= restart_condition.v_restart
+    end
+end
+
+function time_span(tspan, restart_conditions::RestartCondition)
+    return (first(restart_conditions).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
+
+precondition_system!(system, values::Nothing) = system
+
+function precondition_system!(system::OpenBoundarySystem, values)
+    (; pressure_reference_values, boundary_zones_flow_rate) = system.cache
+
+    if !haskey(values, :pressure_reference_values, :boundary_zones_flow_rate)
+        throw(ArgumentError("Missing required fields in `values` for `OpenBoundarySystem`"))
+    end
+
+    foreach_noalloc(pressure_reference_values,
+                    values.pressure_reference_values) do (pressure_model,
+                                                          previous_pressure)
+        if isa(pressure_model, AbstractPressureModel)
+            pressure_model.pressure[] = previous_pressure
+        else
+            if previous_pressure != Inf
+                throw(ArgumentError("Expected `Inf` for non-pressure-model boundary zones"))
+            end
+        end
+    end
+
+    foreach_noalloc(pressure_reference_values,
+                    values.boundary_zones_flow_rate) do (pressure_model,
+                                                         previous_flow_rate)
+        if isa(pressure_model, AbstractPressureModel)
+            pressure_model.flow_rate[] = previous_flow_rate
+        else
+            if previous_flow_rate != Inf
+                throw(ArgumentError("Expected `Inf` for non-pressure-model boundary zones"))
+            end
+        end
+    end
+
+    return system
+end
diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl
index c350511b63..dc991aae95 100644
--- a/src/general/semidiscretization.jl
+++ b/src/general/semidiscretization.jl
@@ -288,7 +288,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_conditions=nothing)
     (; systems) = semi
 
     # Check that all systems have the same eltype
@@ -329,13 +329,7 @@ function semidiscretize(semi, tspan; reset_threads=true)
     end
 
     # Set initial condition
-    foreach_system(semi) do system
-        u0_system = wrap_u(u0_ode, system, semi)
-        v0_system = wrap_v(v0_ode, system, semi)
-
-        write_u0!(u0_system, system)
-        write_v0!(v0_system, system)
-    end
+    set_intial_conditions!(v0_ode, u0_ode, semi, restart_conditions)
 
     # TODO initialize after adapting to the GPU.
     # Requires https://github.com/trixi-framework/PointNeighbors.jl/pull/86.
@@ -371,9 +365,22 @@ function semidiscretize(semi, tspan; reset_threads=true)
     # 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_conditions), semi_new)
 end
 
+function set_intial_conditions!(v0_ode, u0_ode, semi, restart_conditions::Nothing)
+    foreach_system(semi) do system
+        v0_system = wrap_v(v0_ode, system, semi)
+        u0_system = wrap_u(u0_ode, system, semi)
+
+        write_v0!(v0_system, system)
+        write_u0!(u0_system, system)
+    end
+end
+
+time_span(tspan, restart_conditions::Nothing) = tspan
+
 """
     restart_with!(semi, sol)
 
diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl
index b74da7153d..c18416989c 100644
--- a/src/schemes/boundary/open_boundary/system.jl
+++ b/src/schemes/boundary/open_boundary/system.jl
@@ -659,3 +659,48 @@ function interpolate_velocity!(system::OpenBoundarySystem, boundary_zone,
 
     return system
 end
+
+function restart_u(system::OpenBoundarySystem, data)
+    coords_total = zeros(eltype(system), u_nvariables(system),
+                         n_integrated_particles(system))
+    coords_total .= eltype(system)(1e16)
+    system.buffer.active_particle .= false
+
+    coords_active = data.coordinates
+
+    for particle in axes(coords_active, 2)
+        system.buffer.active_particle[particle] = true
+        for dim in 1:ndims(system)
+            coords_total[dim, particle] = coords_active[dim, particle]
+        end
+    end
+
+    update_system_buffer!(system.buffer)
+
+    return coords_total
+end
+
+function restart_v(system::OpenBoundarySystem, data)
+    velocity_total = zeros(eltype(system), v_nvariables(system),
+                           n_integrated_particles(system))
+    velocity_total .= eltype(system)(1e16)
+
+    system.buffer.active_particle .= false
+
+    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.fluid_system,
+                                density_calculator(system), data.pressure, data.density)
+
+    for particle in axes(velocity_active, 2)
+        system.buffer.active_particle[particle] = true
+        for i in 1:axes(velocity_active, 1)
+            velocity_total[i, particle] = velocity_active[i, particle]
+        end
+    end
+
+    update_system_buffer!(system.buffer)
+
+    return coords_total
+end
diff --git a/src/schemes/boundary/wall_boundary/dummy_particles.jl b/src/schemes/boundary/wall_boundary/dummy_particles.jl
index 192eb5fee1..bc26b147fa 100644
--- a/src/schemes/boundary/wall_boundary/dummy_particles.jl
+++ b/src/schemes/boundary/wall_boundary/dummy_particles.jl
@@ -697,3 +697,7 @@ end
 @inline function correction_matrix(system::AbstractBoundarySystem, particle)
     extract_smatrix(system.boundary_model.cache.correction_matrix, system, particle)
 end
+
+@inline function density_calculator(system::WallBoundarySystem{<:BoundaryModelDummyParticles})
+    return system.boundary_model.density_calculator
+end
diff --git a/src/schemes/boundary/wall_boundary/monaghan_kajtar.jl b/src/schemes/boundary/wall_boundary/monaghan_kajtar.jl
index 1df9a38cf9..f5f3b9005a 100644
--- a/src/schemes/boundary/wall_boundary/monaghan_kajtar.jl
+++ b/src/schemes/boundary/wall_boundary/monaghan_kajtar.jl
@@ -119,3 +119,7 @@ end
     # Nothing to do in the update step
     return boundary_model
 end
+
+@inline function density_calculator(system::WallBoundarySystem{<:BoundaryModelMonaghanKajtar})
+    return nothing
+end
diff --git a/src/schemes/boundary/wall_boundary/system.jl b/src/schemes/boundary/wall_boundary/system.jl
index 151787e7aa..ce18344eea 100644
--- a/src/schemes/boundary/wall_boundary/system.jl
+++ b/src/schemes/boundary/wall_boundary/system.jl
@@ -246,6 +246,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 coords_total
+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,
diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl
index 35d14b5e04..0b9387baaf 100644
--- a/src/schemes/fluid/fluid.jl
+++ b/src/schemes/fluid/fluid.jl
@@ -220,6 +220,51 @@ end
     return nothing
 end
 
+function restart_u(system::AbstractFluidSystem, data)
+    coords_total = zeros(eltype(system), u_nvariables(system),
+                         n_integrated_particles(system))
+    coords_total .= eltype(system)(1e16)
+    isnothing(buffer(system)) || system.buffer.active_particle .= false
+
+    coords_active = data.coordinates
+
+    for particle in axes(coords_active, 2)
+        isnothing(buffer(system)) || system.buffer.active_particle[particle] = true
+        for dim in 1:ndims(system)
+            coords_total[dim, particle] = coords_active[dim, particle]
+        end
+    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_total .= eltype(system)(1e16)
+
+    isnothing(buffer(system)) || system.buffer.active_particle .= false
+
+    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)
+        isnothing(buffer(system)) || system.buffer.active_particle[particle] = true
+        for i in 1:axes(velocity_active, 1)
+            velocity_total[i, particle] = velocity_active[i, particle]
+        end
+    end
+
+    update_system_buffer!(system.buffer)
+
+    return coords_total
+end
+
 function system_data(system::AbstractFluidSystem, dv_ode, du_ode, v_ode, u_ode, semi)
     (; mass) = system
 
diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl
index ae1b7735db..18f4886a20 100644
--- a/src/schemes/structure/total_lagrangian_sph/system.jl
+++ b/src/schemes/structure/total_lagrangian_sph/system.jl
@@ -485,6 +485,17 @@ function restart_with!(system::TotalLagrangianSPHSystem, v, u)
     restart_with!(system, system.boundary_model, v, u)
 end
 
+function restart_u(system::AbstractStructureSystem, data)
+    # TODO: Is this correct?
+    # data.coordinates = coords_integrated
+    return data.coordinates[:, each_integrated_particle(system)]
+end
+
+function restart_v(system::AbstractStructureSystem, data)
+    # TODO: Is this correct?
+    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.

From e9750b0b7e2f7235c982de81d642b041e919fc18 Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Tue, 6 Jan 2026 15:19:23 +0100
Subject: [PATCH 31/62] add show

---
 src/general/restart.jl | 30 ++++++++++++++++++++++++++++--
 1 file changed, 28 insertions(+), 2 deletions(-)

diff --git a/src/general/restart.jl b/src/general/restart.jl
index c2a5526296..efbef368c2 100644
--- a/src/general/restart.jl
+++ b/src/general/restart.jl
@@ -4,8 +4,14 @@ struct RestartCondition{V, U}
     t_restart::Real
 end
 
-function RestartCondition(system::AbstractSystem, restart_file; precondition_values=nothing)
-    restart_data = vtk2trixi(restart_file)
+function RestartCondition(system::AbstractSystem, filename; output_directory="out",
+                          precondition_values=nothing)
+                          @autoinfiltrate
+    if !occursin(vtkname(system), basename(splitext(filename)[1]))
+        throw(ArgumentError("Filename '$filename' does not seem to correspond to system of type $(nameof(typeof(system)))."))
+    end
+
+    restart_data = vtk2trixi(joinpath(output_directory, filename))
     v_restart = restart_v(system, restart_data)
     u_restart = restart_u(system, restart_data)
 
@@ -16,6 +22,26 @@ function RestartCondition(system::AbstractSystem, restart_file; precondition_val
     return RestartCondition(v_restart, u_restart, restart_data.time)
 end
 
+function Base.show(io::IO, rc::RestartCondition)
+    @nospecialize rc # reduce precompilation time
+
+    print(io, "RestartCondition{}()")
+end
+
+function Base.show(io::IO, ::MIME"text/plain", rc::RestartCondition)
+    @nospecialize rc # reduce precompilation time
+
+    if get(io, :compact, false)
+        show(io, rc)
+    else
+        summary_header(io, "RestartCondition")
+        summary_line(io, "#particles u", "$(size(rc.u_restart, 2))")
+        summary_line(io, "#particles v", "$(size(rc.v_restart, 2))")
+        summary_line(io, "eltype", "$(eltype(rc.v_restart))")
+        summary_footer(io)
+    end
+end
+
 function set_intial_conditions!(v0_ode, u0_ode, semi, restart_conditions)
     if length(semi.systems) != length(restart_conditions)
         throw(ArgumentError("Number of systems in `semi` does not match number of `restart_conditions`"))

From 695fe6cf07e05a027c1395f779c7dffe4ab028de Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Tue, 6 Jan 2026 15:20:36 +0100
Subject: [PATCH 32/62] remove mass

---
 src/io/read_vtk.jl | 7 ++-----
 1 file changed, 2 insertions(+), 5 deletions(-)

diff --git a/src/io/read_vtk.jl b/src/io/read_vtk.jl
index e30a28d76c..6f0701f560 100644
--- a/src/io/read_vtk.jl
+++ b/src/io/read_vtk.jl
@@ -56,7 +56,7 @@ function vtk2trixi(file; element_type=:default, coordinates_eltype=Float64,
     ndims = first(ReadVTK.get_data(field_data["ndims"]))
     coordinates = convert.(cELTYPE, point_coords[1:ndims, :])
 
-    fields = [:velocity, :density, :pressure, :mass, :particle_spacing]
+    fields = [:velocity, :density, :pressure, :particle_spacing]
     for field in fields
         # Look for any key that contains the field name
         all_keys = keys(point_data)
@@ -64,10 +64,7 @@ function vtk2trixi(file; element_type=:default, coordinates_eltype=Float64,
         if idx !== nothing
             results[field] = convert.(ELTYPE, ReadVTK.get_data(point_data[all_keys[idx]]))
         else
-            # Use zeros as default values when a field is missing
-            results[field] = string(field) in ["mass"] ?
-                             zeros(ELTYPE, size(coordinates, 2)) : zero(coordinates)
-            @info "No '$field' field found in VTK file. Will be set to zero."
+            @info "No '$field' field found in VTK file"
         end
     end
 

From 2d8bce8925f493361585700b8513cffa38ed7fdd Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Tue, 6 Jan 2026 16:00:00 +0100
Subject: [PATCH 33/62] add preconditioning

---
 src/general/restart.jl                       | 41 +++-----------------
 src/io/write_vtk.jl                          |  8 ++++
 src/schemes/boundary/open_boundary/system.jl | 18 +++++++++
 3 files changed, 31 insertions(+), 36 deletions(-)

diff --git a/src/general/restart.jl b/src/general/restart.jl
index 87c8b5b905..e4c85ec0c6 100644
--- a/src/general/restart.jl
+++ b/src/general/restart.jl
@@ -4,19 +4,17 @@ struct RestartCondition{V, U}
     t_restart::Real
 end
 
-function RestartCondition(system::AbstractSystem, filename; output_directory="out",
-                          precondition_values=nothing)
+function RestartCondition(system::AbstractSystem, filename; output_directory="out")
     if !occursin(vtkname(system), basename(splitext(filename)[1]))
         throw(ArgumentError("Filename '$filename' does not seem to correspond to system of type $(nameof(typeof(system)))."))
     end
 
-    restart_data = vtk2trixi(joinpath(output_directory, filename))
+    restart_file = joinpath(output_directory, filename)
+    restart_data = vtk2trixi(restart_file)
     v_restart = restart_v(system, restart_data)
     u_restart = restart_u(system, restart_data)
 
-    if !isnothing(precondition_values)
-        precondition_system!(system, precondition_values)
-    end
+    precondition_system!(system, restart_file)
 
     return RestartCondition(v_restart, u_restart, restart_data.time)
 end
@@ -87,33 +85,4 @@ function write_density_and_pressure!(v_restart, system::EntropicallyDampedSPHSys
     return v_restart
 end
 
-precondition_system!(system, values) = system
-
-function precondition_system!(system::OpenBoundarySystem, values)
-    (; pressure_reference_values, boundary_zones_flow_rate) = system.cache
-
-    if !haskey(values, :pressure_reference_values) &&
-       !haskey(values, :boundary_zones_flow_rate)
-        throw(ArgumentError("Missing required fields in `values` for `OpenBoundarySystem`"))
-    end
-
-    foreach_noalloc(pressure_reference_values,
-                    values.pressure_reference_values) do (pressure_model,
-                                                          pressure_restart)
-        if isa(pressure_model, AbstractPressureModel)
-            pressure_model.pressure[] = pressure_restart
-        else
-            if pressure_restart != Inf
-                throw(ArgumentError("Expected `Inf` for non-pressure-model boundary zones"))
-            end
-        end
-    end
-
-    foreach_noalloc(boundary_zones_flow_rate,
-                    values.boundary_zones_flow_rate) do (flow_rate,
-                                                         flow_rate_restart)
-        flow_rate[] = flow_rate_restart
-    end
-
-    return system
-end
+precondition_system!(system, restart_file) = system
diff --git a/src/io/write_vtk.jl b/src/io/write_vtk.jl
index ff2f2f85eb..14d3d0100e 100644
--- a/src/io/write_vtk.jl
+++ b/src/io/write_vtk.jl
@@ -419,6 +419,14 @@ function write2vtk!(vtk, v, u, t, system::OpenBoundarySystem)
         vtk["Q_total"] = Q_total
     end
 
+    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
+
     return vtk
 end
 
diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl
index 507dcd29f5..1e5040d8c3 100644
--- a/src/schemes/boundary/open_boundary/system.jl
+++ b/src/schemes/boundary/open_boundary/system.jl
@@ -704,3 +704,21 @@ function restart_v(system::OpenBoundarySystem, data)
 
     return v_total
 end
+
+function precondition_system!(system::OpenBoundarySystem, file)
+    if any(pm -> isa(pm, AbstractPressureModel), system.cache.pressure_reference_values)
+        keys_p = [(Symbol(:p, i), "boundary_zone_pressure_$(i)")
+                  for i in eachindex(system.boundary_zones)]
+        keys_Q = [(Symbol(:Q, i), "Q_$(i)") for i in eachindex(system.boundary_zones)]
+        values = vtk2trixi(file; merge(keys_p, keys_Q)...)
+
+        for (i, pressure_model) in enumerate(system.cache.pressure_reference_values)
+            if pressure_model isa AbstractPressureModel
+                pressure_model.pressure[] = values[Symbol(:p, i)]
+                pressure_model.flow_rate[] = values[Symbol(:Q, i)]
+            end
+        end
+    end
+
+    return system
+end

From 616a2dc839bfe52237bded78ecf38b4bdc64ccbd Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Tue, 6 Jan 2026 18:58:48 +0100
Subject: [PATCH 34/62] fix inactive particles

---
 src/general/restart.jl                       |  4 +++
 src/schemes/boundary/open_boundary/system.jl | 28 +++++++++++++++-----
 src/schemes/fluid/fluid.jl                   |  7 -----
 3 files changed, 26 insertions(+), 13 deletions(-)

diff --git a/src/general/restart.jl b/src/general/restart.jl
index e4c85ec0c6..4c5cfdfd09 100644
--- a/src/general/restart.jl
+++ b/src/general/restart.jl
@@ -45,6 +45,8 @@ function set_intial_conditions!(v0_ode, u0_ode, semi, restart_conditions)
     end
 
     foreach_noalloc(semi.systems, restart_conditions) do (system, restart_condition)
+        initialize_before_restart!(system, restart_condition, semi)
+
         v0_system = wrap_v(v0_ode, system, semi)
         u0_system = wrap_u(u0_ode, system, semi)
 
@@ -86,3 +88,5 @@ function write_density_and_pressure!(v_restart, system::EntropicallyDampedSPHSys
 end
 
 precondition_system!(system, restart_file) = system
+
+initialize_before_restart!(system, restart_condition, semi) = system
diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl
index 1e5040d8c3..4498c7e9dc 100644
--- a/src/schemes/boundary/open_boundary/system.jl
+++ b/src/schemes/boundary/open_boundary/system.jl
@@ -134,6 +134,26 @@ function initialize!(system::OpenBoundarySystem, semi)
     return system
 end
 
+function initialize_before_restart!(system::OpenBoundarySystem, restart_condition, semi)
+    set_zero!(system.boundary_zone_indices)
+
+    # 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.
+    for particle in axes(restart_condition.u_restart, 2)
+        particle_coords = current_coords(restart_condition.u_restart, system, particle)
+
+        for (zone_id, boundary_zone) in enumerate(system.boundary_zones)
+            # Check if boundary particle is in the boundary zone
+            if is_in_boundary_zone(boundary_zone, particle_coords)
+                system.boundary_zone_indices[particle] = zone_id
+            end
+        end
+    end
+
+    return system
+end
+
 function create_cache_open_boundary(boundary_model, fluid_system, initial_condition,
                                     calculate_flow_rate, boundary_zones)
     reference_values = map(bz -> bz.reference_values, boundary_zones)
@@ -265,6 +285,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]
@@ -683,9 +705,6 @@ end
 function restart_v(system::OpenBoundarySystem, data)
     v_total = zeros(eltype(system), v_nvariables(system),
                     n_integrated_particles(system))
-    v_total .= eltype(system)(1e16)
-
-    system.buffer.active_particle .= false
 
     v_active = zeros(eltype(system), v_nvariables(system), size(data.velocity, 2))
 
@@ -694,14 +713,11 @@ function restart_v(system::OpenBoundarySystem, data)
                                 density_calculator(system), data.pressure, data.density)
 
     for particle in axes(v_active, 2)
-        system.buffer.active_particle[particle] = true
         for i in axes(v_active, 1)
             v_total[i, particle] = v_active[i, particle]
         end
     end
 
-    update_system_buffer!(system.buffer)
-
     return v_total
 end
 
diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl
index 22a5ad4c03..f3b93d2bb9 100644
--- a/src/schemes/fluid/fluid.jl
+++ b/src/schemes/fluid/fluid.jl
@@ -243,10 +243,6 @@ end
 function restart_v(system::AbstractFluidSystem, data)
     velocity_total = zeros(eltype(system), v_nvariables(system),
                            n_integrated_particles(system))
-    velocity_total .= eltype(system)(1e16)
-
-    isnothing(buffer(system)) || (system.buffer.active_particle .= false)
-
     velocity_active = zeros(eltype(system), v_nvariables(system), size(data.velocity, 2))
 
     velocity_active[1:ndims(system), :] = data.velocity
@@ -254,14 +250,11 @@ function restart_v(system::AbstractFluidSystem, data)
                                 data.pressure, data.density)
 
     for particle in axes(velocity_active, 2)
-        isnothing(buffer(system)) || (system.buffer.active_particle[particle] = true)
         for i in axes(velocity_active, 1)
             velocity_total[i, particle] = velocity_active[i, particle]
         end
     end
 
-    update_system_buffer!(system.buffer)
-
     return velocity_total
 end
 

From a80518cefa8cbd8dfc14042540a2cbe4afd51caa Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Tue, 6 Jan 2026 22:00:39 +0100
Subject: [PATCH 35/62] fix zone IDs

---
 src/general/restart.jl                       | 39 ++++++++++++++++++--
 src/general/semidiscretization.jl            | 22 ++++++++---
 src/io/write_vtk.jl                          |  2 +
 src/schemes/boundary/open_boundary/system.jl | 34 ++++++-----------
 src/schemes/fluid/fluid.jl                   |  7 +++-
 5 files changed, 71 insertions(+), 33 deletions(-)

diff --git a/src/general/restart.jl b/src/general/restart.jl
index 4c5cfdfd09..b5a6d5717e 100644
--- a/src/general/restart.jl
+++ b/src/general/restart.jl
@@ -45,8 +45,6 @@ function set_intial_conditions!(v0_ode, u0_ode, semi, restart_conditions)
     end
 
     foreach_noalloc(semi.systems, restart_conditions) do (system, restart_condition)
-        initialize_before_restart!(system, restart_condition, semi)
-
         v0_system = wrap_v(v0_ode, system, semi)
         u0_system = wrap_u(u0_ode, system, semi)
 
@@ -89,4 +87,39 @@ end
 
 precondition_system!(system, restart_file) = system
 
-initialize_before_restart!(system, restart_condition, semi) = system
+function initialize_neighborhood_searches!(semi, u0_ode, restart_conditions)
+    foreach_system(semi) do system
+        foreach_system(semi) do 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_restart_coordinates(system, u0_ode, semi),
+                                       initial_restart_coordinates(neighbor, u0_ode, semi),
+                                       eachindex_y=each_active_particle(neighbor),
+                                       parallelization_backend=PolyesterBackend())
+        end
+    end
+
+    return semi
+end
+
+function initial_restart_coordinates(system, u0_ode, semi)
+    return 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_conditions)
+    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 dc991aae95..71531aa402 100644
--- a/src/general/semidiscretization.jl
+++ b/src/general/semidiscretization.jl
@@ -333,7 +333,7 @@ function semidiscretize(semi, tspan; reset_threads=true, restart_conditions=noth
 
     # 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_conditions)
 
     if semi.parallelization_backend isa KernelAbstractions.Backend
         # Convert all arrays to the correct array type.
@@ -357,10 +357,7 @@ function semidiscretize(semi, tspan; reset_threads=true, restart_conditions=noth
     end
 
     # Initialize all particle systems
-    foreach_system(semi_new) do system
-        # Initialize this system
-        initialize!(system, semi_new)
-    end
+    initialize!(semi_new, restart_conditions)
 
     # Reset callback flag that will be set by the `UpdateCallback`
     semi_new.update_callback_used[] = false
@@ -381,6 +378,15 @@ end
 
 time_span(tspan, restart_conditions::Nothing) = tspan
 
+function initialize!(semi::Semidiscretization, restart_conditions::Nothing)
+    foreach_system(semi) do system
+        # Initialize this system
+        initialize!(system, semi)
+    end
+
+    return semi
+end
+
 """
     restart_with!(semi, sol)
 
@@ -416,6 +422,10 @@ function restart_with!(semi, sol; reset_threads=true)
     return semi
 end
 
+function initialize_neighborhood_searches!(semi, u0_ode, restart_conditions::Nothing)
+    initialize_neighborhood_searches!(semi)
+end
+
 function initialize_neighborhood_searches!(semi)
     foreach_system(semi) do system
         foreach_system(semi) do neighbor
@@ -1112,7 +1122,7 @@ function check_configuration(system::ImplicitIncompressibleSPHSystem, systems, n
                 neighbor.boundary_model.density_calculator isa PressureBoundaries)
                 time_step_boundary = neighbor.boundary_model.density_calculator.time_step
                 omega_boundary = neighbor.boundary_model.density_calculator.omega
-                if !(time_step==time_step_boundary && omega==omega_boundary)
+                if !(time_step == time_step_boundary && omega == omega_boundary)
                     throw(ArgumentError("`PressureBoundaries` parameters have to be the same as the
                     `ImplicitIncompressibleSPHSystem` parameters"))
                 end
diff --git a/src/io/write_vtk.jl b/src/io/write_vtk.jl
index 14d3d0100e..74b99b5c36 100644
--- a/src/io/write_vtk.jl
+++ b/src/io/write_vtk.jl
@@ -408,6 +408,8 @@ 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 system.cache.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 4498c7e9dc..6270999880 100644
--- a/src/schemes/boundary/open_boundary/system.jl
+++ b/src/schemes/boundary/open_boundary/system.jl
@@ -134,25 +134,8 @@ function initialize!(system::OpenBoundarySystem, semi)
     return system
 end
 
-function initialize_before_restart!(system::OpenBoundarySystem, restart_condition, semi)
-    set_zero!(system.boundary_zone_indices)
-
-    # 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.
-    for particle in axes(restart_condition.u_restart, 2)
-        particle_coords = current_coords(restart_condition.u_restart, system, particle)
-
-        for (zone_id, boundary_zone) in enumerate(system.boundary_zones)
-            # Check if boundary particle is in the boundary zone
-            if is_in_boundary_zone(boundary_zone, particle_coords)
-                system.boundary_zone_indices[particle] = zone_id
-            end
-        end
-    end
-
-    return system
-end
+# Skip during restart, as boundary zone indices are updated in `precondition_system!`
+initialize_restart!(system::OpenBoundarySystem, semi) = system
 
 function create_cache_open_boundary(boundary_model, fluid_system, initial_condition,
                                     calculate_flow_rate, boundary_zones)
@@ -686,17 +669,17 @@ function restart_u(system::OpenBoundarySystem, data)
     coords_total = zeros(eltype(system), u_nvariables(system),
                          n_integrated_particles(system))
     coords_total .= eltype(system)(1e16)
-    system.buffer.active_particle .= false
 
     coords_active = data.coordinates
-
     for particle in axes(coords_active, 2)
-        system.buffer.active_particle[particle] = true
         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
@@ -722,6 +705,13 @@ function restart_v(system::OpenBoundarySystem, data)
 end
 
 function precondition_system!(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; zone_id="zone_id")
+    system.boundary_zone_indices[each_integrated_particle(system)] .= values.zone_id
+
     if any(pm -> isa(pm, AbstractPressureModel), system.cache.pressure_reference_values)
         keys_p = [(Symbol(:p, i), "boundary_zone_pressure_$(i)")
                   for i in eachindex(system.boundary_zones)]
diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl
index f3b93d2bb9..fbdef99f1b 100644
--- a/src/schemes/fluid/fluid.jl
+++ b/src/schemes/fluid/fluid.jl
@@ -224,17 +224,20 @@ function restart_u(system::AbstractFluidSystem, data)
     coords_total = zeros(eltype(system), u_nvariables(system),
                          n_integrated_particles(system))
     coords_total .= eltype(system)(1e16)
-    isnothing(buffer(system)) || (system.buffer.active_particle .= false)
 
     coords_active = data.coordinates
 
     for particle in axes(coords_active, 2)
-        isnothing(buffer(system)) || (system.buffer.active_particle[particle] = true)
         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

From 439cb6b9e0d73c12045faf72da970297fe2a1b77 Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Tue, 6 Jan 2026 22:02:07 +0100
Subject: [PATCH 36/62] add test

---
 test/general/general.jl |  1 +
 test/general/restart.jl | 41 +++++++++++++++++++++++++++++++++++++++++
 2 files changed, 42 insertions(+)
 create mode 100644 test/general/restart.jl

diff --git a/test/general/general.jl b/test/general/general.jl
index 5a75b9078f..e194d8350e 100644
--- a/test/general/general.jl
+++ b/test/general/general.jl
@@ -6,3 +6,4 @@ include("interpolation.jl")
 include("buffer.jl")
 include("util.jl")
 include("custom_quantities.jl")
+include("restart.jl")
diff --git a/test/general/restart.jl b/test/general/restart.jl
new file mode 100644
index 0000000000..c2ddc609c9
--- /dev/null
+++ b/test/general/restart.jl
@@ -0,0 +1,41 @@
+@trixi_testset "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, wall_distance / 2]
+    end_point = [flow_length - 10 * particle_spacing, wall_distance / 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 = RestartCondition(fluid_system, "fluid_1_$iter.vtu")
+    open_boundary_restart = RestartCondition(open_boundary, "open_boundary_1_$iter.vtu")
+    boundary_restart = RestartCondition(boundary_system, "boundary_1_$iter.vtu")
+
+    ode_restart = semidiscretize(semi, (0.3, 0.6);
+                                 restart_conditions=(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=8e-3)
+    @test isapprox(result_full.density, result_restart.density, rtol=8e-4)
+    @test isapprox(result_full.pressure, result_restart.pressure, rtol=7e-2)
+end

From f3049a479e6e389288bad69abe2a83873c494cb8 Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Tue, 6 Jan 2026 22:47:19 +0100
Subject: [PATCH 37/62] make it gpu compatible

---
 src/general/restart.jl | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

diff --git a/src/general/restart.jl b/src/general/restart.jl
index b5a6d5717e..abf8af3bc9 100644
--- a/src/general/restart.jl
+++ b/src/general/restart.jl
@@ -16,7 +16,9 @@ function RestartCondition(system::AbstractSystem, filename; output_directory="ou
 
     precondition_system!(system, restart_file)
 
-    return RestartCondition(v_restart, u_restart, restart_data.time)
+    t_restart = convert(eltype(system), restart_data.time)
+
+    return RestartCondition(v_restart, u_restart, t_restart)
 end
 
 function Base.show(io::IO, rc::RestartCondition)
@@ -48,8 +50,8 @@ function set_intial_conditions!(v0_ode, u0_ode, semi, restart_conditions)
         v0_system = wrap_v(v0_ode, system, semi)
         u0_system = wrap_u(u0_ode, system, semi)
 
-        u0_system .= restart_condition.u_restart
-        v0_system .= restart_condition.v_restart
+        v0_system .= Adapt.adapt(semi.parallelization_backend, restart_condition.v_restart)
+        u0_system .= Adapt.adapt(semi.parallelization_backend, restart_condition.u_restart)
     end
 end
 
@@ -105,7 +107,8 @@ function initialize_neighborhood_searches!(semi, u0_ode, restart_conditions)
 end
 
 function initial_restart_coordinates(system, u0_ode, semi)
-    return wrap_u(u0_ode, system, 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,

From 1c7cd3810b7abac24bb2f6d6bf34bed4aca189f0 Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Tue, 6 Jan 2026 22:50:00 +0100
Subject: [PATCH 38/62] add gpu tests

---
 test/examples/gpu.jl | 56 ++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 56 insertions(+)

diff --git a/test/examples/gpu.jl b/test/examples/gpu.jl
index 0e2655437c..4672bae45c 100644
--- a/test/examples/gpu.jl
+++ b/test/examples/gpu.jl
@@ -645,4 +645,60 @@ 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 = RestartCondition(semi_new.systems[1], "fluid_1_$iter.vtu")
+        open_boundary_restart = RestartCondition(semi_new.systems[2],
+                                                 "open_boundary_1_$iter.vtu")
+        boundary_restart = RestartCondition(semi_new.systems[3], "boundary_1_$iter.vtu")
+
+        ode_restart = semidiscretize(semi_new, (0.3f0, 0.6f0);
+                                     restart_conditions=(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

From c605502510426561473106939e85d242e7045798 Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Wed, 7 Jan 2026 08:44:20 +0100
Subject: [PATCH 39/62] fix tests

---
 src/io/read_vtk.jl | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/src/io/read_vtk.jl b/src/io/read_vtk.jl
index 6f0701f560..f4f4f61781 100644
--- a/src/io/read_vtk.jl
+++ b/src/io/read_vtk.jl
@@ -64,7 +64,10 @@ function vtk2trixi(file; element_type=:default, coordinates_eltype=Float64,
         if idx !== nothing
             results[field] = convert.(ELTYPE, ReadVTK.get_data(point_data[all_keys[idx]]))
         else
-            @info "No '$field' field found in VTK file"
+            # Use zeros as default values when a field is missing
+            results[field] = string(field) in ["velocity"] ?
+                             zero(coordinates) : zeros(ELTYPE, size(coordinates, 2))
+            @info "No '$field' field found in VTK file. Will be set to zero."
         end
     end
 

From 1cada7704c4028c3e78b36e0f4354d84d1ac331b Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Wed, 7 Jan 2026 09:24:23 +0100
Subject: [PATCH 40/62] add checks

---
 src/general/restart.jl                       | 25 +++++++++++++++-----
 src/schemes/boundary/open_boundary/system.jl |  4 ++--
 src/schemes/fluid/fluid.jl                   |  4 ++--
 3 files changed, 23 insertions(+), 10 deletions(-)

diff --git a/src/general/restart.jl b/src/general/restart.jl
index abf8af3bc9..24e387313a 100644
--- a/src/general/restart.jl
+++ b/src/general/restart.jl
@@ -1,7 +1,8 @@
 struct RestartCondition{V, U}
-    v_restart::V
-    u_restart::U
-    t_restart::Real
+    system    :: AbstractSystem
+    v_restart :: V
+    u_restart :: U
+    t_restart :: Real
 end
 
 function RestartCondition(system::AbstractSystem, filename; output_directory="out")
@@ -18,13 +19,13 @@ function RestartCondition(system::AbstractSystem, filename; output_directory="ou
 
     t_restart = convert(eltype(system), restart_data.time)
 
-    return RestartCondition(v_restart, u_restart, t_restart)
+    return RestartCondition(system, v_restart, u_restart, t_restart)
 end
 
 function Base.show(io::IO, rc::RestartCondition)
     @nospecialize rc # reduce precompilation time
 
-    print(io, "RestartCondition{}()")
+    print(io, "RestartCondition{$(nameof(typeof(rc.system)))}()")
 end
 
 function Base.show(io::IO, ::MIME"text/plain", rc::RestartCondition)
@@ -34,18 +35,30 @@ function Base.show(io::IO, ::MIME"text/plain", rc::RestartCondition)
         show(io, rc)
     else
         summary_header(io, "RestartCondition")
+        summary_line(io, "System", "$(nameof(typeof(rc.system)))")
         summary_line(io, "#particles u", "$(size(rc.u_restart, 2))")
         summary_line(io, "#particles v", "$(size(rc.v_restart, 2))")
-        summary_line(io, "eltype", "$(eltype(rc.v_restart))")
+        summary_line(io, "eltype u", "$(eltype(rc.u_restart))")
+        summary_line(io, "eltype v", "$(eltype(rc.v_restart))")
         summary_footer(io)
     end
 end
 
 function set_intial_conditions!(v0_ode, u0_ode, semi, restart_conditions)
+    # Check number of systems
     if length(semi.systems) != length(restart_conditions)
         throw(ArgumentError("Number of systems in `semi` does not match number of `restart_conditions`"))
     end
 
+    # Check that systems match
+    foreach_system(semi) do system
+        system_index = system_indices(system, semi)
+        if !(system == restart_conditions[system_index].system)
+            throw(ArgumentError("System at index $system_index in `semi` does not match system in `restart_conditions`"))
+        end
+    end
+
+    # Set initial conditions
     foreach_noalloc(semi.systems, restart_conditions) do (system, restart_condition)
         v0_system = wrap_v(v0_ode, system, semi)
         u0_system = wrap_u(u0_ode, system, semi)
diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl
index 6270999880..3524e1cfe6 100644
--- a/src/schemes/boundary/open_boundary/system.jl
+++ b/src/schemes/boundary/open_boundary/system.jl
@@ -666,9 +666,9 @@ function interpolate_velocity!(system::OpenBoundarySystem, boundary_zone,
 end
 
 function restart_u(system::OpenBoundarySystem, data)
-    coords_total = zeros(eltype(system), u_nvariables(system),
+    coords_total = zeros(coordinates_eltype(system), u_nvariables(system),
                          n_integrated_particles(system))
-    coords_total .= eltype(system)(1e16)
+    coords_total .= coordinates_eltype(system)(1e16)
 
     coords_active = data.coordinates
     for particle in axes(coords_active, 2)
diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl
index fbdef99f1b..a9fd06f4e6 100644
--- a/src/schemes/fluid/fluid.jl
+++ b/src/schemes/fluid/fluid.jl
@@ -221,9 +221,9 @@ end
 end
 
 function restart_u(system::AbstractFluidSystem, data)
-    coords_total = zeros(eltype(system), u_nvariables(system),
+    coords_total = zeros(coordinates_eltype(system), u_nvariables(system),
                          n_integrated_particles(system))
-    coords_total .= eltype(system)(1e16)
+    coords_total .= coordinates_eltype(system)(1e16)
 
     coords_active = data.coordinates
 

From 100590e5e77097c0ec0618c92341b20659e4f150 Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Wed, 7 Jan 2026 09:26:44 +0100
Subject: [PATCH 41/62] add example

---
 .../restart_poiseuille_flow_2d.jl             | 32 +++++++++++++++++++
 1 file changed, 32 insertions(+)
 create mode 100644 examples/postprocessing/restart_poiseuille_flow_2d.jl

diff --git a/examples/postprocessing/restart_poiseuille_flow_2d.jl b/examples/postprocessing/restart_poiseuille_flow_2d.jl
new file mode 100644
index 0000000000..ff3701f972
--- /dev/null
+++ b/examples/postprocessing/restart_poiseuille_flow_2d.jl
@@ -0,0 +1,32 @@
+# ==========================================================================================
+# Restart Example
+#
+# TODO
+# ==========================================================================================
+using TrixiParticles
+
+# 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)
+
+# 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 = RestartCondition(fluid_system, "fluid_1_$iter.vtu")
+open_boundary_restart = RestartCondition(open_boundary, "open_boundary_1_$iter.vtu")
+boundary_restart = RestartCondition(boundary_system, "boundary_1_$iter.vtu")
+
+ode_restart = semidiscretize(semi, (0.3, 0.6);
+                             restart_conditions=(fluid_restart, open_boundary_restart,
+                                                 boundary_restart))
+
+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)

From 358848df8793a9b043875b0299b6c566588673f6 Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Wed, 7 Jan 2026 09:35:25 +0100
Subject: [PATCH 42/62] add docs

---
 src/general/restart.jl            | 13 +++++++++++++
 src/general/semidiscretization.jl |  5 ++++-
 2 files changed, 17 insertions(+), 1 deletion(-)

diff --git a/src/general/restart.jl b/src/general/restart.jl
index 24e387313a..e49cd9b873 100644
--- a/src/general/restart.jl
+++ b/src/general/restart.jl
@@ -1,3 +1,16 @@
+"""
+    RestartCondition(system::AbstractSystem, filename; output_directory="out")
+
+Create a `RestartCondition` for the given `system` from the specified `filename`.
+
+# Arguments
+- `system`: The system to restart.
+- `filename`: The name of the file from which to restart. This file should be in
+              VTK format and correspond to the specified `system`.
+
+# Keywords
+- `output_directory`: The directory where the restart file is located (default: `"out"`).
+"""
 struct RestartCondition{V, U}
     system    :: AbstractSystem
     v_restart :: V
diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl
index 71531aa402..163f0fc77f 100644
--- a/src/general/semidiscretization.jl
+++ b/src/general/semidiscretization.jl
@@ -253,7 +253,7 @@ end
 @inline foreach_system(f, systems) = foreach_noalloc(f, systems)
 
 """
-    semidiscretize(semi, tspan; reset_threads=true)
+    semidiscretize(semi, tspan; reset_threads=true, restart_conditions=nothing)
 
 Create an `ODEProblem` from the semidiscretization with the specified `tspan`.
 
@@ -262,6 +262,9 @@ Create an `ODEProblem` from the semidiscretization with the specified `tspan`.
 - `tspan`: The time span over which the simulation will be run.
 
 # Keywords
+- `restart_conditions`: A tuple of [`RestartCondition`](@ref) objects specifying
+  from which files to restart each system. The order has to match the order of systems
+  in `semi`. By default, no restart is performed.
 - `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.

From 2b56dc8dc8e4e486e1632c4e5fdeaf8bccb761d5 Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Wed, 7 Jan 2026 09:37:58 +0100
Subject: [PATCH 43/62] fix typo

---
 src/general/restart.jl            | 2 +-
 src/general/semidiscretization.jl | 4 ++--
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/general/restart.jl b/src/general/restart.jl
index e49cd9b873..75de2d8822 100644
--- a/src/general/restart.jl
+++ b/src/general/restart.jl
@@ -57,7 +57,7 @@ function Base.show(io::IO, ::MIME"text/plain", rc::RestartCondition)
     end
 end
 
-function set_intial_conditions!(v0_ode, u0_ode, semi, restart_conditions)
+function set_initial_conditions!(v0_ode, u0_ode, semi, restart_conditions)
     # Check number of systems
     if length(semi.systems) != length(restart_conditions)
         throw(ArgumentError("Number of systems in `semi` does not match number of `restart_conditions`"))
diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl
index 163f0fc77f..35bcf55a29 100644
--- a/src/general/semidiscretization.jl
+++ b/src/general/semidiscretization.jl
@@ -332,7 +332,7 @@ function semidiscretize(semi, tspan; reset_threads=true, restart_conditions=noth
     end
 
     # Set initial condition
-    set_intial_conditions!(v0_ode, u0_ode, semi, restart_conditions)
+    set_initial_conditions!(v0_ode, u0_ode, semi, restart_conditions)
 
     # TODO initialize after adapting to the GPU.
     # Requires https://github.com/trixi-framework/PointNeighbors.jl/pull/86.
@@ -369,7 +369,7 @@ function semidiscretize(semi, tspan; reset_threads=true, restart_conditions=noth
                                time_span(tspan, restart_conditions), semi_new)
 end
 
-function set_intial_conditions!(v0_ode, u0_ode, semi, restart_conditions::Nothing)
+function set_initial_conditions!(v0_ode, u0_ode, semi, restart_conditions::Nothing)
     foreach_system(semi) do system
         v0_system = wrap_v(v0_ode, system, semi)
         u0_system = wrap_u(u0_ode, system, semi)

From 7c087617d516094258948c276cbd6edd8f16d0b8 Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Wed, 7 Jan 2026 12:11:52 +0100
Subject: [PATCH 44/62] implement suggestions

---
 .../boundary/open_boundary/boundary_zones.jl  | 34 +++++++++++--------
 src/schemes/boundary/open_boundary/system.jl  | 16 ++++-----
 2 files changed, 27 insertions(+), 23 deletions(-)

diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl
index 8e42b7088a..75fab302cc 100644
--- a/src/schemes/boundary/open_boundary/boundary_zones.jl
+++ b/src/schemes/boundary/open_boundary/boundary_zones.jl
@@ -79,14 +79,20 @@ There are three ways to specify the actual shape of the boundary zone:
 - `rest_pressure=0.0`: For `BoundaryModelDynamicalPressureZhang`, a rest pressure is required when the pressure is not prescribed.
                        This should match the rest pressure of the fluid system.
                        Per default it is set to zero (assuming a gauge pressure system).
-                       - For `EntropicallyDampedSPHSystem`: Use the initial pressure from the `InitialCondition`
-                       - For `WeaklyCompressibleSPHSystem`: Use the background pressure from the equation of state
-- `sample_points=:default`: Either `:default` to automatically generate sample points on the boundary face (default),
-                            or a matrix of dimensions `(ndims, npoints)` containing sample points
-                            on the boundary face used to compute the volumetric flow rate.
-                            Each sample point represents a discrete area of `particle_spacing^(ndims-1)`.
-                            Therefore, `sample_points` must form a regular grid.
-                            Set to `nothing` to skip sampling.
+    - For `EntropicallyDampedSPHSystem`: Use the initial pressure from the `InitialCondition`
+    - For `WeaklyCompressibleSPHSystem`: Use the background pressure from the equation of state
+- `sample_points=:default`: Controls how the boundary face is sampled to estimate the
+                            volumetric flow rate.
+    - `:default` (default): Automatically generate sampling points
+      on the boundary face using a regular grid with spacing `particle_spacing`.
+    - `sample_points::AbstractMatrix`: Provide your own sampling points
+      as a matrix of size `(ndims, npoints)`, where **each column is
+      one point** on the boundary face (line in 2D, rectangle/face in 3D).
+      The points must lie on the face and form a regular grid.
+    - `nothing`: Disable flow-rate sampling.
+    Note: Each sampling point represents an area element of size
+    `particle_spacing^(ndims-1)`. Therefore, the discretized
+    sampled area is `npoints * particle_spacing^(ndims-1)`.
 
 !!! note "Note"
     The reference values (`reference_velocity`, `reference_pressure`, `reference_density`)
@@ -330,9 +336,9 @@ function create_cache_boundary_zone(initial_condition, boundary_face, face_norma
                                   direction=(-face_normal), n_extrude=1).coordinates
         sample_points_ = convert.(eltype(initial_condition), points)
     else
-        if !(sample_points isa Matrix && size(sample_points, 1) == ndims(initial_condition))
-            throw(ArgumentError("`sample_points` must be a matrix with " *
-                                "`ndims(initial_condition)` rows"))
+        if !(sample_points isa AbstractMatrix &&
+             size(sample_points, 1) == ndims(initial_condition))
+            throw(ArgumentError("`sample_points` must be an `(ndims, npoints)` matrix with $(ndims(initial_condition)) rows"))
         end
         sample_points_ = convert.(eltype(initial_condition), sample_points)
     end
@@ -347,9 +353,9 @@ function create_cache_boundary_zone(initial_condition, boundary_face, face_norma
         face_area = norm(v2 - v1)
     end
 
-    # We only check if the discretized area exceeds the boundary face area.
-    # For 3D boundary zones with complex or non-rectangular flow profiles
-    # (e.g., pipe flow), the cross-sectional area can legitimately be smaller than the boundary face area.
+    # We only warn when the discretized sampled area is larger than the geometric face area.
+    # In 3D, the effective flow cross-section may be smaller than the boundary face, especially for
+    # complex or non-rectangular profiles (e.g. pipe flow), so smaller sampled areas are acceptable.
     if discrete_face_area > (face_area + eps(face_area))
         @warn "The sampled area of the boundary face " *
               "($(discrete_face_area)) is larger than the actual face area ($(face_area)). "
diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl
index 0e09c669c0..08f35b263d 100644
--- a/src/schemes/boundary/open_boundary/system.jl
+++ b/src/schemes/boundary/open_boundary/system.jl
@@ -355,12 +355,8 @@ function calculate_flow_rate!(system::OpenBoundarySystem{<:Any, ELTYPE, NDIMS},
 
             # Compute volumetric flow rate: Q = ∫ v ⋅ n dA
             velocities = reinterpret(reshape, SVector{NDIMS, ELTYPE}, sample_velocity)
-            current_flow_rate = sum(velocities) do velocity
-                vn = dot(velocity, -face_normal)
-                return vn * area_increment
-            end
-
-            flow_rate[] = current_flow_rate
+            flow_rate[] = @inbounds area_increment *
+                                    sum(v -> dot(v, -face_normal), velocities)
         end
     end
 
@@ -624,6 +620,8 @@ function interpolate_velocity!(system::OpenBoundarySystem, boundary_zone,
     set_zero!(shepard_coefficient)
     set_zero!(sample_velocity)
 
+    # Shepard-normalized interpolation:
+    #   v(p) = (Σ_b v_b V_b W_pb) / (Σ_b V_b W_pb)
     foreach_system(semi) do neighbor_system
         v_neighbor = wrap_v(v_ode, neighbor_system, semi)
         u_neighbor = wrap_u(u_ode, neighbor_system, semi)
@@ -639,15 +637,15 @@ function interpolate_velocity!(system::OpenBoundarySystem, boundary_zone,
             shepard_coefficient[point] += volume_b * W_ab
 
             velocity_neighbor = viscous_velocity(v_neighbor, neighbor_system, neighbor)
-            for i in axes(velocity_neighbor, 1)
+            @inbounds for i in axes(velocity_neighbor, 1)
                 sample_velocity[i, point] += velocity_neighbor[i] * volume_b * W_ab
             end
         end
     end
 
     @threaded semi for point in axes(sample_points, 2)
-        if shepard_coefficient[point] > eps()
-            for i in axes(sample_velocity, 1)
+        if shepard_coefficient[point] > eps(eltype(shepard_coefficient))
+            @inbounds for i in axes(sample_velocity, 1)
                 sample_velocity[i, point] /= shepard_coefficient[point]
             end
         end

From 4a6bca19fe3e61bc8d030d88ba115a8beb819abb Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Wed, 7 Jan 2026 12:28:50 +0100
Subject: [PATCH 45/62] implement suggestions

---
 src/io/read_vtk.jl | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/io/read_vtk.jl b/src/io/read_vtk.jl
index 0638c5c21a..37f5cd7d0c 100644
--- a/src/io/read_vtk.jl
+++ b/src/io/read_vtk.jl
@@ -7,8 +7,9 @@ Load VTK file and convert data to an [`InitialCondition`](@ref).
 - `file`: Name of the VTK file to be loaded.
 
 # Keywords
-- `element_type`: Element type for particle fields (`:default` keeps the type
-    stored in the VTK file, otherwise converted to the given type).
+- `element_type`: Element type for particle fields. By default, the type
+                  stored in the VTK file is used.
+                  Otherwise, data is converted to the specified type.
 - `coordinates_eltype`: Element type for particle coordinates (defaults to `Float64`).
 
 !!! warning "Experimental Implementation"

From f8f788c5bfcc3af1d0722cb990e0f512a7e5a852 Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Wed, 7 Jan 2026 14:33:53 +0100
Subject: [PATCH 46/62] implement suggestions

---
 src/io/read_vtk.jl  | 12 +++++++-----
 test/io/read_vtk.jl |  2 +-
 2 files changed, 8 insertions(+), 6 deletions(-)

diff --git a/src/io/read_vtk.jl b/src/io/read_vtk.jl
index 37f5cd7d0c..15618afa72 100644
--- a/src/io/read_vtk.jl
+++ b/src/io/read_vtk.jl
@@ -1,5 +1,5 @@
 """
-	vtk2trixi(file::String; element_type=:default, coordinates_eltype=Float64)
+	vtk2trixi(file::String; element_type=nothing, coordinates_eltype=nothing)
 
 Load VTK file and convert data to an [`InitialCondition`](@ref).
 
@@ -10,7 +10,9 @@ Load VTK file and convert data to an [`InitialCondition`](@ref).
 - `element_type`: Element type for particle fields. By default, the type
                   stored in the VTK file is used.
                   Otherwise, data is converted to the specified type.
-- `coordinates_eltype`: Element type for particle coordinates (defaults to `Float64`).
+- `coordinates_eltype`: Element type for particle coordinates. By default, the type
+                        stored in the VTK file is used.
+                        Otherwise, data is converted to the specified type.
 
 !!! warning "Experimental Implementation"
     This is an experimental feature and may change in any future releases.
@@ -39,7 +41,7 @@ ic = vtk2trixi(joinpath("out", "rectangular.vtu"))
 └──────────────────────────────────────────────────────────────────────────────────────────────────┘
 ```
 """
-function vtk2trixi(file; element_type=:default, coordinates_eltype=Float64)
+function vtk2trixi(file; element_type=nothing, coordinates_eltype=nothing)
     vtk_file = ReadVTK.VTKFile(file)
 
     # Retrieve data fields (e.g., pressure, velocity, ...)
@@ -47,8 +49,8 @@ function vtk2trixi(file; element_type=:default, coordinates_eltype=Float64)
     field_data = ReadVTK.get_field_data(vtk_file)
     point_coords = ReadVTK.get_points(vtk_file)
 
-    cELTYPE = coordinates_eltype
-    ELTYPE = element_type === :default ? eltype(point_coords) : element_type
+    cELTYPE = isnothing(coordinates_eltype) ? eltype(point_coords) : coordinates_eltype
+    ELTYPE = isnothing(element_type) ? eltype(point_coords) : element_type
 
     # Retrieve fields
     ndims = first(ReadVTK.get_data(field_data["ndims"]))
diff --git a/test/io/read_vtk.jl b/test/io/read_vtk.jl
index 1331fdc9b9..cdb350931b 100644
--- a/test/io/read_vtk.jl
+++ b/test/io/read_vtk.jl
@@ -38,7 +38,7 @@
                 @test isapprox(expected_ic_32.density, test_ic.density, rtol=1e-5)
                 @test isapprox(expected_ic_32.pressure, test_ic.pressure, rtol=1e-5)
                 @test eltype(test_ic) === Float32
-                @test eltype(test_ic.coordinates) === Float64
+                @test eltype(test_ic.coordinates) === Float32
             end
 
             @testset verbose=true "Custom Element Type" begin

From 9e4276adc8a1e5a13544cd71ea147a4f5b79b9bf Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Wed, 7 Jan 2026 15:30:49 +0100
Subject: [PATCH 47/62] use strings instead of RestartConditin

---
 .../restart_poiseuille_flow_2d.jl             | 12 ++-
 src/TrixiParticles.jl                         |  2 +-
 src/general/restart.jl                        | 96 +++++--------------
 src/general/semidiscretization.jl             | 27 +++---
 test/examples/gpu.jl                          | 12 +--
 test/general/restart.jl                       | 10 +-
 6 files changed, 57 insertions(+), 102 deletions(-)

diff --git a/examples/postprocessing/restart_poiseuille_flow_2d.jl b/examples/postprocessing/restart_poiseuille_flow_2d.jl
index ff3701f972..fdde06d4b4 100644
--- a/examples/postprocessing/restart_poiseuille_flow_2d.jl
+++ b/examples/postprocessing/restart_poiseuille_flow_2d.jl
@@ -16,13 +16,15 @@ trixi_include(@__MODULE__,
               tspan=(0.0, 0.3), sound_speed_factor=10, particle_spacing=4e-5)
 
 iter = round(Int, 0.3 / 0.02)
-fluid_restart = RestartCondition(fluid_system, "fluid_1_$iter.vtu")
-open_boundary_restart = RestartCondition(open_boundary, "open_boundary_1_$iter.vtu")
-boundary_restart = RestartCondition(boundary_system, "boundary_1_$iter.vtu")
+
+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_conditions=(fluid_restart, open_boundary_restart,
-                                                 boundary_restart))
+                             restart_with=(restart_file_fluid,
+                                           restart_file_open_boundary,
+                                           restart_file_boundary))
 
 saving_callback = SolutionSavingCallback(dt=0.02, prefix="restart")
 
diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl
index 0a5d0d5447..4740e88b65 100644
--- a/src/TrixiParticles.jl
+++ b/src/TrixiParticles.jl
@@ -62,7 +62,7 @@ include("io/io.jl")
 include("general/restart.jl")
 include("visualization/recipes_plots.jl")
 
-export Semidiscretization, semidiscretize, restart_with!, RestartCondition
+export Semidiscretization, semidiscretize, restart_with!
 export InitialCondition
 export WeaklyCompressibleSPHSystem, EntropicallyDampedSPHSystem, TotalLagrangianSPHSystem,
        WallBoundarySystem, DEMSystem, BoundaryDEMSystem, OpenBoundarySystem,
diff --git a/src/general/restart.jl b/src/general/restart.jl
index 75de2d8822..d8488a24f1 100644
--- a/src/general/restart.jl
+++ b/src/general/restart.jl
@@ -1,88 +1,41 @@
-"""
-    RestartCondition(system::AbstractSystem, filename; output_directory="out")
-
-Create a `RestartCondition` for the given `system` from the specified `filename`.
-
-# Arguments
-- `system`: The system to restart.
-- `filename`: The name of the file from which to restart. This file should be in
-              VTK format and correspond to the specified `system`.
-
-# Keywords
-- `output_directory`: The directory where the restart file is located (default: `"out"`).
-"""
-struct RestartCondition{V, U}
-    system    :: AbstractSystem
-    v_restart :: V
-    u_restart :: U
-    t_restart :: Real
-end
-
-function RestartCondition(system::AbstractSystem, filename; output_directory="out")
-    if !occursin(vtkname(system), basename(splitext(filename)[1]))
-        throw(ArgumentError("Filename '$filename' does not seem to correspond to system of type $(nameof(typeof(system)))."))
-    end
-
-    restart_file = joinpath(output_directory, filename)
-    restart_data = vtk2trixi(restart_file)
-    v_restart = restart_v(system, restart_data)
-    u_restart = restart_u(system, restart_data)
-
-    precondition_system!(system, restart_file)
-
-    t_restart = convert(eltype(system), restart_data.time)
-
-    return RestartCondition(system, v_restart, u_restart, t_restart)
-end
-
-function Base.show(io::IO, rc::RestartCondition)
-    @nospecialize rc # reduce precompilation time
-
-    print(io, "RestartCondition{$(nameof(typeof(rc.system)))}()")
-end
-
-function Base.show(io::IO, ::MIME"text/plain", rc::RestartCondition)
-    @nospecialize rc # reduce precompilation time
-
-    if get(io, :compact, false)
-        show(io, rc)
-    else
-        summary_header(io, "RestartCondition")
-        summary_line(io, "System", "$(nameof(typeof(rc.system)))")
-        summary_line(io, "#particles u", "$(size(rc.u_restart, 2))")
-        summary_line(io, "#particles v", "$(size(rc.v_restart, 2))")
-        summary_line(io, "eltype u", "$(eltype(rc.u_restart))")
-        summary_line(io, "eltype v", "$(eltype(rc.v_restart))")
-        summary_footer(io)
-    end
-end
-
-function set_initial_conditions!(v0_ode, u0_ode, semi, restart_conditions)
+function set_initial_conditions!(v0_ode, u0_ode, semi, restart_with::Tuple{Vararg{String}})
     # Check number of systems
-    if length(semi.systems) != length(restart_conditions)
-        throw(ArgumentError("Number of systems in `semi` does not match number of `restart_conditions`"))
+    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)
-        if !(system == restart_conditions[system_index].system)
-            throw(ArgumentError("System at index $system_index in `semi` does not match system in `restart_conditions`"))
+        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_conditions) do (system, restart_condition)
+    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)
 
-        v0_system .= Adapt.adapt(semi.parallelization_backend, restart_condition.v_restart)
-        u0_system .= Adapt.adapt(semi.parallelization_backend, restart_condition.u_restart)
+        restart_data = vtk2trixi(restart_file)
+        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)
+
+        precondition_system!(system, restart_file)
     end
 end
 
-function time_span(tspan, restart_conditions)
-    t_restart = first(restart_conditions).t_restart
+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"
@@ -115,7 +68,8 @@ end
 
 precondition_system!(system, restart_file) = system
 
-function initialize_neighborhood_searches!(semi, u0_ode, restart_conditions)
+function initialize_neighborhood_searches!(semi, u0_ode,
+                                           restart_with::Tuple{Vararg{String}})
     foreach_system(semi) do system
         foreach_system(semi) do neighbor
             # TODO Initialize after adapting to the GPU.
@@ -142,7 +96,7 @@ function initial_restart_coordinates(system::Union{WallBoundarySystem,
     return initial_coordinates(system)
 end
 
-function initialize!(semi::Semidiscretization, restart_conditions)
+function initialize!(semi::Semidiscretization, restart_with::Tuple{Vararg{String}})
     foreach_system(semi) do system
         # Initialize this system
         initialize_restart!(system, semi)
diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl
index 35bcf55a29..0e68818e26 100644
--- a/src/general/semidiscretization.jl
+++ b/src/general/semidiscretization.jl
@@ -253,7 +253,7 @@ end
 @inline foreach_system(f, systems) = foreach_noalloc(f, systems)
 
 """
-    semidiscretize(semi, tspan; reset_threads=true, restart_conditions=nothing)
+    semidiscretize(semi, tspan; reset_threads=true, restart_with=nothing)
 
 Create an `ODEProblem` from the semidiscretization with the specified `tspan`.
 
@@ -262,9 +262,10 @@ Create an `ODEProblem` from the semidiscretization with the specified `tspan`.
 - `tspan`: The time span over which the simulation will be run.
 
 # Keywords
-- `restart_conditions`: A tuple of [`RestartCondition`](@ref) objects specifying
-  from which files to restart each system. The order has to match the order of systems
-  in `semi`. By default, no restart is performed.
+- `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.
@@ -291,7 +292,7 @@ timespan: (0.0, 1.0)
 u0: ([...], [...]) *this line is ignored by filter*
 ```
 """
-function semidiscretize(semi, tspan; reset_threads=true, restart_conditions=nothing)
+function semidiscretize(semi, tspan; reset_threads=true, restart_with=nothing)
     (; systems) = semi
 
     # Check that all systems have the same eltype
@@ -332,11 +333,11 @@ function semidiscretize(semi, tspan; reset_threads=true, restart_conditions=noth
     end
 
     # Set initial condition
-    set_initial_conditions!(v0_ode, u0_ode, semi, restart_conditions)
+    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, u0_ode, restart_conditions)
+    initialize_neighborhood_searches!(semi, u0_ode, restart_with)
 
     if semi.parallelization_backend isa KernelAbstractions.Backend
         # Convert all arrays to the correct array type.
@@ -360,16 +361,16 @@ function semidiscretize(semi, tspan; reset_threads=true, restart_conditions=noth
     end
 
     # Initialize all particle systems
-    initialize!(semi_new, restart_conditions)
+    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,
-                               time_span(tspan, restart_conditions), semi_new)
+                               time_span(tspan, restart_with), semi_new)
 end
 
-function set_initial_conditions!(v0_ode, u0_ode, semi, restart_conditions::Nothing)
+function set_initial_conditions!(v0_ode, u0_ode, semi, restart_with::Nothing)
     foreach_system(semi) do system
         v0_system = wrap_v(v0_ode, system, semi)
         u0_system = wrap_u(u0_ode, system, semi)
@@ -379,9 +380,9 @@ function set_initial_conditions!(v0_ode, u0_ode, semi, restart_conditions::Nothi
     end
 end
 
-time_span(tspan, restart_conditions::Nothing) = tspan
+time_span(tspan, restart_with::Nothing) = tspan
 
-function initialize!(semi::Semidiscretization, restart_conditions::Nothing)
+function initialize!(semi::Semidiscretization, restart_with::Nothing)
     foreach_system(semi) do system
         # Initialize this system
         initialize!(system, semi)
@@ -425,7 +426,7 @@ function restart_with!(semi, sol; reset_threads=true)
     return semi
 end
 
-function initialize_neighborhood_searches!(semi, u0_ode, restart_conditions::Nothing)
+function initialize_neighborhood_searches!(semi, u0_ode, restart_with::Nothing)
     initialize_neighborhood_searches!(semi)
 end
 
diff --git a/test/examples/gpu.jl b/test/examples/gpu.jl
index 4672bae45c..8dcfb422f7 100644
--- a/test/examples/gpu.jl
+++ b/test/examples/gpu.jl
@@ -678,15 +678,13 @@ end
         semi_new = TrixiParticles.Adapt.adapt(Array, sol.prob.p)
 
         iter = round(Int, 0.3 / 0.02)
-        fluid_restart = RestartCondition(semi_new.systems[1], "fluid_1_$iter.vtu")
-        open_boundary_restart = RestartCondition(semi_new.systems[2],
-                                                 "open_boundary_1_$iter.vtu")
-        boundary_restart = RestartCondition(semi_new.systems[3], "boundary_1_$iter.vtu")
+        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_conditions=(fluid_restart,
-                                                         open_boundary_restart,
-                                                         boundary_restart))
+                                     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,
diff --git a/test/general/restart.jl b/test/general/restart.jl
index c2ddc609c9..1af3496707 100644
--- a/test/general/restart.jl
+++ b/test/general/restart.jl
@@ -19,13 +19,13 @@
                   tspan=(0.0, 0.3), sound_speed_factor=10, particle_spacing=4e-5)
 
     iter = round(Int, 0.3 / 0.02)
-    fluid_restart = RestartCondition(fluid_system, "fluid_1_$iter.vtu")
-    open_boundary_restart = RestartCondition(open_boundary, "open_boundary_1_$iter.vtu")
-    boundary_restart = RestartCondition(boundary_system, "boundary_1_$iter.vtu")
+    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_conditions=(fluid_restart, open_boundary_restart,
-                                                     boundary_restart))
+                                 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())

From 7564985ed55bed65d8aa5c59cfd6805c470d3ed6 Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Wed, 7 Jan 2026 16:30:43 +0100
Subject: [PATCH 48/62] add test mixed types

---
 test/io/read_vtk.jl | 14 ++++++++++++++
 1 file changed, 14 insertions(+)

diff --git a/test/io/read_vtk.jl b/test/io/read_vtk.jl
index cdb350931b..ceef5aede4 100644
--- a/test/io/read_vtk.jl
+++ b/test/io/read_vtk.jl
@@ -54,6 +54,20 @@
                 @test eltype(test_ic) === Float32
                 @test eltype(test_ic.coordinates) === Float32
             end
+
+            @testset verbose=true "Custom Element Type Mixed" begin
+                trixi2vtk(expected_ic; filename="tmp_initial_condition_64",
+                          output_directory=tmp_dir)
+                file = joinpath(tmp_dir, "tmp_initial_condition_64.vtu")
+                test_ic = vtk2trixi(file, element_type=Float32, coordinates_eltype=Float64)
+
+                @test isapprox(expected_ic.coordinates, test_ic.coordinates, rtol=1e-5)
+                @test isapprox(expected_ic.velocity, test_ic.velocity, rtol=1e-5)
+                @test isapprox(expected_ic.density, test_ic.density, rtol=1e-5)
+                @test isapprox(expected_ic.pressure, test_ic.pressure, rtol=1e-5)
+                @test eltype(test_ic) === Float32
+                @test eltype(test_ic.coordinates) === Float64
+            end
         end
 
         @testset verbose=true "`AbstractFluidSystem`" begin

From cfba71cec03787d0e36e9d64c4703917fcd8ea50 Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Sat, 10 Jan 2026 14:52:17 +0100
Subject: [PATCH 49/62] add NEWS entry

---
 NEWS.md | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/NEWS.md b/NEWS.md
index a026fa85b3..dac772ece5 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -10,6 +10,8 @@ used in the Julia ecosystem. Notable changes will be documented in this file for
 
 - Keyword argument `n_clamped_particles` of the `TotalLagrangianSPHSystem`
   has been deprecated in favor of a new kwarg `clamped_particles`.
+- Return type of `vtk2trixi` changed to `NamedTuple` including an optional
+  `:initial_condition` field if `create_initial_condition=true` is passed. (#959)
 
 ### Features
 

From 2dcbc76c4cbd480aaa66d5f51149ee8697226f6f Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Sat, 10 Jan 2026 14:52:36 +0100
Subject: [PATCH 50/62] revise

---
 src/io/read_vtk.jl  |  60 +++++++++++------
 test/io/read_vtk.jl | 159 +++++++++++++++++++++++++++++---------------
 2 files changed, 147 insertions(+), 72 deletions(-)

diff --git a/src/io/read_vtk.jl b/src/io/read_vtk.jl
index de095aaf00..2fb3ec4324 100644
--- a/src/io/read_vtk.jl
+++ b/src/io/read_vtk.jl
@@ -1,15 +1,14 @@
 """
-	vtk2trixi(file::String; custom_quantities...)
+	vtk2trixi(file::String; element_type=nothing, coordinates_eltype=nothing,
+              create_initial_condition=true, custom_quantities...)
 
-Load VTK file and convert data to a NamedTuple.
+Read a VTK file and return a `NamedTuple` with keys
+`:coordinates`, `:velocity`, `:density`, `:pressure`, `:particle_spacing`, `:time`,
+plus any requested custom quantities.
+Missing fields are zero-filled; `:particle_spacing` is scalar if constant, otherwise per-particle.
 
 # Arguments
-- `file`:                 Name of the VTK file to be loaded.
-- `custom_quantities...`: Additional custom quantities to be loaded from the VTK file.
-                          Each custom quantity must be explicitly listed in the
-                          `custom_quantities` during the simulation to ensure it is
-                          included in the VTK output and can be successfully loaded.
-                          See [Custom Quantities](@ref custom_quantities) for details.
+- `file`: Name of the VTK file to be loaded.
 
 # Keywords
 - `element_type`: Element type for particle fields. By default, the type
@@ -18,12 +17,20 @@ Load VTK file and convert data to a NamedTuple.
 - `coordinates_eltype`: Element type for particle coordinates. By default, the type
                         stored in the VTK file is used.
                         Otherwise, data is converted to the specified type.
+- `create_initial_condition`: If `true`, an `InitialCondition` object is created
+                              and included in the returned `NamedTuple` under
+                              the key `:initial_condition`. Default is `true`.
+- `custom_quantities...`: Keyword arguments to load additional quantities from the VTK file.
+                          Each keyword becomes a key in the returned `NamedTuple`, with its
+                          string value specifying the VTK field name to read.
+                          Example: `my_data="field_name"` loads VTK field `"field_name"`
+                          as `:my_data` in the result.
 
 !!! warning "Experimental Implementation"
     This is an experimental feature and may change in any future releases.
 
 # Example
-```jldoctest; output = false, filter = r"density = \\[.*\\]|pressure = \\[.*\\]|mass = \\[.*\\]|velocity = \\[.*\\]|coordinates = \\[.*\\]"
+```jldoctest; output = false, filter = r"density = \\[.*\\]|pressure = \\[.*\\]|velocity = \\[.*\\]|coordinates = \\[.*\\]|initial_condition = \\[.*\\]"
 # Create a rectangular shape
 rectangular = RectangularShape(0.1, (10, 10), (0, 0), density=1.5, velocity=(1.0, -2.0),
                                pressure=1000.0)
@@ -37,10 +44,11 @@ data = vtk2trixi(joinpath("out", "rectangular.vtu");
                  my_custom_quantity="my_custom_quantity")
 
 # output
-(particle_spacing = 0.1, density = [...], time = 0.0, pressure = [...], mass = [...], my_custom_quantity = 3.0, velocity = [...], coordinates = [...])
+(particle_spacing = 0.1, density = [...], time = 0.0, pressure = [...], my_custom_quantity = 3.0, velocity = [...], coordinates = [...], initial_condition = [...])
 ```
 """
-function vtk2trixi(file; custom_quantities...)
+function vtk2trixi(file; element_type=nothing, coordinates_eltype=nothing,
+                   create_initial_condition=true, custom_quantities...)
     vtk_file = ReadVTK.VTKFile(file)
 
     # Retrieve data fields (e.g., pressure, velocity, ...)
@@ -58,7 +66,7 @@ function vtk2trixi(file; custom_quantities...)
     ndims = first(ReadVTK.get_data(field_data["ndims"]))
     coordinates = convert.(cELTYPE, point_coords[1:ndims, :])
 
-    fields = [:velocity, :density, :pressure, :mass, :particle_spacing]
+    fields = [:velocity, :density, :pressure, :particle_spacing]
     for field in fields
         # Look for any key that contains the field name
         all_keys = keys(point_data)
@@ -67,8 +75,8 @@ function vtk2trixi(file; custom_quantities...)
             results[field] = convert.(ELTYPE, ReadVTK.get_data(point_data[all_keys[idx]]))
         else
             # Use zeros as default values when a field is missing
-            results[field] = string(field) in ["mass"] ?
-                             zeros(size(coordinates, 2)) : zero(coordinates)
+            results[field] = string(field) in ["velocity"] ?
+                             zero(coordinates) : zeros(ELTYPE, size(coordinates, 2))
             @info "No '$field' field found in VTK file. Will be set to zero."
         end
     end
@@ -78,16 +86,30 @@ function vtk2trixi(file; custom_quantities...)
                                  results[:particle_spacing]
     results[:coordinates] = coordinates
     results[:time] = "time" in keys(field_data) ?
-                     first(ReadVTK.get_data(field_data["time"])) : 0.0
+                     first(ReadVTK.get_data(field_data["time"])) : zero(ELTYPE)
 
-    for (key, quantity) in custom_quantities
+    for (key, quantity_) in custom_quantities
+        quantity = string(quantity_)
         if quantity in keys(point_data)
             results[key] = ReadVTK.get_data(point_data[quantity])
-        end
-        if quantity in keys(field_data)
+        elseif quantity in keys(field_data)
             results[key] = first(ReadVTK.get_data(field_data[quantity]))
+        else
+            throw(ArgumentError("Custom quantity '$quantity' not found in VTK file. " *
+                                "Make sure it was included during the simulation."))
         end
     end
 
-    return NamedTuple(results)
+    results = NamedTuple(results)
+
+    if create_initial_condition
+        ic = InitialCondition(; coordinates=results.coordinates,
+                              particle_spacing=results.particle_spacing,
+                              velocity=results.velocity, density=results.density,
+                              pressure=results.pressure)
+
+        return merge(results, (initial_condition=ic,))
+    else
+        return results
+    end
 end
diff --git a/test/io/read_vtk.jl b/test/io/read_vtk.jl
index c603c597e0..9a3d720553 100644
--- a/test/io/read_vtk.jl
+++ b/test/io/read_vtk.jl
@@ -3,59 +3,118 @@
         coordinates = fill(1.0, 2, 12)
         velocity = fill(2.0, 2, 12)
 
-        expected_ic = InitialCondition(coordinates=coordinates, velocity=velocity,
-                                       density=1000.0, pressure=900.0, mass=50.0)
+        expected_data = InitialCondition(coordinates=coordinates, velocity=velocity,
+                                         density=1000.0, pressure=900.0,
+                                         particle_spacing=0.1)
 
         expected_cq_scalar = 3.0
-        expected_cq_vector = fill(expected_cq_scalar,
-                                  size(expected_ic.coordinates, 2))
+        expected_cq_vector = fill(expected_cq_scalar, nparticles(expected_data))
         scalar_cq_function(system, data, t) = expected_cq_scalar
-        vector_cq_function(system, data,
-                           t) = fill(expected_cq_scalar, nparticles(system))
+        vector_cq_function(system, data, t) = fill(expected_cq_scalar, nparticles(system))
 
         @testset verbose=true "`InitialCondition`" begin
-            @testset verbose=true "Scalar custom quantity" begin
-                trixi2vtk(expected_ic; filename="tmp_initial_condition_scalar",
-                          output_directory=tmp_dir,
-                          cq_scalar=expected_cq_scalar)
+            @testset verbose=true "`Float64`" begin
+                trixi2vtk(expected_data; filename="tmp_initial_condition_64",
+                          output_directory=tmp_dir)
+                file = joinpath(tmp_dir, "tmp_initial_condition_64.vtu")
+                data = vtk2trixi(file)
+
+                @test isapprox(expected_data.coordinates, data.coordinates, rtol=1e-5)
+                @test isapprox(expected_data.velocity, data.velocity, rtol=1e-5)
+                @test isapprox(expected_data.density, data.density, rtol=1e-5)
+                @test isapprox(expected_data.pressure, data.pressure, rtol=1e-5)
+                @test all(key -> eltype(data[key]) === Float64, keys(data))
+                @test eltype(data.coordinates) === Float64
+            end
+
+            @testset verbose=true "`Float32`" begin
+                expected_ic_32 = InitialCondition(;
+                                                  coordinates=convert.(Float32,
+                                                                       coordinates),
+                                                  velocity=convert.(Float32, velocity),
+                                                  density=1000.0f0, pressure=900.0f0,
+                                                  mass=50.0f0, particle_spacing=0.1f0)
+                trixi2vtk(expected_ic_32; filename="tmp_initial_condition_32",
+                          output_directory=tmp_dir)
+                file = joinpath(tmp_dir, "tmp_initial_condition_32.vtu")
+                data = vtk2trixi(file)
+
+                @test isapprox(expected_ic_32.coordinates, data.coordinates, rtol=1e-5)
+                @test isapprox(expected_ic_32.velocity, data.velocity, rtol=1e-5)
+                @test isapprox(expected_ic_32.density, data.density, rtol=1e-5)
+                @test isapprox(expected_ic_32.pressure, data.pressure, rtol=1e-5)
+                @test all(key -> eltype(data[key]) === Float32, keys(data))
+                @test eltype(data.coordinates) === Float32
+            end
+
+            @testset verbose=true "Custom Element Type" begin
+                trixi2vtk(expected_data; filename="tmp_initial_condition_64",
+                          output_directory=tmp_dir)
+                file = joinpath(tmp_dir, "tmp_initial_condition_64.vtu")
+                data = vtk2trixi(file, element_type=Float32, coordinates_eltype=Float32)
+
+                @test isapprox(expected_data.coordinates, data.coordinates, rtol=1e-5)
+                @test isapprox(expected_data.velocity, data.velocity, rtol=1e-5)
+                @test isapprox(expected_data.density, data.density, rtol=1e-5)
+                @test isapprox(expected_data.pressure, data.pressure, rtol=1e-5)
+                @test all(key -> eltype(data[key]) === Float32, keys(data))
+                @test eltype(data.coordinates) === Float32
+            end
+
+            @testset verbose=true "Scalar Custom Quantity" begin
+                trixi2vtk(expected_data; filename="tmp_initial_condition_scalar",
+                          output_directory=tmp_dir, cq_scalar=expected_cq_scalar)
 
                 # Load file containing test data
                 test_data = vtk2trixi(joinpath(tmp_dir, "tmp_initial_condition_scalar.vtu");
-                                      cq_scalar="cq_scalar")
+                                      cq_scalar=:cq_scalar)
 
-                @test isapprox(expected_ic.coordinates, test_data.coordinates, rtol=1e-5)
-                @test isapprox(expected_ic.velocity, test_data.velocity, rtol=1e-5)
-                @test isapprox(expected_ic.density, test_data.density, rtol=1e-5)
-                @test isapprox(expected_ic.pressure, test_data.pressure, rtol=1e-5)
-                @test isapprox(expected_cq_scalar, test_data.cq_scalar,
-                               rtol=1e-5)
+                @test isapprox(expected_data.coordinates, test_data.coordinates, rtol=1e-5)
+                @test isapprox(expected_data.velocity, test_data.velocity, rtol=1e-5)
+                @test isapprox(expected_data.density, test_data.density, rtol=1e-5)
+                @test isapprox(expected_data.pressure, test_data.pressure, rtol=1e-5)
+                @test isapprox(expected_cq_scalar, test_data.cq_scalar, rtol=1e-5)
             end
 
-            @testset verbose=true "Vector custom quantity" begin
-                trixi2vtk(expected_ic; filename="tmp_initial_condition_vector",
+            @testset verbose=true "Vector Custom Quantity" begin
+                trixi2vtk(expected_data; filename="tmp_initial_condition_vector",
                           output_directory=tmp_dir,
                           cq_vector=expected_cq_vector)
 
                 # Load file containing test data
                 test_data = vtk2trixi(joinpath(tmp_dir, "tmp_initial_condition_vector.vtu");
-                                      cq_vector="cq_vector")
+                                      cq_vector=:cq_vector)
 
-                @test isapprox(expected_ic.coordinates, test_data.coordinates, rtol=1e-5)
-                @test isapprox(expected_ic.velocity, test_data.velocity, rtol=1e-5)
-                @test isapprox(expected_ic.density, test_data.density, rtol=1e-5)
-                @test isapprox(expected_ic.pressure, test_data.pressure, rtol=1e-5)
-                @test isapprox(expected_cq_vector, test_data.cq_vector,
-                               rtol=1e-5)
+                @test isapprox(expected_data.coordinates, test_data.coordinates, rtol=1e-5)
+                @test isapprox(expected_data.velocity, test_data.velocity, rtol=1e-5)
+                @test isapprox(expected_data.density, test_data.density, rtol=1e-5)
+                @test isapprox(expected_data.pressure, test_data.pressure, rtol=1e-5)
+                @test isapprox(expected_cq_vector, test_data.cq_vector, rtol=1e-5)
+            end
+
+            @testset verbose=true "Custom Element Type Mixed" begin
+                trixi2vtk(expected_data; filename="tmp_initial_condition_64",
+                          output_directory=tmp_dir)
+                file = joinpath(tmp_dir, "tmp_initial_condition_64.vtu")
+                data = vtk2trixi(file, element_type=Float32, coordinates_eltype=Float64)
+
+                @test isapprox(expected_data.coordinates, data.coordinates, rtol=1e-5)
+                @test isapprox(expected_data.velocity, data.velocity, rtol=1e-5)
+                @test isapprox(expected_data.density, data.density, rtol=1e-5)
+                @test isapprox(expected_data.pressure, data.pressure, rtol=1e-5)
+                @test all(key -> eltype(data[key]) === Float32,
+                          setdiff(keys(data), (:coordinates,)))
+                @test eltype(data.coordinates) === Float64
             end
         end
 
         @testset verbose=true "`AbstractFluidSystem`" begin
-            fluid_system = EntropicallyDampedSPHSystem(expected_ic,
+            fluid_system = EntropicallyDampedSPHSystem(expected_data,
                                                        SchoenbergCubicSplineKernel{2}(),
                                                        1.5, 1.5)
 
             # Overwrite values because we skip the update step
-            fluid_system.cache.density .= expected_ic.density
+            fluid_system.cache.density .= expected_data.density
 
             semi = Semidiscretization(fluid_system)
 
@@ -71,33 +130,30 @@
             vu_ode = (; x)
             dvdu_ode = (; x=(; v_ode=dv_ode, u_ode=du_ode))
 
-            @testset verbose=true "Scalar custom quantity" begin
+            @testset verbose=true "Scalar Custom Quantity" begin
                 trixi2vtk(fluid_system, dvdu_ode, vu_ode, semi, 0.0,
                           nothing; system_name="tmp_file_fluid_scalar",
-                          output_directory=tmp_dir,
-                          iter=1, cq_scalar=scalar_cq_function)
+                          output_directory=tmp_dir, iter=1, cq_scalar=scalar_cq_function)
 
                 # Load file containing test data
                 test_data = vtk2trixi(joinpath(tmp_dir, "tmp_file_fluid_scalar_1.vtu");
-                                      cq_scalar="cq_scalar")
+                                      cq_scalar=:cq_scalar)
 
                 @test isapprox(u, test_data.coordinates, rtol=1e-5)
                 @test isapprox(v, test_data.velocity, rtol=1e-5)
                 @test isapprox(pressure, test_data.pressure, rtol=1e-5)
                 @test isapprox(fluid_system.cache.density, test_data.density, rtol=1e-5)
-                @test isapprox(expected_cq_scalar, test_data.cq_scalar,
-                               rtol=1e-5)
+                @test isapprox(expected_cq_scalar, test_data.cq_scalar, rtol=1e-5)
             end
 
-            @testset verbose=true "Vector custom quantity" begin
+            @testset verbose=true "Vector Custom Quantity" begin
                 trixi2vtk(fluid_system, dvdu_ode, vu_ode, semi, 0.0,
                           nothing; system_name="tmp_file_fluid_vector",
-                          output_directory=tmp_dir,
-                          iter=1, cq_vector=vector_cq_function)
+                          output_directory=tmp_dir, iter=1, cq_vector=vector_cq_function)
 
                 # Load file containing test data
                 test_data = vtk2trixi(joinpath(tmp_dir, "tmp_file_fluid_vector_1.vtu");
-                                      cq_vector="cq_vector")
+                                      cq_vector=:cq_vector)
 
                 @test isapprox(u, test_data.coordinates, rtol=1e-5)
                 @test isapprox(v, test_data.velocity, rtol=1e-5)
@@ -109,17 +165,17 @@
         end
 
         @testset verbose=true "`WallBoundarySystem`" begin
-            boundary_model = BoundaryModelDummyParticles(expected_ic.density,
-                                                         expected_ic.mass,
+            boundary_model = BoundaryModelDummyParticles(expected_data.density,
+                                                         expected_data.mass,
                                                          SummationDensity(),
                                                          SchoenbergCubicSplineKernel{2}(),
                                                          1.5)
 
             # Overwrite values because we skip the update step
-            boundary_model.pressure .= expected_ic.pressure
-            boundary_model.cache.density .= expected_ic.density
+            boundary_model.pressure .= expected_data.pressure
+            boundary_model.cache.density .= expected_data.density
 
-            boundary_system = WallBoundarySystem(expected_ic, boundary_model)
+            boundary_system = WallBoundarySystem(expected_data, boundary_model)
             semi = Semidiscretization(boundary_system)
 
             # Create dummy ODE solutions
@@ -131,15 +187,14 @@
             vu_ode = (; x)
             dvdu_ode = (; x=(; v_ode=dv_ode, u_ode=du_ode))
 
-            @testset verbose=true "Scalar custom quantity" begin
+            @testset verbose=true "Scalar Custom Quantity" begin
                 trixi2vtk(boundary_system, dvdu_ode, vu_ode, semi, 0.0,
                           nothing; system_name="tmp_file_boundary_scalar",
-                          output_directory=tmp_dir,
-                          iter=1, cq_scalar=scalar_cq_function)
+                          output_directory=tmp_dir, iter=1, cq_scalar=scalar_cq_function)
 
                 # Load file containing test data
                 test_data = vtk2trixi(joinpath(tmp_dir, "tmp_file_boundary_scalar_1.vtu");
-                                      cq_scalar="cq_scalar")
+                                      cq_scalar=:cq_scalar)
 
                 @test isapprox(boundary_system.coordinates, test_data.coordinates,
                                rtol=1e-5)
@@ -148,19 +203,17 @@
                                rtol=1e-5)
                 @test isapprox(boundary_model.pressure, test_data.pressure, rtol=1e-5)
                 @test isapprox(boundary_model.cache.density, test_data.density, rtol=1e-5)
-                @test isapprox(expected_cq_scalar, test_data.cq_scalar,
-                               rtol=1e-5)
+                @test isapprox(expected_cq_scalar, test_data.cq_scalar, rtol=1e-5)
             end
 
-            @testset verbose=true "Vector custom quantity" begin
+            @testset verbose=true "Vector Custom Quantity" begin
                 trixi2vtk(boundary_system, dvdu_ode, vu_ode, semi, 0.0,
                           nothing; system_name="tmp_file_boundary_vector",
-                          output_directory=tmp_dir,
-                          iter=1, cq_vector=vector_cq_function)
+                          output_directory=tmp_dir, iter=1, cq_vector=vector_cq_function)
 
                 # Load file containing test data
                 test_data = vtk2trixi(joinpath(tmp_dir, "tmp_file_boundary_vector_1.vtu");
-                                      cq_vector="cq_vector")
+                                      cq_vector=:cq_vector)
 
                 @test isapprox(boundary_system.coordinates, test_data.coordinates,
                                rtol=1e-5)

From 0d587cd5b439e068f148a2589dc670a21b8b48d1 Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Sat, 10 Jan 2026 15:38:02 +0100
Subject: [PATCH 51/62] fix doc test

---
 src/io/read_vtk.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/io/read_vtk.jl b/src/io/read_vtk.jl
index 2fb3ec4324..97d524f8d4 100644
--- a/src/io/read_vtk.jl
+++ b/src/io/read_vtk.jl
@@ -30,7 +30,7 @@ Missing fields are zero-filled; `:particle_spacing` is scalar if constant, other
     This is an experimental feature and may change in any future releases.
 
 # Example
-```jldoctest; output = false, filter = r"density = \\[.*\\]|pressure = \\[.*\\]|velocity = \\[.*\\]|coordinates = \\[.*\\]|initial_condition = \\[.*\\]"
+```jldoctest; output = false, filter = r"density = \\[.*\\]|pressure = \\[.*\\]|velocity = \\[.*\\]|coordinates = \\[.*\\]"
 # Create a rectangular shape
 rectangular = RectangularShape(0.1, (10, 10), (0, 0), density=1.5, velocity=(1.0, -2.0),
                                pressure=1000.0)
@@ -44,7 +44,7 @@ data = vtk2trixi(joinpath("out", "rectangular.vtu");
                  my_custom_quantity="my_custom_quantity")
 
 # output
-(particle_spacing = 0.1, density = [...], time = 0.0, pressure = [...], my_custom_quantity = 3.0, velocity = [...], coordinates = [...], initial_condition = [...])
+(particle_spacing = 0.1, density = [...], time = 0.0, pressure = [...], my_custom_quantity = 3.0, velocity = [...], coordinates = [...], initial_condition = InitialCondition{Float64, Float64}())
 ```
 """
 function vtk2trixi(file; element_type=nothing, coordinates_eltype=nothing,

From 48a79e1c1c78d3dd1ca28935c4ff5430474bcb77 Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Sat, 10 Jan 2026 15:54:16 +0100
Subject: [PATCH 52/62] rm test file

---
 .../boundary/open_boundary/flow_rate.jl       | 64 -------------------
 1 file changed, 64 deletions(-)
 delete mode 100644 test/schemes/boundary/open_boundary/flow_rate.jl

diff --git a/test/schemes/boundary/open_boundary/flow_rate.jl b/test/schemes/boundary/open_boundary/flow_rate.jl
deleted file mode 100644
index 9bbaf0e277..0000000000
--- a/test/schemes/boundary/open_boundary/flow_rate.jl
+++ /dev/null
@@ -1,64 +0,0 @@
-@testset verbose=true "Calculate Flow Rate" begin
-    particle_spacing = 0.01
-
-    # Define a parabolic velocity profile
-    velocity_function(pos) = [4 * (pos[2] - 0.5) * (1.5 - pos[2]), 0.0]
-
-    # Create fluid domain with the specified velocity profile
-    n_particles_y = round(Int, 2 / particle_spacing)
-    fluid = RectangularShape(particle_spacing, (4, n_particles_y), (0.0, 0.0),
-                             density=1000.0, velocity=velocity_function)
-
-    smoothing_length = 1.3 * particle_spacing
-    smoothing_kernel = WendlandC2Kernel{ndims(fluid)}()
-    fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length,
-                                               1.0)
-    fluid_system.cache.density .= fluid.density
-
-    # Use a smaller cross-sectional area to test user-defined area functionality
-    # and to perform interpolation within an embedded domain.
-    # (In a simulation with solid boundaries, wall velocities would also be included.)
-    sample_points = RectangularShape(particle_spacing, (1, round(Int, n_particles_y / 2)),
-                                     (0.0, 0.5), density=1.0).coordinates
-
-    zone = BoundaryZone(; boundary_face=([0.0, 0.0], [0.0, 2.0]),
-                        face_normal=(-1.0, 0.0), open_boundary_layers=10, density=1000.0,
-                        particle_spacing, sample_points=sample_points,
-                        reference_velocity=(pos, t) -> velocity_function(pos))
-
-    open_boundary = OpenBoundarySystem(zone; fluid_system,
-                                       boundary_model=BoundaryModelMirroringTafuni(),
-                                       calculate_flow_rate=true, buffer_size=0)
-
-    semi = Semidiscretization(fluid_system, open_boundary)
-    TrixiParticles.initialize_neighborhood_searches!(semi)
-    TrixiParticles.initialize!(open_boundary, semi)
-
-    # Set up ODE for initial conditions
-    ode = semidiscretize(semi, (0, 1))
-    v_ode, u_ode = ode.u0.x
-
-    boundary_zone = first(open_boundary.boundary_zones)
-
-    @testset verbose=true "Velocity Interpolation" begin
-        TrixiParticles.interpolate_velocity!(open_boundary, boundary_zone,
-                                             v_ode, u_ode, semi)
-
-        coords = reinterpret(reshape, SVector{2, Float64},
-                             boundary_zone.cache.sample_points)
-        vels_analytic = first.(velocity_function.(coords))
-        vels_interpolated = first.(reinterpret(reshape, SVector{2, Float64},
-                                               boundary_zone.cache.sample_velocity))
-
-        @test isapprox(vels_interpolated, vels_analytic, rtol=1e-3)
-    end
-
-    @testset verbose=true "Flow Rate Calculation" begin
-        TrixiParticles.calculate_flow_rate!(open_boundary, v_ode, u_ode, semi)
-
-        Q_analytic = zone.cache.cross_sectional_area * (2 / 3)
-        Q_calculated = first(open_boundary.cache.boundary_zones_flow_rate)[]
-
-        @test isapprox(Q_analytic, Q_calculated; rtol=1e-3)
-    end
-end

From f264c9a6d9d04ea5d56b1af3cbefe0b7cda2df03 Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Sat, 10 Jan 2026 16:00:23 +0100
Subject: [PATCH 53/62] fix merge bugs

---
 src/io/write_vtk.jl                          | 10 ---
 src/schemes/boundary/open_boundary/system.jl | 64 ++++++++++++++++++++
 2 files changed, 64 insertions(+), 10 deletions(-)

diff --git a/src/io/write_vtk.jl b/src/io/write_vtk.jl
index 685ce5b095..6b5b83cd6c 100644
--- a/src/io/write_vtk.jl
+++ b/src/io/write_vtk.jl
@@ -411,16 +411,6 @@ function write2vtk!(vtk, v, u, t, system::OpenBoundarySystem)
     vtk["zone_id"] = [system.boundary_zone_indices[particle]
                       for particle in eachparticle(system)]
 
-    if system.cache.calculate_flow_rate
-        Q_total = zero(eltype(system))
-        for i in eachindex(system.cache.boundary_zones_flow_rate)
-            vtk["Q_$i"] = system.cache.boundary_zones_flow_rate[i][]
-            Q_total += system.cache.boundary_zones_flow_rate[i][]
-        end
-
-        vtk["Q_total"] = Q_total
-    end
-
     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
diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl
index dc3f150cf8..fdae2dd496 100644
--- a/src/schemes/boundary/open_boundary/system.jl
+++ b/src/schemes/boundary/open_boundary/system.jl
@@ -678,3 +678,67 @@ 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 precondition_system!(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; zone_id="zone_id")
+    system.boundary_zone_indices[each_integrated_particle(system)] .= values.zone_id
+
+    if any(pm -> isa(pm, AbstractPressureModel), system.cache.pressure_reference_values)
+        keys_p = [(Symbol(:p, i), "boundary_zone_pressure_$(i)")
+                  for i in eachindex(system.boundary_zones)]
+        keys_Q = [(Symbol(:Q, i), "Q_$(i)") for i in eachindex(system.boundary_zones)]
+        values = vtk2trixi(file; merge(keys_p, keys_Q)...)
+
+        for (i, pressure_model) in enumerate(system.cache.pressure_reference_values)
+            if pressure_model isa AbstractPressureModel
+                pressure_model.pressure[] = values[Symbol(:p, i)]
+                pressure_model.flow_rate[] = values[Symbol(:Q, i)]
+            end
+        end
+    end
+
+    return system
+end

From df63263aad3fb501689e74075f7b62d586912b2d Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Sat, 10 Jan 2026 17:04:29 +0100
Subject: [PATCH 54/62] fix preconditioning

---
 src/schemes/boundary/open_boundary/system.jl | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl
index fdae2dd496..53509faf89 100644
--- a/src/schemes/boundary/open_boundary/system.jl
+++ b/src/schemes/boundary/open_boundary/system.jl
@@ -727,10 +727,12 @@ function precondition_system!(system::OpenBoundarySystem, file)
     system.boundary_zone_indices[each_integrated_particle(system)] .= values.zone_id
 
     if any(pm -> isa(pm, AbstractPressureModel), system.cache.pressure_reference_values)
+        pm_zone_ids = findall(pm -> isa(pm, AbstractPressureModel),
+                              system.cache.pressure_reference_values)
         keys_p = [(Symbol(:p, i), "boundary_zone_pressure_$(i)")
-                  for i in eachindex(system.boundary_zones)]
+                  for i in eachindex(system.boundary_zones)[pm_zone_ids]]
         keys_Q = [(Symbol(:Q, i), "Q_$(i)") for i in eachindex(system.boundary_zones)]
-        values = vtk2trixi(file; merge(keys_p, keys_Q)...)
+        values = vtk2trixi(file; vcat(keys_p, keys_Q)...)
 
         for (i, pressure_model) in enumerate(system.cache.pressure_reference_values)
             if pressure_model isa AbstractPressureModel

From 47bc4c2fce41ecea99f39dd4abfcbb57690b63e2 Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Tue, 13 Jan 2026 08:32:36 +0100
Subject: [PATCH 55/62] don't read IC

---
 src/general/restart.jl                       | 2 +-
 src/schemes/boundary/open_boundary/system.jl | 4 ++--
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/general/restart.jl b/src/general/restart.jl
index d8488a24f1..5ee092ad5c 100644
--- a/src/general/restart.jl
+++ b/src/general/restart.jl
@@ -22,7 +22,7 @@ function set_initial_conditions!(v0_ode, u0_ode, semi, restart_with::Tuple{Varar
         v0_system = wrap_v(v0_ode, system, semi)
         u0_system = wrap_u(u0_ode, system, semi)
 
-        restart_data = vtk2trixi(restart_file)
+        restart_data = vtk2trixi(restart_file; create_initial_condition=false)
         v_restart = restart_v(system, restart_data)
         u_restart = restart_u(system, restart_data)
 
diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl
index 53509faf89..c6adac30c1 100644
--- a/src/schemes/boundary/open_boundary/system.jl
+++ b/src/schemes/boundary/open_boundary/system.jl
@@ -723,7 +723,7 @@ function precondition_system!(system::OpenBoundarySystem, file)
     # 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; zone_id="zone_id")
+    values = vtk2trixi(file; create_initial_condition=false, zone_id="zone_id")
     system.boundary_zone_indices[each_integrated_particle(system)] .= values.zone_id
 
     if any(pm -> isa(pm, AbstractPressureModel), system.cache.pressure_reference_values)
@@ -732,7 +732,7 @@ function precondition_system!(system::OpenBoundarySystem, file)
         keys_p = [(Symbol(:p, i), "boundary_zone_pressure_$(i)")
                   for i in eachindex(system.boundary_zones)[pm_zone_ids]]
         keys_Q = [(Symbol(:Q, i), "Q_$(i)") for i in eachindex(system.boundary_zones)]
-        values = vtk2trixi(file; vcat(keys_p, keys_Q)...)
+        values = vtk2trixi(file; create_initial_condition=false, vcat(keys_p, keys_Q)...)
 
         for (i, pressure_model) in enumerate(system.cache.pressure_reference_values)
             if pressure_model isa AbstractPressureModel

From ae5496a284a978c52c49d60a954868d76a94dd2c Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Thu, 15 Jan 2026 09:24:33 +0100
Subject: [PATCH 56/62] fix merge bug

---
 src/general/restart.jl | 27 +++++++++++++++++++--------
 1 file changed, 19 insertions(+), 8 deletions(-)

diff --git a/src/general/restart.jl b/src/general/restart.jl
index 5ee092ad5c..e5016ba847 100644
--- a/src/general/restart.jl
+++ b/src/general/restart.jl
@@ -72,20 +72,31 @@ function initialize_neighborhood_searches!(semi, u0_ode,
                                            restart_with::Tuple{Vararg{String}})
     foreach_system(semi) do system
         foreach_system(semi) do 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_restart_coordinates(system, u0_ode, semi),
-                                       initial_restart_coordinates(neighbor, u0_ode, semi),
-                                       eachindex_y=each_active_particle(neighbor),
-                                       parallelization_backend=PolyesterBackend())
+            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))

From a1fc021a3cbabd6ed9e90958bb16ea52fe7658fd Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Tue, 20 Jan 2026 13:26:04 +0100
Subject: [PATCH 57/62] fix NEWS entry

---
 NEWS.md | 13 ++++++++-----
 1 file changed, 8 insertions(+), 5 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index afecdb755b..cad5c3ac24 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -4,6 +4,13 @@ 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
+
+### API Changes
+
+- Return type of `vtk2trixi` changed to `NamedTuple` including an optional
+  `:initial_condition` field if `create_initial_condition=true` is passed. (#959)
+
 ## Version 0.4.3
 
 ### API Changes
@@ -27,7 +34,7 @@ used in the Julia ecosystem. Notable changes will be documented in this file for
 
 - Added GPU and FP32 support for DEM (#979).
 
-  
+
 ### Performance
 - Improved GPU performance with shifting up to a factor of 10x (#974, #993).
 
@@ -39,8 +46,6 @@ used in the Julia ecosystem. Notable changes will be documented in this file for
 
 - Keyword argument `n_clamped_particles` of the `TotalLagrangianSPHSystem`
   has been deprecated in favor of a new kwarg `clamped_particles`.
-- Return type of `vtk2trixi` changed to `NamedTuple` including an optional
-  `:initial_condition` field if `create_initial_condition=true` is passed. (#959)
 
 ### Features
 
@@ -355,5 +360,3 @@ Features:
 #### TLSPH
 
 An implementation of TLSPH (Total Lagrangian Smoothed Particle Hydrodynamics) for solid bodies enabling FSI (Fluid Structure Interactions).
-
-

From 8c50353f57bd3fddd50d298455f402658139307d Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Thu, 22 Jan 2026 08:51:26 +0100
Subject: [PATCH 58/62] fix merge bug

---
 NEWS.md | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/NEWS.md b/NEWS.md
index cad5c3ac24..f24457e7e3 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -11,6 +11,12 @@ used in the Julia ecosystem. Notable changes will be documented in this file for
 - Return type of `vtk2trixi` changed to `NamedTuple` including an optional
   `:initial_condition` field if `create_initial_condition=true` is passed. (#959)
 
+## Version 0.4.4
+
+### Features
+
+- Added `StateEquationAdaptiveCole` an adaptive sound speed version of the Cole state equation (#875)
+
 ## Version 0.4.3
 
 ### API Changes
@@ -34,7 +40,6 @@ used in the Julia ecosystem. Notable changes will be documented in this file for
 
 - Added GPU and FP32 support for DEM (#979).
 
-
 ### Performance
 - Improved GPU performance with shifting up to a factor of 10x (#974, #993).
 

From 4b3b91177a793b948e5db6ff903e6d23352758fe Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Thu, 22 Jan 2026 08:57:18 +0100
Subject: [PATCH 59/62] fix formatting

---
 src/general/semidiscretization.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl
index f09231c9e7..85b5e12285 100644
--- a/src/general/semidiscretization.jl
+++ b/src/general/semidiscretization.jl
@@ -1207,7 +1207,7 @@ function check_configuration(system::ImplicitIncompressibleSPHSystem, systems, n
                 neighbor.boundary_model.density_calculator isa PressureBoundaries)
                 time_step_boundary = neighbor.boundary_model.density_calculator.time_step
                 omega_boundary = neighbor.boundary_model.density_calculator.omega
-                if !(time_step == time_step_boundary && omega == omega_boundary)
+                if !(time_step==time_step_boundary && omega==omega_boundary)
                     throw(ArgumentError("`PressureBoundaries` parameters have to be the same as the
                     `ImplicitIncompressibleSPHSystem` parameters"))
                 end

From 55635fca26ca776040ee50013e2574180b7af008 Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Thu, 22 Jan 2026 09:09:49 +0100
Subject: [PATCH 60/62] better name

---
 src/general/restart.jl                               | 4 ++--
 src/schemes/boundary/open_boundary/system.jl         | 4 ++--
 src/schemes/structure/total_lagrangian_sph/system.jl | 5 ++---
 3 files changed, 6 insertions(+), 7 deletions(-)

diff --git a/src/general/restart.jl b/src/general/restart.jl
index e5016ba847..0f0f80dcc6 100644
--- a/src/general/restart.jl
+++ b/src/general/restart.jl
@@ -29,7 +29,7 @@ function set_initial_conditions!(v0_ode, u0_ode, semi, restart_with::Tuple{Varar
         v0_system .= Adapt.adapt(semi.parallelization_backend, v_restart)
         u0_system .= Adapt.adapt(semi.parallelization_backend, u_restart)
 
-        precondition_system!(system, restart_file)
+        restore_previous_state!(system, restart_file)
     end
 end
 
@@ -66,7 +66,7 @@ function write_density_and_pressure!(v_restart, system::EntropicallyDampedSPHSys
     return v_restart
 end
 
-precondition_system!(system, restart_file) = system
+restore_previous_state!(system, restart_file) = system
 
 function initialize_neighborhood_searches!(semi, u0_ode,
                                            restart_with::Tuple{Vararg{String}})
diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl
index c6adac30c1..b63c0f56eb 100644
--- a/src/schemes/boundary/open_boundary/system.jl
+++ b/src/schemes/boundary/open_boundary/system.jl
@@ -145,7 +145,7 @@ function initialize!(system::OpenBoundarySystem, semi)
     return system
 end
 
-# Skip during restart, as boundary zone indices are updated in `precondition_system!`
+# 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,
@@ -718,7 +718,7 @@ function restart_v(system::OpenBoundarySystem, data)
     return v_total
 end
 
-function precondition_system!(system::OpenBoundarySystem, file)
+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.
diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl
index 123730ee99..2208ac53d8 100644
--- a/src/schemes/structure/total_lagrangian_sph/system.jl
+++ b/src/schemes/structure/total_lagrangian_sph/system.jl
@@ -573,13 +573,12 @@ function restart_with!(system::TotalLagrangianSPHSystem, v, u)
 end
 
 function restart_u(system::AbstractStructureSystem, data)
-    # TODO: Is this correct?
-    # data.coordinates = coords_integrated
+    # 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)
-    # TODO: Is this correct?
+    # The integration array requires only the particles that need to be integrated
     return data.velocity[:, each_integrated_particle(system)]
 end
 

From f02d72406fb5d8412635753f5e4013db4ddc5cf1 Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Thu, 22 Jan 2026 09:36:34 +0100
Subject: [PATCH 61/62] add more tests

---
 test/general/restart.jl | 110 ++++++++++++++++++++++++++--------------
 1 file changed, 71 insertions(+), 39 deletions(-)

diff --git a/test/general/restart.jl b/test/general/restart.jl
index 1af3496707..6bc9427547 100644
--- a/test/general/restart.jl
+++ b/test/general/restart.jl
@@ -1,41 +1,73 @@
 @trixi_testset "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, wall_distance / 2]
-    end_point = [flow_length - 10 * particle_spacing, wall_distance / 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=8e-3)
-    @test isapprox(result_full.density, result_restart.density, rtol=8e-4)
-    @test isapprox(result_full.pressure, result_restart.pressure, rtol=7e-2)
+    @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, wall_distance / 2]
+        end_point = [flow_length - 10 * particle_spacing, wall_distance / 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=8e-3)
+        @test isapprox(result_full.density, result_restart.density, rtol=8e-4)
+        @test isapprox(result_full.pressure, result_restart.pressure, rtol=7e-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"),
+                      reference_pressure_out=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

From 6221cbc173bed0e3a99ff9394331d6d6f9badd4a Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.nehe@web.de>
Date: Fri, 17 Apr 2026 14:52:15 +0200
Subject: [PATCH 62/62] revise example

---
 .../postprocessing/restart_poiseuille_flow_2d.jl  | 15 ++++++---------
 1 file changed, 6 insertions(+), 9 deletions(-)

diff --git a/examples/postprocessing/restart_poiseuille_flow_2d.jl b/examples/postprocessing/restart_poiseuille_flow_2d.jl
index fdde06d4b4..0332a1c31d 100644
--- a/examples/postprocessing/restart_poiseuille_flow_2d.jl
+++ b/examples/postprocessing/restart_poiseuille_flow_2d.jl
@@ -1,21 +1,18 @@
 # ==========================================================================================
-# Restart Example
+# Restart Example: Poiseuille Flow 2D
 #
-# TODO
+# 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
 
-# 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)
-
-# 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)
+# 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")