Hotfix for Monaghan-Kajtar repulsive boundary particles (#354)

* Hotfix for Monaghan-Kajtar repulsive boundary particles

* Remove singularity in repulsive force

* Update comment

* Add tests for repulsive particles

src/schemes/boundary/boundary.jl

@@ -1,3 +1,4 @@
 include("dummy_particles/dummy_particles.jl")
-include("monaghan_kajtar/monaghan_kajtar.jl")
 include("system.jl")
+# Monaghan-Kajtar repulsive boundary particles require the `BoundarySPHSystem`
+# and the `TotalLagrangianSPHSystem` and are therefore included later.

src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl

@@ -77,15 +77,27 @@ function Base.show(io::IO, model::BoundaryModelMonaghanKajtar)
 end
 
 @inline function pressure_acceleration(particle_system,
-                                       neighbor_system::Union{BoundarySystem{<:BoundaryModelMonaghanKajtar},
-                                                              SolidSystem{<:BoundaryModelMonaghanKajtar}},
+                                       neighbor_system::Union{BoundarySPHSystem{<:BoundaryModelMonaghanKajtar},
+                                                              TotalLagrangianSPHSystem{<:BoundaryModelMonaghanKajtar}},
                                        neighbor, m_a, m_b, p_a, p_b, rho_a, rho_b,
                                        pos_diff, distance, grad_kernel,
                                        pressure_correction, correction)
     (; K, beta, boundary_particle_spacing) = neighbor_system.boundary_model
 
+    # This is `distance - boundary_particle_spacing` in the paper. This factor makes
+    # the force grow infinitely close to the boundary, with a singularity where
+    # `distance` = `boundary_particle_spacing`. However, when the time step is large
+    # enough for a particle to end up behind the singularity in a time integration stage,
+    # the force will switch sign and become smaller again.
+    #
+    # In order to avoid this, we clip the force at a "large" value, large enough to prevent
+    # penetration when a reasonable `K` is used, but small enough to not cause instabilites
+    # or super small time steps.
+    distance_from_singularity = max(0.01 * boundary_particle_spacing,
+                                    distance - boundary_particle_spacing)
+
     return K / beta^(ndims(particle_system) - 1) * pos_diff /
-           (distance * (distance - boundary_particle_spacing)) *
+           (distance * distance_from_singularity) *
            boundary_kernel(distance, particle_system.smoothing_length)
 end
 

src/schemes/schemes.jl

@@ -3,6 +3,9 @@
 include("fluid/fluid.jl")
 include("boundary/boundary.jl")
 include("solid/total_lagrangian_sph/total_lagrangian_sph.jl")
+# Monaghan-Kajtar repulsive boundary particles require the `BoundarySPHSystem`
+# and the `TotalLagrangianSPHSystem`.
+include("boundary/monaghan_kajtar/monaghan_kajtar.jl")
 
 # Include rhs for all schemes
 include("fluid/weakly_compressible_sph/rhs.jl")

test/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl

@@ -1,9 +1,69 @@
 
-@testset verbose=true "Dummy Particles" begin
-    @testset "show" begin
+@testset verbose=true "Monghan-Kajtar Repulsive Particles" begin
+    @testset "`show`" begin
         boundary_model = BoundaryModelMonaghanKajtar(10.0, 3.0, 0.1, [1.0])
 
         show_compact = "BoundaryModelMonaghanKajtar(10.0, 3.0, NoViscosity)"
         @test repr(boundary_model) == show_compact
     end
+
+    @testset "RHS" begin
+        particle_spacing = 0.1
+
+        # The state equation is only needed to unpack `sound_speed`, so we can mock
+        # it by using a `NamedTuple`.
+        state_equation = (; sound_speed=0.0)
+        smoothing_kernel = SchoenbergCubicSplineKernel{2}()
+        smoothing_length = 1.2particle_spacing
+        search_radius = TrixiParticles.compact_support(smoothing_kernel, smoothing_length)
+
+        # 3x3 fluid particles to the left of a 1x3 vertical wall, with the rightmost
+        # fluid particle one `particle_spacing` away from the boundary (which is at x=0)
+        fluid = rectangular_patch(particle_spacing, (3, 3), perturbation_factor=0.0,
+                                  perturbation_factor_position=0.0,
+                                  offset=(-1.5particle_spacing, 0.0))
+        fluid_system = WeaklyCompressibleSPHSystem(fluid, ContinuityDensity(),
+                                                   state_equation, smoothing_kernel,
+                                                   smoothing_length)
+
+        # Use double spacing for the boundary (exactly the opposite of what we would do
+        # in a simulation) to test that forces grow infinitely when a fluid particle
+        # comes too close to the boundary, independent of the boundary particle spacing.
+        boundary = rectangular_patch(2particle_spacing, (1, 3), perturbation_factor=0.0,
+                                     perturbation_factor_position=0.0)
+        K = 1.0
+        spacing_ratio = 0.5
+        boundary_model = BoundaryModelMonaghanKajtar(K, spacing_ratio, 2particle_spacing,
+                                                     boundary.mass)
+        boundary_system = BoundarySPHSystem(boundary, boundary_model)
+
+        # Density is integrated with `ContinuityDensity`
+        v = vcat(fluid.velocity, fluid.density')
+        u = fluid.coordinates
+
+        v_neighbor = zeros(0, TrixiParticles.nparticles(boundary_system))
+        u_neighbor = boundary.coordinates
+
+        nhs = TrixiParticles.TrivialNeighborhoodSearch{2}(search_radius,
+                                                          TrixiParticles.eachparticle(boundary_system))
+
+        # Result
+        dv = zero(fluid.velocity)
+        TrixiParticles.interact!(dv, v, u, v_neighbor, u_neighbor, nhs, fluid_system,
+                                 boundary_system)
+
+        # Due to the symmetric setup, all particles will only be accelerated horizontally
+        @test isapprox(dv[2, :], zeros(TrixiParticles.nparticles(fluid_system)), atol=1e-14)
+
+        # For the leftmost column of fluid particles, the boundary particles are outside the
+        # compact support of the kernel.
+        @test iszero(dv[:, [1, 4, 7]])
+
+        # The rightmost column of fluid particles should experience strong accelerations
+        # towards the left.
+        @test all(dv[1, [3, 6, 9]] .< -300)
+
+        # The middle column of fluid particles should experience weaker accelerations
+        @test isapprox(dv[1, [2, 5, 8]], [-26.052449, -95.162888, -26.052449])
+    end
 end