From 3a63138e46230c1f76008db9f13c2c7d54641cf5 Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Tue, 17 Feb 2026 12:09:49 +0100
Subject: [PATCH 1/7] Implement velocity averaging for TLSPH

---
 src/callbacks/split_integration.jl            | 53 +++++++++++++
 src/callbacks/update.jl                       | 11 ++-
 src/general/semidiscretization.jl             |  2 +-
 src/preprocessing/particle_packing/system.jl  |  2 +-
 src/schemes/boundary/open_boundary/system.jl  |  2 +-
 .../boundary/wall_boundary/dummy_particles.jl |  6 +-
 src/schemes/fluid/shifting_techniques.jl      |  4 +-
 .../structure/total_lagrangian_sph/system.jl  | 40 ++++++++--
 .../total_lagrangian_sph.jl                   |  1 +
 .../velocity_averaging.jl                     | 79 +++++++++++++++++++
 10 files changed, 187 insertions(+), 13 deletions(-)
 create mode 100644 src/schemes/structure/total_lagrangian_sph/velocity_averaging.jl

diff --git a/src/callbacks/split_integration.jl b/src/callbacks/split_integration.jl
index 2bb9dfc9e2..60ef940eeb 100644
--- a/src/callbacks/split_integration.jl
+++ b/src/callbacks/split_integration.jl
@@ -122,6 +122,15 @@ function initialize_split_integration!(cb, u, t, integrator)
     ode_split = DynamicalODEProblem(kick_split!, drift_split!, v0_ode_split, u0_ode_split,
                                     tspan, p)
 
+    # Add lightweight callback to (potentially) update the averaged velocity
+    # during the split integration.
+    callback = UpdateAveragedVelocityCallback()
+    if haskey(kwargs, :callback)
+        kwargs[:callback] = CallbackSet(kwargs[:callback], callback)
+    else
+        kwargs[:callback] = callback
+    end
+
     # Create the split integrator.
     # We need the timer here to keep the output clean because this will call `kick!` once.
     @trixi_timeit timer() "split integration" begin
@@ -356,3 +365,47 @@ function Base.show(io::IO, ::MIME"text/plain",
         summary_box(io, "SplitIntegrationCallback", setup)
     end
 end
+
+# === Non-public callback for updating the averaged velocity ===
+# When no split integration is used, this is done from the `UpdateCallback`.
+# With split integration, we use this lightweight callback to avoid updating the systems.
+function UpdateAveragedVelocityCallback()
+    # The first one is the `condition`, the second the `affect!`
+    return DiscreteCallback(update_averaged_velocity_callback!,
+                            update_averaged_velocity_callback!,
+                            initialize=(initialize_averaged_velocity_callback!),
+                            save_positions=(false, false))
+end
+
+# `initialize`
+function initialize_averaged_velocity_callback!(cb, vu_ode, t, integrator)
+    v_ode, u_ode = vu_ode.x
+    semi = integrator.p
+
+    foreach_system(semi) do system
+        initialize_averaged_velocity!(system, v_ode, semi, t)
+    end
+
+    return cb
+end
+
+# `condition`
+function update_averaged_velocity_callback!(u, t, integrator)
+    return true
+end
+
+# `affect!`
+function update_averaged_velocity_callback!(integrator)
+    t_new = integrator.t
+    semi = integrator.p
+    v_ode, u_ode = integrator.u.x
+
+    foreach_system(semi) do system
+        compute_averaged_velocity!(system, v_ode, semi, t_new)
+    end
+
+    # Tell OrdinaryDiffEq that `integrator.u` has not been modified
+    u_modified!(integrator, false)
+
+    return integrator
+end
diff --git a/src/callbacks/update.jl b/src/callbacks/update.jl
index ab76dfacc2..b5e8b66cdb 100644
--- a/src/callbacks/update.jl
+++ b/src/callbacks/update.jl
@@ -52,12 +52,17 @@ function initial_update!(cb, u, t, integrator)
     initial_update!(cb.affect!, u, t, integrator)
 end
 
-function initial_update!(cb::UpdateCallback, u, t, integrator)
+function initial_update!(cb::UpdateCallback, vu_ode, t, integrator)
+    v_ode, u_ode = vu_ode.x
     semi = integrator.p
 
     # Tell the semidiscretization that the `UpdateCallback` is used
     semi.update_callback_used[] = true
 
+    foreach_system(semi) do system
+        initialize_averaged_velocity!(system, v_ode, semi, t)
+    end
+
     return cb(integrator)
 end
 
@@ -105,6 +110,10 @@ function (update_callback!::UpdateCallback)(integrator)
             particle_shifting_from_callback!(u_ode, shifting_technique(system), system,
                                              v_ode, semi, integrator)
         end
+
+        foreach_system(semi) do system
+            compute_averaged_velocity!(system, v_ode, semi, t)
+        end
     end
 
     return integrator
diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl
index eff642b2cf..5498fb9b72 100644
--- a/src/general/semidiscretization.jl
+++ b/src/general/semidiscretization.jl
@@ -694,7 +694,7 @@ 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) && !semi.update_callback_used[]
+        if requires_update_callback(system, semi) && !semi.update_callback_used[]
             system_name = system |> typeof |> nameof
             throw(ArgumentError("`UpdateCallback` is required for `$system_name`"))
         end
diff --git a/src/preprocessing/particle_packing/system.jl b/src/preprocessing/particle_packing/system.jl
index d6439741ac..65c5160cb2 100644
--- a/src/preprocessing/particle_packing/system.jl
+++ b/src/preprocessing/particle_packing/system.jl
@@ -214,7 +214,7 @@ end
     return ndims(system)
 end
 
-@inline requires_update_callback(system::ParticlePackingSystem) = true
+@inline requires_update_callback(system::ParticlePackingSystem, semi) = true
 
 function write2vtk!(vtk, v, u, t, system::ParticlePackingSystem)
     vtk["velocity"] = [advection_velocity(v, system, particle)
diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl
index d74a69b34b..087f258565 100644
--- a/src/schemes/boundary/open_boundary/system.jl
+++ b/src/schemes/boundary/open_boundary/system.jl
@@ -244,7 +244,7 @@ end
 @inline buffer(system::OpenBoundarySystem) = system.buffer
 
 # The `UpdateCallback` is required to update particle positions between time steps
-@inline requires_update_callback(system::OpenBoundarySystem) = true
+@inline requires_update_callback(system::OpenBoundarySystem, semi) = true
 
 function smoothing_length(system::OpenBoundarySystem, particle)
     return system.smoothing_length
diff --git a/src/schemes/boundary/wall_boundary/dummy_particles.jl b/src/schemes/boundary/wall_boundary/dummy_particles.jl
index 9f5d8926a1..c88658eb28 100644
--- a/src/schemes/boundary/wall_boundary/dummy_particles.jl
+++ b/src/schemes/boundary/wall_boundary/dummy_particles.jl
@@ -681,8 +681,10 @@ end
     (; volume, wall_velocity) = cache
 
     # Prescribed velocity of the boundary particle.
-    # This velocity is zero when not using moving boundaries.
-    v_boundary = current_velocity(v, system, particle)
+    # This velocity is zero when not using moving boundaries or TLSPH.
+    # If not using TLSPH with velocity averaging, this function simply
+    # forwards to `current_velocity`.
+    v_boundary = velocity_for_viscosity(v, system, particle)
 
     for dim in eachindex(v_boundary)
         # The second term is the precalculated smoothed velocity field of the fluid
diff --git a/src/schemes/fluid/shifting_techniques.jl b/src/schemes/fluid/shifting_techniques.jl
index ceaff12a21..3aecb4ea45 100644
--- a/src/schemes/fluid/shifting_techniques.jl
+++ b/src/schemes/fluid/shifting_techniques.jl
@@ -5,7 +5,9 @@ abstract type AbstractShiftingTechnique end
 
 # WARNING: Be careful if defining this function for a specific system type.
 # The version for a specific system type will override this generic version.
-requires_update_callback(system) = requires_update_callback(shifting_technique(system))
+function requires_update_callback(system, semi)
+    return requires_update_callback(shifting_technique(system))
+end
 requires_update_callback(::Nothing) = false
 requires_update_callback(::AbstractShiftingTechnique) = false
 
diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl
index 2d2792cd21..f8e13560c7 100644
--- a/src/schemes/structure/total_lagrangian_sph/system.jl
+++ b/src/schemes/structure/total_lagrangian_sph/system.jl
@@ -7,7 +7,8 @@
                              acceleration=ntuple(_ -> 0.0, NDIMS),
                              penalty_force=nothing, viscosity=nothing,
                              source_terms=nothing, boundary_model=nothing,
-                             self_interaction_nhs=:default)
+                             self_interaction_nhs=:default,
+                             velocity_averaging=nothing)
 
 System for particles of an elastic structure.
 
@@ -61,6 +62,9 @@ See [Total Lagrangian SPH](@ref tlsph) for more details on the method.
                     `transpose_backend` of the [`PrecomputedNeighborhoodSearch`](@ref)
                     is set to `true` when running on GPUs and `false` on CPUs.
                     Alternatively, a user-defined neighborhood search can be passed here.
+- `velocity_averaging`: Velocity averaging technique to be applied on the velocity field
+                    to obtain an averaged velocity to be used in fluid-structure interaction.
+                    See [`VelocityAveraging`](@ref) for details.
 
 !!! note
     If specifying the clamped particles manually (via `n_clamped_particles`),
@@ -83,7 +87,7 @@ See [Total Lagrangian SPH](@ref tlsph) for more details on the method.
     where `beam` and `clamped_particles` are of type [`InitialCondition`](@ref).
 """
 struct TotalLagrangianSPHSystem{BM, NDIMS, ELTYPE <: Real, IC, ARRAY1D, ARRAY2D, ARRAY3D,
-                                YM, PR, LL, LM, K, PF, V, ST, M, IM, NHS,
+                                YM, PR, LL, LM, K, PF, V, ST, M, IM, NHS, VA,
                                 C} <: AbstractStructureSystem{NDIMS}
     initial_condition   :: IC
     initial_coordinates :: ARRAY2D # Array{ELTYPE, 2}: [dimension, particle]
@@ -109,6 +113,7 @@ struct TotalLagrangianSPHSystem{BM, NDIMS, ELTYPE <: Real, IC, ARRAY1D, ARRAY2D,
     clamped_particles_motion :: M
     clamped_particles_moving :: IM
     self_interaction_nhs     :: NHS
+    velocity_averaging       :: VA
     cache                    :: C
 end
 
@@ -121,7 +126,8 @@ function TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing
                                                       ndims(smoothing_kernel)),
                                   penalty_force=nothing, viscosity=nothing,
                                   source_terms=nothing, boundary_model=nothing,
-                                  self_interaction_nhs=:default)
+                                  self_interaction_nhs=:default,
+                                  velocity_averaging=nothing)
     NDIMS = ndims(initial_condition)
     ELTYPE = eltype(initial_condition)
     n_particles = nparticles(initial_condition)
@@ -186,7 +192,8 @@ function TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing
     initialize_prescribed_motion!(clamped_particles_motion, initial_condition_sorted,
                                   n_clamped_particles)
 
-    cache = create_cache_tlsph(clamped_particles_motion, initial_condition_sorted)
+    cache = (; create_cache_tlsph(clamped_particles_motion, initial_condition_sorted)...,
+             create_cache_tlsph(velocity_averaging, initial_condition_sorted)...)
 
     return TotalLagrangianSPHSystem(initial_condition_sorted, initial_coordinates,
                                     current_coordinates, mass, correction_matrix,
@@ -197,7 +204,7 @@ function TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing
                                     smoothing_length, acceleration_, boundary_model,
                                     penalty_force, viscosity, source_terms,
                                     clamped_particles_motion, ismoving,
-                                    self_interaction_nhs, cache)
+                                    self_interaction_nhs, velocity_averaging, cache)
 end
 
 # Initialize self-interaction neighborhood search if not provided by the user
@@ -259,7 +266,8 @@ function initialize_self_interaction_nhs(system::TotalLagrangianSPHSystem,
                                     system.viscosity, system.source_terms,
                                     system.clamped_particles_motion,
                                     system.clamped_particles_moving,
-                                    self_interaction_nhs, system.cache)
+                                    self_interaction_nhs, system.velocity_averaging,
+                                    system.cache)
 end
 
 extract_periodic_box(::Nothing) = nothing
@@ -274,6 +282,22 @@ function create_cache_tlsph(::PrescribedMotion, initial_condition)
     return (; velocity, acceleration)
 end
 
+function create_cache_tlsph(::VelocityAveraging, initial_condition)
+    averaged_velocity = zero(initial_condition.velocity)
+    t_last_averaging = Ref(zero(eltype(initial_condition)))
+
+    return (; averaged_velocity, t_last_averaging)
+end
+
+@inline function requires_update_callback(system::TotalLagrangianSPHSystem, semi)
+    # Velocity averaging used and `integrate_tlsph == false` means that split integration
+    # is used and the averaged velocity is updated there.
+    # Velocity averaging used, `integrate_tlsph == true`, and no fluid system in the semi
+    # means we are inside the split integration
+    return !isnothing(system.velocity_averaging) && semi.integrate_tlsph[] &&
+        any(system -> system isa AbstractFluidSystem, semi.systems)
+end
+
 @inline function Base.eltype(::TotalLagrangianSPHSystem{<:Any, <:Any, ELTYPE}) where {ELTYPE}
     return ELTYPE
 end
@@ -384,6 +408,10 @@ end
     return poisson_ratio[particle]
 end
 
+@inline function velocity_averaging(system::TotalLagrangianSPHSystem)
+    return system.velocity_averaging
+end
+
 function initialize!(system::TotalLagrangianSPHSystem, semi)
     (; correction_matrix) = system
 
diff --git a/src/schemes/structure/total_lagrangian_sph/total_lagrangian_sph.jl b/src/schemes/structure/total_lagrangian_sph/total_lagrangian_sph.jl
index 168d6b634e..a9b859b661 100644
--- a/src/schemes/structure/total_lagrangian_sph/total_lagrangian_sph.jl
+++ b/src/schemes/structure/total_lagrangian_sph/total_lagrangian_sph.jl
@@ -1,5 +1,6 @@
 # Penalty force needs to be included first, so that `dv_penalty_force` is available
 # in the closure of `foreach_point_neighbor`.
 include("penalty_force.jl")
+include("velocity_averaging.jl")
 include("system.jl")
 include("viscosity.jl")
diff --git a/src/schemes/structure/total_lagrangian_sph/velocity_averaging.jl b/src/schemes/structure/total_lagrangian_sph/velocity_averaging.jl
new file mode 100644
index 0000000000..f72da49887
--- /dev/null
+++ b/src/schemes/structure/total_lagrangian_sph/velocity_averaging.jl
@@ -0,0 +1,79 @@
+struct VelocityAveraging{ELTYPE}
+    time_constant::ELTYPE
+end
+
+# No velocity averaging for a system by default
+@inline function velocity_averaging(system)
+    return nothing
+end
+
+@inline function velocity_for_viscosity(v, system, particle)
+    return velocity_for_viscosity(v, velocity_averaging(system), system, particle)
+end
+
+@inline function velocity_for_viscosity(v, ::Nothing, system, particle)
+    return current_velocity(v, system, particle)
+end
+
+@inline function velocity_for_viscosity(v, ::VelocityAveraging, system, particle)
+    return extract_svector(system.cache.average_velocity, system, particle)
+end
+
+initialize_averaged_velocity!(system, v_ode, semi, t) = system
+
+function initialize_averaged_velocity!(system, v_ode, semi, t)
+    initialize_averaged_velocity!(system, velocity_averaging(system), v_ode, semi, t)
+end
+
+initialize_averaged_velocity!(system, velocity_averaging, v_ode, semi, t) = system
+
+function initialize_averaged_velocity!(system, ::VelocityAveraging, v_ode, semi, t)
+    v = wrap_v(v_ode, system, semi)
+
+    # Make sure the averaged velocity is zero for non-integrated particles
+    set_zero!(system.cache.averaged_velocity)
+
+    @threaded semi for particle in each_integrated_particle(system)
+        v_particle = current_velocity(v, system, particle)
+
+        for i in eachindex(v_particle)
+            system.cache.averaged_velocity[i, particle] = v_particle[i]
+        end
+    end
+    system.cache.t_last_averaging[] = t
+
+    return system
+end
+
+compute_averaged_velocity!(system, v_ode, semi, t_new) = system
+
+function compute_averaged_velocity!(system, v_ode, semi, t_new)
+    compute_averaged_velocity!(system, velocity_averaging(system), v_ode, semi, t_new)
+end
+
+compute_averaged_velocity!(system, velocity_averaging, v_ode, semi, t_new) = system
+
+function compute_averaged_velocity!(system, velocity_averaging::VelocityAveraging,
+                                    v_ode, semi, t_new)
+    time_constant = velocity_averaging.time_constant
+    averaged_velocity = system.cache.averaged_velocity
+    dt = t_new - system.cache.t_last_averaging[]
+    system.cache.t_last_averaging[] = t_new
+
+    @trixi_timeit timer() "compute averaged velocity" begin
+        v = wrap_v(v_ode, system, semi)
+
+        # Compute an exponential moving average of the velocity with the given time constant
+        alpha = 1 - exp(-dt / time_constant)
+
+        @threaded semi for particle in each_integrated_particle(system)
+            v_particle = current_velocity(v, system, particle)
+
+            for i in eachindex(v_particle)
+                averaged_velocity[i, particle] = (1 - alpha) * averaged_velocity[i, particle] + alpha * v_particle[i]
+            end
+        end
+    end
+
+    return system
+end

From 5ae0f38b8b8f964449b377abe0be41609457dfe9 Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Mon, 23 Feb 2026 18:51:07 +0100
Subject: [PATCH 2/7] Fix averaged velocity

---
 src/callbacks/split_integration.jl            | 20 +++++++++----------
 .../velocity_averaging.jl                     | 10 +++++-----
 2 files changed, 15 insertions(+), 15 deletions(-)

diff --git a/src/callbacks/split_integration.jl b/src/callbacks/split_integration.jl
index 60ef940eeb..71320f8a53 100644
--- a/src/callbacks/split_integration.jl
+++ b/src/callbacks/split_integration.jl
@@ -52,7 +52,16 @@ callback = SplitIntegrationCallback(RDPK3SpFSAL49(), abstol=1e-6, reltol=1e-4)
 ```
 """
 function SplitIntegrationCallback(alg; kwargs...)
-    split_integration_callback = SplitIntegrationCallback(nothing, alg, kwargs)
+    # Add lightweight callback to (potentially) update the averaged velocity
+    # during the split integration.
+    if haskey(kwargs, :callback)
+        # Note that `CallbackSet`s can be nested
+        kwargs = (; kwargs..., callback=CallbackSet(values(kwargs).callback,
+                                                    UpdateAveragedVelocityCallback()))
+    else
+        kwargs = (; kwargs..., callback=UpdateAveragedVelocityCallback())
+    end
+    split_integration_callback = SplitIntegrationCallback(nothing, alg, pairs(kwargs))
 
     # The first one is the `condition`, the second the `affect!`
     return DiscreteCallback(split_integration_callback, split_integration_callback,
@@ -122,15 +131,6 @@ function initialize_split_integration!(cb, u, t, integrator)
     ode_split = DynamicalODEProblem(kick_split!, drift_split!, v0_ode_split, u0_ode_split,
                                     tspan, p)
 
-    # Add lightweight callback to (potentially) update the averaged velocity
-    # during the split integration.
-    callback = UpdateAveragedVelocityCallback()
-    if haskey(kwargs, :callback)
-        kwargs[:callback] = CallbackSet(kwargs[:callback], callback)
-    else
-        kwargs[:callback] = callback
-    end
-
     # Create the split integrator.
     # We need the timer here to keep the output clean because this will call `kick!` once.
     @trixi_timeit timer() "split integration" begin
diff --git a/src/schemes/structure/total_lagrangian_sph/velocity_averaging.jl b/src/schemes/structure/total_lagrangian_sph/velocity_averaging.jl
index f72da49887..79cb99f01f 100644
--- a/src/schemes/structure/total_lagrangian_sph/velocity_averaging.jl
+++ b/src/schemes/structure/total_lagrangian_sph/velocity_averaging.jl
@@ -16,10 +16,12 @@ end
 end
 
 @inline function velocity_for_viscosity(v, ::VelocityAveraging, system, particle)
-    return extract_svector(system.cache.average_velocity, system, particle)
-end
+    if particle <= system.n_integrated_particles
+        return extract_svector(system.cache.averaged_velocity, system, particle)
+    end
 
-initialize_averaged_velocity!(system, v_ode, semi, t) = system
+    return current_clamped_velocity(v, system, system.clamped_particles_motion, particle)
+end
 
 function initialize_averaged_velocity!(system, v_ode, semi, t)
     initialize_averaged_velocity!(system, velocity_averaging(system), v_ode, semi, t)
@@ -45,8 +47,6 @@ function initialize_averaged_velocity!(system, ::VelocityAveraging, v_ode, semi,
     return system
 end
 
-compute_averaged_velocity!(system, v_ode, semi, t_new) = system
-
 function compute_averaged_velocity!(system, v_ode, semi, t_new)
     compute_averaged_velocity!(system, velocity_averaging(system), v_ode, semi, t_new)
 end

From cae086c8fedc44bea49500e443d952e539f0e988 Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Tue, 3 Mar 2026 15:02:52 +0100
Subject: [PATCH 3/7] Fix

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

diff --git a/src/callbacks/split_integration.jl b/src/callbacks/split_integration.jl
index 71320f8a53..e62bf7c993 100644
--- a/src/callbacks/split_integration.jl
+++ b/src/callbacks/split_integration.jl
@@ -380,7 +380,7 @@ end
 # `initialize`
 function initialize_averaged_velocity_callback!(cb, vu_ode, t, integrator)
     v_ode, u_ode = vu_ode.x
-    semi = integrator.p
+    semi = integrator.p.semi_split
 
     foreach_system(semi) do system
         initialize_averaged_velocity!(system, v_ode, semi, t)
@@ -397,7 +397,7 @@ end
 # `affect!`
 function update_averaged_velocity_callback!(integrator)
     t_new = integrator.t
-    semi = integrator.p
+    semi = integrator.p.semi_split
     v_ode, u_ode = integrator.u.x
 
     foreach_system(semi) do system

From 2a80fe575503d0f406361c2fc648d7b1aa88855f Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Tue, 3 Mar 2026 16:14:49 +0100
Subject: [PATCH 4/7] Add docs

---
 docs/src/systems/total_lagrangian_sph.md      | 23 +++++++++++
 src/TrixiParticles.jl                         |  2 +-
 src/callbacks/split_integration.jl            |  7 ++++
 .../structure/total_lagrangian_sph/system.jl  |  2 +-
 .../velocity_averaging.jl                     | 38 +++++++++++++++++++
 5 files changed, 70 insertions(+), 2 deletions(-)

diff --git a/docs/src/systems/total_lagrangian_sph.md b/docs/src/systems/total_lagrangian_sph.md
index 86167d164b..923c23436f 100644
--- a/docs/src/systems/total_lagrangian_sph.md
+++ b/docs/src/systems/total_lagrangian_sph.md
@@ -129,3 +129,26 @@ density are present, instabilities in the fluid can be induced by the structure.
 In these cases, artificial viscosity is effective at stabilizing the fluid close to the
 structure, and we recommend using it in combination with penalty force to both
 prevent hourglass modes and stabilize the fluid close to the fluid-structure interface.
+
+
+## [Velocity Averaging](@id velocity_averaging)
+
+In challenging FSI cases with very stiff structures, sometimes the two techniques
+above are not sufficient to prevent instabilities.
+High-frequency noise in the structure velocity can trigger instabilities or spurious
+pressure waves (due to aliasing) in the fluid.
+Another stabilization technique is an exponential moving average (EMA) of the structure
+velocity used **only** for the fluid-structure viscous coupling (no-slip boundary condition).
+
+The averaged velocity ``\bar v`` is updated (by the [`UpdateCallback`](@ref)
+or in every sub-step of the [`SplitIntegrationCallback`](@ref) if it is used) as
+```math
+\bar v^{n+1} = (1-\alpha)\,\bar v^n + \alpha\, v^{n+1}, \qquad
+\alpha = 1 - \exp(-\Delta t/\tau),
+```
+where ``\tau`` is the time constant.
+
+```@autodocs
+Modules = [TrixiParticles]
+Pages = [joinpath("schemes", "structure", "total_lagrangian_sph", "velocity_averaging.jl")]
+```
diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl
index 5530a79a86..2a0c15afea 100644
--- a/src/TrixiParticles.jl
+++ b/src/TrixiParticles.jl
@@ -76,7 +76,7 @@ export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback,
 export ContinuityDensity, SummationDensity
 export PenaltyForceGanzenmueller, TransportVelocityAdami, ParticleShiftingTechnique,
        ParticleShiftingTechniqueSun2017, ConsistentShiftingSun2019,
-       ContinuityEquationTermSun2019, MomentumEquationTermSun2019
+       ContinuityEquationTermSun2019, MomentumEquationTermSun2019, VelocityAveraging
 export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel,
        SchoenbergQuinticSplineKernel, GaussianKernel, WendlandC2Kernel, WendlandC4Kernel,
        WendlandC6Kernel, SpikyKernel, Poly6Kernel, LaguerreGaussKernel
diff --git a/src/callbacks/split_integration.jl b/src/callbacks/split_integration.jl
index e62bf7c993..b4bf1f3730 100644
--- a/src/callbacks/split_integration.jl
+++ b/src/callbacks/split_integration.jl
@@ -48,6 +48,7 @@ callback = SplitIntegrationCallback(RDPK3SpFSAL49(), abstol=1e-6, reltol=1e-4)
 │ alg: …………………………………………………………………… RDPK3SpFSAL49                                                    │
 │ abstol: …………………………………………………………… 1.0e-6                                                           │
 │ reltol: …………………………………………………………… 0.0001                                                           │
+│ callback: ……………………………………………………… UpdateAveragedVelocityCallback                                   │
 └──────────────────────────────────────────────────────────────────────────────────────────────────┘
 ```
 """
@@ -409,3 +410,9 @@ function update_averaged_velocity_callback!(integrator)
 
     return integrator
 end
+
+function Base.show(io::IO,
+                   cb::DiscreteCallback{typeof(update_averaged_velocity_callback!)})
+    @nospecialize cb # reduce precompilation time
+    print(io, "UpdateAveragedVelocityCallback")
+end
diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl
index f8e13560c7..a9919e4c26 100644
--- a/src/schemes/structure/total_lagrangian_sph/system.jl
+++ b/src/schemes/structure/total_lagrangian_sph/system.jl
@@ -64,7 +64,7 @@ See [Total Lagrangian SPH](@ref tlsph) for more details on the method.
                     Alternatively, a user-defined neighborhood search can be passed here.
 - `velocity_averaging`: Velocity averaging technique to be applied on the velocity field
                     to obtain an averaged velocity to be used in fluid-structure interaction.
-                    See [`VelocityAveraging`](@ref) for details.
+                    See [the docs](@ref velocity_averaging) for details.
 
 !!! note
     If specifying the clamped particles manually (via `n_clamped_particles`),
diff --git a/src/schemes/structure/total_lagrangian_sph/velocity_averaging.jl b/src/schemes/structure/total_lagrangian_sph/velocity_averaging.jl
index 79cb99f01f..4144c14a47 100644
--- a/src/schemes/structure/total_lagrangian_sph/velocity_averaging.jl
+++ b/src/schemes/structure/total_lagrangian_sph/velocity_averaging.jl
@@ -1,5 +1,43 @@
+@doc raw"""
+    VelocityAveraging(; time_constant)
+
+Exponential moving average (EMA) of the structure velocity used **only** for the
+fluid-structure viscous coupling (no-slip boundary condition) of a
+[`TotalLagrangianSPHSystem`](@ref).
+
+This is a stabilization option that is typically only needed for challenging FSI
+cases with very stiff structures, where high-frequency noise in the structure velocity
+can trigger instabilities or spurious pressure waves (due to aliasing) in the fluid.
+
+See [the docs](@ref velocity_averaging) for more details.
+
+!!! note "Callback Required"
+    Velocity averaging requires either an [`UpdateCallback`](@ref)
+    or a [`SplitIntegrationCallback`](@ref) to be used in the simulation.
+    In the typical use case of stabilizing a challenging FSI simulation, a
+    [`SplitIntegrationCallback`](@ref) will usually be used anyway for performance reasons.
+
+# Keywords
+- `time_constant`: Time constant ``\tau`` of the EMA (same units as simulation time).
+                   Smaller ``\tau`` smooths only very fast oscillations, which is often
+                   sufficient for stabilizing FSI simulations; larger ``\tau``
+                   can also suppress slower oscillations that are causing spurious
+                   pressure waves without destabilizing the simulation.
+
+# Examples
+```jldoctest; output=false
+VelocityAveraging(time_constant=1e-5)
+
+# output
+VelocityAveraging{Float64}(1.0e-5)
+```
+"""
 struct VelocityAveraging{ELTYPE}
     time_constant::ELTYPE
+
+    function VelocityAveraging(; time_constant)
+        return new{typeof(time_constant)}(time_constant)
+    end
 end
 
 # No velocity averaging for a system by default

From e12592a944f9ea28397e961492fca891d4054186 Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Tue, 3 Mar 2026 16:17:19 +0100
Subject: [PATCH 5/7] Reformat

---
 src/callbacks/split_integration.jl                         | 5 +++--
 src/schemes/structure/total_lagrangian_sph/system.jl       | 2 +-
 .../structure/total_lagrangian_sph/velocity_averaging.jl   | 7 +++++--
 3 files changed, 9 insertions(+), 5 deletions(-)

diff --git a/src/callbacks/split_integration.jl b/src/callbacks/split_integration.jl
index b4bf1f3730..b3b9c822fe 100644
--- a/src/callbacks/split_integration.jl
+++ b/src/callbacks/split_integration.jl
@@ -57,8 +57,9 @@ function SplitIntegrationCallback(alg; kwargs...)
     # during the split integration.
     if haskey(kwargs, :callback)
         # Note that `CallbackSet`s can be nested
-        kwargs = (; kwargs..., callback=CallbackSet(values(kwargs).callback,
-                                                    UpdateAveragedVelocityCallback()))
+        kwargs = (; kwargs...,
+                  callback=CallbackSet(values(kwargs).callback,
+                                       UpdateAveragedVelocityCallback()))
     else
         kwargs = (; kwargs..., callback=UpdateAveragedVelocityCallback())
     end
diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl
index a9919e4c26..05c7c2da7d 100644
--- a/src/schemes/structure/total_lagrangian_sph/system.jl
+++ b/src/schemes/structure/total_lagrangian_sph/system.jl
@@ -295,7 +295,7 @@ end
     # Velocity averaging used, `integrate_tlsph == true`, and no fluid system in the semi
     # means we are inside the split integration
     return !isnothing(system.velocity_averaging) && semi.integrate_tlsph[] &&
-        any(system -> system isa AbstractFluidSystem, semi.systems)
+           any(system -> system isa AbstractFluidSystem, semi.systems)
 end
 
 @inline function Base.eltype(::TotalLagrangianSPHSystem{<:Any, <:Any, ELTYPE}) where {ELTYPE}
diff --git a/src/schemes/structure/total_lagrangian_sph/velocity_averaging.jl b/src/schemes/structure/total_lagrangian_sph/velocity_averaging.jl
index 4144c14a47..a496ac1ead 100644
--- a/src/schemes/structure/total_lagrangian_sph/velocity_averaging.jl
+++ b/src/schemes/structure/total_lagrangian_sph/velocity_averaging.jl
@@ -105,10 +105,13 @@ function compute_averaged_velocity!(system, velocity_averaging::VelocityAveragin
         alpha = 1 - exp(-dt / time_constant)
 
         @threaded semi for particle in each_integrated_particle(system)
-            v_particle = current_velocity(v, system, particle)
+            v_particle = @inbounds current_velocity(v, system, particle)
 
             for i in eachindex(v_particle)
-                averaged_velocity[i, particle] = (1 - alpha) * averaged_velocity[i, particle] + alpha * v_particle[i]
+                @inbounds averaged_velocity[i,
+                                            particle] = (1 - alpha) *
+                                                        averaged_velocity[i, particle] +
+                                                        alpha * v_particle[i]
             end
         end
     end

From b74d8d014de79b6bba882da2ef9960209dd860d2 Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Tue, 3 Mar 2026 16:54:01 +0100
Subject: [PATCH 6/7] Add test

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

diff --git a/test/examples/examples.jl b/test/examples/examples.jl
index c7d9fadb92..55c0707387 100644
--- a/test/examples/examples.jl
+++ b/test/examples/examples.jl
@@ -143,6 +143,8 @@
                                              plate_position=(0.2, 0.0),
                                              tspan=(0.0, 0.2),
                                              E=1e7, # Stiffer plate
+                                             # Test velocity averaging as well
+                                             velocity_averaging=VelocityAveraging(time_constant=1e-4),
                                              maxiters=400,
                                              extra_callback=split_integration) [
                 r"\[ Info: To create the self-interaction neighborhood search.*\n"
@@ -156,6 +158,15 @@
                 @test count_rhs_allocations(sol, semi) == 0
             end
 
+            # Verify that velocity averaging is actually used
+            system = sol.prob.p.systems[3]
+            @test system isa TotalLagrangianSPHSystem
+            # Test that `current_velocity` is not the same as `velocity_for_viscosity`.
+            # `current_velocity` should fail because it tries to access `v`, for which
+            # we passed `nothing`.
+            @test_throws "no method" TrixiParticles.current_velocity(nothing, system, 1)
+            @test !iszero(TrixiParticles.velocity_for_viscosity(nothing, system, 1))
+
             # Now use split integration and verify that it is actually used for TLSPH
             # by using a time step that is too large and verifying that it is crashing.
             split_integration = SplitIntegrationCallback(RDPK3SpFSAL35(), adaptive=false,

From f66e10267a32366ab340987348d9ababfe316a5d Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Tue, 3 Mar 2026 16:55:23 +0100
Subject: [PATCH 7/7] Fix tests

---
 docs/src/systems/total_lagrangian_sph.md | 4 ++--
 examples/fsi/dam_break_plate_2d.jl       | 3 ++-
 2 files changed, 4 insertions(+), 3 deletions(-)

diff --git a/docs/src/systems/total_lagrangian_sph.md b/docs/src/systems/total_lagrangian_sph.md
index 923c23436f..d9495c08e6 100644
--- a/docs/src/systems/total_lagrangian_sph.md
+++ b/docs/src/systems/total_lagrangian_sph.md
@@ -133,8 +133,8 @@ prevent hourglass modes and stabilize the fluid close to the fluid-structure int
 
 ## [Velocity Averaging](@id velocity_averaging)
 
-In challenging FSI cases with very stiff structures, sometimes the two techniques
-above are not sufficient to prevent instabilities.
+In challenging FSI cases with very stiff structures, the two techniques above might not be
+sufficient to prevent instabilities.
 High-frequency noise in the structure velocity can trigger instabilities or spurious
 pressure waves (due to aliasing) in the fluid.
 Another stabilization technique is an exponential moving average (EMA) of the structure
diff --git a/examples/fsi/dam_break_plate_2d.jl b/examples/fsi/dam_break_plate_2d.jl
index 44d86c5317..176a3efd68 100644
--- a/examples/fsi/dam_break_plate_2d.jl
+++ b/examples/fsi/dam_break_plate_2d.jl
@@ -129,7 +129,8 @@ structure_system = TotalLagrangianSPHSystem(structure,
                                             E, nu, boundary_model=boundary_model_structure,
                                             clamped_particles=1:nparticles(clamped_particles),
                                             acceleration=(0.0, -gravity),
-                                            penalty_force=PenaltyForceGanzenmueller(alpha=0.01))
+                                            penalty_force=PenaltyForceGanzenmueller(alpha=0.01),
+                                            velocity_averaging=nothing)
 
 # ==========================================================================================
 # ==== Simulation