Merge branch 'main' into fix-visualization

src/schemes/boundary/rhs.jl

@@ -12,6 +12,7 @@ function interact!(dv, v_particle_system, u_particle_system,
                    particle_system::BoundarySPHSystem{<:BoundaryModelDummyParticles{ContinuityDensity}},
                    neighbor_system::FluidSystem)
     (; boundary_model) = particle_system
+    fluid_density_calculator = neighbor_system.density_calculator
 
     system_coords = current_coordinates(u_particle_system, particle_system)
     neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system)
@@ -20,16 +21,33 @@ function interact!(dv, v_particle_system, u_particle_system,
     for_particle_neighbor(particle_system, neighbor_system,
                           system_coords, neighbor_coords,
                           neighborhood_search) do particle, neighbor, pos_diff, distance
-        # Continuity equation
-        vdiff = current_velocity(v_particle_system, particle_system, particle) -
-                current_velocity(v_neighbor_system, neighbor_system, neighbor)
-
-        # For boundary particles, the velocity is not integrated.
-        # Therefore, the density is stored in the first dimension of `dv`.
-        dv[1, particle] += sum(neighbor_system.mass[neighbor] * vdiff .*
-                               smoothing_kernel_grad(boundary_model, pos_diff,
-                                                     distance))
+        m_b = hydrodynamic_mass(neighbor_system, neighbor)
+
+        rho_a = particle_density(v_particle_system, particle_system, particle)
+        rho_b = particle_density(v_neighbor_system, neighbor_system, neighbor)
+
+        v_diff = current_velocity(v_particle_system, particle_system, particle) -
+                 current_velocity(v_neighbor_system, neighbor_system, neighbor)
+
+        grad_kernel = smoothing_kernel_grad(boundary_model, pos_diff, distance)
+
+        continuity_equation!(dv, fluid_density_calculator, m_b, rho_a, rho_b, v_diff,
+                             grad_kernel, particle)
     end
 
     return dv
 end
+
+# This is the derivative of the density summation, which is compatible with the
+# `SummationDensity` pressure acceleration.
+# Energy preservation tests will fail with the other formulation.
+function continuity_equation!(dv, fluid_density_calculator::SummationDensity,
+                              m_b, rho_a, rho_b, v_diff, grad_kernel, particle)
+    dv[end, particle] += m_b * dot(v_diff, grad_kernel)
+end
+
+# This is identical to the continuity equation of the fluid
+function continuity_equation!(dv, fluid_density_calculator::ContinuityDensity,
+                              m_b, rho_a, rho_b, v_diff, grad_kernel, particle)
+    dv[end, particle] += rho_a / rho_b * m_b * dot(v_diff, grad_kernel)
+end

src/schemes/solid/total_lagrangian_sph/rhs.jl

@@ -100,9 +100,11 @@ function interact!(dv, v_particle_system, u_particle_system,
         p_a = particle_pressure(v_particle_system, particle_system, particle)
         p_b = particle_pressure(v_neighbor_system, neighbor_system, neighbor)
 
-        # Boundary forces
-        # Note: neighbor and particle pressure are switched in this call
-        #       and `pressure_correction` is set to `1.0` (no correction)
+        # Particle and neighbor (and corresponding systems and all corresponding quantities)
+        # are switched in this call.
+        # This way, we obtain the exact same force as for the fluid-solid interaction,
+        # but with a flipped sign (because `pos_diff` is flipped compared to fluid-solid).
+        # `pressure_correction` is set to `1.0` (no correction).
         dv_boundary = pressure_acceleration(1.0, m_a, p_b, p_a,
                                             rho_b, rho_a, pos_diff, distance, grad_kernel,
                                             neighbor_system, neighbor, particle_system,
@@ -111,61 +113,44 @@ function interact!(dv, v_particle_system, u_particle_system,
 
         for i in 1:ndims(particle_system)
             # Multiply `dv` (acceleration on fluid particle b) by the mass of
-            # particle b to obtain the force.
+            # particle b to obtain the same force as for the fluid-solid interaction.
             # Divide by the material mass of particle a to obtain the acceleration
             # of solid particle a.
-            dv[i, particle] += dv_particle[i] * neighbor_system.mass[neighbor] /
-                               particle_system.mass[particle]
+            dv[i, particle] += dv_particle[i] * m_b / particle_system.mass[particle]
         end
 
         continuity_equation!(dv, v_particle_system, v_neighbor_system,
                              particle, neighbor, pos_diff, distance,
-                             particle_system, neighbor_system)
+                             m_b, rho_a, rho_b,
+                             particle_system, neighbor_system, grad_kernel)
     end
 
     return dv
 end
 
 @inline function continuity_equation!(dv, v_particle_system, v_neighbor_system,
                                       particle, neighbor, pos_diff, distance,
+                                      m_b, rho_a, rho_b,
                                       particle_system::TotalLagrangianSPHSystem,
-                                      neighbor_system::WeaklyCompressibleSPHSystem)
+                                      neighbor_system::WeaklyCompressibleSPHSystem,
+                                      grad_kernel)
     return dv
 end
 
 @inline function continuity_equation!(dv, v_particle_system, v_neighbor_system,
                                       particle, neighbor, pos_diff, distance,
-                                      particle_system::TotalLagrangianSPHSystem{<:BoundaryModelDummyParticles},
-                                      neighbor_system::WeaklyCompressibleSPHSystem)
-    (; density_calculator) = particle_system.boundary_model
-
-    continuity_equation!(dv, density_calculator,
-                         v_particle_system, v_neighbor_system,
-                         particle, neighbor, pos_diff, distance,
-                         particle_system, neighbor_system)
-end
-
-@inline function continuity_equation!(dv, density_calculator,
-                                      u_particle_system, u_neighbor_system,
-                                      particle, neighbor, pos_diff, distance,
-                                      particle_system::TotalLagrangianSPHSystem,
-                                      neighbor_system::WeaklyCompressibleSPHSystem)
-    return dv
-end
-
-@inline function continuity_equation!(dv, ::ContinuityDensity,
-                                      v_particle_system, v_neighbor_system,
-                                      particle, neighbor, pos_diff, distance,
-                                      particle_system::TotalLagrangianSPHSystem,
-                                      neighbor_system::WeaklyCompressibleSPHSystem)
-    vdiff = current_velocity(v_particle_system, particle_system, particle) -
-            current_velocity(v_neighbor_system, neighbor_system, neighbor)
-
-    dv[end, particle] += sum(neighbor_system.mass[neighbor] * vdiff .*
-                             smoothing_kernel_grad(neighbor_system, pos_diff,
-                                                   distance))
-
-    return dv
+                                      m_b, rho_a, rho_b,
+                                      particle_system::TotalLagrangianSPHSystem{<:BoundaryModelDummyParticles{ContinuityDensity}},
+                                      neighbor_system::WeaklyCompressibleSPHSystem,
+                                      grad_kernel)
+    fluid_density_calculator = neighbor_system.density_calculator
+
+    v_diff = current_velocity(v_particle_system, particle_system, particle) -
+             current_velocity(v_neighbor_system, neighbor_system, neighbor)
+
+    # Call the dummy BC version of the continuity equation
+    continuity_equation!(dv, fluid_density_calculator, m_b, rho_a, rho_b, v_diff,
+                         grad_kernel, particle)
 end
 
 # Solid-boundary interaction

test/schemes/boundary/dummy_particles/dummy_particles.jl

@@ -1,4 +1,3 @@
-
 @testset verbose=true "Dummy Particles" begin
     @testset "show" begin
         boundary_model = BoundaryModelDummyParticles([1000.0], [1.0],
@@ -136,3 +135,5 @@
         end
     end
 end
+
+include("rhs.jl")