Skip to content

Commit

Permalink
New boundary density calculators PressureZeroing and `PressureMirro…
Browse files Browse the repository at this point in the history 
…ring` (#223)

* Add `PressureZeroing` and `PressureMirroring` boundary methods

* Fix bug where damping was applied to the boundary density

* Reformat code

* Add docs for new boundary methods

* Fix typo

* Implement suggestions

* Use correct pressure acceleration formulation for `ContinuityDensity`

* Fix typo

---------

Co-authored-by: Niklas Neher <[email protected]>
  • Loading branch information
@efaulhaber @LasNikas
efaulhaber and LasNikas authored Oct 27, 2023
1 parent bdcb03e commit cb7a6d7
Show file tree
Hide file tree
Showing 3 changed files with 204 additions and 33 deletions.
3 changes: 2 additions & 1 deletion src/TrixiParticles.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
222 changes: 196 additions & 26 deletions src/schemes/boundary/dummy_particles/dummy_particles.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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)

Expand All @@ -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()
Expand All @@ -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)
Expand Down Expand Up @@ -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]
Expand All @@ -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
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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
12 changes: 6 additions & 6 deletions src/schemes/fluid/weakly_compressible_sph/rhs.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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,
Expand Down

0 comments on commit cb7a6d7

Please sign in to comment.