IISPH algorithm:
This example shows the IISPH method introduced by Ihmsen et al. [ICS+14]:
- perform neighborhood search
- compute particle densities $\rho_i$
- compute diagonal matrix elements $a_{ii}$
- compute viscosity (XSPH)
- integrate velocities only considering non-pressure forces
- compute source term
- loop
- compute pressure accelerations
- determine $\mathbf{A} \mathbf{p}$
- update pressure values
- time integration with pressure forces
IISPH solves the pressure Poisson equation (PPE):
$$\Delta t \nabla^2 p = \frac{\rho_0 - \rho^*}{\Delta t},$$
where $\rho^*$ is the predicted density after applying all non-pressure forces.
Discretizing the PPE using the SPH formulation yields a linear system
$$\mathbf{A} \mathbf{p} = \mathbf{s}$$
for the unknown pressure values $\mathbf{p}$ and the source term $\mathbf{s}$.
1. Neighborhood search
The neighborhood search is performed using spatial hashing. In each step all particles are added to a spatial grid using a cell size that equals the support radius of the SPH kernel function. Hence, the neighbor particles of a particle in cell (i,j) lie in one of the 9 neighboring cells: (i±1, j±1).
2. SPH density computation
Using the SPH formulation the density is determined as:
$$\rho_i = \sum_j \frac{m_j}{\rho_j} \rho_j W_{ij} = \sum_j m_j W_{ij}.$$
Hence, the density only depends on the masses and positions of the particles.
3. Compute diagonal matrix elements $a_{ii}$
The diagonal elements of the system matrix $\mathbf{A}$ are determined as
$$a_{ii} = -\frac{\Delta t}{\rho_i^2} \left (\sum_j \left \| m_j \nabla W_{ij}\right \|^2 + \left \| \sum_j m_j \nabla W_{ij}\right \|^2 \right ).$$
4. Viscosity (XSPH)
Artificial viscosity is introduced by smoothing the velocity field as:
$$\mathbf {a}^{\text{visco}}_i = \frac{\alpha}{\Delta t} \sum_j \frac{m_j}{\rho_j} (\mathbf v_j - \mathbf v_i) W_{ij}.$$
5. Integrate velocities only considering non-pressure forces
The particles velocities are updated using an Euler method:
$$\begin{align*}
\mathbf v^* &= \mathbf v(t) + \Delta t \mathbf a^\text{np}(t) \\
\end{align*}$$
where $\mathbf a^\text{np}$ are the accelerations corresponding to the sum of non-pressure forces (e.g. viscosity).
6. Compute source term
Determine the source term on the right hand side of the linear system:
$$s_i = \frac{\rho_0 - \rho_i^*}{\Delta t},$$
where
$$\rho_i^* = \rho_i + \Delta t \frac{D \rho_i}{D t} = \rho_i + \Delta t \sum_j m_j (\mathbf{v}^*_i - \mathbf{v}^*_j) \cdot \nabla W_{ij}$$
7. Solver loop
The PPE is solved using relaxed Jacobi iterations with a relaxation factor of $\omega = 0.5$ until the density error is smaller than a given tolerance or a maximum number of iterations is reached.
8. Compute pressure accelerations
The pressure acceleration of a particle is determined by a symmetric SPH formulation to preserves linear and angular momentum:
$$\mathbf a^\text{p}_i = - \sum_j m_j \left ( \frac{p_i}{\rho_i^2} + \frac{p_j}{\rho_j^2}\right ) \nabla W_{ij}.$$
9. Determine $\mathbf{A} \mathbf{p}$
Compute the product of the system matrix $\mathbf{A}$ and the current pressure vector:
$$(\mathbf{A} \mathbf{p})_i = \Delta t \sum_j m_j \left (\mathbf{a}^\text{p}_i - \mathbf{a}^\text{p}_j \right ) \cdot \nabla W_{ij}.$$
10. Update pressure values
In each solver iteration we update the pressure value for each particle $i$ as
$$p_i := p_i + \frac{\omega}{a_{ii}} \left (s_i - (\mathbf{A} \mathbf{p})_i \right ).$$
11. Time integration with pressure forces
Finally, the particles are advected using a symplectic Euler method:
$$\begin{align*}
\mathbf v(t + \Delta t) &= \mathbf v(t) + \Delta t \mathbf a^\text{p}(t) \\
\mathbf x(t + \Delta t) &= \mathbf x(t) + \Delta t \mathbf v(t + \Delta t),
\end{align*}$$
where $\mathbf a^\text{p}$ are the accelerations corresponding to the sum of pressure forces.
Boundary handling
In order to implement a rigid-fluid coupling the SPH equation for the density computation is extended by a second sum over the contributing neighboring boundary particles $k$ [AIA+12]:
$$\rho_i = \sum_j m_j W_{ij} + \sum_k \Psi_k W_{ik}.$$
Each boundary particle $k$ has a pseudo-mass $\Psi$ which considers the density of the boundary sampling points:
$$\Psi_k = \rho_0 V_k = \rho_0 \frac{m_k}{\rho_k} = \rho_0 \frac{m_k}{\sum_l m_k W_{kl}} = \rho_0 \frac{1}{\sum_l W_{kl}}$$
where $l$ denotes the boundary particle neighbors of particle $k$.
References
- [ICS+14] Markus Ihmsen, Jens Cornelis, Barbara Solenthaler, Christopher Horvath, Matthias Teschner. Implicit incompressible SPH. IEEE Transactions on Visualization and Computer Graphics, 20(3), 2014.
- [AIA+12] Nadir Akinci, Markus Ihmsen, Gizem Akinci, Barbara Solenthaler, and Matthias Teschner, "Versatile rigid-fluid coupling for incompressible SPH", ACM Transactions on Graphics 31(4), 2012
|