Rework pressure_acceleration (again) (#334)

* proof of concept * finalize * fix correction bug * remove inter particle averaged pressure * implement suggestions * fix bugs * get rid of dispatching * add TLSP to Monaghan `pressure_acceleration` * implement suggestions * fix monaghan interaction * implement suggestions * add tests * remove unused line * apply formatter * implement suggestions * fix bug --------- Co-authored-by: Erik Faulhaber <[email protected]>
trixi-framework · Jan 25, 2024 · 882d1c6 · 882d1c6
1 parent 898d8d8
commit 882d1c6
Show file tree

Hide file tree

Showing 13 changed files with 267 additions and 212 deletions.
diff --git a/src/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl
@@ -274,65 +274,6 @@ 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)
-
-    # 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)
-
-    # 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)
-end
-
 @inline function particle_density(v, model::BoundaryModelDummyParticles, system, particle)
     return particle_density(v, model.density_calculator, model, particle)
 end

diff --git a/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl b/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl
@@ -76,19 +76,17 @@ 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}
-    (; K, beta, boundary_particle_spacing) = boundary_model
-
-    distance = norm(pos_diff)
-    return K / beta^(NDIMS - 1) * pos_diff /
+@inline function pressure_acceleration(particle_system,
+                                       neighbor_system::Union{BoundarySystem{<:BoundaryModelMonaghanKajtar},
+                                                              SolidSystem{<:BoundaryModelMonaghanKajtar}},
+                                       neighbor, m_a, m_b, p_a, p_b, rho_a, rho_b,
+                                       pos_diff, distance, grad_kernel,
+                                       pressure_correction, correction)
+    (; K, beta, boundary_particle_spacing) = neighbor_system.boundary_model
+
+    return K / beta^(ndims(particle_system) - 1) * pos_diff /
            (distance * (distance - boundary_particle_spacing)) *
-           boundary_kernel(distance, smoothing_length)
+           boundary_kernel(distance, particle_system.smoothing_length)
 end
 
 @fastpow @inline function boundary_kernel(r, h)

diff --git a/src/schemes/boundary/system.jl b/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
diff --git a/src/schemes/fluid/fluid.jl b/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")
diff --git a/src/schemes/fluid/pressure_acceleration.jl b/src/schemes/fluid/pressure_acceleration.jl
@@ -0,0 +1,116 @@
+# 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 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
+
+# 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 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
+
+# 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 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 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
+
+function choose_pressure_acceleration_formulation(pressure_acceleration,
+                                                  density_calculator, NDIMS, ELTYPE,
+                                                  correction)
+    if correction isa KernelCorrection ||
+       correction isa GradientCorrection ||
+       correction isa BlendedGradientCorrection ||
+       correction isa MixedKernelGradientCorrection
+        if isempty(methods(pressure_acceleration,
+                           (ELTYPE, ELTYPE, ELTYPE, ELTYPE, ELTYPE, ELTYPE,
+                            SVector{NDIMS, ELTYPE}, SVector{NDIMS, ELTYPE})))
+            throw(ArgumentError("when a correction with an asymmetric kernel gradient is " *
+                                "used, the passed pressure acceleration formulation must " *
+                                "provide a version with the arguments " *
+                                "`m_a, m_b, rho_a, rho_b, p_a, p_b, W_a, W_b`"))
+        end
+    else
+        if isempty(methods(pressure_acceleration,
+                           (ELTYPE, ELTYPE, ELTYPE, ELTYPE, ELTYPE, ELTYPE,
+                            SVector{NDIMS, ELTYPE})))
+            throw(ArgumentError("when not using a correction with an asymmetric kernel " *
+                                "gradient, the passed pressure acceleration formulation must " *
+                                "provide a version with the arguments " *
+                                "`m_a, m_b, rho_a, rho_b, p_a, p_b, W_a`, " *
+                                "using the symmetry of the kernel gradient"))
+        end
+    end
+
+    return pressure_acceleration
+end
+
+function choose_pressure_acceleration_formulation(pressure_acceleration::Nothing,
+                                                  density_calculator::SummationDensity,
+                                                  NDIMS, ELTYPE,
+                                                  correction)
+
+    # Choose the pressure acceleration formulation corresponding to the density calculator.
+    return pressure_acceleration_summation_density
+end
+
+function choose_pressure_acceleration_formulation(pressure_acceleration::Nothing,
+                                                  density_calculator::ContinuityDensity,
+                                                  NDIMS, ELTYPE,
+                                                  correction)
+
+    # Choose the pressure acceleration formulation corresponding to the density calculator.
+    return pressure_acceleration_continuity_density
+end
+
+# Formulation using symmetric gradient formulation for corrections not depending on local neighborhood.
+@inline function pressure_acceleration(particle_system, neighbor_system, neighbor,
+                                       m_a, m_b, p_a, p_b, rho_a, rho_b, pos_diff,
+                                       distance, W_a, pressure_correction,
+                                       correction)
+    (; pressure_acceleration_formulation) = particle_system
+
+    # Without correction or with `AkinciFreeSurfaceCorrection`, the kernel gradient is
+    # symmetric, so call the symmetric version of the pressure acceleration formulation.
+    return pressure_acceleration_formulation(m_a, m_b, rho_a, rho_b, p_a, p_b, W_a) *
+           pressure_correction
+end
+
+# Formulation using asymmetric gradient formulation for corrections depending on local neighborhood.
+@inline function pressure_acceleration(particle_system, neighbor_system, neighbor,
+                                       m_a, m_b, p_a, p_b, rho_a, rho_b, pos_diff,
+                                       distance, W_a, pressure_correction,
+                                       correction::Union{KernelCorrection,
+                                                         GradientCorrection,
+                                                         BlendedGradientCorrection,
+                                                         MixedKernelGradientCorrection})
+    (; pressure_acceleration_formulation) = 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 version of the pressure acceleration formulation.
+    return pressure_acceleration_formulation(m_a, m_b, rho_a, rho_b, p_a, p_b, W_a, W_b) *
+           pressure_correction
+end