Make dv_shifting, surface_tension_force and adhesion_force inplace (#1124)

* Make `dv_shifting`, `surface_tension_force` and `adhesion_force` inplace

* Move free surface correction into the function

* Fix

---------

Co-authored-by: Sven Berger <berger.sven@gmail.com>
Co-authored-by: Niklas Neher <73897120+LasNikas@users.noreply.github.com>

Author: 3 people

src/io/write_vtk.jl

@@ -350,18 +350,13 @@ function write2vtk!(vtk, v, u, t, system::AbstractFluidSystem)
             rho_b = current_density(v, system, neighbor)
             grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle)
 
-            surface_tension[1:ndims(system),
-                            particle] .+= surface_tension_force(surface_tension_a,
-                                                                surface_tension_b,
-                                                                system,
-                                                                system,
-                                                                particle,
-                                                                neighbor,
-                                                                pos_diff,
-                                                                distance,
-                                                                rho_a,
-                                                                rho_b,
-                                                                grad_kernel)
+            dv_surface_tension = Ref(zero(pos_diff))
+            surface_tension_force!(dv_surface_tension,
+                                   surface_tension_a, surface_tension_b,
+                                   system, system, particle, neighbor,
+                                   pos_diff, distance, rho_a, rho_b, grad_kernel, 1)
+
+            surface_tension[1:ndims(system), particle] .+= dv_surface_tension[]
         end
         vtk["surface_tension"] = surface_tension
 

src/schemes/fluid/entropically_damped_sph/rhs.jl

@@ -52,26 +52,27 @@ function interact!(dv, v_particle_system, u_particle_system,
                                                sound_speed, m_a, m_b, rho_a, rho_b,
                                                grad_kernel)
 
-        # Extra terms in the momentum equation when using a shifting technique
-        dv_tvf = @inbounds dv_shifting(shifting_technique(particle_system),
-                                       particle_system, neighbor_system,
-                                       v_particle_system, v_neighbor_system,
-                                       particle, neighbor, m_a, m_b, rho_a, rho_b,
-                                       pos_diff, distance, grad_kernel, correction)
-
-        dv_surface_tension = surface_tension_force(surface_tension_a, surface_tension_b,
-                                                   particle_system, neighbor_system,
-                                                   particle, neighbor, pos_diff, distance,
-                                                   rho_a, rho_b, grad_kernel)
+        dv_particle = Ref(dv_pressure + dv_viscosity_)
 
-        dv_adhesion = adhesion_force(surface_tension_a, particle_system, neighbor_system,
-                                     particle, neighbor, pos_diff, distance)
-
-        dv_particle = dv_pressure + dv_viscosity_ + dv_tvf + dv_surface_tension +
-                      dv_adhesion
+        # Extra terms in the momentum equation when using a shifting technique
+        @inbounds dv_shifting!(dv_particle, shifting_technique(particle_system),
+                               particle_system, neighbor_system,
+                               v_particle_system, v_neighbor_system,
+                               particle, neighbor, m_a, m_b, rho_a, rho_b,
+                               pos_diff, distance, grad_kernel, correction)
+
+        @inbounds surface_tension_force!(dv_particle, surface_tension_a,
+                                         surface_tension_b,
+                                         particle_system, neighbor_system,
+                                         particle, neighbor, pos_diff, distance,
+                                         rho_a, rho_b, grad_kernel, 1)
+
+        @inbounds adhesion_force!(dv_particle, surface_tension_a, particle_system,
+                                  neighbor_system,
+                                  particle, neighbor, pos_diff, distance)
 
         for i in 1:ndims(particle_system)
-            @inbounds dv[i, particle] += dv_particle[i]
+            @inbounds dv[i, particle] += dv_particle[][i]
         end
 
         v_diff = current_velocity(v_particle_system, particle_system, particle) -

src/schemes/fluid/shifting_techniques.jl

@@ -41,11 +41,11 @@ function update_shifting!(system, shifting, v, u, v_ode, u_ode, semi)
 end
 
 # Additional term in the momentum equation due to the shifting technique
-@inline function dv_shifting(shifting, system, neighbor_system,
-                             v_system, v_neighbor_system, particle, neighbor,
-                             m_a, m_b, rho_a, rho_b, pos_diff, distance,
-                             grad_kernel, correction)
-    return zero(grad_kernel)
+@inline function dv_shifting!(dv_particle, shifting, system, neighbor_system,
+                              v_system, v_neighbor_system, particle, neighbor,
+                              m_a, m_b, rho_a, rho_b, pos_diff, distance,
+                              grad_kernel, correction)
+    return dv_particle
 end
 
 # Additional term(s) in the continuity equation due to the shifting technique
@@ -324,31 +324,35 @@ See [`ParticleShiftingTechnique`](@ref).
 struct MomentumEquationTermSun2019 end
 
 # Additional term in the momentum equation due to the shifting technique
-@propagate_inbounds function dv_shifting(shifting::ParticleShiftingTechnique, system,
-                                         neighbor_system,
-                                         v_system, v_neighbor_system, particle, neighbor,
-                                         m_a, m_b, rho_a, rho_b, pos_diff, distance,
-                                         grad_kernel, correction)
-    return dv_shifting(shifting.momentum_equation_term, system, neighbor_system,
-                       v_system, v_neighbor_system, particle, neighbor,
-                       m_a, m_b, rho_a, rho_b, pos_diff, distance,
-                       grad_kernel, correction)
+@propagate_inbounds function dv_shifting!(dv_particle,
+                                          shifting::ParticleShiftingTechnique,
+                                          system, neighbor_system,
+                                          v_system, v_neighbor_system, particle, neighbor,
+                                          m_a, m_b, rho_a, rho_b, pos_diff, distance,
+                                          grad_kernel, correction)
+    return dv_shifting!(dv_particle, shifting.momentum_equation_term, system,
+                        neighbor_system, v_system, v_neighbor_system,
+                        particle, neighbor, m_a, m_b, rho_a, rho_b,
+                        pos_diff, distance, grad_kernel, correction)
 end
 
-@propagate_inbounds function dv_shifting(::MomentumEquationTermSun2019,
-                                         system, neighbor_system,
-                                         v_system, v_neighbor_system,
-                                         particle, neighbor, m_a, m_b, rho_a, rho_b,
-                                         pos_diff, distance, grad_kernel, correction)
+@propagate_inbounds function dv_shifting!(dv_particle, ::MomentumEquationTermSun2019,
+                                          system, neighbor_system,
+                                          v_system, v_neighbor_system,
+                                          particle, neighbor, m_a, m_b, rho_a, rho_b,
+                                          pos_diff, distance, grad_kernel, correction)
     delta_v_a = delta_v(system, particle)
     delta_v_b = delta_v(neighbor_system, neighbor)
 
     v_a = current_velocity(v_system, system, particle)
     v_b = current_velocity(v_neighbor_system, neighbor_system, neighbor)
 
     tensor_product = v_a * delta_v_a' + v_b * delta_v_b'
-    return m_b / rho_b *
-           (tensor_product * grad_kernel + v_a * dot(delta_v_a - delta_v_b, grad_kernel))
+    dv_particle[] += m_b / rho_b *
+                     (tensor_product * grad_kernel +
+                      v_a * dot(delta_v_a - delta_v_b, grad_kernel))
+
+    return dv_particle
 end
 
 # `ParticleShiftingTechnique{<:Any, <:Any, true}` means `modify_continuity_equation=true`
@@ -568,10 +572,11 @@ struct TransportVelocityAdami{modify_continuity_equation, T <: Real} <:
     end
 end
 
-@propagate_inbounds function dv_shifting(::TransportVelocityAdami, system, neighbor_system,
-                                         v_system, v_neighbor_system, particle, neighbor,
-                                         m_a, m_b, rho_a, rho_b, pos_diff, distance,
-                                         grad_kernel, correction)
+@propagate_inbounds function dv_shifting!(dv_particle, ::TransportVelocityAdami,
+                                          system, neighbor_system,
+                                          v_system, v_neighbor_system, particle, neighbor,
+                                          m_a, m_b, rho_a, rho_b, pos_diff, distance,
+                                          grad_kernel, correction)
     v_a = current_velocity(v_system, system, particle)
     delta_v_a = delta_v(system, particle)
 
@@ -588,9 +593,11 @@ end
     #   m_b * (A_a + A_b) / (ρ_a * ρ_b) * ∇W_ab.
     # In order to obtain this, we pass `p_a = A_a` and `p_b = A_b` to the
     # `pressure_acceleration` function.
-    return pressure_acceleration(system, neighbor_system, particle, neighbor,
-                                 m_a, m_b, A_a, A_b, rho_a, rho_b, pos_diff,
-                                 distance, grad_kernel, correction)
+    dv_particle[] += pressure_acceleration(system, neighbor_system, particle, neighbor,
+                                           m_a, m_b, A_a, A_b, rho_a, rho_b,
+                                           pos_diff, distance, grad_kernel, correction)
+
+    return dv_particle
 end
 
 # The function above misuses the pressure acceleration function by passing a Matrix as `p_a`.