diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 583dab898..1b1c88c5e 100644 --- a/src/TrixiParticles.jl +++ b/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 diff --git a/src/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl index 747b55b45..6f8283dfd 100644 --- a/src/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/src/schemes/boundary/dummy_particles/dummy_particles.jl @@ -33,10 +33,13 @@ 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 @@ -44,8 +47,14 @@ We provide three options to compute the boundary density and pressure, determine 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 diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index ef5c2bbf3..b01aab34b 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/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,