Speed up continuity equation

Author: efaulhaber

src/schemes/fluid/fluid.jl

@@ -136,21 +136,38 @@ end
     return dv
 end
 
-# This formulation was chosen to be consistent with the used pressure_acceleration formulations
-@propagate_inbounds function continuity_equation!(dv, density_calculator::ContinuityDensity,
+# For compatibility with the old continuity equation function
+@propagate_inbounds function continuity_equation!(dv::AbstractArray,
+                                                  density_calculator::ContinuityDensity,
                                                   particle_system::AbstractFluidSystem,
                                                   neighbor_system,
                                                   v_particle_system, v_neighbor_system,
                                                   particle, neighbor, pos_diff, distance,
                                                   m_b, rho_a, rho_b, grad_kernel)
+    drho_particle = Ref(zero(rho_a))
     vdiff = current_velocity(v_particle_system, particle_system, particle) -
             current_velocity(v_neighbor_system, neighbor_system, neighbor)
 
+    continuity_equation!(drho_particle, density_calculator, particle_system,
+                         neighbor_system, v_particle_system, v_neighbor_system,
+                         particle, neighbor, pos_diff, distance,
+                         m_b, rho_a, rho_b, vdiff, grad_kernel)
+    dv[end, particle] += drho_particle[]
+end
+
+# This formulation was chosen to be consistent with the used pressure_acceleration formulations
+@propagate_inbounds function continuity_equation!(drho_particle,
+                                                  density_calculator::ContinuityDensity,
+                                                  particle_system::AbstractFluidSystem,
+                                                  neighbor_system,
+                                                  v_particle_system, v_neighbor_system,
+                                                  particle, neighbor, pos_diff, distance,
+                                                  m_b, rho_a, rho_b, vdiff, grad_kernel)
     vdiff += continuity_equation_shifting_term(shifting_technique(particle_system),
                                                particle_system, neighbor_system,
                                                particle, neighbor, rho_a, rho_b)
 
-    dv[end, particle] += rho_a / rho_b * m_b * dot(vdiff, grad_kernel)
+    drho_particle[] += div_fast(rho_a, rho_b) * m_b * dot(vdiff, grad_kernel)
 
     # Artificial density diffusion should only be applied to systems representing a fluid
     # with the same physical properties i.e. density and viscosity.

src/schemes/fluid/weakly_compressible_sph/rhs.jl

@@ -55,6 +55,7 @@ function interact!(dv, v_particle_system, u_particle_system,
             v_b, rho_b = @inbounds velocity_and_density(v_neighbor_system, neighbor_system,
                                                         neighbor)
             rho_mean = (rho_a + rho_b) / 2
+            vdiff = v_a - v_b
 
             # The following call is equivalent to
             #     `p_b = current_pressure(v_neighbor_system, neighbor_system, neighbor)`
@@ -107,15 +108,17 @@ function interact!(dv, v_particle_system, u_particle_system,
 
             # TODO If variable smoothing_length is used, this should use the neighbor smoothing length
             # Propagate `@inbounds` to the continuity equation, which accesses particle data
-            @inbounds continuity_equation!(dv, density_calculator, particle_system,
+            @inbounds continuity_equation!(drho_particle, density_calculator, particle_system,
                                            neighbor_system, v_particle_system,
                                            v_neighbor_system, particle, neighbor,
-                                           pos_diff, distance, m_b, rho_a, rho_b, grad_kernel)
+                                           pos_diff, distance, m_b, rho_a, rho_b, vdiff,
+                                           grad_kernel)
         end
 
         for i in eachindex(dv_particle[])
             @inbounds dv[i, particle] += dv_particle[][i]
         end
+        @inbounds dv[end, particle] += drho_particle[]
     end
 
     return dv