From 960319e5a0374e3db66b499d7b8791e441b813c3 Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Thu, 2 Apr 2026 17:58:55 +0200
Subject: [PATCH 1/7] Make `continuity_equation!` and `density_diffusion!` work
 on Ref values

---
 src/schemes/boundary/open_boundary/rhs.jl     |  9 +++---
 src/schemes/boundary/wall_boundary/rhs.jl     | 27 ++++++++++-------
 .../fluid/entropically_damped_sph/rhs.jl      | 14 ++++++---
 src/schemes/fluid/fluid.jl                    | 29 ++++++++++---------
 src/schemes/fluid/shifting_techniques.jl      | 28 ++++++++++--------
 .../density_diffusion.jl                      | 19 ++++++------
 .../fluid/weakly_compressible_sph/rhs.jl      | 16 +++++++---
 src/schemes/structure/structure.jl            | 27 ++++++++++-------
 8 files changed, 101 insertions(+), 68 deletions(-)

diff --git a/src/schemes/boundary/open_boundary/rhs.jl b/src/schemes/boundary/open_boundary/rhs.jl
index 186a9e7b7d..136912f2af 100644
--- a/src/schemes/boundary/open_boundary/rhs.jl
+++ b/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
diff --git a/src/schemes/boundary/wall_boundary/rhs.jl b/src/schemes/boundary/wall_boundary/rhs.jl
index e0553764ef..05ee411fb1 100644
--- a/src/schemes/boundary/wall_boundary/rhs.jl
+++ b/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
diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl
index cea18484d4..7e906c27fa 100644
--- a/src/schemes/fluid/entropically_damped_sph/rhs.jl
+++ b/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 dv[end, particle] += drho_particle[]
     end
 
     return dv
diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl
index 8725190506..dd2251c2c3 100644
--- a/src/schemes/fluid/fluid.jl
+++ b/src/schemes/fluid/fluid.jl
@@ -128,39 +128,42 @@ 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
 
 function calculate_dt(v_ode, u_ode, cfl_number, system::AbstractFluidSystem, semi)
diff --git a/src/schemes/fluid/shifting_techniques.jl b/src/schemes/fluid/shifting_techniques.jl
index ceaff12a21..e5c595155b 100644
--- a/src/schemes/fluid/shifting_techniques.jl
+++ b/src/schemes/fluid/shifting_techniques.jl
@@ -49,10 +49,10 @@ end
 end
 
 # Additional term(s) in the continuity equation due to the shifting technique
-@inline function continuity_equation_shifting_term(shifting, particle_system,
+@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 +352,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 +363,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 +379,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,7 +602,8 @@ 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,
@@ -607,7 +611,7 @@ end
     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,
diff --git a/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl b/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl
index ff8d8e079e..d3180b90b1 100644
--- a/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl
+++ b/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
@@ -228,13 +227,13 @@ end
 
     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
+    drho_particle[] += delta * smoothing_length_avg * 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
diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl
index acb0f06894..c95ff84498 100644
--- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl
+++ b/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 dv[end, particle] += drho_particle[]
     end
     # Debug example
     # periodic_box = neighborhood_search.periodic_box
diff --git a/src/schemes/structure/structure.jl b/src/schemes/structure/structure.jl
index 1ac304d79f..51652e9cca 100644
--- a/src/schemes/structure/structure.jl
+++ b/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,36 @@ 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)
+
+        dv[end, particle] += drho_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)
 
-    continuity_equation!(dv, density_calculator(neighbor_system), m_b, rho_a, rho_b, v_diff,
-                         grad_kernel, particle)
+    return drho_particle
 end

From 63d3eec90bd3a8aeb3ccb2f4dfa058361c09c7d4 Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Tue, 7 Apr 2026 15:45:45 +0200
Subject: [PATCH 2/7] Fix writing drho to dv for summation density

---
 src/schemes/fluid/entropically_damped_sph/rhs.jl |  2 +-
 src/schemes/fluid/fluid.jl                       | 11 +++++++++++
 src/schemes/fluid/weakly_compressible_sph/rhs.jl |  2 +-
 src/schemes/structure/structure.jl               | 11 ++++++++++-
 4 files changed, 23 insertions(+), 3 deletions(-)

diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl
index 7e906c27fa..77183b4a52 100644
--- a/src/schemes/fluid/entropically_damped_sph/rhs.jl
+++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl
@@ -92,7 +92,7 @@ function interact!(dv, v_particle_system, u_particle_system,
                                        pos_diff, distance, m_b, rho_a, rho_b,
                                        v_a, v_b, grad_kernel)
 
-        @inbounds dv[end, particle] += drho_particle[]
+        @inbounds write_drho_particle!(dv, density_calculator, drho_particle, particle)
     end
 
     return dv
diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl
index dd2251c2c3..ffd7df4ed8 100644
--- a/src/schemes/fluid/fluid.jl
+++ b/src/schemes/fluid/fluid.jl
@@ -166,6 +166,17 @@ 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)
     (; viscosity, acceleration, surface_tension) = system
 
diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl
index c95ff84498..154f433f10 100644
--- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl
+++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl
@@ -108,7 +108,7 @@ function interact!(dv, v_particle_system, u_particle_system,
                                        particle, neighbor, pos_diff, distance,
                                        m_b, rho_a, rho_b, v_a, v_b, grad_kernel)
 
-        @inbounds dv[end, particle] += drho_particle[]
+        @inbounds write_drho_particle!(dv, density_calculator, drho_particle, particle)
     end
     # Debug example
     # periodic_box = neighborhood_search.periodic_box
diff --git a/src/schemes/structure/structure.jl b/src/schemes/structure/structure.jl
index 51652e9cca..b4939232ee 100644
--- a/src/schemes/structure/structure.jl
+++ b/src/schemes/structure/structure.jl
@@ -87,7 +87,7 @@ function interact_structure_fluid!(dv, v_particle_system,
                              m_b, rho_a, rho_b, v_a, v_b,
                              particle_system, neighbor_system, grad_kernel)
 
-        dv[end, particle] += drho_particle[]
+        @inbounds write_drho_particle!(dv, particle_system, drho_particle, particle)
     end
 
     return dv
@@ -114,3 +114,12 @@ end
 
     return drho_particle
 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

From 9887d1e6d69889c331caf6670afe1fc2a5f09ebd Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Tue, 7 Apr 2026 18:37:03 +0200
Subject: [PATCH 3/7] Implement copilot suggestions

---
 src/schemes/fluid/shifting_techniques.jl | 3 ++-
 src/schemes/structure/structure.jl       | 4 ++++
 2 files changed, 6 insertions(+), 1 deletion(-)

diff --git a/src/schemes/fluid/shifting_techniques.jl b/src/schemes/fluid/shifting_techniques.jl
index e5c595155b..f1fdc03c63 100644
--- a/src/schemes/fluid/shifting_techniques.jl
+++ b/src/schemes/fluid/shifting_techniques.jl
@@ -48,7 +48,8 @@ end
     return zero(grad_kernel)
 end
 
-# Additional term(s) in the continuity equation due to the shifting technique
+# 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)
diff --git a/src/schemes/structure/structure.jl b/src/schemes/structure/structure.jl
index b4939232ee..038389ce3c 100644
--- a/src/schemes/structure/structure.jl
+++ b/src/schemes/structure/structure.jl
@@ -115,6 +115,10 @@ end
     return drho_particle
 end
 
+@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}}},

From ffc212086d196da9d724a7e5f3148c7a3c64dd35 Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Tue, 7 Apr 2026 19:18:10 +0200
Subject: [PATCH 4/7] Fix density diffusion smoothing length

---
 .../fluid/weakly_compressible_sph/density_diffusion.jl      | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl b/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl
index d3180b90b1..aa3d23da87 100644
--- a/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl
+++ b/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl
@@ -225,10 +225,8 @@ 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
-    drho_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

From 89d40f299e550028fc9603c6a06f41cdbb152dbe Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Tue, 14 Apr 2026 11:20:53 +0200
Subject: [PATCH 5/7] Fix missing suggestion from #1043

---
 src/schemes/boundary/open_boundary/rhs.jl |  2 +-
 src/schemes/fluid/fluid.jl                | 11 -----------
 2 files changed, 1 insertion(+), 12 deletions(-)

diff --git a/src/schemes/boundary/open_boundary/rhs.jl b/src/schemes/boundary/open_boundary/rhs.jl
index e3d5b055e8..debaf085f4 100644
--- a/src/schemes/boundary/open_boundary/rhs.jl
+++ b/src/schemes/boundary/open_boundary/rhs.jl
@@ -83,7 +83,7 @@ function interact!(dv, v_particle_system, u_particle_system,
 
         # Propagate `@inbounds` to the continuity equation, which accesses particle data
         drho_particle = Ref(zero(rho_a))
-        @inbounds continuity_equation!(drho_particle,
+        @inbounds continuity_equation!(drho_particle, ContinuityDensity(),
                                        particle_system, neighbor_system,
                                        v_particle_system, v_neighbor_system,
                                        particle, neighbor, pos_diff, distance,
diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl
index 8d8ac38ead..45bb9c6be2 100644
--- a/src/schemes/fluid/fluid.jl
+++ b/src/schemes/fluid/fluid.jl
@@ -145,17 +145,6 @@ end
                                                   particle, neighbor, pos_diff, distance,
                                                   m_b, rho_a, rho_b, v_a, v_b,
                                                   grad_kernel)
-    return continuity_equation!(drho_particle, 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)
-end
-
-@propagate_inbounds function continuity_equation!(drho_particle,
-                                                  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)
     v_diff = v_a - v_b
 
     v_diff = continuity_equation_shifting_term(v_diff,

From 6bc3034dcb6e1f37d960158e662b41ad269343da Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Tue, 14 Apr 2026 11:22:49 +0200
Subject: [PATCH 6/7] Fix

---
 src/schemes/boundary/open_boundary/rhs.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/schemes/boundary/open_boundary/rhs.jl b/src/schemes/boundary/open_boundary/rhs.jl
index debaf085f4..1ecb7414a4 100644
--- a/src/schemes/boundary/open_boundary/rhs.jl
+++ b/src/schemes/boundary/open_boundary/rhs.jl
@@ -88,6 +88,7 @@ function interact!(dv, v_particle_system, u_particle_system,
                                        v_particle_system, v_neighbor_system,
                                        particle, neighbor, pos_diff, distance,
                                        m_b, rho_a, rho_b, v_a, v_b, grad_kernel)
+        dv[end, particle] += drho_particle[]
 
         # Open boundary pressure evolution matches the corresponding fluid system:
         # - EDAC: Compute pressure evolution like the fluid system

From 5f36a062e61b77ccd14d257165634fb6f079bd8c Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Tue, 14 Apr 2026 13:01:23 +0200
Subject: [PATCH 7/7] Fix `continuity_equation!` functions

---
 src/schemes/boundary/open_boundary/rhs.jl        |  3 +--
 src/schemes/boundary/wall_boundary/rhs.jl        | 10 ++++++----
 src/schemes/fluid/entropically_damped_sph/rhs.jl |  9 ++++-----
 src/schemes/fluid/fluid.jl                       | 14 ++++++++++----
 src/schemes/fluid/weakly_compressible_sph/rhs.jl |  1 -
 src/schemes/structure/structure.jl               | 15 ++++++---------
 6 files changed, 27 insertions(+), 25 deletions(-)

diff --git a/src/schemes/boundary/open_boundary/rhs.jl b/src/schemes/boundary/open_boundary/rhs.jl
index 1ecb7414a4..74a47faa55 100644
--- a/src/schemes/boundary/open_boundary/rhs.jl
+++ b/src/schemes/boundary/open_boundary/rhs.jl
@@ -83,9 +83,8 @@ function interact!(dv, v_particle_system, u_particle_system,
 
         # Propagate `@inbounds` to the continuity equation, which accesses particle data
         drho_particle = Ref(zero(rho_a))
-        @inbounds continuity_equation!(drho_particle, ContinuityDensity(),
+        @inbounds continuity_equation!(drho_particle,
                                        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)
         dv[end, particle] += drho_particle[]
diff --git a/src/schemes/boundary/wall_boundary/rhs.jl b/src/schemes/boundary/wall_boundary/rhs.jl
index f6bc17f1f2..63811938e5 100644
--- a/src/schemes/boundary/wall_boundary/rhs.jl
+++ b/src/schemes/boundary/wall_boundary/rhs.jl
@@ -61,16 +61,18 @@ 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!(drho_particle, ::SummationDensity,
-                              m_b, rho_a, rho_b, v_a, v_b, grad_kernel, particle)
+@propagate_inbounds 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!(drho_particle, ::ContinuityDensity,
-                              m_b, rho_a, rho_b, v_a, v_b, grad_kernel, particle)
+@propagate_inbounds 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
diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl
index f1e7dc0cd5..142dd2ccf1 100644
--- a/src/schemes/fluid/entropically_damped_sph/rhs.jl
+++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl
@@ -105,11 +105,10 @@ 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!(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 continuity_equation!(drho_particle, density_calculator,
+                                       particle_system, 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
diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl
index 45bb9c6be2..732fa9b82a 100644
--- a/src/schemes/fluid/fluid.jl
+++ b/src/schemes/fluid/fluid.jl
@@ -130,7 +130,6 @@ end
 # With 'SummationDensity', density is calculated in wcsph/system.jl:compute_density!
 @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, v_a, v_b, grad_kernel)
     return drho_particle
@@ -141,10 +140,17 @@ end
                                                   ::ContinuityDensity,
                                                   particle_system::AbstractFluidSystem,
                                                   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)
+                                                  m_b, rho_a, rho_b, v_a, v_b, grad_kernel)
+    continuity_equation!(drho_particle, particle_system, neighbor_system,
+                         particle, neighbor, pos_diff, distance,
+                         m_b, rho_a, rho_b, v_a, v_b, grad_kernel)
+end
+
+@propagate_inbounds function continuity_equation!(drho_particle,
+                                                  particle_system, neighbor_system,
+                                                  particle, neighbor, pos_diff, distance,
+                                                  m_b, rho_a, rho_b, v_a, v_b, grad_kernel)
     v_diff = v_a - v_b
 
     v_diff = continuity_equation_shifting_term(v_diff,
diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl
index bac52e78b2..27c3ff8dbf 100644
--- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl
+++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl
@@ -123,7 +123,6 @@ function interact!(dv, v_particle_system, u_particle_system,
         # Propagate `@inbounds` to the continuity equation, which accesses particle data
         @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)
 
diff --git a/src/schemes/structure/structure.jl b/src/schemes/structure/structure.jl
index 5a4bfc963d..2a7f98d3c1 100644
--- a/src/schemes/structure/structure.jl
+++ b/src/schemes/structure/structure.jl
@@ -94,10 +94,9 @@ function interact_structure_fluid!(dv, v_particle_system, u_particle_system,
         accumulate_structure_fluid_pair!(dv, dv_particle[], particle_system, particle, m_b)
 
         drho_particle = Ref(zero(rho_a))
-        continuity_equation!(drho_particle,
+        continuity_equation!(drho_particle, particle_system, neighbor_system,
                              particle, neighbor, pos_diff, distance,
-                             m_b, rho_a, rho_b, v_a, v_b,
-                             particle_system, neighbor_system, grad_kernel)
+                             m_b, rho_a, rho_b, v_a, v_b, grad_kernel)
 
         @inbounds write_drho_particle!(dv, particle_system, drho_particle, particle)
     end
@@ -106,21 +105,19 @@ function interact_structure_fluid!(dv, v_particle_system, u_particle_system,
 end
 
 @inline function continuity_equation!(drho_particle,
-                                      particle, neighbor, pos_diff, distance,
-                                      m_b, rho_a, rho_b, v_a, v_b,
                                       particle_system::AbstractStructureSystem,
                                       neighbor_system::AbstractFluidSystem,
-                                      grad_kernel)
+                                      particle, neighbor, pos_diff, distance,
+                                      m_b, rho_a, rho_b, v_a, v_b, grad_kernel)
     return drho_particle
 end
 
 @inline function continuity_equation!(drho_particle,
-                                      particle, neighbor, pos_diff, distance,
-                                      m_b, rho_a, rho_b, v_a, v_b,
                                       particle_system::Union{RigidBodySystem{<:BoundaryModelDummyParticles{ContinuityDensity}},
                                                              TotalLagrangianSPHSystem{<:BoundaryModelDummyParticles{ContinuityDensity}}},
                                       neighbor_system::AbstractFluidSystem,
-                                      grad_kernel)
+                                      particle, neighbor, pos_diff, distance,
+                                      m_b, rho_a, rho_b, v_a, v_b, grad_kernel)
     continuity_equation!(drho_particle, density_calculator(neighbor_system),
                          m_b, rho_a, rho_b, v_a, v_b, grad_kernel, particle)