Rework pressure_acceleration (again)

src/schemes/boundary/dummy_particles/dummy_particles.jl

@@ -274,63 +274,44 @@ function initial_boundary_pressure(initial_density, ::PressureZeroing, ::Nothing
     return zero(initial_density)
 end
 
-@inline function pressure_acceleration_bnd(pressure_correction, m_b, p_a, p_b,
-                                           rho_a, rho_b, pos_diff, distance,
-                                           smoothing_length, grad_kernel,
-                                           particle_system, neighbor, neighbor_system,
-                                           boundary_model::BoundaryModelDummyParticles{<:PressureMirroring},
-                                           fluid_density_calculator,
-                                           correction)
+# ==== Pressure acceleration without correction
+@inline function pressure_acceleration(pressure_correction, m_a, m_b, p_a, p_b,
+                                       rho_a, rho_b, pos_diff, distance,
+                                       smoothing_length, W_a,
+                                       boundary_model::BoundaryModelDummyParticles{<:PressureMirroring},
+                                       pressure_gradient)
 
     # Use `p_a` as pressure for both particles with `PressureMirroring`
-    return pressure_acceleration_symmetric(pressure_correction, m_b, p_a, p_a, rho_a, rho_b,
-                                           grad_kernel, fluid_density_calculator)
-end
-
-@inline function pressure_acceleration_bnd(pressure_correction, m_b, p_a, p_b,
-                                           rho_a, rho_b, pos_diff, distance,
-                                           smoothing_length, grad_kernel,
-                                           particle_system, neighbor, neighbor_system,
-                                           boundary_model::BoundaryModelDummyParticles,
-                                           fluid_density_calculator,
-                                           correction)
-    return pressure_acceleration_symmetric(pressure_correction, m_b, p_a, p_b, rho_a, rho_b,
-                                           grad_kernel, fluid_density_calculator)
-end
-
-@inline function pressure_acceleration_bnd(pressure_correction, m_b, p_a, p_b,
-                                           rho_a, rho_b, pos_diff, distance,
-                                           smoothing_length, W_a,
-                                           particle_system, neighbor, neighbor_system,
-                                           boundary_model::BoundaryModelDummyParticles{<:PressureMirroring},
-                                           fluid_density_calculator,
-                                           correction::Union{KernelCorrection,
-                                                             GradientCorrection,
-                                                             BlendedGradientCorrection,
-                                                             MixedKernelGradientCorrection})
-    W_b = smoothing_kernel_grad(neighbor_system, -pos_diff, distance, neighbor)
+    return pressure_gradient(m_a, m_b, rho_a, rho_b, p_a, p_a, W_a)
+end
+
+@inline function pressure_acceleration(pressure_correction, m_a, m_b, p_a, p_b,
+                                       rho_a, rho_b, pos_diff, distance,
+                                       smoothing_length, W_a,
+                                       boundary_model::BoundaryModelDummyParticles,
+                                       pressure_gradient)
+    return pressure_gradient(m_a, m_b, rho_a, rho_b, p_a, p_b, W_a)
+end
+
+# ==== Pressure acceleration with correction
+@inline function pressure_acceleration(pressure_correction, m_a, m_b, p_a, p_b,
+                                       rho_a, rho_b, pos_diff, distance,
+                                       smoothing_length, W_a, W_b,
+                                       boundary_model::BoundaryModelDummyParticles{<:PressureMirroring},
+                                       pressure_gradient)
 
     # Use `p_a` as pressure for both particles with `PressureMirroring`
-    return pressure_acceleration_asymmetric(pressure_correction, m_b, p_a, p_a,
-                                            rho_a, rho_b, W_a, W_b,
-                                            fluid_density_calculator)
-end
-
-@inline function pressure_acceleration_bnd(pressure_correction, m_b, p_a, p_b,
-                                           rho_a, rho_b, pos_diff, distance,
-                                           smoothing_length, W_a,
-                                           particle_system, neighbor, neighbor_system,
-                                           boundary_model::BoundaryModelDummyParticles,
-                                           fluid_density_calculator,
-                                           correction::Union{KernelCorrection,
-                                                             GradientCorrection,
-                                                             BlendedGradientCorrection,
-                                                             MixedKernelGradientCorrection})
-    W_b = smoothing_kernel_grad(neighbor_system, -pos_diff, distance, neighbor)
-
-    return pressure_acceleration_asymmetric(pressure_correction, m_b, p_a, p_b,
-                                            rho_a, rho_b, W_a, W_b,
-                                            fluid_density_calculator)
+    return pressure_gradient(m_a, m_b, rho_a, rho_b, p_a, p_a, W_a, W_b) *
+           pressure_correction
+end
+
+@inline function pressure_acceleration(pressure_correction, m_a, m_b, p_a, p_b,
+                                       rho_a, rho_b, pos_diff, distance,
+                                       smoothing_length, W_a, W_b,
+                                       boundary_model::BoundaryModelDummyParticles,
+                                       pressure_gradient)
+    return pressure_gradient(m_a, m_b, rho_a, rho_b, p_a, p_b, W_a, W_b) *
+           pressure_correction
 end
 
 @inline function particle_density(v, model::BoundaryModelDummyParticles, system, particle)

src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl

@@ -76,16 +76,13 @@ function Base.show(io::IO, model::BoundaryModelMonaghanKajtar)
     print(io, ")")
 end
 
-@inline function pressure_acceleration_bnd(pressure_correction, m_b, p_a, p_b,
-                                           rho_a, rho_b, pos_diff::SVector{NDIMS}, distance,
-                                           smoothing_length, grad_kernel,
-                                           particle_system, neighbor, neighbor_system,
-                                           boundary_model::BoundaryModelMonaghanKajtar,
-                                           fluid_density_calculator,
-                                           correction) where {NDIMS}
+@inline function pressure_acceleration(pressure_correction, m_a, m_b, p_a, p_b,
+                                       rho_a, rho_b, pos_diff::SVector{NDIMS}, distance,
+                                       smoothing_length, grad_kernel,
+                                       boundary_model::BoundaryModelMonaghanKajtar,
+                                       pressure_gradient) where {NDIMS}
     (; K, beta, boundary_particle_spacing) = boundary_model
 
-    distance = norm(pos_diff)
     return K / beta^(NDIMS - 1) * pos_diff /
            (distance * (distance - boundary_particle_spacing)) *
            boundary_kernel(distance, smoothing_length)

src/schemes/boundary/system.jl

@@ -264,19 +264,3 @@ end
 function viscosity_model(system::BoundarySPHSystem)
     return system.boundary_model.viscosity
 end
-
-@inline function pressure_acceleration(pressure_correction, m_b, p_a, p_b,
-                                       rho_a, rho_b, pos_diff, distance, grad_kernel,
-                                       particle_system, neighbor,
-                                       neighbor_system::BoundarySPHSystem,
-                                       density_calculator, correction)
-    (; boundary_model) = neighbor_system
-    (; smoothing_length) = particle_system
-
-    return pressure_acceleration_bnd(pressure_correction, m_b, p_a, p_b,
-                                     rho_a, rho_b, pos_diff, distance,
-                                     smoothing_length, grad_kernel,
-                                     particle_system, neighbor, neighbor_system,
-                                     boundary_model,
-                                     density_calculator, correction)
-end

src/schemes/fluid/fluid.jl

@@ -15,6 +15,7 @@ end
 
 @inline viscosity_model(system::FluidSystem) = system.viscosity
 
+include("pressure_acceleration.jl")
 include("viscosity.jl")
 include("weakly_compressible_sph/weakly_compressible_sph.jl")
 include("entropically_damped_sph/entropically_damped_sph.jl")

src/schemes/fluid/pressure_acceleration.jl

@@ -0,0 +1,155 @@
+# As shown in "Variational and momentum preservation aspects of Smooth Particle Hydrodynamic
+# formulations" by Bonet and Lok (1999), for a consistent formulation this form has to be
+# used with `SummationDensity`.
+# This can also be seen in the tests for total energy conservation, which fail with the
+# other `pressure_acceleration` form.
+# We assume symmetry of the kernel gradient in this formulation. See below for the
+# asymmetric version.
+@inline function symmetric_pressure_acceleration_summation_density(m_a, m_b, rho_a, rho_b,
+                                                                   p_a, p_b, W_a)
+    return -m_b * (p_a / rho_a^2 + p_b / rho_b^2) * W_a
+end
+
+# As shown in "Variational and momentum preservation aspects of Smooth Particle Hydrodynamic
+# formulations" by Bonet and Lok (1999), for a consistent formulation this form has to be
+# used with `ContinuityDensity` with the formulation `\rho_a * \sum m_b / \rho_b ...`.
+# This can also be seen in the tests for total energy conservation, which fail with the
+# other `pressure_acceleration` form.
+# We assume symmetry of the kernel gradient in this formulation. See below for the
+# asymmetric version.
+@inline function symmetric_pressure_acceleration_continuity_density(m_a, m_b, rho_a, rho_b,
+                                                                    p_a, p_b, W_a)
+    return -m_b * (p_a + p_b) / (rho_a * rho_b) * W_a
+end
+
+# Same as above, but not assuming symmetry of the kernel gradient. To be used with
+# corrections that do not produce a symmetric kernel gradient.
+@inline function asymmetric_pressure_acceleration_summation_density(m_a, m_b, rho_a, rho_b,
+                                                                    p_a, p_b, W_a, W_b)
+    return -m_b * (p_a / rho_a^2 * W_a - p_b / rho_b^2 * W_b)
+end
+
+# Same as above, but not assuming symmetry of the kernel gradient. To be used with
+# corrections that do not produce a symmetric kernel gradient.
+@inline function asymmetric_pressure_acceleration_continuity_density(m_a, m_b, rho_a, rho_b,
+                                                                     p_a, p_b, W_a, W_b)
+    return -m_b / (rho_a * rho_b) * (p_a * W_a - p_b * W_b)
+end
+
+@inline function interparticle_averaged_pressure_acceleration(m_a, m_b, rho_a, rho_b,
+                                                              p_a, p_b, W_a)
+    volume_a = m_a / rho_a
+    volume_b = m_b / rho_b
+    volume_term = (volume_a^2 + volume_b^2) / m_a
+
+    pressure_tilde = (rho_b * p_a + rho_a * p_b) / (rho_a + rho_b)
+
+    return -volume_term * pressure_tilde * W_a
+end
+
+function set_pressure_acceleration(pressure_acceleration,
+                                   density_calculator, correction)
+    return pressure_acceleration
+end
+
+function set_pressure_acceleration(pressure_acceleration::Nothing,
+                                   density_calculator::SummationDensity,
+                                   correction::Nothing)
+    return symmetric_pressure_acceleration_summation_density
+end
+
+function set_pressure_acceleration(pressure_acceleration::Nothing,
+                                   density_calculator::ContinuityDensity,
+                                   correction::Nothing)
+    return symmetric_pressure_acceleration_continuity_density
+end
+
+function set_pressure_acceleration(pressure_acceleration::Nothing,
+                                   density_calculator::SummationDensity,
+                                   correction::Union{KernelCorrection,
+                                                     GradientCorrection,
+                                                     BlendedGradientCorrection,
+                                                     MixedKernelGradientCorrection})
+    return asymmetric_pressure_acceleration_summation_density
+end
+
+function set_pressure_acceleration(pressure_acceleration::Nothing,
+                                   density_calculator::ContinuityDensity,
+                                   correction::Union{KernelCorrection,
+                                                     GradientCorrection,
+                                                     BlendedGradientCorrection,
+                                                     MixedKernelGradientCorrection})
+    return asymmetric_pressure_acceleration_continuity_density
+end
+
+# ==== Fluid System
+
+# No correction
+@inline function pressure_acceleration(pressure_correction, m_a, m_b, p_a, p_b,
+                                       rho_a, rho_b, pos_diff, distance,
+                                       W_a, particle_system, neighbor_system::FluidSystem,
+                                       correction)
+    (; pressure_gradient) = neighbor_system
+
+    # Without correction, the kernel gradient is symmetric, so call the symmetric
+    # pressure acceleration formulation corresponding to the density calculator.
+    return pressure_gradient(m_a, m_b, rho_a, rho_b, p_a, p_b, W_a)
+end
+
+# Correction
+@inline function pressure_acceleration(pressure_correction, m_a, m_b, p_a, p_b,
+                                       rho_a, rho_b, pos_diff, distance,
+                                       W_a, particle_system, neighbor_system::FluidSystem,
+                                       correction::Union{KernelCorrection,
+                                                         GradientCorrection,
+                                                         BlendedGradientCorrection,
+                                                         MixedKernelGradientCorrection})
+    (; pressure_gradient, smoothing_kernel, smoothing_length) = neighbor_system
+
+    W_b = kernel_grad(smoothing_kernel, -pos_diff, distance, smoothing_length)
+
+    # With correction, the kernel gradient is not necessarily symmetric, so call the
+    # asymmetric pressure acceleration formulation corresponding to the density calculator.
+    return pressure_gradient(m_a, m_b, rho_a, rho_b, p_a, p_b, W_a, W_b) *
+           pressure_correction
+end
+
+# ==== Boundary and Solid System
+
+# No correction
+@inline function pressure_acceleration(pressure_correction, m_a, m_b, p_a, p_b,
+                                       rho_a, rho_b, pos_diff, distance, W_a,
+                                       particle_system,
+                                       neighbor_system::Union{BoundarySystem, SolidSystem},
+                                       correction)
+    (; boundary_model) = neighbor_system
+    (; smoothing_length, pressure_gradient) = particle_system
+
+    # Without correction, the kernel gradient is symmetric, so call the symmetric
+    # pressure acceleration formulation corresponding to the density calculator.
+    return pressure_acceleration(pressure_correction, m_a, m_b, p_a, p_b,
+                                 rho_a, rho_b, pos_diff, distance,
+                                 smoothing_length, W_a,
+                                 boundary_model, pressure_gradient)
+end
+
+# Correction
+@inline function pressure_acceleration(pressure_correction, m_a, m_b, p_a, p_b,
+                                       rho_a, rho_b, pos_diff, distance, W_a,
+                                       particle_system,
+                                       neighbor_system::Union{BoundarySystem, SolidSystem},
+                                       correction::Union{KernelCorrection,
+                                                         GradientCorrection,
+                                                         BlendedGradientCorrection,
+                                                         MixedKernelGradientCorrection})
+    (; boundary_model) = neighbor_system
+    (; smoothing_length) = particle_system
+    W_b = smoothing_kernel_grad(neighbor_system, -pos_diff, distance, neighbor)
+
+    # With correction, the kernel gradient is not necessarily symmetric, so call the
+    # asymmetric pressure acceleration formulation corresponding to the density calculator.
+    return pressure_acceleration(pressure_correction, m_a, m_b, p_a, p_b,
+                                 rho_a, rho_b, pos_diff, distance,
+                                 smoothing_length, W_a, W_b,
+                                 boundary_model, pressure_gradient)
+end