diff --git a/src/callbacks/density_reinit.jl b/src/callbacks/density_reinit.jl
index ded47149f9..79f19c79d0 100644
--- a/src/callbacks/density_reinit.jl
+++ b/src/callbacks/density_reinit.jl
@@ -6,7 +6,9 @@ Callback to reinitialize the density field when using [`ContinuityDensity`](@ref
 # Keywords
 - `interval=0`:              Reinitialize the density every `interval` time steps.
 - `dt`:                      Reinitialize the density in regular intervals of `dt` in terms
-                             of integration time.
+                             of integration time. This callback does not add extra time
+                             steps / `tstops`; instead, reinitialization is triggered at
+                             the first solver step after each `dt` interval has elapsed.
 - `reinit_initial_solution`: Reinitialize the initial solution (default=false)
 """
 mutable struct DensityReinitializationCallback{I}
@@ -102,4 +104,6 @@ function (reinit_callback::DensityReinitializationCallback)(integrator)
     @trixi_timeit timer() "reinit density" reinit_density!(vu_ode, semi)
 
     reinit_callback.last_t = integrator.t
+
+    u_modified!(integrator, true)
 end
diff --git a/src/callbacks/solution_saving.jl b/src/callbacks/solution_saving.jl
index a694b3f5ae..07c1eca306 100644
--- a/src/callbacks/solution_saving.jl
+++ b/src/callbacks/solution_saving.jl
@@ -104,7 +104,9 @@ function SolutionSavingCallback(; interval::Integer=0, dt=0.0,
                                                -1, Ref("UnknownVersion"))
 
     if length(save_times) > 0
-        return PresetTimeCallback(save_times, solution_callback)
+        return PresetTimeCallback(save_times, solution_callback,
+                                  save_positions=(false, false),
+                                  initialize=(initialize_save_cb!))
     elseif dt > 0
         # Add a `tstop` every `dt`, and save the final solution
         return PeriodicCallback(solution_callback, dt,
diff --git a/src/callbacks/steady_state_reached.jl b/src/callbacks/steady_state_reached.jl
index dfbb8ec75f..6390308abc 100644
--- a/src/callbacks/steady_state_reached.jl
+++ b/src/callbacks/steady_state_reached.jl
@@ -50,25 +50,20 @@ function SteadyStateReachedCallback(; interval::Integer=0, dt=0.0,
     end
 end
 
-# `affect!` (`PeriodicCallback`)
+# `affect!` (`PeriodicCallback` and `DiscreteCallback`)
 function (cb::SteadyStateReachedCallback)(integrator)
-    steady_state_condition!(cb, integrator) || return nothing
+    if steady_state_condition!(cb, integrator)
+        print_summary(integrator)
 
-    print_summary(integrator)
-
-    terminate!(integrator)
-end
-
-# `affect!` (`DiscreteCallback`)
-function (cb::SteadyStateReachedCallback{Int})(integrator)
-    print_summary(integrator)
-
-    terminate!(integrator)
+        terminate!(integrator)
+    end
 end
 
 # `condition` (`DiscreteCallback`)
 function (steady_state_callback::SteadyStateReachedCallback)(vu_ode, t, integrator)
-    return steady_state_condition!(steady_state_callback, integrator)
+    (; interval) = steady_state_callback
+
+    return condition_integrator_interval(integrator, interval)
 end
 
 @inline function steady_state_condition!(cb, integrator)
diff --git a/src/general/initial_condition.jl b/src/general/initial_condition.jl
index f69f5b78d1..b78ae8013f 100644
--- a/src/general/initial_condition.jl
+++ b/src/general/initial_condition.jl
@@ -368,8 +368,8 @@ function InitialCondition(sol::ODESolution, system, semi; use_final_velocity=fal
 
     v_ode, u_ode = sol.u[end].x
 
+    v = wrap_v(v_ode, system, semi)
     u = wrap_u(u_ode, system, semi)
-    v = wrap_u(v_ode, system, semi)
 
     # Check if particles come too close especially when the surface exhibits large curvature
     too_close = find_too_close_particles(u, min_particle_distance)
diff --git a/src/general/interpolation.jl b/src/general/interpolation.jl
index fa77dab0e0..a7f01ef419 100644
--- a/src/general/interpolation.jl
+++ b/src/general/interpolation.jl
@@ -555,8 +555,9 @@ end
     ref_id = system_indices(ref_system, semi)
     ref_smoothing_kernel = ref_system.smoothing_kernel
 
-    # If we don't cut at the boundary, we only need to iterate over the reference system
-    systems = cut_off_bnd ? semi : (ref_system,)
+    # If we don't cut at the boundary and don't need the wall velocity of the boundary,
+    # we only need to iterate over the reference system.
+    systems = (cut_off_bnd || include_wall_velocity) ? semi : (ref_system,)
 
     foreach_system(systems) do neighbor_system
         system_id = system_indices(neighbor_system, semi)
@@ -605,8 +606,9 @@ end
     end
 
     @threaded parallelization_backend for point in axes(point_coords, 2)
-        if other_density[point] > computed_density[point] ||
-           computed_density[point] < eps()
+        set_nan = computed_density[point] < eps() ||
+                  (cut_off_bnd && other_density[point] > computed_density[point])
+        if set_nan
             # Return NaN values that can be filtered out in ParaView
             computed_density[point] = NaN
             neighbor_count[point] = 0
diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl
index f84849f096..bc535f18ba 100644
--- a/src/schemes/fluid/weakly_compressible_sph/system.jl
+++ b/src/schemes/fluid/weakly_compressible_sph/system.jl
@@ -392,6 +392,13 @@ end
 
 function reinit_density!(system::WeaklyCompressibleSPHSystem, v, u,
                          v_ode, u_ode, semi)
+    (; density_calculator) = system
+
+    reinit_density!(system, density_calculator, v, u, v_ode, u_ode, semi)
+end
+
+function reinit_density!(system::WeaklyCompressibleSPHSystem, ::ContinuityDensity, v, u,
+                         v_ode, u_ode, semi)
     # Compute density with `SummationDensity` and store the result in `v`,
     # overwriting the previous integrated density.
     summation_density!(system, semi, u, u_ode, v[end, :])
@@ -407,6 +414,11 @@ function reinit_density!(system::WeaklyCompressibleSPHSystem, v, u,
     return system
 end
 
+function reinit_density!(system::WeaklyCompressibleSPHSystem, ::SummationDensity, v, u,
+                         v_ode, u_ode, semi)
+    return system
+end
+
 function reinit_density!(system, v, u, v_ode, u_ode, semi)
     return system
 end