Extend shared rigid contact runtime model

docs/src/systems/rigid_body.md

@@ -17,6 +17,37 @@ Rigid contact is configured through the contact model. This is separate from the
 boundary model used for fluid-structure interaction; see
 [Boundary Models](@ref boundary_models) for that part of the rigid-body setup.
 
+`RigidContactModel` currently defines a normal spring-dashpot contact law with the
+parameters `normal_stiffness`, `normal_damping`, and `contact_distance`.
+
+The current implementation uses the same model for rigid-wall and rigid-rigid contact:
+
+- rigid-wall contact groups penetrating wall neighbors into a small number of contact
+  manifolds per rigid particle and applies one normal contact force per manifold,
+- rigid-rigid contact evaluates direct pairwise normal contact forces between rigid
+  particles,
+- and both paths are currently normal-only, i.e. there are no tangential/frictional
+  forces or contact-history terms yet.
+
+Here, a contact manifold is a discrete approximation of one locally smooth contact
+patch. A rigid particle touching a flat wall will usually produce one manifold,
+while corners or edges can produce several.
+
+The number of cached rigid-wall manifolds per rigid particle is controlled by the
+`RigidBodySystem(...; max_manifolds=8)` keyword argument. If more wall-contact
+patches are detected than cached manifold slots are available, the implementation
+falls back to the best-matching existing manifold for that particle.
+
+`contact_distance` defines when contact starts. If `contact_distance == 0`, the
+particle spacing of the `RigidBodySystem` is used.
+
+If no `contact_model` is specified for a rigid body, rigid-wall and rigid-rigid contact
+for that system are disabled.
+
+For output and postprocessing, rigid bodies also expose the diagnostics
+`contact_count` and `max_contact_penetration`. They are available through rigid-body
+system data and VTK output.
+
 ```@autodocs
 Modules = [TrixiParticles]
 Pages = [joinpath("schemes", "structure", "rigid_body", "contact_models.jl")]

src/general/abstract_system.jl

@@ -179,7 +179,11 @@ function update_boundary_interpolation!(system, v, u, v_ode, u_ode, semi, t)
     return system
 end
 
-function update_final!(system, v, u, v_ode, u_ode, semi, t)
+function update_final!(system, v, u, v_ode, u_ode, semi, t; kwargs...)
+    return system
+end
+
+function reset_interaction_caches!(system)
     return system
 end
 

src/general/semidiscretization.jl

@@ -529,6 +529,16 @@ function update_systems_and_nhs(v_ode, u_ode, semi, t)
     end
 end
 
+# Some systems accumulate pairwise interaction state outside `dv_ode`. Reset that state once
+# at the beginning of every explicitly assembled interaction pass.
+function reset_interaction_caches!(semi::Union{NamedTuple, Semidiscretization})
+    foreach_system(semi) do system
+        reset_interaction_caches!(system)
+    end
+
+    return semi
+end
+
 # The `SplitIntegrationCallback` overwrites `semi_wrap` to use a different
 # semidiscretization for wrapping arrays.
 # TODO `semi` is not used yet, but will be used when the source terms API is modified
@@ -668,6 +678,8 @@ end
 end
 
 function system_interaction!(dv_ode, v_ode, u_ode, semi)
+    reset_interaction_caches!(semi)
+
     # Call `interact!` for each pair of systems
     foreach_system(semi) do system
         foreach_system(semi) do neighbor
@@ -700,6 +712,8 @@ end
 # Function barrier to make benchmarking interactions easier.
 # One can benchmark, e.g. the fluid-fluid interaction, with:
 # dv_ode, du_ode = copy(sol.u[end]).x; v_ode, u_ode = copy(sol.u[end]).x;
+# For manual multi-pair interaction assembly, call `reset_interaction_caches!(semi)` once
+# before the first direct `interact!` call.
 # @btime TrixiParticles.interact!($dv_ode, $v_ode, $u_ode, $fluid_system, $fluid_system, $semi);
 @inline function interact!(dv_ode, v_ode, u_ode, system, neighbor, semi; timer_str="")
     dv = wrap_v(dv_ode, system, semi)

src/io/write_vtk.jl

@@ -433,6 +433,8 @@ function write2vtk!(vtk, v, u, t, system::RigidBodySystem)
     vtk["resultant_torque"] = [system.resultant_torque[]]
     vtk["angular_acceleration_force"] = [system.angular_acceleration_force[]]
     vtk["gyroscopic_acceleration"] = [system.gyroscopic_acceleration[]]
+    vtk["contact_count"] = [system.cache.contact_count[]]
+    vtk["max_contact_penetration"] = [system.cache.max_contact_penetration[]]
 
     write2vtk!(vtk, v, u, t, system.boundary_model, system)
 end

src/schemes/fluid/entropically_damped_sph/system.jl

@@ -306,7 +306,8 @@ 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)
+function update_final!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, semi, t;
+                       kwargs...)
     (; surface_tension) = system
 
     # Surface normal of neighbor and boundary needs to have been calculated already

src/schemes/fluid/weakly_compressible_sph/system.jl

@@ -338,7 +338,8 @@ function update_pressure!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_od
     return system
 end
 
-function update_final!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t)
+function update_final!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t;
+                       kwargs...)
     (; surface_tension) = system
 
     # Surface normal of neighbor and boundary needs to have been calculated already

src/schemes/structure/rigid_body/contact_models.jl

@@ -5,10 +5,10 @@ abstract type AbstractRigidContactModel end
                       normal_damping=0.0,
                       contact_distance=0.0)
 
-Basic rigid contact model stored on a rigid body.
-It is currently used for both rigid-wall and rigid-rigid contact.
-The contact force consists of a linear normal spring-dashpot contribution only.
-If `contact_distance == 0`, the particle spacing of the `RigidBodySystem` will be used as contact distance.
+Shared rigid-contact model used by the active rigid-wall and rigid-rigid contact paths.
+The current contact force consists of a linear normal spring-dashpot contribution only.
+If `contact_distance == 0`, the particle spacing of the `RigidBodySystem` will be used
+as contact distance.
 
 !!! warning "Experimental implementation"
     This is an experimental feature and may change in future releases.

src/schemes/structure/rigid_body/rhs.jl

@@ -108,6 +108,10 @@ function interact!(dv, v_particle_system, u_particle_system,
 
     NDIMS = ndims(particle_system)
     ELTYPE = eltype(particle_system)
+    set_zero!(particle_system.cache.contact_count_per_particle)
+    set_zero!(particle_system.cache.max_contact_penetration_per_particle)
+    contact_count_per_particle = particle_system.cache.contact_count_per_particle
+    max_contact_penetration_per_particle = particle_system.cache.max_contact_penetration_per_particle
 
     # First gather all penetrating wall neighbors into a small set of contact manifolds per
     # rigid particle. This avoids applying one noisy contact force per wall particle.
@@ -128,6 +132,8 @@ function interact!(dv, v_particle_system, u_particle_system,
     # Apply one force contribution per manifold using the averaged normal, penetration, and
     # wall velocity stored in the cache.
     @threaded semi for particle in each_integrated_particle(particle_system)
+        local_contact_count = 0
+        local_max_contact_penetration = zero(ELTYPE)
         n_manifolds = particle_system.cache.contact_manifold_count[particle]
         if n_manifolds > 0
             v_particle = current_velocity(v_particle_system, particle_system, particle)
@@ -160,7 +166,6 @@ function interact!(dv, v_particle_system, u_particle_system,
                         relative_velocity = v_particle - v_boundary
                         normal_velocity = dot(relative_velocity, normal)
 
-                        # Only the normal spring-dashpot part is kept in the basic collision.
                         elastic_force = contact_model.normal_stiffness *
                                         penetration_effective
                         damping_force = -contact_model.normal_damping * normal_velocity
@@ -174,13 +179,24 @@ function interact!(dv, v_particle_system, u_particle_system,
                                 particle_system.force_per_particle[dim,
                                                                    particle] += interaction_force[dim]
                             end
+
+                            local_contact_count += 1
+                            local_max_contact_penetration = max(local_max_contact_penetration,
+                                                                penetration_effective)
                         end
                     end
                 end
             end
         end
+
+        contact_count_per_particle[particle] = local_contact_count
+        max_contact_penetration_per_particle[particle] = local_max_contact_penetration
     end
 
+    particle_system.cache.contact_count[] += sum(contact_count_per_particle)
+    particle_system.cache.max_contact_penetration[] = max(particle_system.cache.max_contact_penetration[],
+                                                          maximum(max_contact_penetration_per_particle))
+
     return dv
 end
 
@@ -331,6 +347,14 @@ function interact!(dv, v_particle_system, u_particle_system,
     ELTYPE = eltype(particle_system)
     system_coords = current_coordinates(u_particle_system, particle_system)
     neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system)
+    set_zero!(particle_system.cache.contact_count_per_particle)
+    set_zero!(particle_system.cache.max_contact_penetration_per_particle)
+    contact_count_per_particle = particle_system.cache.contact_count_per_particle
+    max_contact_penetration_per_particle = particle_system.cache.max_contact_penetration_per_particle
+    pair_normal_stiffness = (contact_model.normal_stiffness +
+                             neighbor_contact_model.normal_stiffness) / 2
+    pair_normal_damping = (contact_model.normal_damping +
+                           neighbor_contact_model.normal_damping) / 2
 
     foreach_point_neighbor(particle_system, neighbor_system, system_coords, neighbor_coords,
                            semi;
@@ -354,10 +378,8 @@ function interact!(dv, v_particle_system, u_particle_system,
         relative_velocity = particle_velocity - neighbor_velocity
         normal_velocity = dot(relative_velocity, normal)
 
-        elastic_force = (contact_model.normal_stiffness +
-                         neighbor_contact_model.normal_stiffness) / 2 * penetration
-        damping_force = -(contact_model.normal_damping +
-                          neighbor_contact_model.normal_damping) / 2 * normal_velocity
+        elastic_force = pair_normal_stiffness * penetration
+        damping_force = -pair_normal_damping * normal_velocity
         normal_force_magnitude = max(elastic_force + damping_force, zero(ELTYPE))
         normal_force_magnitude <= 0 && return dv
 
@@ -366,8 +388,18 @@ function interact!(dv, v_particle_system, u_particle_system,
         for dim in 1:ndims(particle_system)
             particle_system.force_per_particle[dim, particle] += interaction_force[dim]
         end
+
+        # `foreach_point_neighbor` parallelizes the outer loop over `points`, i.e. particles.
+        # This makes these per-particle reductions race-free under the regular backends.
+        contact_count_per_particle[particle] += 1
+        max_contact_penetration_per_particle[particle] = max(max_contact_penetration_per_particle[particle],
+                                                             penetration)
     end
 
+    particle_system.cache.contact_count[] += sum(contact_count_per_particle)
+    particle_system.cache.max_contact_penetration[] = max(particle_system.cache.max_contact_penetration[],
+                                                          maximum(max_contact_penetration_per_particle))
+
     return dv
 end
 

src/schemes/structure/rigid_body/system.jl

@@ -119,7 +119,8 @@ function RigidBodySystem(initial_condition; boundary_model=nothing,
         inverse_inertia = Ref(zero(SMatrix{3, 3, ELTYPE, 9}))
     end
 
-    cache = (;
+    cache = (; contact_count=Ref(0),
+             max_contact_penetration=Ref(zero(ELTYPE)),
              create_cache_contact_manifold(contact_model_, Val(NDIMS), ELTYPE,
                                            nparticles(initial_condition),
                                            max_manifolds_)...,
@@ -149,15 +150,17 @@ function create_cache_contact_manifold(::Nothing, ::Val{NDIMS}, ELTYPE,
     return (;)
 end
 
-# Allocate per-particle manifold scratch arrays for rigid-wall contact.
+# Allocate per-particle scratch arrays for rigid contact.
 #
-# The cache shape is `[dimension, manifold, particle]` for vector-valued sums and
+# The manifold cache shape is `[dimension, manifold, particle]` for vector-valued sums and
 # `[manifold, particle]` for scalar sums. It is rebuilt for each rigid-wall system pair in
-# the RHS and therefore acts purely as transient manifold assembly storage; the persistent
-# collision effect is the force accumulated in `force_per_particle`.
+# the RHS and therefore acts purely as transient manifold assembly storage. The per-particle
+# diagnostic scratch is reused for both rigid-wall and rigid-rigid contact reductions.
 function create_cache_contact_manifold(contact_model, ::Val{NDIMS}, ELTYPE,
                                        n_particles, max_manifolds) where {NDIMS}
-    return (; contact_manifold_count=zeros(Int, n_particles),
+    return (; contact_count_per_particle=zeros(Int, n_particles),
+            max_contact_penetration_per_particle=zeros(ELTYPE, n_particles),
+            contact_manifold_count=zeros(Int, n_particles),
             contact_manifold_weight_sum=zeros(ELTYPE, max_manifolds, n_particles),
             contact_manifold_penetration_sum=zeros(ELTYPE, max_manifolds, n_particles),
             contact_manifold_normal_sum=zeros(ELTYPE, NDIMS, max_manifolds, n_particles),
@@ -374,6 +377,14 @@ function update_boundary_interpolation!(boundary_model, system::RigidBodySystem,
     return system
 end
 
+function reset_interaction_caches!(system::RigidBodySystem)
+    set_zero!(system.force_per_particle)
+    system.cache.contact_count[] = 0
+    system.cache.max_contact_penetration[] = zero(eltype(system))
+
+    return system
+end
+
 function update_final!(system::RigidBodySystem, v, u, v_ode, u_ode, semi, t)
     system_coords = current_coordinates(u, system)
     system_velocity = current_velocity(v, system)
@@ -387,10 +398,6 @@ function update_final!(system::RigidBodySystem, v, u, v_ode, u_ode, semi, t)
                                                         center_of_mass_velocity;
                                                         relative_coordinates=system.relative_coordinates)
 
-    # Reset interaction caches before RHS assembly so pairwise rigid-fluid forces can
-    # accumulate from scratch and non-RHS update paths do not expose stale resultants.
-    set_zero!(system.force_per_particle)
-
     system.center_of_mass[] = center_of_mass
     system.center_of_mass_velocity[] = center_of_mass_velocity
     system.inertia[] = rotational_kinematics.inertia
@@ -555,6 +562,8 @@ function system_data(system::RigidBodySystem, dv_ode, du_ode, v_ode, u_ode, semi
     resultant_torque = system.resultant_torque[]
     angular_acceleration_force = system.angular_acceleration_force[]
     gyroscopic_acceleration = system.gyroscopic_acceleration[]
+    contact_count = system.cache.contact_count[]
+    max_contact_penetration = system.cache.max_contact_penetration[]
     relative_coordinates = system.relative_coordinates
 
     return (; coordinates, velocity, mass=system.mass,
@@ -564,6 +573,7 @@ function system_data(system::RigidBodySystem, dv_ode, du_ode, v_ode, u_ode, semi
             angular_velocity,
             resultant_force, resultant_torque,
             angular_acceleration_force, gyroscopic_acceleration,
+            contact_count, max_contact_penetration,
             density, pressure, acceleration)
 end
 
@@ -573,6 +583,7 @@ function available_data(::RigidBodySystem)
             :center_of_mass, :center_of_mass_velocity,
             :angular_velocity, :resultant_force, :resultant_torque,
             :angular_acceleration_force, :gyroscopic_acceleration,
+            :contact_count, :max_contact_penetration,
             :density, :pressure, :acceleration)
 end
 

test/io/read_vtk.jl

@@ -199,6 +199,10 @@
                   rigid_system.angular_acceleration_force[]
             @test only(Array(TrixiParticles.ReadVTK.get_data(field_data["gyroscopic_acceleration"]))) ==
                   rigid_system.gyroscopic_acceleration[]
+            @test only(Array(TrixiParticles.ReadVTK.get_data(field_data["contact_count"]))) ==
+                  rigid_system.cache.contact_count[]
+            @test only(Array(TrixiParticles.ReadVTK.get_data(field_data["max_contact_penetration"]))) ==
+                  rigid_system.cache.max_contact_penetration[]
         end
     end
 end