WCSPH algorithm:
This example shows the WCSPH method introduced by Becker and Teschner [BT07].
 perform neighborhood search
 compute particle densities $\rho_i$
 compute particle pressure values $p_i$ using an equation of state
 compute pressure forces
 compute viscosity (XSPH)
 time integration
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. Equation of State
Pressure is determined locally by the density:
$$p_i = \frac{\kappa \rho_0}{\gamma} \left ( \left ( \frac{\rho_i}{\rho_0} \right )^\gamma  1 \right ),$$
where he constants $\kappa, \gamma$ define the stiffness.
4. Pressure force
The pressure force of a particle is determined by a symmetric SPH formulation to preserves linear and angular momentum:
$$\mathbf f_i =  m_i \sum_j m_j \left ( \frac{p_i}{\rho_i^2} + \frac{p_j}{\rho_j^2}\right ) \nabla W_{ij}.$$
5. Viscosity (XSPH)
Artificial viscosity is introduced by smoothing the velocity field as:
$$\hat {\mathbf v}_i = \mathbf v_i + \alpha \sum_j \frac{m_j}{\rho_j} (\mathbf v_j  \mathbf v_i) W_{ij}.$$
6. Time integration
Finally, the particles are advected using a symplectic Euler method:
$$\begin{align*}
\mathbf v(t + \Delta t) &= \mathbf v(t) + \Delta t \mathbf a(t) \\
\mathbf x(t + \Delta t) &= \mathbf x(t) + \Delta t \mathbf v(t + \Delta t),
\end{align*}$$
where $\mathbf a$ are the accelerations corresponding to the sum of nonpressure and pressure forces.
Boundary handling
In order to implement a rigidfluid 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 pseudomass $\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
 [BT07] Markus Becker and Matthias Teschner. Weakly compressible SPH for free surface flows. In Proceedings of ACM SIGGRAPH/Eurographics Symposium on Computer Animation, 2007. Eurographics Association.
 [AIA+12] Nadir Akinci, Markus Ihmsen, Gizem Akinci, Barbara Solenthaler, and Matthias Teschner, "Versatile rigidfluid coupling for incompressible SPH", ACM Transactions on Graphics 31(4), 2012
