New boundary density calculators PressureZeroing and PressureMirroring

src/TrixiParticles.jl

@@ -44,7 +44,8 @@ export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel,
        SchoenbergQuinticSplineKernel
 export StateEquationIdealGas, StateEquationCole
 export ArtificialViscosityMonaghan, ViscosityAdami
-export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureExtrapolation
+export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureExtrapolation,
+       PressureMirroring, PressureZeroing
 export BoundaryMovement
 export GridNeighborhoodSearch
 export examples_dir, trixi_include

src/schemes/boundary/dummy_particles/dummy_particles.jl

@@ -33,19 +33,28 @@ f_{ab} = m_a m_b \left( \frac{p_a}{\rho_a^2} + \frac{p_b}{\rho_b^2} \right) \nab
 The quantities to be defined here are the density ``\rho_b`` and pressure ``p_b``
 of the boundary particle ``b``.
 
-We provide three options to compute the boundary density and pressure, determined by the `density_calculator`:
-1. With [`SummationDensity`](@ref), the density is calculated by summation over the neighboring particles,
+We provide five options to compute the boundary density and pressure, determined by the `density_calculator`:
+1. With [`AdamiPressureExtrapolation`](@ref), the pressure is extrapolated from the pressure of the
+   fluid according to (Adami et al., 2012), and the density is obtained by applying the inverse of the state equation.
+   This option usually yields the best results of the options listed here.
+2. With [`SummationDensity`](@ref), the density is calculated by summation over the neighboring particles,
    and the pressure is computed from the density with the state equation.
-2. With [`ContinuityDensity`](@ref), the density is integrated from the continuity equation,
+3. With [`ContinuityDensity`](@ref), the density is integrated from the continuity equation,
    and the pressure is computed from the density with the state equation.
    Note that this causes a gap between fluid and boundary where the boundary is initialized
    without any contact to the fluid. This is due to overestimation of the boundary density
    as soon as the fluid comes in contact with boundary particles that initially did not have
    contact to the fluid.
    Therefore, in dam break simulations, there is a visible "step", even though the boundary is supposed to be flat.
    See also [dual.sphysics.org/faq/#Q_13](https://dual.sphysics.org/faq/#Q_13).
-3. With [`AdamiPressureExtrapolation`](@ref), the pressure is extrapolated from the pressure of the
-   fluid according to (Adami et al., 2012), and the density is obtained by applying the inverse of the state equation.
+4. With [`PressureZeroing`](@ref), the density is set to the reference density and the pressure
+   is computed from the density with the state equation.
+   This option is not recommended. The other options yield significantly better results.
+5. With [`PressureMirroring`](@ref), the density is set to the reference density. The pressure
+   is not used. Instead, the fluid pressure is mirrored as boundary pressure in the
+   momentum equation.
+   This option is not recommended due to stability issues. See [`PressureMirroring`](@ref)
+   for more details.
 
 ## References:
 - S. Adami, X. Y. Hu, N. A. Adams.
@@ -79,7 +88,8 @@ struct BoundaryModelDummyParticles{ELTYPE <: Real, SE, DC, K, V, C}
                                          density_calculator, smoothing_kernel,
                                          smoothing_length; viscosity=NoViscosity(),
                                          state_equation=nothing)
-        pressure = similar(initial_density)
+        pressure = initial_boundary_pressure(initial_density, density_calculator,
+                                             state_equation)
 
         n_particles = length(initial_density)
 
@@ -103,22 +113,6 @@ function Base.show(io::IO, model::BoundaryModelDummyParticles)
     print(io, ")")
 end
 
-@inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system,
-                                       v_particle_system, boundary_particle,
-                                       boundary_system,
-                                       v_boundary_system,
-                                       boundary_model::BoundaryModelDummyParticles, rho_a,
-                                       rho_b, pos_diff, distance, grad_kernel,
-                                       density_calculator)
-    rho_a = particle_density(v_particle_system, particle_system, particle)
-    rho_b = particle_density(v_boundary_system, boundary_system, boundary_particle)
-
-    return -m_b *
-           (particle_system.pressure[particle] / rho_a^2 +
-            boundary_model.pressure[boundary_particle] / rho_b^2) *
-           grad_kernel
-end
-
 @doc raw"""
     AdamiPressureExtrapolation()
 
@@ -140,7 +134,172 @@ where the sum is over all fluid particles, ``\rho_f`` and ``p_f`` denote the den
 """
 struct AdamiPressureExtrapolation end
 
-function create_cache(initial_density, ::SummationDensity)
+@doc raw"""
+    PressureMirroring()
+
+Instead of calculating density and pressure for each boundary particle, we modify the
+momentum equation,
+```math
+\frac{\mathrm{d}v_a}{\mathrm{d}t} = -\sum_b m_b \left( \frac{p_a}{\rho_a^2} + \frac{p_b}{\rho_b^2} \right) \nabla_a W_{ab}
+```
+to replace the unknown density $\rho_b$ if $b$ is a boundary particle by the reference density
+and the unknown pressure $p_b$ if $b$ is a boundary particle by the pressure $p_a$ of the
+interacting fluid particle.
+The momentum equation therefore becomes
+```math
+\frac{\mathrm{d}v_a}{\mathrm{d}t} = -\sum_f m_f \left( \frac{p_a}{\rho_a^2} + \frac{p_f}{\rho_f^2} \right) \nabla_a W_{af}
+-\sum_b m_b \left( \frac{p_a}{\rho_a^2} + \frac{p_a}{\rho_0^2} \right) \nabla_a W_{ab},
+```
+where the first sum is over all fluid particles and the second over all boundary particles.
+
+This approach was first mentioned by Akinci et al. (2012) and written down in this form
+by Band et al. (2018).
+
+!! note
+    This boundary model requires high viscosity for stability with WCSPH.
+    It also produces significantly worse results than [`AdamiPressureExtrapolation`](@ref)
+    and is not more efficient because smaller time steps are required due to more noise
+    in the pressure.
+    We added this model only for research purposes and for comparison with
+    [SPlisHSPlasH](https://github.com/InteractiveComputerGraphics/SPlisHSPlasH).
+
+## References:
+- Nadir Akinci, Markus Ihmsen, Gizem Akinci, Barbara Solenthaler, and Matthias Teschner.
+  "Versatile Rigid-Fluid Coupling for Incompressible SPH."
+  In: ACM Transactions on Graphics 31, 4 (2012), pages 1–8.
+  [doi: 10.1145/2185520.2185558](https://doi.org/10.1145/2185520.2185558)
+- Stefan Band, Christoph Gissler, Andreas Peer, and Matthias Teschner.
+  "MLS Pressure Boundaries for Divergence-Free and Viscous SPH Fluids."
+  In: Computers & Graphics 76 (2018), pages 37–46.
+  [doi: 10.1016/j.cag.2018.08.001](https://doi.org/10.1016/j.cag.2018.08.001)
+"""
+struct PressureMirroring end
+
+@doc raw"""
+    PressureZeroing()
+
+This is the simplest way to implement dummy boundary particles.
+The density of each particle is set to the reference density and the pressure to the
+reference pressure (the corresponding pressure to the reference density by the state equation).
+
+!! note
+    This boundary model produces significantly worse results than all other models and
+    is only included for research purposes.
+"""
+struct PressureZeroing end
+
+# For most density calculators, the pressure is updated in every step
+initial_boundary_pressure(initial_density, density_calculator, _) = similar(initial_density)
+# Pressure mirroring does not use the pressure, so we set it to zero for the visualization
+initial_boundary_pressure(initial_density, ::PressureMirroring, _) = zero(initial_density)
+
+# For pressure zeroing, set the pressure to the reference pressure (zero with free surfaces)
+function initial_boundary_pressure(initial_density, ::PressureZeroing, state_equation)
+    return state_equation.(initial_density)
+end
+
+# With EDAC, just use zero pressure
+function initial_boundary_pressure(initial_density, ::PressureZeroing, ::Nothing)
+    return zero(initial_density)
+end
+
+@inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system,
+                                       v_particle_system, boundary_particle,
+                                       boundary_system,
+                                       v_boundary_system,
+                                       boundary_model::BoundaryModelDummyParticles,
+                                       rho_a, rho_b, pos_diff, distance, grad_kernel,
+                                       fluid_density_calculator)
+    (; density_calculator) = boundary_model
+
+    pressure_acceleration(pressure_correction, m_b, particle, particle_system,
+                          v_particle_system, boundary_particle,
+                          boundary_system, v_boundary_system,
+                          boundary_model, density_calculator,
+                          rho_a, rho_b, pos_diff, distance, grad_kernel,
+                          fluid_density_calculator)
+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.
+@inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system,
+                                       v_particle_system, boundary_particle,
+                                       boundary_system,
+                                       v_boundary_system,
+                                       boundary_model::BoundaryModelDummyParticles,
+                                       boundary_density_calculator,
+                                       rho_a, rho_b, pos_diff, distance, grad_kernel,
+                                       fluid_density_calculator::ContinuityDensity)
+    rho_a = particle_density(v_particle_system, particle_system, particle)
+    rho_b = particle_density(v_boundary_system, boundary_system, boundary_particle)
+
+    return -m_b *
+           (particle_system.pressure[particle] + boundary_model.pressure[boundary_particle]) /
+           (rho_a * rho_b) * grad_kernel
+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 SummationDensity.
+@inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system,
+                                       v_particle_system, boundary_particle,
+                                       boundary_system,
+                                       v_boundary_system,
+                                       boundary_model::BoundaryModelDummyParticles,
+                                       boundary_density_calculator,
+                                       rho_a, rho_b, pos_diff, distance, grad_kernel,
+                                       fluid_density_calculator::SummationDensity)
+    rho_a = particle_density(v_particle_system, particle_system, particle)
+    rho_b = particle_density(v_boundary_system, boundary_system, boundary_particle)
+
+    return -m_b *
+           (particle_system.pressure[particle] / rho_a^2 +
+            boundary_model.pressure[boundary_particle] / rho_b^2) *
+           grad_kernel
+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.
+@inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system,
+                                       v_particle_system, boundary_particle,
+                                       boundary_system,
+                                       v_boundary_system,
+                                       boundary_model::BoundaryModelDummyParticles,
+                                       ::PressureMirroring,
+                                       rho_a, rho_b, pos_diff, distance, grad_kernel,
+                                       fluid_density_calculator::ContinuityDensity)
+    rho_a = particle_density(v_particle_system, particle_system, particle)
+    rho_b = particle_density(v_boundary_system, boundary_system, boundary_particle)
+
+    return -m_b *
+           (particle_system.pressure[particle] + particle_system.pressure[particle]) /
+           (rho_a * rho_b) * grad_kernel
+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 SummationDensity.
+@inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system,
+                                       v_particle_system, boundary_particle,
+                                       boundary_system,
+                                       v_boundary_system,
+                                       boundary_model::BoundaryModelDummyParticles,
+                                       ::PressureMirroring,
+                                       rho_a, rho_b, pos_diff, distance, grad_kernel,
+                                       fluid_density_calculator::SummationDensity)
+    rho_a = particle_density(v_particle_system, particle_system, particle)
+    rho_b = particle_density(v_boundary_system, boundary_system, boundary_particle)
+
+    return -m_b *
+           (particle_system.pressure[particle] / rho_a^2 +
+            particle_system.pressure[particle] / rho_b^2) *
+           grad_kernel
+end
+
+function create_cache(initial_density,
+                      ::Union{SummationDensity, PressureMirroring, PressureZeroing})
     density = copy(initial_density)
 
     return (; density)
@@ -185,7 +344,10 @@ end
 end
 
 # Note that the other density calculators are dispatched in `density_calculators.jl`
-@inline function particle_density(v, ::AdamiPressureExtrapolation, boundary_model, particle)
+@inline function particle_density(v,
+                                  ::Union{AdamiPressureExtrapolation, PressureMirroring,
+                                          PressureZeroing},
+                                  boundary_model, particle)
     (; cache) = boundary_model
 
     return cache.density[particle]
@@ -206,9 +368,10 @@ end
 end
 
 function compute_density!(boundary_model,
-                          ::Union{ContinuityDensity, AdamiPressureExtrapolation},
+                          ::Union{ContinuityDensity, AdamiPressureExtrapolation,
+                                  PressureMirroring, PressureZeroing},
                           system, system_index, v, u, v_ode, u_ode, semi)
-    # No density update for `ContinuityDensity`.
+    # No density update for `ContinuityDensity`, `PressureMirroring` and `PressureZeroing`.
     # For `AdamiPressureExtrapolation`, the density is updated in `compute_pressure!`.
     return boundary_model
 end
@@ -293,6 +456,12 @@ function compute_pressure!(boundary_model, ::AdamiPressureExtrapolation,
     end
 end
 
+function compute_pressure!(boundary_model, ::Union{PressureMirroring, PressureZeroing},
+                           system, system_index, v, u, v_ode, u_ode, semi)
+    # No pressure update needed with `PressureMirroring` and `PressureZeroing`.
+    return boundary_model
+end
+
 @inline function adami_pressure_extrapolation!(boundary_model, system,
                                                neighbor_system::FluidSystem,
                                                system_coords, neighbor_coords,
@@ -382,5 +551,6 @@ end
 
 @inline function inverse_state_equation!(density, state_equation::Nothing, pressure,
                                          particle)
+    # The density is constant when using EDAC
     return density
 end

src/schemes/fluid/weakly_compressible_sph/rhs.jl

@@ -71,9 +71,9 @@ function interact!(dv, v_particle_system, u_particle_system,
     return dv
 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.
+# 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.
 @inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system,
                                        v_particle_system, neighbor,
                                        neighbor_system::WeaklyCompressibleSPHSystem,
@@ -85,9 +85,9 @@ end
            pressure_correction
 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 SummationDensity.
+# 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.
 @inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system,
                                        v_particle_system, neighbor,
                                        neighbor_system::WeaklyCompressibleSPHSystem,