Implement velocity averaging for TLSPH (#1069)

* Implement velocity averaging for TLSPH

* Fix averaged velocity

* Fix

* Add docs

* Reformat

* Add test

* Fix tests

* Add unittests

Co-authored-by: Copilot <copilot@github.com>

* Fix duplicate velocity averaging with `UpdateCallback` and split integration

* Implement suggestions

* Reformat

* Fix tests

* Fix tests

* Reformat

* Implement suggestions

---------

Co-authored-by: Sven Berger <berger.sven@gmail.com>
Co-authored-by: Copilot <copilot@github.com>
Co-authored-by: Niklas Neher <73897120+LasNikas@users.noreply.github.com>

Author: 4 people

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 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 in the fluid due to aliasing.
+Another stabilization technique is to use an exponential moving average (EMA) of the structure
+velocity used **only** for the fluid-structure viscous coupling (i.e., the no-slip boundary condition).
+
+The averaged velocity ``\bar v`` is updated (by the [`UpdateCallback`](@ref)
+or at 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")]
+```

examples/fsi/dam_break_plate_2d.jl

@@ -138,7 +138,8 @@ structure_system = TotalLagrangianSPHSystem(structure;
                                             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

src/TrixiParticles.jl

@@ -79,7 +79,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

src/callbacks/split_integration.jl

@@ -82,14 +82,24 @@ callback = SplitIntegrationCallback(CarpenterKennedy2N54(williamson_condition=fa
 │ stage_coupling: ……………………………………… true                                                             │
 │ predict_positions: ……………………………… true                                                             │
 │ dt: ……………………………………………………………………… 1.0                                                              │
-│ callback: ……………………………………………………… StepsizeCallback(is_constant=true, cfl_number=1.6)               │
+│ callback: ……………………………………………………… CallbackSet{Tuple{}, Tuple{Discr…pdateAveragedVelocityCallback)) │
 └──────────────────────────────────────────────────────────────────────────────────────────────────┘
 ```
 """
 function SplitIntegrationCallback(alg; stage_coupling=false, predict_positions=true,
                                   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(alg, stage_coupling,
-                                                          predict_positions, kwargs)
+                                                          predict_positions, pairs(kwargs))
 
     # The first one is the `condition`, the second the `affect!`
     return DiscreteCallback(split_integration_callback, split_integration_callback,
@@ -522,3 +532,53 @@ 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.semi
+
+    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.semi
+    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
+
+function Base.show(io::IO,
+                   cb::DiscreteCallback{typeof(update_averaged_velocity_callback!)})
+    @nospecialize cb # reduce precompilation time
+    print(io, "UpdateAveragedVelocityCallback")
+end

src/callbacks/update.jl

@@ -52,12 +52,21 @@ 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.semi
 
     # Tell the semidiscretization that the `UpdateCallback` is used
     semi.update_callback_used[] = true
 
+    # If TLSPH is not integrated, the averaged velocity will be initialized in the
+    # split integration.
+    if semi.integrate_tlsph[]
+        foreach_system(semi) do system
+            initialize_averaged_velocity!(system, v_ode, semi, t)
+        end
+    end
+
     return cb(integrator)
 end
 
@@ -105,6 +114,14 @@ function (update_callback!::UpdateCallback)(integrator)
             particle_shifting_from_callback!(u_ode, shifting_technique(system), system,
                                              v_ode, semi, integrator)
         end
+
+        # If TLSPH is not integrated, the averaged velocity will be updated in the
+        # split integration.
+        if semi.integrate_tlsph[]
+            foreach_system(semi) do system
+                compute_averaged_velocity!(system, v_ode, semi, t)
+            end
+        end
     end
 
     return integrator

src/general/semidiscretization.jl

@@ -764,7 +764,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

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)

src/schemes/boundary/open_boundary/system.jl

@@ -245,7 +245,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

src/schemes/boundary/wall_boundary/dummy_particles.jl

@@ -731,8 +731,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)
         # Compute the no-slip wall velocity using the formula v_w = 2 * v_a - ṽ_fluid,

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