Implement velocity averaging for TLSPH

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

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

src/general/semidiscretization.jl

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

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

@@ -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

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
 

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
 

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")

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)
+    if particle <= system.n_integrated_particles
+        return extract_svector(system.cache.averaged_velocity, system, particle)
+    end
+
+    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)
+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
+
+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