Rename compute_smoothed_velocity! → interpolate_fluid_velocity! and fix dispatch

NEWS.md

@@ -25,6 +25,7 @@ used in the Julia ecosystem. Notable changes will be documented in this file for
 
 - Fixed the periodic array of cylinders example file (#975).
 - A `StepsizeCallback` can now be used with open boundaries (#1074).
+- Fixed a bug with no-slip boundary conditions when using any viscosity model other than `ViscosityAdami` (#1089).
 
 ### Documentation
 

src/schemes/boundary/wall_boundary/dummy_particles.jl

@@ -277,12 +277,13 @@ function create_cache_model(viscosity, n_particles, n_dims)
     return (; wall_velocity)
 end
 
-@inline reset_cache!(cache, viscosity) = set_zero!(cache.volume)
+@inline reset_cache!(cache, viscosity::Nothing) = set_zero!(cache.volume)
 
-function reset_cache!(cache, viscosity::ViscosityAdami)
+function reset_cache!(cache, viscosity)
     (; volume, wall_velocity) = cache
 
     set_zero!(volume)
+    # Reset the accumulated fluid velocity interpolation used in `interpolate_fluid_velocity!`
     set_zero!(wall_velocity)
 
     return cache
@@ -461,7 +462,7 @@ function compute_pressure!(boundary_model,
 
     set_zero!(pressure)
 
-    # Set `volume` to zero. For `ViscosityAdami` the `wall_velocity` is also set to zero.
+    # Set `volume` to zero. When using a viscosity model, the `wall_velocity` is also set to zero.
     reset_cache!(cache, viscosity)
 
     system_coords = current_coordinates(u, system)
@@ -625,8 +626,8 @@ end
     pressure[particle] += sum_pressures * kernel_weight
     cache.volume[particle] += kernel_weight
 
-    compute_smoothed_velocity!(cache, viscosity, neighbor_system, v_neighbor_system,
-                               kernel_weight, particle, neighbor)
+    interpolate_fluid_velocity!(cache, viscosity, neighbor_system, v_neighbor_system,
+                                kernel_weight, particle, neighbor)
 end
 
 @inline function dynamic_pressure(boundary_density_calculator, v,
@@ -664,15 +665,22 @@ end
            dot(normal_velocity, normal_velocity) / 2
 end
 
-@inline function compute_smoothed_velocity!(cache, viscosity, neighbor_system,
-                                            v_neighbor_system, kernel_weight,
-                                            particle, neighbor)
+# No-op when no viscosity model is set
+@inline function interpolate_fluid_velocity!(cache, viscosity::Nothing, neighbor_system,
+                                             v_neighbor_system, kernel_weight,
+                                             particle, neighbor)
     return cache
 end
 
-@propagate_inbounds function compute_smoothed_velocity!(cache, viscosity::ViscosityAdami,
-                                                        neighbor_system, v_neighbor_system,
-                                                        kernel_weight, particle, neighbor)
+# For any viscosity model, interpolate the fluid velocity to the boundary particle location.
+# This is used later in `compute_wall_velocity!` to impose a no-slip boundary condition
+# by computing the wall velocity as v_w = 2 * v_boundary - ṽ_fluid, where ṽ_fluid is the
+# kernel-weighted interpolation of the surrounding fluid velocities:
+#   ṽ_fluid = (Σ_b v_b W_ab) / (Σ_b W_ab)
+# The denominator Σ_b W_ab is stored as `volume` and accumulated elsewhere.
+@propagate_inbounds function interpolate_fluid_velocity!(cache, viscosity,
+                                                         neighbor_system, v_neighbor_system,
+                                                         kernel_weight, particle, neighbor)
     v_b = current_velocity(v_neighbor_system, neighbor_system, neighbor)
 
     for dim in eachindex(v_b)
@@ -696,7 +704,9 @@ end
     v_boundary = current_velocity(v, system, particle)
 
     for dim in eachindex(v_boundary)
-        # The second term is the precalculated smoothed velocity field of the fluid
+        # Compute the no-slip wall velocity using the formula v_w = 2 * v_a - ṽ_fluid,
+        # where ṽ_fluid = wall_velocity / volume is the kernel-weighted fluid velocity
+        # interpolated to the boundary particle location by `interpolate_fluid_velocity!`.
         new_velocity = @inbounds 2 * v_boundary[dim] -
                                  wall_velocity[dim, particle] / volume[particle]
         @inbounds wall_velocity[dim, particle] = new_velocity