Make continuity_equation! and density_diffusion! work on Ref values

src/schemes/boundary/open_boundary/rhs.jl

@@ -62,10 +62,11 @@ function interact!(dv, v_particle_system, u_particle_system,
         # Continuity equation
         @inbounds dv[end, particle] += rho_a / rho_b * m_b * dot(v_diff, grad_kernel)
 
-        density_diffusion!(dv, density_diffusion(particle_system),
-                           v_particle_system, particle, neighbor,
-                           pos_diff, distance, m_b, rho_a, rho_b,
-                           particle_system, grad_kernel)
+        drho_particle = Ref(zero(rho_a))
+        density_diffusion!(drho_particle, density_diffusion(particle_system),
+                           particle_system, particle, neighbor,
+                           pos_diff, distance, m_b, rho_a, rho_b, grad_kernel)
+        @inbounds dv[end, particle] += drho_particle[]
 
         # Open boundary pressure evolution matches the corresponding fluid system:
         # - EDAC: Compute pressure evolution like the fluid system

src/schemes/boundary/wall_boundary/rhs.jl

@@ -27,13 +27,16 @@ function interact!(dv, v_particle_system, u_particle_system,
         rho_a = current_density(v_particle_system, particle_system, particle)
         rho_b = current_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)
+        v_a = current_velocity(v_particle_system, particle_system, particle)
+        v_b = current_velocity(v_neighbor_system, neighbor_system, neighbor)
 
         grad_kernel = smoothing_kernel_grad(boundary_model, pos_diff, distance, particle)
 
-        continuity_equation!(dv, density_calculator(neighbor_system),
-                             m_b, rho_a, rho_b, v_diff, grad_kernel, particle)
+        drho_particle = Ref(zero(rho_a))
+        continuity_equation!(drho_particle, density_calculator(neighbor_system),
+                             m_b, rho_a, rho_b, v_a, v_b, grad_kernel, particle)
+
+        dv[end, particle] += drho_particle[]
     end
 
     return dv
@@ -42,13 +45,17 @@ 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)
+function continuity_equation!(drho_particle, ::SummationDensity,
+                              m_b, rho_a, rho_b, v_a, v_b, grad_kernel, particle)
+    drho_particle[] += m_b * dot(v_a - v_b, grad_kernel)
+
+    return drho_particle
 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)
+function continuity_equation!(drho_particle, ::ContinuityDensity,
+                              m_b, rho_a, rho_b, v_a, v_b, grad_kernel, particle)
+    drho_particle[] += rho_a / rho_b * m_b * dot(v_a - v_b, grad_kernel)
+
+    return drho_particle
 end

src/schemes/fluid/entropically_damped_sph/rhs.jl

@@ -74,19 +74,25 @@ function interact!(dv, v_particle_system, u_particle_system,
             @inbounds dv[i, particle] += dv_particle[i]
         end
 
-        v_diff = current_velocity(v_particle_system, particle_system, particle) -
-                 current_velocity(v_neighbor_system, neighbor_system, neighbor)
+        v_a = current_velocity(v_particle_system, particle_system, particle)
+        v_b = current_velocity(v_neighbor_system, neighbor_system, neighbor)
+        v_diff = v_a - v_b
 
         pressure_evolution!(dv, particle_system, neighbor_system, v_diff, grad_kernel,
                             particle, neighbor, pos_diff, distance,
                             sound_speed, m_a, m_b, p_a, p_b, rho_a, rho_b, nu_edac)
 
+        drho_particle = Ref(zero(rho_a))
+
         # 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,
+                                       v_a, v_b, grad_kernel)
+
+        @inbounds write_drho_particle!(dv, density_calculator, drho_particle, particle)
     end
 
     return dv

src/schemes/fluid/fluid.jl

@@ -128,39 +128,53 @@ function compute_density!(system, u, u_ode, semi, ::SummationDensity)
 end
 
 # With 'SummationDensity', density is calculated in wcsph/system.jl:compute_density!
-@inline function continuity_equation!(dv, density_calculator::SummationDensity,
+@inline function continuity_equation!(drho_particle, ::SummationDensity,
                                       particle_system, neighbor_system,
                                       v_particle_system, v_neighbor_system,
                                       particle, neighbor, pos_diff, distance,
-                                      m_b, rho_a, rho_b, grad_kernel)
-    return dv
+                                      m_b, rho_a, rho_b, v_a, v_b, grad_kernel)
+    return drho_particle
 end
 
 # This formulation was chosen to be consistent with the used pressure_acceleration formulations
-@propagate_inbounds function continuity_equation!(dv, density_calculator::ContinuityDensity,
+@propagate_inbounds function continuity_equation!(drho_particle,
+                                                  ::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)
-    vdiff = current_velocity(v_particle_system, particle_system, particle) -
-            current_velocity(v_neighbor_system, neighbor_system, neighbor)
+                                                  m_b, rho_a, rho_b, v_a, v_b,
+                                                  grad_kernel)
+    v_diff = v_a - v_b
 
-    vdiff += continuity_equation_shifting_term(shifting_technique(particle_system),
+    v_diff = continuity_equation_shifting_term(v_diff,
+                                               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[] += rho_a / rho_b * m_b * dot(v_diff, grad_kernel)
 
     # Artificial density diffusion should only be applied to systems representing a fluid
     # with the same physical properties i.e. density and viscosity.
     # TODO: shouldn't be applied to particles on the interface (depends on PR #539)
     if particle_system === neighbor_system
-        density_diffusion!(dv, density_diffusion(particle_system),
-                           v_particle_system, particle, neighbor,
-                           pos_diff, distance, m_b, rho_a, rho_b, particle_system,
-                           grad_kernel)
+        density_diffusion!(drho_particle, density_diffusion(particle_system),
+                           particle_system, particle, neighbor,
+                           pos_diff, distance, m_b, rho_a, rho_b, grad_kernel)
     end
+
+    return drho_particle
+end
+
+@inline function write_drho_particle!(dv, density_calculator, drho_particle, particle)
+    return dv
+end
+
+@propagate_inbounds function write_drho_particle!(dv, ::ContinuityDensity,
+                                                  drho_particle, particle)
+    dv[end, particle] += drho_particle[]
+
+    return dv
 end
 
 function calculate_dt(v_ode, u_ode, cfl_number, system::AbstractFluidSystem, semi)

src/schemes/fluid/shifting_techniques.jl

@@ -48,11 +48,12 @@ end
     return zero(grad_kernel)
 end
 
-# Additional term(s) in the continuity equation due to the shifting technique
-@inline function continuity_equation_shifting_term(shifting, particle_system,
+# Add additional term(s) in the continuity equation due to the shifting technique
+# and return the modified term.
+@inline function continuity_equation_shifting_term(v_diff, shifting, particle_system,
                                                    neighbor_system,
                                                    particle, neighbor, rho_a, rho_b)
-    return zero(SVector{ndims(particle_system), eltype(particle_system)})
+    return v_diff
 end
 
 @doc raw"""
@@ -352,7 +353,8 @@ end
 end
 
 # `ParticleShiftingTechnique{<:Any, <:Any, true}` means `modify_continuity_equation=true`
-@propagate_inbounds function continuity_equation_shifting_term(shifting::ParticleShiftingTechnique{<:Any,
+@propagate_inbounds function continuity_equation_shifting_term(v_diff,
+                                                               shifting::ParticleShiftingTechnique{<:Any,
                                                                                                    <:Any,
                                                                                                    true},
                                                                system, neighbor_system,
@@ -362,9 +364,10 @@ end
     delta_v_b = delta_v(neighbor_system, neighbor)
     delta_v_diff = delta_v_a - delta_v_b
 
-    second_term = second_continuity_equation_term(shifting.second_continuity_equation_term,
-                                                  delta_v_a, delta_v_b, rho_a, rho_b)
-    return delta_v_diff + second_term
+    shifting_term = second_continuity_equation_term(delta_v_diff,
+                                                    shifting.second_continuity_equation_term,
+                                                    delta_v_a, delta_v_b, rho_a, rho_b)
+    return v_diff + shifting_term
 end
 
 """
@@ -377,15 +380,16 @@ See [`ParticleShiftingTechnique`](@ref).
 """
 struct ContinuityEquationTermSun2019 end
 
-@propagate_inbounds function second_continuity_equation_term(::ContinuityEquationTermSun2019,
+@propagate_inbounds function second_continuity_equation_term(v_diff,
+                                                             ::ContinuityEquationTermSun2019,
                                                              delta_v_a, delta_v_b,
                                                              rho_a, rho_b)
-    return delta_v_a + rho_b / rho_a * delta_v_b
+    return v_diff + delta_v_a + rho_b / rho_a * delta_v_b
 end
 
-@inline function second_continuity_equation_term(second_continuity_equation_term,
+@inline function second_continuity_equation_term(v_diff, second_continuity_equation_term,
                                                  delta_v_a, delta_v_b, rho_a, rho_b)
-    return zero(delta_v_a)
+    return v_diff
 end
 
 # `ParticleShiftingTechnique{<:Any, true}` means `update_everystage=true`
@@ -599,15 +603,16 @@ end
     return pressure_acceleration_continuity_density(m_a, m_b, rho_a, rho_b, p_a, p_b, W_a)
 end
 
-@propagate_inbounds function continuity_equation_shifting_term(::TransportVelocityAdami{true},
+@propagate_inbounds function continuity_equation_shifting_term(v_diff,
+                                                               ::TransportVelocityAdami{true},
                                                                particle_system,
                                                                neighbor_system,
                                                                particle, neighbor,
                                                                rho_a, rho_b)
     delta_v_diff = delta_v(particle_system, particle) -
                    delta_v(neighbor_system, neighbor)
 
-    return delta_v_diff
+    return v_diff + delta_v_diff
 end
 
 function update_shifting!(system, shifting::TransportVelocityAdami, v, u, v_ode,

src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl

@@ -206,13 +206,12 @@ function update!(density_diffusion::DensityDiffusionAntuono, v, u, system, semi)
     return density_diffusion
 end
 
-@propagate_inbounds function density_diffusion!(dv,
+@propagate_inbounds function density_diffusion!(drho_particle,
                                                 density_diffusion::AbstractDensityDiffusion,
-                                                v_particle_system, particle, neighbor,
-                                                pos_diff, distance, m_b, rho_a, rho_b,
                                                 particle_system::Union{AbstractFluidSystem,
                                                                        OpenBoundarySystem{<:BoundaryModelDynamicalPressureZhang}},
-                                                grad_kernel)
+                                                particle, neighbor, pos_diff, distance,
+                                                m_b, rho_a, rho_b, grad_kernel)
     # Density diffusion terms are all zero for distance zero.
     # See `src/general/smoothing_kernels.jl` for more details.
     distance^2 < eps(initial_smoothing_length(particle_system)^2) && return
@@ -226,15 +225,13 @@ end
                                 particle_system, particle, neighbor)
     density_diffusion_term = dot(psi, grad_kernel) * volume_b
 
-    smoothing_length_avg = (smoothing_length(particle_system, particle) +
-                            smoothing_length(particle_system, neighbor)) / 2
-    dv[end, particle] += delta * smoothing_length_avg * sound_speed *
-                         density_diffusion_term
+    h = smoothing_length(particle_system, particle)
+    drho_particle[] += delta * h * sound_speed * density_diffusion_term
 end
 
 # Density diffusion `nothing` or interaction other than fluid-fluid
-@inline function density_diffusion!(dv, density_diffusion, v_particle_system, particle,
-                                    neighbor, pos_diff, distance, m_b, rho_a, rho_b,
-                                    particle_system, grad_kernel)
-    return dv
+@inline function density_diffusion!(drho_particle, density_diffusion, particle_system,
+                                    particle, neighbor, pos_diff, distance,
+                                    m_b, rho_a, rho_b, grad_kernel)
+    return drho_particle
 end

src/schemes/fluid/weakly_compressible_sph/rhs.jl

@@ -33,6 +33,9 @@ function interact!(dv, v_particle_system, u_particle_system,
         rho_b = @inbounds current_density(v_neighbor_system, neighbor_system, neighbor)
         rho_mean = (rho_a + rho_b) / 2
 
+        v_a = @inbounds current_velocity(v_particle_system, particle_system, particle)
+        v_b = @inbounds current_velocity(v_neighbor_system, neighbor_system, neighbor)
+
         # Determine correction factors.
         # This can be ignored, as these are all 1 when no correction is used.
         (viscosity_correction, pressure_correction,
@@ -95,12 +98,17 @@ function interact!(dv, v_particle_system, u_particle_system,
             # debug_array[i, particle] += dv_pressure[i]
         end
 
+        drho_particle = Ref(zero(rho_a))
+
         # 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,
-                                       neighbor_system, v_particle_system,
-                                       v_neighbor_system, particle, neighbor,
-                                       pos_diff, distance, m_b, rho_a, rho_b, grad_kernel)
+        @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, v_a, v_b, grad_kernel)
+
+        @inbounds write_drho_particle!(dv, density_calculator, drho_particle, particle)
     end
     # Debug example
     # periodic_box = neighborhood_search.periodic_box

src/schemes/structure/structure.jl

@@ -46,6 +46,9 @@ function interact_structure_fluid!(dv, v_particle_system,
         rho_a = current_density(v_particle_system, particle_system, particle)
         rho_b = current_density(v_neighbor_system, neighbor_system, neighbor)
 
+        v_a = current_velocity(v_particle_system, particle_system, particle)
+        v_b = current_velocity(v_neighbor_system, neighbor_system, neighbor)
+
         surface_tension = surface_tension_model(neighbor_system)
 
         # In fluid-structure interaction, use the "hydrodynamic mass" of the structure particles
@@ -78,34 +81,49 @@ function interact_structure_fluid!(dv, v_particle_system,
 
         accumulate_structure_fluid_pair!(dv, dv_fs, particle_system, particle, m_b)
 
-        continuity_equation!(dv, v_particle_system, v_neighbor_system,
+        drho_particle = Ref(zero(rho_a))
+        continuity_equation!(drho_particle,
                              particle, neighbor, pos_diff, distance,
-                             m_b, rho_a, rho_b,
+                             m_b, rho_a, rho_b, v_a, v_b,
                              particle_system, neighbor_system, grad_kernel)
+
+        @inbounds write_drho_particle!(dv, particle_system, drho_particle, particle)
     end
 
     return dv
 end
 
-@inline function continuity_equation!(dv, v_particle_system, v_neighbor_system,
+@inline function continuity_equation!(drho_particle,
                                       particle, neighbor, pos_diff, distance,
-                                      m_b, rho_a, rho_b,
+                                      m_b, rho_a, rho_b, v_a, v_b,
                                       particle_system::AbstractStructureSystem,
                                       neighbor_system::AbstractFluidSystem,
                                       grad_kernel)
-    return dv
+    return drho_particle
 end
 
-@inline function continuity_equation!(dv, v_particle_system, v_neighbor_system,
+@inline function continuity_equation!(drho_particle,
                                       particle, neighbor, pos_diff, distance,
-                                      m_b, rho_a, rho_b,
+                                      m_b, rho_a, rho_b, v_a, v_b,
                                       particle_system::Union{RigidBodySystem{<:BoundaryModelDummyParticles{ContinuityDensity}},
                                                              TotalLagrangianSPHSystem{<:BoundaryModelDummyParticles{ContinuityDensity}}},
                                       neighbor_system::AbstractFluidSystem,
                                       grad_kernel)
-    v_diff = current_velocity(v_particle_system, particle_system, particle) -
-             current_velocity(v_neighbor_system, neighbor_system, neighbor)
+    continuity_equation!(drho_particle, density_calculator(neighbor_system),
+                         m_b, rho_a, rho_b, v_a, v_b, grad_kernel, particle)
+
+    return drho_particle
+end
 
-    continuity_equation!(dv, density_calculator(neighbor_system), m_b, rho_a, rho_b, v_diff,
-                         grad_kernel, particle)
+@inline function write_drho_particle!(dv, ::AbstractSystem, drho_particle, particle)
+    return dv
+end
+
+@propagate_inbounds function write_drho_particle!(dv,
+                                                  ::Union{RigidBodySystem{<:BoundaryModelDummyParticles{ContinuityDensity}},
+                                                          TotalLagrangianSPHSystem{<:BoundaryModelDummyParticles{ContinuityDensity}}},
+                                                  drho_particle, particle)
+    dv[end, particle] += drho_particle[]
+
+    return dv
 end