Predictive-Corrective Incompressible SPH

B. Solenthaler ∗ R. Pajarola †
University of Zurich University of Zurich

Figure 1: Three examples produced with our incompressible simulation: (Left) 2M particles splashing against the simulation boundaries.
(Center) Close-up view of a wave tank. (Right) A fluid represented by 700k particles colliding with cylinder obstacles.

Abstract process and thus renders particle methods less attractive for high
quality and photorealistic water animations. In the context of
We present a novel, incompressible fluid simulation method based Smoothed Particle Hydrodynamics (SPH), two different strategies
on the Lagrangian Smoothed Particle Hydrodynamics (SPH) model. have been pursued to model incompressibility. First, the weakly
In our method, incompressibility is enforced by using a prediction- compressible SPH (WCSPH) method has been used where pressure
correction scheme to determine the particle pressures. For this, is modeled using a stiff equation of state (EOS), and second, incom-
the information about density fluctuations is actively propagated pressibility has been achieved by solving a pressure Poisson equa-
through the fluid and pressure values are updated until the targeted tion. Although both methods satisfy incompressibility, the compu-
density is satisfied. With this approach, we avoid the computational tational expenses of simulating high resolution fluid animations are
expenses of solving a pressure Poisson equation, while still being too large for practical use.
able to use large time steps in the simulation. The achieved results
show that our predictive-corrective incompressible SPH (PCISPH) In the standard SPH and WCSPH model the particle pressures are
method clearly outperforms the commonly used weakly compress- determined by an EOS. The characteristics of this equation and
ible SPH (WCSPH) model by more than an order of magnitude the stiffness parameter determine the speed of the acoustic waves
while the computations are in good agreement with the WCSPH in a medium. The EOS-based SPH with low stiffness accord-
results. ing to [Desbrun and Cani 1996] was used in a series of papers to
simulate water [Müller et al. 2003; Adams et al. 2007], multiple
CR Categories: I.3.5 [Computer Graphics]: Computational Ge- fluids [Müller et al. 2005; Solenthaler and Pajarola 2008], fluid-
ometry and Object Modeling—Physically based modeling; I.3.7 solid coupling [Müller et al. 2004b; Lenaerts et al. 2008], melting
[Computer Graphics]: Three-Dimensional Graphics and Realism— solids [Müller et al. 2004a; Keiser et al. 2005; Solenthaler et al.
Animation. 2007], and fluid control [Thürey et al. 2006]. In contrast to the
standard SPH formulation, WCSPH uses a stiff EOS [Monaghan
2005; Becker and Teschner 2007; Becker et al. 2009] resulting
Keywords: fluid simulation, SPH, incompressibility
in acoustic waves traveling closer to their real speed through the
medium. Typically, the stiffness value is chosen so large that the
1 Introduction and Previous Work density fluctuations do not exceed 1%. The required stiffness value
to achieve this, however, is difficult or even impossible to deter-
Enforcing incompressibility in fully particle-based fluid simula- mine before running the simulation. Consequently, an animator
tions represents the most expensive part of the whole simulation cannot get around extensive testing and parameter tuning. Another
∗ e-mail:
drawback is that WCSPH imposes a severe time step restriction as
[email protected] the stiffness of the fluid usually dominates the Courant-Friedrichs-
† e-mail: [email protected] Levy (CFL) condition. Thus the computational cost increases with
decreasing compressibility – since higher stiffness requires smaller
time steps, making it infeasible to simulate high resolution fluids
within reasonable time.
Rather than simulating acoustic waves, incompressibility in La-
grangian methods can be enforced by solving a pressure projection
similar to Eulerian methods (e.g. [Enright et al. 2002]). These in-
compressible SPH (ISPH) methods first integrate the velocity field
in time without enforcing incompressibility. Then, either the inter-
mediate velocity field [Cummins and Rudman 1999], the resulting
variation in particle density [Shao 2006], or both [J. Liu and Oka
2005; Hu and Adams 2007; Losasso et al. 2008] are projected onto Algorithm 1 SPH / WCSPH
a divergence-free space to satisfy incompressibility through a pres- 1 while animating do
sure Poisson equation. With these ISPH methods density fluctua- 2 for all i do
tions of 1% to 3% have been reported. A problem with these meth- 3 find neighborhoods Ni (t)
ods, however, is the complexity to formulate and solve the equation 4 for all i do
system on unstructured particle configurations. Although ISPH al- 5 compute density ρi (t)
lows larger time steps than WCSPH, the computational cost per 6 compute pressure pi (t)
physics step is much higher. A Poisson solver was also used in [Pre- 7 for all i do
moze et al. 2003] for the particle method Moving-Particle Semi- 8 compute forces Fp,v,g,ext (t)
Implicit (MPS), increasing the cost per physics time step enor- 9 for all i do
mously. In contrast to the fully Lagrangian models, [Zhu and Brid- 10 compute new velocity vi (t + 1)
son 2005] propose to use an auxiliary background grid to simplify 11 compute new position xi (t + 1)
the equation system to a sparse set of linear equations which can be
efficiently solved. A similar hybrid solver for vorticity confinement
is presented in [Selle et al. 2005]. In [Losasso et al. 2008], a two- Algorithm 2 PCISPH
way coupled level set method with an SPH solver is introduced to 1 while animating do
simulate dense and diffuse water volumes. They demonstrate how 2 for all i do
to enforce incompressibility and target the particle number density 3 find neighborhoods Ni (t)
with a single Poisson solve. 4 for all i do
5 compute forces Fv,g,ext (t)
In this paper, we propose a novel, fully Lagrangian, incompressible 6 initialize pressure p(t) = 0.0
SPH method featuring the advantages of both WCSPH and ISPH in 7 initialize pressure force Fp (t) = 0.0
one model, namely low computational cost per physics update and 8 while (ρ∗err (t + 1) > η) || (iter < minIterations) do
large time steps. Our method makes use of a prediction-correction 9 for all i do
scheme which propagates the estimated density values through the 10 predict velocity vi∗ (t + 1)
fluid and updates the pressures in such a way that incompressibility 11 predict position x∗i (t + 1)
is achieved. The propagation stops as soon as a previously user- 12 for all i do
defined density variation limit is reached for each individual parti- 13 predict density ρ∗i (t + 1)
cle. We will show in this paper that our new predictive-corrective 14 predict density variation ρ∗err (t + 1)
incompressible SPH (PCISPH) method outperforms WCSPH by 15 update pressure pi (t)+= f (ρ∗err (t + 1))
more than an order of magnitude while the computations are in 16 for all i do
good agreement with the WCSPH results. The efficiency of our 17 compute pressure force Fp (t)
method enables an animator to produce high-resolution fluid ani- 18 for all i do
mations within reasonable time without compressibility artifacts. 19 compute new velocity vi (t + 1)
20 compute new position xi (t + 1)
2 PCISPH Model

2.1 Basic SPH / WCSPH Algorithm equations and given by

X pi pj
In SPH, liquids are approximated by artificial, slightly compress- Fpi = −m2 ( 2 + 2 )∇W (xij , h).
ible fluids. The basic SPH method is summarized in Algorithm 1. j
ρ i ρ j
In each physics update, the local neighborhood Ni of each particle
i is found and then used to evaluate the density, pressure, and the In our implementation, we use the viscous force and weighting ker-
resulting forces acting on each particle [Monaghan 1992]. The den- nels presented in [Monaghan 1992].
sity ρi of particle i at location xi can be found by summing up the
weighted contributions of the neighboring particles j 2.2 PCISPH Algorithm
X
ρi = m W (xij , h), (1) To avoid the time step restriction of WCSPH we propose to
j use a prediction-correction scheme based on the SPH algorithm
(PCISPH). In our method, the velocities and positions are tem-
where m is the particle mass (we assume that all particles have porarily forwarded in time and the new particle densities are es-
equal masses), W is the weighting kernel with smoothing length h, timated. Then, for each particle, the predicted variation from the
and xij = xi − xj . The pressure pi of a particle is then derived reference density is computed and used to update the pressure val-
from the EOS according to [Batchelor 1967] ues, which in turn enter the recomputation of the pressure forces.
Similar to a Jacobi iteration for linear systems, this process is it-
kρ0 ρi γ erated until it converges, i.e. until all particle density fluctuations
pi = (( ) − 1), are smaller than a user-defined threshold η (for example 1%). Note
γ ρ0
that this is a nonlinear problem since we include collision handling
where k is a stiffness parameter and ρ0 the reference density. [Des- and updated kernel values in our iteration process. As a final step,
brun and Cani 1996] use a γ of 1 and a small value for k, while in the velocities and positions of the next physics update step are com-
WCSPH (e.g. [Monaghan 2005]) γ is set to 7 and k is chosen so puted. The PCISPH method is illustrated in Algorithm 2.
that the speed of sound is large enough to keep the density fluctu-
ations small (∼1%). Note that the CFL condition [Courant et al. 2.3 Pressure Derivation
1967] requires smaller time steps for stiffer fluids which increases
the overall computation cost tremendously when simulating water. One of the main difficulties is to derive the pressure change from the
The pressure force field is directly derived from the Navier-Stokes predicted density variation (line 15 of Algorithm 2). This pressure
update is executed in each iteration, reducing the density fluctua- and the position of j changes by
tion of the particle. The aim is to find a pressure p which changes
the particle positions in such a way that the predicted density corre- 2p̃i
∆xj|i = ∆t2 m ∇Wij . (6)
sponds to the reference density. Over the course of this section, a set ρ20
of approximations will be made to derive a simple update rule for
the pressure (Equations 8 to 10). Although the approximations in- Note that we only consider the effect of the central particle i here,
crease the number of convergence iterations which are needed until i.e. ∆xj = ∆xj|i . Equation 5 and Equation 6 can now be inserted
the desired density fluctuation limit is reached, they keep the final into Equation 2 resulting in
pressure update rule simple and thus efficient to compute. “ 2p̃i X X
∆ρi (t) = m −∆t2 m 2 ∇Wij · ∇Wij −
For a given kernel smoothing length h, the density at a point in ρ0 j j
time t + 1 is computed using the SPH density summation equation X 2p̃i ”
analogously to Equation 1 (∇Wij · ∆t2 m 2 ∇Wij )
X j
ρ0
ρi (t + 1) = m W (xi (t + 1) − xj (t + 1)) 2p̃ i
“ X X
j = ∆t2 m2 2 − ∇Wij · ∇Wij −
X ρ0 j j
= m W (xi (t) + ∆xi (t) − xj (t) − ∆xj (t)) X ”
j (∇Wij · ∇Wij )
X j
= m W (dij (t) + ∆dij (t))
j
After solving for p̃i we get
where dij (t) = xi (t) − xj (t), and ∆dij (t) = ∆xi (t) − ∆xj (t).
Assuming that ∆dij is relatively small, the first order Taylor ap- ∆ρi (t)
p̃i = P P P (7)
proximation can be applied to the term W (dij (t) + ∆dij (t)) re- β(− j ∇Wij · j ∇Wij − j (∇Wij · ∇Wij ))
sulting in
X where β is
ρi (t + 1) = m W (dij (t)) + ∇W (dij (t)) · ∆dij (t) 2
β = ∆t2 m2 .
j ρ20
X
= m W (xi (t) − xj (t)) + The meaning of Equation 7 is that a pressure p̃i is needed to achieve
j a change in density of ∆ρi (t). As we know the predicted density
error ρ∗erri = ρ∗i − ρ0 of a particle, we can thus reverse that error
X
m ∇W (xi (t) − xj (t)) · (∆xi (t) − ∆xj (t))
j by applying a pressure of
= ρi (t) + ∆ρi (t).
−ρ∗erri
In this equation, the term ∆ρi (t) is unknown and, as we show later, p̃i = P P P .
β(− j ∇Wij · j ∇Wij − j (∇Wij · ∇Wij ))
a function of p which we are looking for. After reformulation and
using Wij = W (xi (t) − xj (t)) we get This formula shows problems in situations where i is suffering from
X particle deficiency in the neighborhood resulting in falsified values.
∆ρi (t) = m ∇Wij · (∆xi (t) − ∆xj (t)) To circumvent that problem, we precompute a single scaling factor
j δ according to the following formula which is evaluated for a pro-
“X X ”
= m ∇Wij ∆xi (t) − ∇Wij ∆xj (t) totype particle with a filled neighborhood. The resulting value is
then used for all particles. Finally, we end up with the following
“ j j
equations which are used in the PCISPH method
X X ”
= m ∆xi (t) ∇Wij − ∇Wij ∆xj (t) (2)
j j −1
δ= P P P (8)
∆x can be derived from the time integration scheme (Leap-Frog). β(− j ∇Wij · j ∇W ij − j (∇Wij · ∇Wij ))
Neglecting all forces but the pressure force we get
p
and
∆xi = ∆t .2 Fi
(3) p̃i = δρ∗erri . (9)
m
Since we repeat the prediction-correction step as long as the incom-
If we make the simplistic assumption that neighbors have equal pressibility condition is not yet satisfied, the correction pressures of
pressures p˜i and that the density corresponds to the rest density ρ0 the individual iterations are accumulated as indicated on line 15 of
(according to the incompressibility condition), this results in Algorithm 2
X p̃i p̃i 2p̃i X pi += p̃i . (10)
Fpi = −m2 ( 2 + 2 )∇Wij = −m2 2 ∇Wij . (4)
j
ρ 0 ρ 0 ρ0 j
2.4 Implementation
Inserting Equation 4 into Equation 3 we get
2.4.1 Neighborhood Approximation
2p̃i X
∆xi = −∆t2 m 2 ∇Wij . (5)
ρ0 j Before predicting the density ρ∗i (t + 1) of a particle (line 13 of
Algorithm 2), the neighborhood should be recomputed using the
Due to the pressure pi of particle i the position of a neighboring predicted positions x∗ (t + 1). However, for efficiency reasons we
particle changes by ∆xj|i . As the pressure forces are symmetric, reuse the current neighbors Ni (t) at time t and only recompute the
particle j gets the following contribution from i distances and the kernel values. This approximation leads to small
p̃i p̃i 2p̃i errors in the density and pressure estimates. In the case of density
Fpj|i = m2 ( + 2 )∇Wij = m2 2 ∇Wij , overestimation the final real densities show lower fluctuations than
ρ20 ρ0 ρ0
the requested threshold η. In the opposite case – density underes-
timation – the correction loop might be aborted prematurely. Such
situations are not yet handled in the current implementation but can
be avoided by using sufficiently small time steps, or by recomputing
the neighborhoods in these particular situations.

2.4.2 Information Propagation

To limit temporal fluctuations in the resulting pressure field we

found it advantageous to employ a minimum number of iterations
in the pressure update loop. This gives the particles enough time to
propagate information about predicted particle locations. We found
a minimum of 3 iterations generally sufficient to achieve a low level
of pressure fluctuations. (a) Average number of convergence iterations over time. After
8s of simulated real time, an average of 3.49 is reached.
3 Results 16 η
at t=0.2s
at t=0.4s
14
at t=0.6s
3.1 Performance Comparison at t=0.8s
12 at t=1.0s

We set up a test scene (Figure 2) to compare the simulation times 10

ρ_err [%]
and visual results of both the commonly used WCSPH and our new 8
PCISPH method. The performance measurements and simulation
6
data are summarized in Table 1. All timings are given for an Intel
Core2 2.66 GHz CPU. 4

2
We executed different simulation runs with varying particle resolu-
tions (10k and 100k) and varying error threshold η (1% and 0.1%) 0
1 2 3 4 5 6 7
which defines the maximally allowed density fluctuation from the iterations
reference density. The 10k and 100k examples have corresponding (b) Several convergence examples at different points in time t.
scene setups but different fluid discretizations, meaning that a par-
ticle in the 10k example represents a larger fluid volume than one
in the 100k example. Since in SPH a particle always needs to have Figure 3: Convergence statistics of the 100k particles simulation
around 30-40 neighbors, the support radius has to be increased with shown in Table 1 and Figure 2.
increasing particle volume, which in turn influences the time step
size. The time step is set according to a CFL condition where the
force terms, the stiffness parameter k, and the viscous term are in- The peaks indicate particle collisions with the ground and the side
volved [Monaghan 1992]. While in WCSPH the time step is domi- walls as in such situations larger density errors are predicted. Fig-
nated by k, it has no influence in PCISPH and can be omitted. Thus, ure 3(b) shows several examples of the convergence within a single
for low viscosity fluids, the time step in PCISPH is dominated by physics update step. It can be seen that the density error is approx-
the force terms, allowing significantly larger time steps than those imately halved after the first iteration and continuously reduced in
used in WCSPH. The stiffness parameter k of WCSPH was de- the following iterations until the error drops below η. In our expe-
termined by testing in such a way that η was satisfied, which was rience this algorithm proved to be very robust and we did not en-
k = 7 · 104 for η = 1%, and k = 6 · 106 for η = 0.1%. In contrast, counter any divergence problems. However, it is likely that certain
PCISPH does not have to cope with finding an appropriate stiffness particle configurations exist that might show such problems.
value since the desired η can be specified directly. In the case of
η=1%, the time step of PCISPH is determined to be in fact 35 times 3.3 Visual Result
larger than the one of WCSPH. With a smaller η the difference
is even larger, η = 0.1% leads to an increase of the time step for The physical behavior and visual results of WCSPH and PCISPH
PCISPH by a factor of 151. While in WCSPH the computation time are compared in Figure 2. It can be seen that the PCISPH compu-
per simulation step stays more or less constant, it varies in PCISPH tations are in full agreement with the WCSPH results with only
since the time per simulation step depends on the number of ex- very minor detail differences. The comparison of WCSPH and
ecuted convergence iterations. Therefore, we compare the overall PCISPH using a simulation time constraint of 298min is shown in
computation time of WCSPH and PCISPH over the entire simu- Figure 4. While with PCISPH a resolution of 100k particles can
lated time period. Although the cost per physics time step is higher be simulated within the given time (see the corresponding entry in
with PCISPH than WCSPH, the overall speed-up over WCSPH still Table 1), with WCSPH the resolution has to be reduced to 17k par-
reaches a factor of 15 and 16 for η = 1% and 55 for η = 0.1%, ticles. The lower resolution leads to less surface details and notably
respectively. damped fluid movement. A higher resolution example computed
with PCISPH is shown in Figure 1 (center and right) and Figure 5
3.2 Convergence Analysis where a wave generator agitates a water body consisting of 700k
particles to interact with cylindrical obstacles in a tank. In Figure 1
In the previously described test scenes, the average number of con- (left) and Figure 6, 2M particles are used to simulate the collapsing
vergence iterations executed per physics step is between 3.24 and column example with PCISPH. In both of these simulations, a η of
4.46. Note that the particle resolution has no effect on the average 1% is enforced which eliminates compression artifacts and enables
number of iterations. For the simulation run with 100k particles realistic wave breaking and splashing behavior. In all examples, the
the average number of iterations is plotted over time in Figure 3(a). surface of the fluid is reconstructed and rendered with the raytracing
The end time of 8s corresponds to the simulated real time. approach presented in [Solenthaler et al. 2007].
Model η [%] #p k ∆t [s] ∆t ratio [s] avgIterations tsim [min] speed-up
WCSPH 1.0 10k 7·104 3.78e-5 - - 142.05 -
PCISPH 1.0 10k - 0.0013 35 3.24 9.37 15.2
WCSPH 1.0 100k 7·104 1.78e-5 - - 4941.5 -
PCISPH 1.0 100k - 0.00062 35 3.49 297.7 16.6
WCSPH 0.1 10k 6·106 4.08e-6 - - 1327.66 -
PCISPH 0.1 10k - 0.00062 151.96 4.46 23.97 55.39

Table 1: Comparison of WCSPH and PCISPH. The stiffness value k of WCSPH is chosen so that the density fluctuation percentage is below
η, and the time step size is determined according to the CFL condition. With our PCISPH method and η=1%, a speed-up of a factor of 15
and 16 over WCSPH is reached. By restricting the error to η=0.1%, PCISPH reduces the computation time by a factor of 55.

Figure 2: Side-by-side comparison of a fluid discretized by 100k particles and simulated with WCSPH (upper row) and PCISPH (lower row),
respectively. The computations correspond to the statistics given in Table 1 for the 100k particles simulation.

proximation which can lead to underestimated density errors abort-

ing the convergence loop prematurely as we have discussed in Sec-
tion 2.4.1. This problem can be addressed by detecting such sit-
uations and adapting the time step size or recomputing the neigh-
bors in this particular simulation step. Besides that, our current im-
plementation does not yet account for the particle deficiency near
boundaries. In these situations, density values are falsified and
compression artifacts can occur. The inclusion of ghost particles
in the density computation or the use of the corrected SPH formu-
lation as described in [Becker et al. 2009] can solve this problem.
Figure 4: Comparison of WCSPH (left, 17k particles) and PCISPH
(right, 100k particles) with equal computation times.
References

4 Conclusion A DAMS , B., PAULY, M., K EISER , R., AND G UIBAS , L. J. 2007.
Adaptively sampled particle fluids. ACM Trans. Graph. 26, 3,
We proposed a novel incompressible SPH solver which combines 48–54.
the advantages of both WCSPH and ISPH in one model, namely BATCHELOR , G. 1967. An Introduction to Fluid Dynamics. Cam-
low computational cost per physics update and large time steps. bridge University Press.
Our method includes a convergence loop which is executed in each
physics update step consisting of a prediction and correction itera- B ECKER , M., AND T ESCHNER , M. 2007. Weakly compressible
tion. In each convergence iteration, the new particle positions and SPH for free surface flows. In Symposium on Computer Anima-
their densities are predicted and the variations from the reference tion, 209–217.
density are computed. We have derived a formulation which relates
the density fluctuation and the pressure, to reduce the density errors B ECKER , M., T ESSENDORF, H., AND T ESCHNER , M. 2009. Di-
and to approach incompressibility. With this method, we gained a rect forcing for Lagrangian rigid-fluid coupling. IEEE Transac-
speed-up of more than an order of magnitude over the commonly tions on Visualization and Computer Graphics 15, 3, 493–503.
used WCSPH method and we showed that the simulation results are
in agreement with WCSPH. C OURANT, R., F RIEDRICHS , K., AND L EWY, H. 1967. On the
partial difference equations of mathematical physics. IBM J. 11,
One issue of the current implementation is the neighborhood ap- 215–234.
 Figure 5: Wave breaking and splashing in a wave tank simulated with the proposed incompressible PCISPH method.

L OSASSO , F., TALTON , J., K WATRA , J., AND F EDKIW, R. 2008.

Two-way coupled SPH and particle level set fluid simulation.
IEEE TVCG 14, 4, 797–804.
M ONAGHAN , J. 1992. Smoothed particle hydrodynamics. Annu.
Rev. Astron. Physics 30, 543.
M ONAGHAN , J. 2005. Smoothed particle hydrodynamics. Rep.
Prog. Phys. 68, 1703–1759.
M ÜLLER , M., C HARYPAR , D., AND G ROSS , M. 2003. Particle-
based fluid simulation for interactive applications. In Symposium
on Computer Animation, 154–159.
M ÜLLER , M., K EISER , R., N EALEN , A., PAULY, M., G ROSS ,
M., AND A LEXA , M. 2004. Point based animation of elastic,
plastic and melting objects. In Symposium on Computer Anima-
tion, 141–151.
M ÜLLER , M., S CHIRM , S., T ESCHNER , M., H EIDELBERGER ,
B., AND G ROSS , M. 2004. Interaction of fluids with deformable
solids. Journal of Computer Animation and Virtual Worlds 15,
Figure 6: PCISPH simulation of 2M particles interacting with
3-4, 159–171.
cylinder obstacles.
M ÜLLER , M., S OLENTHALER , B., K EISER , R., AND G ROSS , M.
2005. Particle-based fluid-fluid interaction. In Symposium on
Computer Animation, 237–244.
C UMMINS , S. J., AND RUDMAN , M. 1999. An SPH projection
method. J. Comput. Phys. 152, 2, 584–607. P REMOZE , S., TASDIZEN , T., B IGLER , J., L EFOHN , A., AND
W HITAKER , R. T. 2003. Particle-based simulation of fluids. In
D ESBRUN , M., AND C ANI , M.-P. 1996. Smoothed particles: A Proceedings of Eurographics, 401–410.
new paradigm for animating highly deformable bodies. In Eu-
rographics Workshop on Computer Animation and Simulation, S ELLE , A., R ASMUSSEN , N., AND F EDKIW, R. 2005. A vortex
61–76. particle method for smoke, water and explosions. ACM Trans.
Graph. 24, 3, 910–914.
E NRIGHT, D., M ARSCHNER , S., AND F EDKIW, R. 2002. Ani- S HAO , S. 2006. Incompressible SPH simulation of wave breaking
mation and rendering of complex water surfaces. ACM Trans. and overtopping with turbulence modelling. Int. J. Numer. Meth.
Graph. 21, 3, 736–744. Fluids 50, 597–621.

H U , X. Y., AND A DAMS , N. A. 2007. An incompressible multi- S OLENTHALER , B., AND PAJAROLA , R. 2008. Density contrast
phase SPH method. J. Comput. Phys. 227, 1, 264–278. SPH interfaces. In Symposium on Computer Animation, 211–
218.
J. L IU , S. K., AND O KA , Y. 2005. A hybrid particle-mesh method S OLENTHALER , B., S CHL ÄFLI , J., AND PAJAROLA , R. 2007.
for viscous, incompressible, multiphase flows. J. Comput. Phys. A unified particle model for fluid-solid interactions. Journal of
202, 1, 65–93. Computer Animation and Virtual Worlds 18, 1, 69–82.
K EISER , R., A DAMS , B., G ASSER , D., BAZZI , P., D UTRE , P., T H ÜREY, N., K EISER , R., PAULY, M., AND R ÜDE , U. 2006.
AND G ROSS , M. 2005. A unified Lagrangian approach to solid- Detail-preserving fluid control. In Symposium on Computer An-
fluid animation. In Proceedings of Eurographics Symposium on imation, 7–15.
Point-Based Graphics, 125–133.
Z HU , Y., AND B RIDSON , R. 2005. Animating sand as a fluid.
ACM Trans. Graph. 24, 3, 965–972.
L ENAERTS , T., A DAMS , B., AND D UTR É , P. 2008. Porous flow in
particle-based fluid simulations. ACM Trans. Graph. 27, 3, 1–8.