diff --git a/examples/n_body/n_body_system.jl b/examples/n_body/n_body_system.jl
index 06e760fbe..9493f2333 100644
--- a/examples/n_body/n_body_system.jl
+++ b/examples/n_body/n_body_system.jl
@@ -1,8 +1,8 @@
 using TrixiParticles
 using LinearAlgebra
 
-struct NBodySystem{NDIMS, ELTYPE <: Real} <: TrixiParticles.System{NDIMS}
-    initial_condition :: InitialCondition{ELTYPE}
+struct NBodySystem{NDIMS, ELTYPE <: Real, IC} <: TrixiParticles.System{NDIMS}
+    initial_condition :: IC
     mass              :: Array{ELTYPE, 1} # [particle]
     G                 :: ELTYPE
     buffer            :: Nothing
@@ -11,13 +11,13 @@ struct NBodySystem{NDIMS, ELTYPE <: Real} <: TrixiParticles.System{NDIMS}
         mass = copy(initial_condition.mass)
 
         new{size(initial_condition.coordinates, 1),
-            eltype(mass)}(initial_condition, mass, G, nothing)
+            eltype(mass), typeof(initial_condition)}(initial_condition, mass, G, nothing)
     end
 end
 
 TrixiParticles.timer_name(::NBodySystem) = "nbody"
 
-@inline Base.eltype(system::NBodySystem) = eltype(system.initial_condition.coordinates)
+@inline Base.eltype(system::NBodySystem{NDIMS, ELTYPE}) where {NDIMS, ELTYPE} = ELTYPE
 
 function TrixiParticles.write_u0!(u0, system::NBodySystem)
     u0 .= system.initial_condition.coordinates
diff --git a/src/callbacks/density_reinit.jl b/src/callbacks/density_reinit.jl
index 50f54e286..ded47149f 100644
--- a/src/callbacks/density_reinit.jl
+++ b/src/callbacks/density_reinit.jl
@@ -69,7 +69,7 @@ function initialize_reinit_cb!(cb::DensityReinitializationCallback, u, t, integr
         # Update systems to compute quantities like density and pressure.
         semi = integrator.p
         v_ode, u_ode = u.x
-        update_systems_and_nhs(v_ode, u_ode, semi, t; update_from_callback=true)
+        update_systems_and_nhs(v_ode, u_ode, semi, t)
 
         # Apply the callback.
         cb(integrator)
diff --git a/src/callbacks/particle_shifting.jl b/src/callbacks/particle_shifting.jl
index 5ff735fcb..921370f3b 100644
--- a/src/callbacks/particle_shifting.jl
+++ b/src/callbacks/particle_shifting.jl
@@ -39,8 +39,7 @@ function particle_shifting!(integrator)
         # still have the values from the last stage of the previous step if not updated here.
         @trixi_timeit timer() "update systems and nhs" begin
             # Don't create sub-timers here to avoid cluttering the timer output
-            @notimeit timer() update_systems_and_nhs(v_ode, u_ode, semi, t;
-                                                     update_from_callback=true)
+            @notimeit timer() update_systems_and_nhs(v_ode, u_ode, semi, t)
         end
 
         @trixi_timeit timer() "particle shifting" foreach_system(semi) do system
diff --git a/src/callbacks/post_process.jl b/src/callbacks/post_process.jl
index 876401316..36e16ce7f 100644
--- a/src/callbacks/post_process.jl
+++ b/src/callbacks/post_process.jl
@@ -240,8 +240,7 @@ function (pp::PostprocessCallback)(integrator)
         # still have the values from the last stage of the previous step if not updated here.
         @trixi_timeit timer() "update systems and nhs" begin
             # Don't create sub-timers here to avoid cluttering the timer output
-            @notimeit timer() update_systems_and_nhs(v_ode, u_ode, semi, t;
-                                                     update_from_callback=true)
+            @notimeit timer() update_systems_and_nhs(v_ode, u_ode, semi, t)
         end
 
         foreach_system(semi) do system
diff --git a/src/callbacks/update.jl b/src/callbacks/update.jl
index 7cc479d25..b98985220 100644
--- a/src/callbacks/update.jl
+++ b/src/callbacks/update.jl
@@ -57,7 +57,7 @@ function initial_update!(cb::UpdateCallback, u, t, integrator)
 
     # Tell systems that `UpdateCallback` is used
     foreach_system(semi) do system
-        update_callback_used!(system)
+        set_callback_flag!(system, true)
     end
 
     return cb(integrator)
@@ -78,7 +78,7 @@ function (update_callback!::UpdateCallback)(integrator)
 
     # Update quantities that are stored in the systems. These quantities (e.g. pressure)
     # still have the values from the last stage of the previous step if not updated here.
-    update_systems_and_nhs(v_ode, u_ode, semi, t; update_from_callback=true)
+    update_systems_and_nhs(v_ode, u_ode, semi, t)
 
     # Update open boundaries first, since particles might be activated or deactivated
     @trixi_timeit timer() "update open boundary" foreach_system(semi) do system
@@ -89,6 +89,7 @@ function (update_callback!::UpdateCallback)(integrator)
         update_particle_packing(system, v_ode, u_ode, semi, integrator)
     end
 
+    # This is only used by the particle packing system and should be removed in the future
     @trixi_timeit timer() "update TVF" foreach_system(semi) do system
         update_transport_velocity!(system, v_ode, semi)
     end
@@ -141,5 +142,3 @@ function Base.show(io::IO, ::MIME"text/plain",
         summary_box(io, "UpdateCallback", setup)
     end
 end
-
-update_callback_used!(system) = system
diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl
index 3e033e447..cdc8cd34e 100644
--- a/src/general/semidiscretization.jl
+++ b/src/general/semidiscretization.jl
@@ -338,8 +338,8 @@ function semidiscretize(semi, tspan; reset_threads=true)
         # Initialize this system
         initialize!(system, semi_new)
 
-        # Only for systems requiring a mandatory callback
-        reset_callback_flag!(system)
+        # Only for systems requiring the use of the `UpdateCallback`
+        set_callback_flag!(system, false)
     end
 
     return DynamicalODEProblem(kick!, drift!, v0_ode, u0_ode, tspan, semi_new)
@@ -373,8 +373,8 @@ function restart_with!(semi, sol; reset_threads=true)
 
         restart_with!(system, v, u)
 
-        # Only for systems requiring a mandatory callback
-        reset_callback_flag!(system)
+        # Only for systems requiring the use of the `UpdateCallback`
+        set_callback_flag!(system, false)
     end
 
     return semi
@@ -464,17 +464,33 @@ function drift!(du_ode, v_ode, u_ode, semi, t)
 end
 
 @inline function add_velocity!(du, v, particle, system)
+    # Generic fallback for all systems that don't define this function
     for i in 1:ndims(system)
-        du[i, particle] = v[i, particle]
+        @inbounds du[i, particle] = v[i, particle]
     end
 
     return du
 end
 
+# Solid wall boundary system doesn't integrate the particle positions
 @inline add_velocity!(du, v, particle, system::BoundarySPHSystem) = du
 
+@inline function add_velocity!(du, v, particle, system::FluidSystem)
+    # This is zero unless a transport velocity is used
+    delta_v_ = delta_v(system, particle)
+
+    for i in 1:ndims(system)
+        @inbounds du[i, particle] = v[i, particle] + delta_v_[i]
+    end
+
+    return du
+end
+
 function kick!(dv_ode, v_ode, u_ode, semi, t)
     @trixi_timeit timer() "kick!" begin
+        # Check that the `UpdateCallback` is used if required
+        check_update_callback(semi)
+
         @trixi_timeit timer() "reset ∂v/∂t" set_zero!(dv_ode)
 
         @trixi_timeit timer() "update systems and nhs" update_systems_and_nhs(v_ode, u_ode,
@@ -490,8 +506,9 @@ function kick!(dv_ode, v_ode, u_ode, semi, t)
     return dv_ode
 end
 
-# Update the systems and neighborhood searches (NHS) for a simulation before calling `interact!` to compute forces
-function update_systems_and_nhs(v_ode, u_ode, semi, t; update_from_callback=false)
+# Update the systems and neighborhood searches (NHS) for a simulation
+# before calling `interact!` to compute forces.
+function update_systems_and_nhs(v_ode, u_ode, semi, t)
     # First update step before updating the NHS
     # (for example for writing the current coordinates in the solid system)
     foreach_system(semi) do system
@@ -523,12 +540,21 @@ function update_systems_and_nhs(v_ode, u_ode, semi, t; update_from_callback=fals
         update_pressure!(system, v, u, v_ode, u_ode, semi, t)
     end
 
+    # This update depends on the computed quantities of the fluid system and therefore
+    # needs to be after `update_quantities!`.
+    foreach_system(semi) do system
+        v = wrap_v(v_ode, system, semi)
+        u = wrap_u(u_ode, system, semi)
+
+        update_boundary_interpolation!(system, v, u, v_ode, u_ode, semi, t)
+    end
+
     # Final update step for all remaining systems
     foreach_system(semi) do system
         v = wrap_v(v_ode, system, semi)
         u = wrap_u(u_ode, system, semi)
 
-        update_final!(system, v, u, v_ode, u_ode, semi, t; update_from_callback)
+        update_final!(system, v, u, v_ode, u_ode, semi, t)
     end
 end
 
@@ -858,6 +884,16 @@ function update!(neighborhood_search, x, y, semi; points_moving=(true, false),
                            parallelization_backend=semi.parallelization_backend)
 end
 
+function check_update_callback(semi)
+    foreach_system(semi) do system
+        # This check will be optimized away if the system does not require the callback
+        if requires_update_callback(system) && !update_callback_used(system)
+            system_name = system |> typeof |> nameof
+            throw(ArgumentError("`UpdateCallback` is required for `$system_name`"))
+        end
+    end
+end
+
 function check_configuration(systems,
                              nhs::Union{Nothing, PointNeighbors.AbstractNeighborhoodSearch})
     foreach_system(systems) do system
diff --git a/src/general/system.jl b/src/general/system.jl
index 69141a8ff..c75761771 100644
--- a/src/general/system.jl
+++ b/src/general/system.jl
@@ -126,7 +126,7 @@ end
                                  system_correction(system), system, particle)
 end
 
-# System update orders. This can be dispatched if needed.
+# System updates do nothing by default, but can be dispatched if needed
 function update_positions!(system, v, u, v_ode, u_ode, semi, t)
     return system
 end
@@ -139,21 +139,27 @@ function update_pressure!(system, v, u, v_ode, u_ode, semi, t)
     return system
 end
 
-function update_final!(system, v, u, v_ode, u_ode, semi, t; update_from_callback=false)
+function update_boundary_interpolation!(system, v, u, v_ode, u_ode, semi, t)
     return system
 end
 
-# Only for systems requiring a mandatory callback
-reset_callback_flag!(system) = system
+function update_final!(system, v, u, v_ode, u_ode, semi, t)
+    return system
+end
+
+# Only for systems requiring the use of the `UpdateCallback`
+@inline requires_update_callback(system) = false
+@inline update_callback_used(system) = false
+@inline set_callback_flag!(system, value) = system
 
-initial_smoothing_length(system) = smoothing_length(system, nothing)
+@inline initial_smoothing_length(system) = smoothing_length(system, nothing)
 
-function smoothing_length(system, particle)
+@inline function smoothing_length(system, particle)
     return system.smoothing_length
 end
 
-system_smoothing_kernel(system) = system.smoothing_kernel
-system_correction(system) = nothing
+@inline system_smoothing_kernel(system) = system.smoothing_kernel
+@inline system_correction(system) = nothing
 
 @inline particle_spacing(system, particle) = system.initial_condition.particle_spacing
 
diff --git a/src/io/write_vtk.jl b/src/io/write_vtk.jl
index a2b40f31e..b93e25f32 100644
--- a/src/io/write_vtk.jl
+++ b/src/io/write_vtk.jl
@@ -58,8 +58,7 @@ function trixi2vtk(vu_ode, semi, t; iter=nothing, output_directory="out", prefix
     # still have the values from the last stage of the previous step if not updated here.
     @trixi_timeit timer() "update systems" begin
         # Don't create sub-timers here to avoid cluttering the timer output
-        @notimeit timer() update_systems_and_nhs(v_ode, u_ode, semi, t;
-                                                 update_from_callback=true)
+        @notimeit timer() update_systems_and_nhs(v_ode, u_ode, semi, t)
     end
 
     filenames = system_names(systems)
diff --git a/src/preprocessing/particle_packing/system.jl b/src/preprocessing/particle_packing/system.jl
index 1318f6354..fb97501e4 100644
--- a/src/preprocessing/particle_packing/system.jl
+++ b/src/preprocessing/particle_packing/system.jl
@@ -215,14 +215,15 @@ end
     return ndims(system)
 end
 
-function reset_callback_flag!(system::ParticlePackingSystem)
-    system.update_callback_used[] = false
+@inline requires_update_callback(system::ParticlePackingSystem) = true
+@inline update_callback_used(system::ParticlePackingSystem) = system.update_callback_used[]
+
+function set_callback_flag!(system::ParticlePackingSystem, value)
+    system.update_callback_used[] = value
 
     return system
 end
 
-update_callback_used!(system::ParticlePackingSystem) = system.update_callback_used[] = true
-
 function write2vtk!(vtk, v, u, t, system::ParticlePackingSystem; write_meta_data=true)
     vtk["velocity"] = [advection_velocity(v, system, particle)
                        for particle in active_particles(system)]
@@ -289,15 +290,6 @@ function update_position!(u, system::ParticlePackingSystem, semi)
     return u
 end
 
-function update_final!(system::ParticlePackingSystem, v, u, v_ode, u_ode, semi, t;
-                       update_from_callback=false)
-    if !update_from_callback && !(system.update_callback_used[])
-        throw(ArgumentError("`UpdateCallback` is required when using `ParticlePackingSystem`"))
-    end
-
-    return system
-end
-
 # Skip for systems without `SignedDistanceField`
 constrain_particles_onto_surface!(u, system::ParticlePackingSystem{Nothing}, semi) = u
 
@@ -398,6 +390,10 @@ end
     return system
 end
 
+@inline function update_transport_velocity!(system, v_ode, semi)
+    return system
+end
+
 # Skip for fixed systems
 @inline add_velocity!(du, v, particle, system::ParticlePackingSystem{<:Any, true}) = du
 
diff --git a/src/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl
index 1d30a5274..f0051bdb3 100644
--- a/src/schemes/boundary/dummy_particles/dummy_particles.jl
+++ b/src/schemes/boundary/dummy_particles/dummy_particles.jl
@@ -311,17 +311,16 @@ end
                                   system, v, u, v_ode, u_ode, semi)
     (; correction, density_calculator) = boundary_model
 
-    compute_correction_values!(system, correction, u, v_ode, u_ode, semi)
+    compute_pressure!(boundary_model, density_calculator, system, v, u, v_ode, u_ode, semi)
 
+    # These are only computed when using corrections
+    compute_correction_values!(system, correction, u, v_ode, u_ode, semi)
     compute_gradient_correction_matrix!(correction, boundary_model, system, u, v_ode, u_ode,
                                         semi)
-
     # `kernel_correct_density!` only performed for `SummationDensity`
     kernel_correct_density!(boundary_model, v, u, v_ode, u_ode, semi, correction,
                             density_calculator)
 
-    compute_pressure!(boundary_model, density_calculator, system, v, u, v_ode, u_ode, semi)
-
     return boundary_model
 end
 
diff --git a/src/schemes/boundary/open_boundary/method_of_characteristics.jl b/src/schemes/boundary/open_boundary/method_of_characteristics.jl
index addc29210..a37a9cd60 100644
--- a/src/schemes/boundary/open_boundary/method_of_characteristics.jl
+++ b/src/schemes/boundary/open_boundary/method_of_characteristics.jl
@@ -89,7 +89,8 @@ end
 end
 
 # Called from semidiscretization
-function update_final!(system, ::BoundaryModelLastiwka, v, u, v_ode, u_ode, semi, t)
+function update_boundary_model!(system, ::BoundaryModelLastiwka, v, u, v_ode, u_ode,
+                                semi, t)
     @trixi_timeit timer() "evaluate characteristics" begin
         evaluate_characteristics!(system, v, u, v_ode, u_ode, semi, t)
     end
diff --git a/src/schemes/boundary/open_boundary/mirroring.jl b/src/schemes/boundary/open_boundary/mirroring.jl
index 3d9fde0af..45a98ad7c 100644
--- a/src/schemes/boundary/open_boundary/mirroring.jl
+++ b/src/schemes/boundary/open_boundary/mirroring.jl
@@ -125,7 +125,7 @@ function update_boundary_quantities!(system, boundary_model::BoundaryModelTafuni
     end
 end
 
-update_final!(system, ::BoundaryModelTafuni, v, u, v_ode, u_ode, semi, t) = system
+update_boundary_model!(system, ::BoundaryModelTafuni, v, u, v_ode, u_ode, semi, t) = system
 
 function extrapolate_values!(system,
                              mirror_method::Union{FirstOrderMirroring, SimpleMirroring},
diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl
index 19bf85792..1d35bf726 100644
--- a/src/schemes/boundary/open_boundary/system.jl
+++ b/src/schemes/boundary/open_boundary/system.jl
@@ -226,14 +226,15 @@ end
     return ELTYPE
 end
 
-function reset_callback_flag!(system::OpenBoundarySPHSystem)
-    system.update_callback_used[] = false
+@inline requires_update_callback(system::OpenBoundarySPHSystem) = true
+@inline update_callback_used(system::OpenBoundarySPHSystem) = system.update_callback_used[]
+
+function set_callback_flag!(system::OpenBoundarySPHSystem, value)
+    system.update_callback_used[] = value
 
     return system
 end
 
-update_callback_used!(system::OpenBoundarySPHSystem) = system.update_callback_used[] = true
-
 function corresponding_fluid_system(system::OpenBoundarySPHSystem, semi)
     return system.fluid_system
 end
@@ -252,13 +253,9 @@ end
     return system.pressure
 end
 
-function update_final!(system::OpenBoundarySPHSystem, v, u, v_ode, u_ode, semi, t;
-                       update_from_callback=false)
-    if !update_from_callback && !(system.update_callback_used[])
-        throw(ArgumentError("`UpdateCallback` is required when using `OpenBoundarySPHSystem`"))
-    end
-
-    update_final!(system, system.boundary_model, v, u, v_ode, u_ode, semi, t)
+function update_boundary_interpolation!(system::OpenBoundarySPHSystem, v, u, v_ode, u_ode,
+                                        semi, t)
+    update_boundary_model!(system, system.boundary_model, v, u, v_ode, u_ode, semi, t)
 end
 
 # This function is called by the `UpdateCallback`, as the integrator array might be modified
diff --git a/src/schemes/boundary/system.jl b/src/schemes/boundary/system.jl
index bc216983b..1e3a52655 100644
--- a/src/schemes/boundary/system.jl
+++ b/src/schemes/boundary/system.jl
@@ -358,12 +358,13 @@ function update_quantities!(system::BoundarySPHSystem, v, u, v_ode, u_ode, semi,
 end
 
 # This update depends on the computed quantities of the fluid system and therefore
-# has to be in `update_final!` after `update_quantities!`.
-function update_final!(system::BoundarySPHSystem, v, u, v_ode, u_ode, semi, t;
-                       update_from_callback=false)
+# has to be in `update_boundary_interpolation!` after `update_quantities!`.
+function update_boundary_interpolation!(system::BoundarySPHSystem, v, u, v_ode, u_ode,
+                                        semi, t)
     (; boundary_model) = system
 
-    # Note that `update_pressure!(::BoundarySPHSystem, ...)` is empty
+    # Note that `update_pressure!(::BoundarySPHSystem, ...)` is empty,
+    # so no pressure is updated in the previous update steps.
     update_pressure!(boundary_model, system, v, u, v_ode, u_ode, semi)
 
     return system
diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl
index 06e34272b..7a0d1ded5 100644
--- a/src/schemes/fluid/entropically_damped_sph/rhs.jl
+++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl
@@ -50,10 +50,12 @@ function interact!(dv, v_particle_system, u_particle_system,
                                      sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel)
 
         # Add convection term (only when using `TransportVelocityAdami`)
-        dv_convection = momentum_convection(particle_system, neighbor_system,
-                                            v_particle_system, v_neighbor_system,
-                                            rho_a, rho_b, m_a, m_b,
-                                            particle, neighbor, grad_kernel)
+        dv_tvf = dv_transport_velocity(transport_velocity(particle_system),
+                                       particle_system, neighbor_system,
+                                       particle, neighbor,
+                                       v_particle_system, v_neighbor_system,
+                                       m_a, m_b, rho_a, rho_b, pos_diff, distance,
+                                       grad_kernel, correction)
 
         dv_surface_tension = surface_tension_force(surface_tension_a, surface_tension_b,
                                                    particle_system, neighbor_system,
@@ -64,9 +66,9 @@ function interact!(dv, v_particle_system, u_particle_system,
                                      particle, neighbor, pos_diff, distance)
 
         for i in 1:ndims(particle_system)
-            dv[i,
-               particle] += dv_pressure[i] + dv_viscosity_[i] + dv_convection[i] +
-                            dv_surface_tension[i] + dv_adhesion[i]
+            @inbounds dv[i,
+                         particle] += dv_pressure[i] + dv_viscosity_[i] + dv_tvf[i] +
+                                      dv_surface_tension[i] + dv_adhesion[i]
         end
 
         v_diff = current_velocity(v_particle_system, particle_system, particle) -
@@ -76,10 +78,6 @@ function interact!(dv, v_particle_system, u_particle_system,
                             particle, neighbor, pos_diff, distance,
                             sound_speed, m_a, m_b, p_a, p_b, rho_a, rho_b)
 
-        # Apply the transport velocity (only when using a transport velocity)
-        transport_velocity!(dv, particle_system, neighbor_system, particle, neighbor,
-                            m_a, m_b, grad_kernel)
-
         continuity_equation!(dv, density_calculator, v_diff, particle, m_b, rho_a, rho_b,
                              particle_system, grad_kernel)
     end
diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl
index 044104916..f9875e684 100644
--- a/src/schemes/fluid/entropically_damped_sph/system.jl
+++ b/src/schemes/fluid/entropically_damped_sph/system.jl
@@ -130,7 +130,7 @@ function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel,
     nu_edac = (alpha * smoothing_length * sound_speed) / 8
 
     cache = (; create_cache_density(initial_condition, density_calculator)...,
-             create_cache_tvf(Val(:edac), initial_condition, transport_velocity)...,
+             create_cache_tvf_edac(initial_condition, transport_velocity)...,
              create_cache_surface_normal(surface_normal_method, ELTYPE, NDIMS,
                                          n_particles)...,
              create_cache_surface_tension(surface_tension, ELTYPE, NDIMS,
@@ -208,16 +208,6 @@ function Base.show(io::IO, ::MIME"text/plain", system::EntropicallyDampedSPHSyst
     end
 end
 
-create_cache_tvf(::Val{:edac}, initial_condition, ::Nothing) = (;)
-
-function create_cache_tvf(::Val{:edac}, initial_condition, ::TransportVelocityAdami)
-    pressure_average = copy(initial_condition.pressure)
-    neighbor_counter = Vector{Int}(undef, nparticles(initial_condition))
-    update_callback_used = Ref(false)
-
-    return (; pressure_average, neighbor_counter, update_callback_used)
-end
-
 @inline function Base.eltype(::EntropicallyDampedSPHSystem{<:Any, ELTYPE}) where {ELTYPE}
     return ELTYPE
 end
@@ -226,12 +216,12 @@ end
     return v_nvariables(system, system.density_calculator)
 end
 
-@inline function v_nvariables(system::EntropicallyDampedSPHSystem, density_calculator)
-    return ndims(system) * factor_tvf(system) + 1
+@inline function v_nvariables(system::EntropicallyDampedSPHSystem, ::SummationDensity)
+    return ndims(system) + 1
 end
 
 @inline function v_nvariables(system::EntropicallyDampedSPHSystem, ::ContinuityDensity)
-    return ndims(system) * factor_tvf(system) + 2
+    return ndims(system) + 2
 end
 
 system_correction(system::EntropicallyDampedSPHSystem) = system.correction
@@ -260,10 +250,12 @@ end
 
 @inline system_sound_speed(system::EntropicallyDampedSPHSystem) = system.sound_speed
 
+@inline transport_velocity(system::EntropicallyDampedSPHSystem) = system.transport_velocity
+
 @inline average_pressure(system, particle) = zero(eltype(system))
 
 @inline function average_pressure(system::EntropicallyDampedSPHSystem, particle)
-    average_pressure(system, system.transport_velocity, particle)
+    average_pressure(system, transport_velocity(system), particle)
 end
 
 @inline function average_pressure(system, ::TransportVelocityAdami, particle)
@@ -303,18 +295,14 @@ function update_pressure!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_od
     compute_surface_delta_function!(system, system.surface_tension, semi)
 end
 
-function update_final!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, semi, t;
-                       update_from_callback=false)
+function update_final!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, semi, t)
     (; surface_tension) = system
 
     # Surface normal of neighbor and boundary needs to have been calculated already
     compute_curvature!(system, surface_tension, v, u, v_ode, u_ode, semi, t)
     compute_stress_tensors!(system, surface_tension, v, u, v_ode, u_ode, semi, t)
-    update_average_pressure!(system, system.transport_velocity, v_ode, u_ode, semi)
-
-    # Check that TVF is only used together with `UpdateCallback`
-    check_tvf_configuration(system, system.transport_velocity, v, u, v_ode, u_ode, semi, t;
-                            update_from_callback)
+    update_average_pressure!(system, transport_velocity(system), v_ode, u_ode, semi)
+    update_tvf!(system, transport_velocity(system), v, u, v_ode, u_ode, semi, t)
 end
 
 function update_average_pressure!(system, ::Nothing, v_ode, u_ode, semi)
@@ -347,8 +335,7 @@ function update_average_pressure!(system, ::TransportVelocityAdami, v_ode, u_ode
                                points=each_moving_particle(system)) do particle, neighbor,
                                                                        pos_diff, distance
             pressure_average[particle] += current_pressure(v_neighbor_system,
-                                                           neighbor_system,
-                                                           neighbor)
+                                                           neighbor_system, neighbor)
             neighbor_counter[particle] += 1
         end
     end
diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl
index 308217170..d3e6adf97 100644
--- a/src/schemes/fluid/fluid.jl
+++ b/src/schemes/fluid/fluid.jl
@@ -85,7 +85,6 @@ function write_v0!(v0, system::FluidSystem)
     copyto!(v0, indices, system.initial_condition.velocity, indices)
 
     write_v0!(v0, system, system.density_calculator)
-    write_v0!(v0, system, system.transport_velocity)
 
     return v0
 end
@@ -199,25 +198,3 @@ include("surface_tension.jl")
 include("surface_normal_sph.jl")
 include("weakly_compressible_sph/weakly_compressible_sph.jl")
 include("entropically_damped_sph/entropically_damped_sph.jl")
-
-@inline function add_velocity!(du, v, particle,
-                               system::Union{EntropicallyDampedSPHSystem,
-                                             WeaklyCompressibleSPHSystem})
-    add_velocity!(du, v, particle, system, system.transport_velocity)
-end
-
-@inline function momentum_convection(system, neighbor_system,
-                                     v_particle_system, v_neighbor_system, rho_a, rho_b,
-                                     m_a, m_b, particle, neighbor, grad_kernel)
-    return zero(grad_kernel)
-end
-
-@inline function momentum_convection(system,
-                                     neighbor_system::Union{EntropicallyDampedSPHSystem,
-                                                            WeaklyCompressibleSPHSystem},
-                                     v_particle_system, v_neighbor_system, rho_a, rho_b,
-                                     m_a, m_b, particle, neighbor, grad_kernel)
-    momentum_convection(system, neighbor_system, system.transport_velocity,
-                        v_particle_system, v_neighbor_system, rho_a, rho_b,
-                        m_a, m_b, particle, neighbor, grad_kernel)
-end
diff --git a/src/schemes/fluid/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl
index a69af10e4..bbf40bde2 100644
--- a/src/schemes/fluid/transport_velocity.jl
+++ b/src/schemes/fluid/transport_velocity.jl
@@ -1,179 +1,168 @@
 """
     TransportVelocityAdami(background_pressure::Real)
 
-Transport Velocity Formulation (TVF) to suppress pairing and tensile instability.
+Transport Velocity Formulation (TVF) by [Adami et al. (2013)](@cite Adami2013)
+to suppress pairing and tensile instability.
 See [TVF](@ref transport_velocity_formulation) for more details of the method.
 
 # Arguments
 - `background_pressure`: Background pressure. Suggested is a background pressure which is
                          on the order of the reference pressure.
-
-!!! note
-    There is no need for an artificial viscosity to suppress tensile instability when using `TransportVelocityAdami`.
-    Thus, it is highly recommended to use [`ViscosityAdami`](@ref) as viscosity model,
-    since [`ArtificialViscosityMonaghan`](@ref) leads to bad results.
 """
 struct TransportVelocityAdami{T <: Real}
     background_pressure::T
 end
 
-# Calculate `v_nvariables` appropriately
-@inline factor_tvf(system::FluidSystem) = factor_tvf(system, system.transport_velocity)
-@inline factor_tvf(system, ::Nothing) = 1
-@inline factor_tvf(system, ::TransportVelocityAdami) = 2
-
-@inline update_transport_velocity!(system, v_ode, semi) = system
+# No TVF for a system by default
+@inline transport_velocity(system) = nothing
 
-@inline function update_transport_velocity!(system::FluidSystem, v_ode, semi)
-    update_transport_velocity!(system, v_ode, semi, system.transport_velocity)
-end
+create_cache_tvf_wcsph(initial_condition, ::Nothing) = (;)
 
-@inline update_transport_velocity!(system, v_ode, semi, ::Nothing) = system
-
-@inline function update_transport_velocity!(system, v_ode, semi, ::TransportVelocityAdami)
-    v = wrap_v(v_ode, system, semi)
-    @threaded semi for particle in each_moving_particle(system)
-        for i in 1:ndims(system)
-            v[ndims(system) + i, particle] = v[i, particle]
-        end
-    end
+function create_cache_tvf_wcsph(initial_condition, ::TransportVelocityAdami)
+    delta_v = zeros(eltype(initial_condition), ndims(initial_condition),
+                    nparticles(initial_condition))
 
-    return system
+    return (; delta_v)
 end
 
-function write_v0!(v0, system::FluidSystem, ::TransportVelocityAdami)
-    (; initial_condition) = system
+create_cache_tvf_edac(initial_condition, ::Nothing) = (;)
 
-    # Write particle velocities
-    for dim in 1:ndims(system)
-        v0[ndims(system) + dim, :] = initial_condition.velocity[dim, :]
-    end
+function create_cache_tvf_edac(initial_condition, ::TransportVelocityAdami)
+    delta_v = zeros(eltype(initial_condition), ndims(initial_condition),
+                    nparticles(initial_condition))
+    pressure_average = copy(initial_condition.pressure)
+    neighbor_counter = Vector{Int}(undef, nparticles(initial_condition))
 
-    return v0
+    return (; delta_v, pressure_average, neighbor_counter)
 end
 
-# Add momentum velocity.
-@inline function add_velocity!(du, v, particle, system, ::Nothing)
-    for i in 1:ndims(system)
-        du[i, particle] = v[i, particle]
-    end
-
-    return du
+# `δv` is the correction to the particle velocity due to the TVF.
+# Particles are advected with the velocity `v + δv`.
+@propagate_inbounds function delta_v(system, particle)
+    return delta_v(system, transport_velocity(system), particle)
 end
 
-# Add advection velocity.
-@inline function add_velocity!(du, v, particle, system, ::TransportVelocityAdami)
-    for i in 1:ndims(system)
-        du[i, particle] = v[ndims(system) + i, particle]
-    end
-
-    return du
+@propagate_inbounds function delta_v(system, ::TransportVelocityAdami, particle)
+    return extract_svector(system.cache.delta_v, system, particle)
 end
 
-@inline function advection_velocity(v, system, particle)
-    return SVector(ntuple(@inline(dim->v[ndims(system) + dim, particle]), ndims(system)))
+# Zero when no TVF is used
+@inline function delta_v(system, transport_velocity, particle)
+    return zero(SVector{ndims(system), eltype(system)})
 end
 
-@inline function momentum_convection(system, neighbor_system, ::Nothing,
-                                     v_particle_system, v_neighbor_system, rho_a, rho_b,
-                                     m_a, m_b, particle, neighbor, grad_kernel)
+@inline function dv_transport_velocity(::Nothing, system, neighbor_system,
+                                       particle, neighbor, v_system, v_neighbor_system,
+                                       m_a, m_b, rho_a, rho_b, pos_diff, distance,
+                                       grad_kernel, correction)
     return zero(grad_kernel)
 end
 
-@inline function momentum_convection(system, neighbor_system, ::TransportVelocityAdami,
-                                     v_particle_system, v_neighbor_system, rho_a, rho_b,
-                                     m_a, m_b, particle, neighbor, grad_kernel)
-    volume_a = m_a / rho_a
-    volume_b = m_b / rho_b
-    volume_term = (volume_a^2 + volume_b^2) / m_a
-
-    momentum_velocity_a = current_velocity(v_particle_system, system, particle)
-    advection_velocity_a = advection_velocity(v_particle_system, system, particle)
-
-    momentum_velocity_b = current_velocity(v_neighbor_system, neighbor_system, neighbor)
-    advection_velocity_b = advection_velocity(v_neighbor_system, neighbor_system, neighbor)
-
-    A_a = rho_a * momentum_velocity_a * (advection_velocity_a - momentum_velocity_a)'
-    A_b = rho_b * momentum_velocity_b * (advection_velocity_b - momentum_velocity_b)'
-
-    return volume_term * ((A_a + A_b) / 2) * grad_kernel
-end
-
-@inline transport_velocity!(dv, system, neighbor_system, particle, neighbor, m_a, m_b,
-                            grad_kernel) = dv
-
-@inline function transport_velocity!(dv, system::FluidSystem, neighbor_system,
-                                     particle, neighbor, m_a, m_b, grad_kernel)
-    transport_velocity!(dv, system, neighbor_system, system.transport_velocity,
-                        particle, neighbor, m_a, m_b, grad_kernel)
-end
-
-@inline function transport_velocity!(dv, system, neighbor_system, ::Nothing,
-                                     particle, neighbor, m_a, m_b, grad_kernel)
-    return dv
-end
-
-@inline function transport_velocity!(dv, system, neighbor_system, ::TransportVelocityAdami,
-                                     particle, neighbor, m_a, m_b, grad_kernel)
-    (; transport_velocity) = system
-    (; background_pressure) = transport_velocity
-
-    NDIMS = ndims(system)
+@inline function dv_transport_velocity(::TransportVelocityAdami, system, neighbor_system,
+                                       particle, neighbor, v_system, v_neighbor_system,
+                                       m_a, m_b, rho_a, rho_b, pos_diff, distance,
+                                       grad_kernel, correction)
+    v_a = current_velocity(v_system, system, particle)
+    delta_v_a = delta_v(system, particle)
 
-    # The TVF is based on the assumption that the pressure gradient is only accurately
-    # computed when the particle distribution is isotropic.
-    # That means, the force contribution vanishes only if the particle distribution is
-    # isotropic AND the field being differentiated by the kernel gradient is spatially constant.
-    # So we must guarantee a constant field and therefore the reference density is used
-    # instead of the locally computed one.
-    # TODO:
-    # volume_a = particle_spacing(system, particle)^ndims(system)
-    # volume_b = particle_spacing(neighbor_system, neighbor)^ndims(neighbor_system)
-    volume_a = m_a / system.initial_condition.density[particle]
-    volume_b = m_b / neighbor_system.initial_condition.density[neighbor]
+    v_b = current_velocity(v_neighbor_system, neighbor_system, neighbor)
+    delta_v_b = delta_v(neighbor_system, neighbor)
 
-    volume_term = (volume_a^2 + volume_b^2) / m_a
+    A_a = rho_a * v_a * delta_v_a'
+    A_b = rho_b * v_b * delta_v_b'
 
-    for dim in 1:NDIMS
-        dv[NDIMS + dim, particle] -= volume_term * background_pressure * grad_kernel[dim]
-    end
-
-    return dv
+    # The following term depends on the pressure acceleration formulation.
+    # See the large comment below. In the original paper (Adami et al., 2013), this is
+    #   (V_a^2 + V_b^2) / m_a * ((A_a + A_b) / 2) * ∇W_ab.
+    # With the most common pressure acceleration formulation, this is
+    #   m_b * (A_a + A_b) / (ρ_a * ρ_b) * ∇W_ab.
+    # In order to obtain this, we pass `p_a = A_a` and `p_b = A_b` to the
+    # `pressure_acceleration` function.
+    return pressure_acceleration(system, neighbor_system, particle, neighbor,
+                                 m_a, m_b, A_a, A_b, rho_a, rho_b, pos_diff,
+                                 distance, grad_kernel, correction)
 end
 
-function reset_callback_flag!(system::FluidSystem)
-    reset_callback_flag!(system, system.transport_velocity)
-end
-
-reset_callback_flag!(system, ::Nothing) = system
-
-function reset_callback_flag!(system::FluidSystem, ::TransportVelocityAdami)
-    system.cache.update_callback_used[] = false
-
+function update_tvf!(system, transport_velocity, v, u, v_ode, u_ode, semi, t)
     return system
 end
 
-function update_callback_used!(system::FluidSystem)
-    update_callback_used!(system, system.transport_velocity)
-end
-
-update_callback_used!(system, ::Nothing) = system
-
-function update_callback_used!(system, transport_velocity)
-    system.cache.update_callback_used[] = true
-end
-
-function check_tvf_configuration(system::FluidSystem, ::Nothing,
-                                 v, u, v_ode, u_ode, semi, t; update_from_callback=false)
-    return system
-end
+function update_tvf!(system, transport_velocity::TransportVelocityAdami, v, u, v_ode,
+                     u_ode, semi, t)
+    (; cache, correction) = system
+    (; delta_v) = cache
+    (; background_pressure) = transport_velocity
 
-function check_tvf_configuration(system::FluidSystem, ::TransportVelocityAdami,
-                                 v, u, v_ode, u_ode, semi, t; update_from_callback=false)
-    # The `UpdateCallback` will set `system.cache.update_callback_used[]` to `true`
-    # in the initialization. However, other callbacks might update the system first, hence `update_from_callback`.
-    if !update_from_callback && !(system.cache.update_callback_used[])
-        throw(ArgumentError("`UpdateCallback` is required when using `TransportVelocityAdami`"))
+    sound_speed = system_sound_speed(system)
+
+    set_zero!(delta_v)
+
+    foreach_system(semi) do neighbor_system
+        v_neighbor = wrap_v(v_ode, neighbor_system, semi)
+        u_neighbor = wrap_u(u_ode, neighbor_system, semi)
+
+        system_coords = current_coordinates(u, system)
+        neighbor_coords = current_coordinates(u_neighbor, neighbor_system)
+
+        foreach_point_neighbor(system, neighbor_system, system_coords, neighbor_coords,
+                               semi;
+                               points=each_moving_particle(system)) do particle, neighbor,
+                                                                       pos_diff, distance
+            m_a = @inbounds hydrodynamic_mass(system, particle)
+            m_b = @inbounds hydrodynamic_mass(neighbor_system, neighbor)
+
+            rho_a = @inbounds current_density(v, system, particle)
+            rho_b = @inbounds current_density(v_neighbor, neighbor_system, neighbor)
+
+            h = smoothing_length(system, particle)
+
+            grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle)
+
+            # In the original paper (Adami et al., 2013), the transport velocity is applied
+            # as follows:
+            #   v_{1/2} = v_0 + Δt/2 * a,
+            # where a is the regular SPH acceleration term (pressure, viscosity, etc.).
+            #   r_1 = r_0 + Δt * (v_{1/2},
+            # where ̃v_{1/2} = v_{1/2} + Δt/2 * p_0 / m_a * \sum_b[ (V_a^2 + V_b^2) * ∇W_ab ]
+            # is the transport velocity.
+            # We call δv_{1/2} = ̃v_{1/2} - v_{1/2} the shifting velocity.
+            # We will call δv_{1/2} the shifting velocity, which is given by
+            #   δv = -Δt/2 * p_0 / m_a * \sum_b[ (V_a^2 + V_b^2) * ∇W_ab ],
+            # where p_0 is the background pressure, V_a = m_a / ρ_a, V_b = m_b / ρ_b.
+            # This term depends on the pressure acceleration formulation.
+            # In Zhang et al. (2017), the pressure acceleration term
+            #   m_b * (p_a / ρ_a^2 + p_b / ρ_b^2) * ∇W_ab
+            # is used. They consequently changed the shifting velocity to
+            #   δv = -Δt/2 * p_0 * \sum_b[ m_b * (1 / ρ_a^2 + 1 / ρ_b^2) * ∇W_ab ].
+            # We therefore use the function `pressure_acceleration` to compute the
+            # shifting velocity according to the used pressure acceleration formulation.
+            # In most cases, this will be
+            #   δv = -Δt/2 * p_0 * \sum_b[ m_b * (1 + 1) / (ρ_a * ρ_b) * ∇W_ab ].
+            #
+            # In these papers, the shifting velocity is scaled by the time step Δt.
+            # We generally want the spatial discretization to be independent of the time step.
+            # Scaling the shifting velocity by the time step would lead to less shifting
+            # when very small time steps are used for testing/debugging purposes.
+            # This is especially problematic in TrixiParticles.jl, as the time step can vary
+            # significantly between different time integration methods (low vs high order).
+            # In order to eliminate the time step from the shifting velocity, we apply the
+            # CFL condition used in Adami et al. (2013):
+            #   Δt <= 0.25 * h / c,
+            # where h is the smoothing length and c is the sound speed.
+            # Applying this equation as equality yields the shifting velocity
+            #   δv = -p_0 / 8 * h / c * \sum_b[ m_b * (1 + 1) / (ρ_a * ρ_b) * ∇W_ab ].
+            # The last part is achieved by passing `p_a = 1` and `p_b = 1` to the
+            # `pressure_acceleration` function.
+            delta_v_ = background_pressure / 8 * h / sound_speed *
+                       pressure_acceleration(system, neighbor_system, particle, neighbor,
+                                             m_a, m_b, 1, 1, rho_a, rho_b, pos_diff,
+                                             distance, grad_kernel, correction)
+
+            # Write into the buffer
+            for i in eachindex(delta_v_)
+                @inbounds delta_v[i, particle] += delta_v_[i]
+            end
+        end
     end
 
     return system
diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl
index a1987a8c0..af088e0a1 100644
--- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl
+++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl
@@ -5,8 +5,9 @@
 function interact!(dv, v_particle_system, u_particle_system,
                    v_neighbor_system, u_neighbor_system,
                    particle_system::WeaklyCompressibleSPHSystem, neighbor_system, semi)
-    (; density_calculator, state_equation, correction) = particle_system
-    (; sound_speed) = state_equation
+    (; density_calculator, correction) = particle_system
+
+    sound_speed = system_sound_speed(particle_system)
 
     surface_tension_a = surface_tension_model(particle_system)
     surface_tension_b = surface_tension_model(neighbor_system)
@@ -70,10 +71,12 @@ function interact!(dv, v_particle_system, u_particle_system,
                                                grad_kernel)
 
         # Add convection term (only when using `TransportVelocityAdami`)
-        dv_convection = momentum_convection(particle_system, neighbor_system,
-                                            v_particle_system, v_neighbor_system,
-                                            rho_a, rho_b, m_a, m_b,
-                                            particle, neighbor, grad_kernel)
+        dv_tvf = dv_transport_velocity(transport_velocity(particle_system),
+                                       particle_system, neighbor_system,
+                                       particle, neighbor,
+                                       v_particle_system, v_neighbor_system,
+                                       m_a, m_b, rho_a, rho_b, pos_diff, distance,
+                                       grad_kernel, correction)
 
         dv_surface_tension = surface_tension_correction *
                              surface_tension_force(surface_tension_a, surface_tension_b,
@@ -86,17 +89,12 @@ function interact!(dv, v_particle_system, u_particle_system,
 
         for i in 1:ndims(particle_system)
             @inbounds dv[i,
-                         particle] += dv_pressure[i] + dv_viscosity_[i] +
-                                      dv_convection[i] + dv_surface_tension[i] +
-                                      dv_adhesion[i]
+                         particle] += dv_pressure[i] + dv_viscosity_[i] + dv_tvf[i] +
+                                      dv_surface_tension[i] + dv_adhesion[i]
             # Debug example
             # debug_array[i, particle] += dv_pressure[i]
         end
 
-        # Apply the transport velocity (only when using a transport velocity)
-        transport_velocity!(dv, particle_system, neighbor_system, particle, neighbor,
-                            m_a, m_b, grad_kernel)
-
         # TODO If variable smoothing_length is used, this should use the neighbor smoothing length
         # Propagate `@inbounds` to the continuity equation, which accesses particle data
         @inbounds continuity_equation!(dv, density_calculator, v_particle_system,
diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl
index 690803960..bc9d541ea 100644
--- a/src/schemes/fluid/weakly_compressible_sph/system.jl
+++ b/src/schemes/fluid/weakly_compressible_sph/system.jl
@@ -147,7 +147,7 @@ function WeaklyCompressibleSPHSystem(initial_condition,
                                           n_particles)...,
              create_cache_refinement(initial_condition, particle_refinement,
                                      smoothing_length)...,
-             create_cache_tvf(Val(:wcsph), initial_condition, transport_velocity)...,
+             create_cache_tvf_wcsph(initial_condition, transport_velocity)...,
              color=Int(color_value))
 
     # If the `reference_density_spacing` is set calculate the `ideal_neighbor_count`
@@ -221,14 +221,6 @@ function Base.show(io::IO, ::MIME"text/plain", system::WeaklyCompressibleSPHSyst
     end
 end
 
-create_cache_tvf(::Val{:wcsph}, initial_condition, ::Nothing) = (;)
-
-function create_cache_tvf(::Val{:wcsph}, initial_condition, ::TransportVelocityAdami)
-    update_callback_used = Ref(false)
-
-    return (; update_callback_used)
-end
-
 @inline function Base.eltype(::WeaklyCompressibleSPHSystem{<:Any, ELTYPE}) where {ELTYPE}
     return ELTYPE
 end
@@ -237,12 +229,12 @@ end
     return v_nvariables(system, system.density_calculator)
 end
 
-@inline function v_nvariables(system::WeaklyCompressibleSPHSystem, density_calculator)
-    return ndims(system) * factor_tvf(system)
+@inline function v_nvariables(system::WeaklyCompressibleSPHSystem, ::SummationDensity)
+    return ndims(system)
 end
 
 @inline function v_nvariables(system::WeaklyCompressibleSPHSystem, ::ContinuityDensity)
-    return ndims(system) * factor_tvf(system) + 1
+    return ndims(system) + 1
 end
 
 system_correction(system::WeaklyCompressibleSPHSystem) = system.correction
@@ -286,6 +278,8 @@ end
 
 @inline system_sound_speed(system::WeaklyCompressibleSPHSystem) = system.state_equation.sound_speed
 
+@inline transport_velocity(system::WeaklyCompressibleSPHSystem) = system.transport_velocity
+
 function update_quantities!(system::WeaklyCompressibleSPHSystem, v, u,
                             v_ode, u_ode, semi, t)
     (; density_calculator, density_diffusion, correction) = system
@@ -301,30 +295,28 @@ end
 function update_pressure!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t)
     (; density_calculator, correction, surface_normal_method, surface_tension) = system
 
-    compute_correction_values!(system, correction, u, v_ode, u_ode, semi)
+    compute_pressure!(system, v, semi)
 
+    # These are only computed when using corrections
+    compute_correction_values!(system, correction, u, v_ode, u_ode, semi)
     compute_gradient_correction_matrix!(correction, system, u, v_ode, u_ode, semi)
-
     # `kernel_correct_density!` only performed for `SummationDensity`
     kernel_correct_density!(system, v, u, v_ode, u_ode, semi, correction,
                             density_calculator)
-    compute_pressure!(system, v, semi)
+
+    # These are only computed when using surface tension
     compute_surface_normal!(system, surface_normal_method, v, u, v_ode, u_ode, semi, t)
     compute_surface_delta_function!(system, surface_tension, semi)
     return system
 end
 
-function update_final!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t;
-                       update_from_callback=false)
+function update_final!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t)
     (; surface_tension) = system
 
     # Surface normal of neighbor and boundary needs to have been calculated already
     compute_curvature!(system, surface_tension, v, u, v_ode, u_ode, semi, t)
     compute_stress_tensors!(system, surface_tension, v, u, v_ode, u_ode, semi, t)
-
-    # Check that TVF is only used together with `UpdateCallback`
-    check_tvf_configuration(system, system.transport_velocity, v, u, v_ode, u_ode, semi, t;
-                            update_from_callback)
+    update_tvf!(system, transport_velocity(system), v, u, v_ode, u_ode, semi, t)
 end
 
 function kernel_correct_density!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode,
diff --git a/src/schemes/solid/total_lagrangian_sph/system.jl b/src/schemes/solid/total_lagrangian_sph/system.jl
index acf269d7b..b465ca2d7 100644
--- a/src/schemes/solid/total_lagrangian_sph/system.jl
+++ b/src/schemes/solid/total_lagrangian_sph/system.jl
@@ -261,8 +261,8 @@ function update_quantities!(system::TotalLagrangianSPHSystem, v, u, v_ode, u_ode
     return system
 end
 
-function update_final!(system::TotalLagrangianSPHSystem, v, u, v_ode, u_ode, semi, t;
-                       update_from_callback=false)
+function update_boundary_interpolation!(system::TotalLagrangianSPHSystem, v, u,
+                                        v_ode, u_ode, semi, t)
     (; boundary_model) = system
 
     # Only update boundary model
diff --git a/test/schemes/fluid/fluid.jl b/test/schemes/fluid/fluid.jl
index 570c3bc86..664325a1a 100644
--- a/test/schemes/fluid/fluid.jl
+++ b/test/schemes/fluid/fluid.jl
@@ -2,6 +2,5 @@ include("weakly_compressible_sph/weakly_compressible_sph.jl")
 include("rhs.jl")
 include("pressure_acceleration.jl")
 include("surface_normal_sph.jl")
-include("transport_velocity.jl")
 include("surface_tension.jl")
 include("viscosity.jl")
diff --git a/test/schemes/fluid/transport_velocity.jl b/test/schemes/fluid/transport_velocity.jl
deleted file mode 100644
index bfa71e125..000000000
--- a/test/schemes/fluid/transport_velocity.jl
+++ /dev/null
@@ -1,35 +0,0 @@
-@trixi_testset "Transport Velocity Formulation" begin
-    particle_spacing = 0.1
-    smoothing_kernel = SchoenbergCubicSplineKernel{2}()
-    smoothing_length = 1.2particle_spacing
-
-    fluid = rectangular_patch(particle_spacing, (3, 3))
-
-    v0_tvf = zeros(5, nparticles(fluid))
-
-    system_tvf = EntropicallyDampedSPHSystem(fluid, smoothing_kernel,
-                                             transport_velocity=TransportVelocityAdami(0.0),
-                                             smoothing_length, 0.0)
-    system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length, 0.0)
-
-    @testset "Number of Variables" begin
-        @test TrixiParticles.v_nvariables(system_tvf) == 5
-        @test TrixiParticles.v_nvariables(system) == 3
-    end
-
-    @testset "write_v0!" begin
-        TrixiParticles.write_v0!(v0_tvf, system_tvf)
-
-        @test isapprox(vcat(fluid.velocity, fluid.velocity, fluid.pressure'), v0_tvf)
-    end
-
-    @testset "Update" begin
-        semi = Semidiscretization(system_tvf)
-        fill!(v0_tvf, 1.5)
-        v0_tvf[1:2, :] .= 2.5
-
-        TrixiParticles.update_transport_velocity!(system_tvf, vec(v0_tvf), semi)
-
-        @test isapprox(fill(2.5, (4, nparticles(system_tvf))), v0_tvf[1:4, :])
-    end
-end
diff --git a/test/systems/edac_system.jl b/test/systems/edac_system.jl
index 424680abf..95549ed81 100644
--- a/test/systems/edac_system.jl
+++ b/test/systems/edac_system.jl
@@ -228,7 +228,7 @@
         TrixiParticles.initialize_neighborhood_searches!(semi)
 
         u_ode = vec(fluid.coordinates)
-        v_ode = vec(vcat(fluid.velocity, fluid.velocity, fluid.pressure'))
+        v_ode = vec(vcat(fluid.velocity, fluid.pressure'))
 
         TrixiParticles.update_average_pressure!(system, system.transport_velocity, v_ode,
                                                 u_ode, semi)