Make viscosity inplace and avoid duplicate velocity loads

src/general/interpolation.jl

@@ -593,7 +593,9 @@ end
                 other_density[point] += m_b * W_ab
 
                 if include_wall_velocity
-                    velocity_neighbor = viscous_velocity(v, neighbor_system, neighbor)
+                    velocity_neighbor_ = current_velocity(v, neighbor_system, neighbor)
+                    velocity_neighbor = viscous_velocity(v, neighbor_system, neighbor,
+                                                         velocity_neighbor_)
                     for i in axes(velocity_neighbor, 1)
                         cache.velocity[i, point] += velocity_neighbor[i] * volume_b * W_ab
                     end

src/schemes/boundary/open_boundary/rhs.jl

@@ -41,6 +41,9 @@ function interact!(dv, v_particle_system, u_particle_system,
         rho_a = @inbounds current_density(v_particle_system, particle_system, particle)
         rho_b = @inbounds current_density(v_neighbor_system, neighbor_system, neighbor)
 
+        v_a = @inbounds current_velocity(v_particle_system, particle_system, particle)
+        v_b = @inbounds current_velocity(v_neighbor_system, neighbor_system, neighbor)
+
         m_a = @inbounds hydrodynamic_mass(particle_system, particle)
         m_b = @inbounds hydrodynamic_mass(neighbor_system, neighbor)
 
@@ -58,15 +61,17 @@ function interact!(dv, v_particle_system, u_particle_system,
         dv_pressure_boundary = 2 * p_boundary * (m_b / (rho_a * rho_b)) * grad_kernel
 
         # Propagate `@inbounds` to the viscosity function, which accesses particle data
-        dv_viscosity_ = @inbounds dv_viscosity(viscosity_model(fluid_system,
-                                                               neighbor_system),
-                                               particle_system, neighbor_system,
-                                               v_particle_system, v_neighbor_system,
-                                               particle, neighbor, pos_diff, distance,
-                                               sound_speed, m_a, m_b, rho_a, rho_b,
-                                               grad_kernel)
-
-        dv_particle = dv_pressure + dv_viscosity_ + dv_pressure_boundary
+        dv_viscosity_ = Ref(zero(pos_diff))
+        @inbounds dv_viscosity!(dv_viscosity_,
+                                viscosity_model(fluid_system,
+                                                neighbor_system),
+                                particle_system, neighbor_system,
+                                v_particle_system, v_neighbor_system,
+                                particle, neighbor, pos_diff, distance,
+                                sound_speed, m_a, m_b, rho_a, rho_b,
+                                v_a, v_b, grad_kernel)
+
+        dv_particle = dv_pressure + dv_viscosity_[] + dv_pressure_boundary
 
         for i in 1:ndims(particle_system)
             @inbounds dv[i, particle] += dv_particle[i]

src/schemes/boundary/open_boundary/system.jl

@@ -654,8 +654,10 @@ function interpolate_velocity!(system::OpenBoundarySystem, boundary_zone,
             W_ab = kernel(smoothing_kernel, distance, smoothing_length)
             @inbounds shepard_coefficient[point] += volume_b * W_ab
 
+            velocity_neighbor_ = @inbounds current_velocity(v_neighbor, neighbor_system,
+                                                            neighbor)
             velocity_neighbor = @inbounds viscous_velocity(v_neighbor, neighbor_system,
-                                                           neighbor)
+                                                           neighbor, velocity_neighbor_)
             for i in axes(velocity_neighbor, 1)
                 @inbounds sample_velocity[i,
                                           point] += velocity_neighbor[i] * volume_b * W_ab

src/schemes/boundary/wall_boundary/system.jl

@@ -145,15 +145,20 @@ end
     return zero(SVector{ndims(system), eltype(system)})
 end
 
-@inline function viscous_velocity(v, system::WallBoundarySystem, particle)
-    return viscous_velocity(v, system.boundary_model.viscosity, system, particle)
+@propagate_inbounds function viscous_velocity(v, system::WallBoundarySystem,
+                                              particle, v_particle)
+    return viscous_velocity(v, system.boundary_model.viscosity, system,
+                            particle, v_particle)
 end
 
-@inline function viscous_velocity(v, ::Nothing, system, particle)
-    return current_velocity(v, system, particle)
+@inline function viscous_velocity(v, ::Nothing, system, particle, v_particle)
+    # Regular particle velocity is used for the viscosity calculation by default
+    return v_particle
 end
 
-@inline function viscous_velocity(v, viscosity, system, particle)
+@propagate_inbounds function viscous_velocity(v, viscosity, system, particle, v_particle)
+    # Wall velocity in the viscosity calculation contains the physical wall velocity
+    # and an interpolated velocity when a wall viscosity (no-slip BC) is used.
     return extract_svector(system.boundary_model.cache.wall_velocity, system, particle)
 end
 

src/schemes/fluid/entropically_damped_sph/rhs.jl

@@ -42,6 +42,9 @@ function interact!(dv, v_particle_system, u_particle_system,
         rho_a = @inbounds current_density(v_particle_system, particle_system, particle)
         rho_b = @inbounds current_density(v_neighbor_system, neighbor_system, neighbor)
 
+        v_a = @inbounds current_velocity(v_particle_system, particle_system, particle)
+        v_b = @inbounds current_velocity(v_neighbor_system, neighbor_system, neighbor)
+
         p_a = @inbounds current_pressure(v_particle_system, particle_system, particle)
         p_b = @inbounds current_pressure(v_neighbor_system, neighbor_system, neighbor)
 
@@ -62,13 +65,12 @@ function interact!(dv, v_particle_system, u_particle_system,
                                             rho_b, pos_diff, distance, grad_kernel,
                                             correction)
 
-        dv_viscosity_ = @inbounds dv_viscosity(particle_system, neighbor_system,
-                                               v_particle_system, v_neighbor_system,
-                                               particle, neighbor, pos_diff, distance,
-                                               sound_speed, m_a, m_b, rho_a, rho_b,
-                                               grad_kernel)
-
-        dv_particle = Ref(dv_pressure + dv_viscosity_)
+        dv_particle = Ref(dv_pressure)
+        @inbounds dv_viscosity!(dv_particle, particle_system, neighbor_system,
+                                v_particle_system, v_neighbor_system,
+                                particle, neighbor, pos_diff, distance,
+                                sound_speed, m_a, m_b, rho_a, rho_b,
+                                v_a, v_b, grad_kernel)
 
         # Extra terms in the momentum equation when using a shifting technique
         @inbounds dv_shifting!(dv_particle, shifting_technique(particle_system),

src/schemes/fluid/implicit_incompressible_sph/rhs.jl

@@ -41,6 +41,9 @@ function interact!(dv, v_particle_system, u_particle_system,
         rho_a = @inbounds current_density(v_particle_system, particle_system, particle)
         rho_b = @inbounds current_density(v_neighbor_system, neighbor_system, neighbor)
 
+        v_a = @inbounds current_velocity(v_particle_system, particle_system, particle)
+        v_b = @inbounds current_velocity(v_neighbor_system, neighbor_system, neighbor)
+
         m_a = @inbounds hydrodynamic_mass(particle_system, particle)
         m_b = @inbounds hydrodynamic_mass(neighbor_system, neighbor)
 
@@ -62,14 +65,15 @@ function interact!(dv, v_particle_system, u_particle_system,
                                             distance, grad_kernel, nothing)
 
         # Propagate `@inbounds` to the viscosity function, which accesses particle data
-        dv_viscosity_ = @inbounds dv_viscosity(particle_system, neighbor_system,
-                                               v_particle_system, v_neighbor_system,
-                                               particle, neighbor, pos_diff, distance,
-                                               sound_speed, m_a, m_b, rho_a, rho_b,
-                                               grad_kernel)
+        dv_viscosity_ = Ref(zero(pos_diff))
+        @inbounds dv_viscosity!(dv_viscosity_, particle_system, neighbor_system,
+                                v_particle_system, v_neighbor_system,
+                                particle, neighbor, pos_diff, distance,
+                                sound_speed, m_a, m_b, rho_a, rho_b,
+                                v_a, v_b, grad_kernel)
 
         for i in 1:ndims(particle_system)
-            @inbounds dv[i, particle] += dv_pressure[i] + dv_viscosity_[i]
+            @inbounds dv[i, particle] += dv_pressure[i] + dv_viscosity_[][i]
         end
     end
     return dv

src/schemes/fluid/implicit_incompressible_sph/system.jl

@@ -281,16 +281,21 @@ function calculate_predicted_velocity_and_d_ii_values!(system::ImplicitIncompres
             rho_a = @inbounds current_density(v_particle_system, system, particle)
             rho_b = @inbounds current_density(v_neighbor_system, neighbor_system, neighbor)
 
+            v_a = @inbounds current_velocity(v_particle_system, system, particle)
+            v_b = @inbounds current_velocity(v_neighbor_system, neighbor_system, neighbor)
+
             grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle)
 
-            dv_viscosity_ = @inbounds dv_viscosity(system, neighbor_system,
-                                                   v_particle_system, v_neighbor_system,
-                                                   particle, neighbor, pos_diff, distance,
-                                                   sound_speed, m_a, m_b, rho_a, rho_b,
-                                                   grad_kernel)
+            dv_viscosity_ = Ref(zero(pos_diff))
+            @inbounds dv_viscosity!(dv_viscosity_, system, neighbor_system,
+                                    v_particle_system, v_neighbor_system,
+                                    particle, neighbor, pos_diff, distance,
+                                    sound_speed, m_a, m_b, rho_a, rho_b,
+                                    v_a, v_b, grad_kernel)
             # Add all other non-pressure forces
             for i in 1:ndims(system)
-                @inbounds advection_velocity[i, particle] += time_step * dv_viscosity_[i]
+                @inbounds advection_velocity[i,
+                                             particle] += time_step * dv_viscosity_[][i]
             end
             # Calculate d_ii with eq. 9 in Ihmsen et al. (2013)
             for i in 1:ndims(system)