To appear in ACM TOG 32(4). Position Based Fluids Miles Macklin Matthias M ?uller   NVIDIA Abstract In fluid simulation, enforcing incompressibility is crucial for real- ism; it is also computationally expensive. Recent work has im- proved efficiency, but still requires time-steps that are impractical for real-time applications. In this work we present an iterative den- sity solver integrated into the Position Based Dynamics framework (PBD). By formulating and solving a set of positional constraints that enforce constant density, our method allows similar incom- pressibility and convergence to modern smoothed particle hydro- dynamic (SPH) solvers, but inherits the stability of the geometric, position based dynamics method, allowing large time steps suit- able for real-time applications. We incorporate an artificial pressure term that improves particle distribution, creates surface tension, and lowers the neighborhood requirements of traditional SPH. Finally, we address the issue of energy loss by applying vorticity confine- ment as a velocity post process. CR Categories: I.3.5 [Computer Graphics]: Computational Geometry and Object ModelingÑPhysically based modeling I.3.7 [Computer Graphics]: Three-Dimensional Graphics and RealismÑ Animation; Keywords: fluid simulation, SPH, PCISPH, constraint fluids, po- sition based dynamics 1 Introduction Fluids, in particular liquids such as water, are responsible for many visually rich phenomena, and simulating them has been an area of long-standing interest and challenge in computer graphics. There are a variety of techniques available, but here we focus on particle methods, which are popular for their simplicity and flexibility. Smoothed Particle Hydrodynamics (SPH) [Monaghan 1992][1994] is a well known particle based method for fluid simulation. It has many attractive properties: mass-conservation, Lagrangian dis- cretization (particularly useful in games where the simulation do- main is not necessarily known in advance), and conceptual simplic- ity. However, SPH is sensitive to density fluctuations from neigh- borhood deficiencies, and enforcing incompressibility is costly due to the unstructured nature of the model. Recent work has im- proved efficiency by an order of magnitude [Solenthaler and Pa- jarola 2009], but small time steps remain a requirement, limiting real-time applications. e-mail:mmacklin@nvidia.com   e-mail:matthiasm@nvidia.com (a) Real-time rendered fluid surface using ellipsoid splatting (b) Underlying simulation particles Figure 1: Bunny taking a bath. 128k particles, 2 sub-steps, 3 den- sity iterations per frame, average simulation time per frame 10ms. For interactive environments, robustness is a key issue: the simula- tion must handle degenerate situations gracefully. SPH algorithms often become unstable if particles do not have enough neighbors for accurate density estimates. The typical solution is to try to avoid these situations by taking sufficiently small time steps, or by using sufficiently many particles, at the cost of increased computation. In this paper, we show how incompressible flow can be simulated inside the Position Based Dynamics (PBD) framework [M ?uller et al. 2007]. We choose PBD for its unconditionally stable time integration and robustness, which has made it popular with game developers and film makers. By addressing the issue of particle deficiency at free surfaces, and handling large density errors, our method allows users to trade incompressibility for performance, while remaining stable. 2 Related Work M ?uller [2003] showed that SPH can be used for interactive fluid simulation with viscosity and surface tension, by using a low stiff- ness equation of state (EOS). However to maintain incompressibil- ity, standard SPH or weakly compressible SPH (WCSPH) [Becker and Teschner 2007] require stiff equations, resulting in large forces that limit the time-step size. Predictive-corrective incompressible SPH (PCISPH) [Solenthaler and Pajarola 2009] relaxes this time- step restriction by using an iterative Jacobi-style method that accu- mulates pressure changes and applies forces until convergence. It has the advantage of not requiring a user-specified stiffness value and of amortizing the cost of neighbor finding over many density 1 To appear in ACM TOG 32(4). corrections. Bodin et al [2012] achieve uniform density fluid by posing incom- pressibility as a system of velocity constraints. They construct a linear complementarity problem using linearized constraint func- tions, which are solved using Gauss-Seidel iteration. In contrast, our method (and PCISPH) attempts to solve the non-linear problem by operating on particles directly, and re-evaluating constraint error and gradients each Jacobi iteration. Hybrid methods, such as Fluid Implicit-Particle (FLIP) [Brackbill and Ruppel 1986] use a grid for the pressure solve and transfer ve- locity changes back to particles. FLIP was later extended to in- compressible flow with free surfaces by Zhu and Bridson [2005]. Raveendran et al. [2011] use a coarse grid to solve for an approx- imately divergence free velocity field before an adaptive SPH up- date. Clavet et al. [2005] also use a position based approach to simu- late viscoelastic fluids. However, because the time step appears in various places of their position projections, their approach is only conditionally stable as in regular explicit integration. Position Based Dynamics [M ?uller et al. 2007] provides a method for simulating dynamics in games based on Verlet integration. It solves a system of non-linear constraints using Gauss-Seidel itera- tion by updating particle positions directly. By eschewing forces, and deriving momentum changes implicitly from the position up- dates, the typical instabilities associated with explicit methods are avoided. 3 Enforcing Incompressibility Following [Bodin et al. 2012] the scalar density constraint on a par- ticle with index i and its set of neighboring particles p1 , á á á , pn is defined using an equation of state: C(p1 , ..., pn ) = i 0 1, (1) where i is given by the standard SPH density estimator i = j m j W (pi p j , h). (2) We treat all particles as having equal mass and will drop this term from subsequent equations. In our implementation we use the Poly6 kernel for density estimation, and the Spiky kernel for gradi- ent calculation, as in [M ?uller et al. 2003]. Now we give some background on the position based dynamics method and then show how to incorporate our density constraint. PBD aims to find a particle position correction that satisfies the constraint: C(p + 0 (3) This is found by a series of Newton steps along the constraint gra- dient: (4) C(p + C(p) + C T = 0 (5) C(p) + C T = 0. (6) [Monaghan 1992] gives the SPH recipe for the gradient of a func- tion defined on the particles. Applying this, the gradient of the Algorithm 1 Simulation Loop 1: for all particles i do 2: apply forces vi vi + fext (xi ) 3: predict position x i xi + vi 4: end for 5: for all particles i do 6: find neighboring particles Ni (x i ) 7: end for 8: while it er < sol verIt erat ions do 9: for all particles i do 10: calculate constraint error C and scaling factor s 11: end for 12: for all particles i do 13: calculate i 14: perform collision detection and response 15: end for 16: for all particles i do 17: update position x i x + i 18: end for 19: end while 20: for all particles i do 21: update velocity vi 1 x xi 23: update position xi x 24: end for constraint function (1) with respect to the it h particle is given by: pi C = 1 0 pi W (pi p j , h) (7) Plugging this into (6) to find and substituting into (4) we have i = j 1 0 pi W (pi p j , h), (8) where the scaling factor is s = C(p1 , ..., pn ) j 1 0 p j W (pi p j , h) 2 . (9) Because the constraint function (1) is non-linear, with a vanish- ing gradient at the smoothing kernel boundary, the denominator in equation (9) causes instability when particles are close to separat- ing. As in PCISPH this can be solved by pre-computing a conser- vative corrective scale based on a reference particle configuration with a filled neighborhood. Alternatively, constraint force mixing (CFM) [Smith 2006] can be used to regularize the constraint. The idea behind CFM is to soften the constraint by mixing in some of the constraint force back into the constraint function, in the case of PBD this changes (6) to C(p + C(p) + C T + = 0. (10) Where is a user specified relaxation parameter that is constant over the simulation. The scaling factor is now s = C(p1 , ..., pn ) j 1 0 p j W (pi p j , h) 2 + , (11) and the final position update i j on particle i due to a neighbor j i j = 1 0 si + s j (pi p j , h). (12) 2 To appear in ACM TOG 32(4). Figure 2: Armadillo Splash, Top: particle clumping due to neigh- bor deficiencies, Bottom: with artificial pressure term, note the im- proved particle distribution and surface tension. 4 Tensile Instability A common problem in SPH simulations is particle clustering or clumping caused by negative pressures when a particle has only a few neighbors and is unable to satisfy the rest density (Figure 2). This may be avoided by clamping pressures to be non-negative, but at the cost of reduced particle cohesion. Clavet et al. [2005] use a second Õnear pressureÕ term, while Alduan and Otaduy [2011] use discrete element (DEM) forces [Bell et al. 2005] to push apart particles closer than half the smoothing kernel width. Schechter and Bridson [2012] place ghost particles around the free surface to ensure consistent density estimates. We follow the approach of [Monaghan 2000] which adds an arti- ficial pressure specified in terms of the smoothing kernel itself as scorr = k W (pi p j , h) W ( h) n , (13) where is a point some fixed distance inside the smoothing kernel radius and k is a small positive constant. | = 0.1h á á á 0.3h, k = 0.1 and n = 4 work well. We then include this term in the particle position update as i j = si + s j + scorr (pi p j , h). (14) This positive pressure term keeps particle density slightly lower than the rest density. Consequently, particles pull their neighbors inwards causing surface tension-like effects similar to the ones de- scribed in [Clavet et al. 2005]. We note that this effect is a non- physical artifact of the anti-clustering term and requires a trade off between clustering errors and surface tension strength. Without clustering problems our algorithm is free from the rule of thumb that in SPH a particle must have 30-40 neighbors at all times, improving efficiency. 5 Vorticity Confinement and Viscosity Position based methods introduce additional damping which is of- ten undesirable. Fedkiw et al. [2001] introduced vorticity confine- ment to computer graphics to address numerical dissipation in the simulation of smoke, which was later extended to energy conserv- ing fluid simulation in [Lentine et al. 2011]. In Bubbles Alive, Hong et al. [2008] show how vorticity confinement can be used in a hy- brid setup where by vorticity is transferred from a grid to the SPH particles to introduce turbulent motion. We optionally use vorticity confinement to replace lost energy (Fig- ure 5). This requires first calculating the vorticity at a particleÕs lo- cation, for which we use the estimator given in [Monaghan 1992]: i = v = j vi j p j W (pi p j , h) (15) where vi j = v j vi . Once we have the vorticity we calculate a corrective force using the location vector N = | | with = |i fvort icit yi = (N i ) . (16) Unlike [Hong et al. 2008] we do not use normalized as this would increase vorticity indiscriminately. Instead we use the unnormal- ized value, which only adds vorticity where it already exists, as in [Fedkiw et al. 2001]. In addition, we apply XSPH viscosity [Schechter and Bridson 2012], which is important for coherent motion. The parameter c is typically chosen to be 0.01 in our simulations: vnewi = vi + c j vi j á W (pi p j , h) (17) 6 Algorithm Our simulation loop is outlined in Algorithm 1. It is similar to the original Position Based Dynamics update except that each con- straint is solved independently in a Jacobi fashion, rather than through sequential Gauss-Seidel iteration. We perform collision de- tection against solids as part of the constraint solving loop. We recompute particle neighborhoods once per-step and re- calculate distance and constraint values each solver iteration. This optimization can lead to density underestimates, for example if a particle separates from the initial set of neighbors. In PCISPH this can cause serious problems, once a particle becomes isolated, each iteration makes its pressure increasingly negative. If it then comes back into contact on a subsequent iteration, large erroneous pres- sure forces are applied. Our algorithm considers only the current particle positions (not accumulated pressure), so this does not oc- cur. 7 Rendering Real-time fluid surfacing is performed using a GPU based ellipsoid splatting technique. Particle anisotropy is first computed using the method of Yu and Turk [2013], and the surface is reconstructed using a method based on the screen-space filtering presented in [van der Laan et al. 2009]. 3 To appear in ACM TOG 32(4). Figure 3: Dropping a liquid bunny into a pool of water (80k parti- cles). (a) Average Density (b) Maximum Density Figure 4: Density over the bunny drop simulation. Our algorithm maintains compressibility similar to PCISPH at time-steps more than twice as large. Color key: Blue, rest density. Red, PCISPH. Green, our method. 8 Results We tested our algorithm by dropping a liquid bunny into a pool of water (Figure 3) and compared our results with a PCISPH imple- mentation. For this scenario PCISPH is not stable with less than 10 sub-steps per frame ( = 0.0016s). In contrast our algorithm is stable with a single step ( = 0.016s). To compare compressibility we run PCISPH with 10 sub-steps and 4 pressure iterations, and our algorithm with 4 sub-steps and 10 iterations per sub-step, so that each performs 40 pressure iterations per-frame in total. The point of this comparison is to show that our method can achieve comparable results with larger time-steps, allowing us to amortize the per-step costs of grid construction and neighbor finding over more density iterations. Our results are in good accordance, and a plot of density over the simulation confirms that the level of compression is similar despite the larger time-step for our method (Figure 4). Tables 1 and 2 sum- marize the performance of our algorithm in a selection of scenarios. Because we are interested in real-time applications with predictable performance, we set the number of iterations to a fixed value (typ- ically 2-4) rather than solving for a specific error threshold. How- ever, we also show the convergence of our method over multiple iterations in Figure 6. We implemented our algorithm in CUDA and ran our simulations on an NVIDIA GTX 680. Each stage of our algorithm is fully par- allelizable so we are able to take advantage of parallel architec- tures such as GPUs. For neighbor finding we use the method of [Green 2008]. We also perform particle-solid collision detection on the GPU where we use signed distance fields [Bridson et al. 2006] stored as volume textures. 9 Limitations and Future Work Occasionally particle stacking along boundaries can occur due to incorrect density estimates when particles are in contact with solids. Recent work by Akinci et al. [2012] would help address this issue. Jacobi methods only propagate information (in our case position corrections) between a particleÕs immediate neighbors each itera- tion. This can lead to slow convergence as the number of particles increases. More sophisticated parallel solvers such as red-black or multi-scale schemes such as [Solenthaler and Gross 2011] should help improve convergence speed. Because our artificial pressure term is dependent on the spatial res- olution and time-step it can be difficult to adjust parameters inde- pendently. Decoupling these parameters and making anti-clustering independent from surface tension effects would be important future work. Position based dynamics is popular for simulating deformable ob- jects such as cloth. We have prototyped two-way interaction be- tween position based cloth and fluid with promising results. Table 1: Performance results for several examples. A frame time of 16ms is used in all cases. Scene particles steps/frame iters/step time/step [ms] Armadillo Splash 128k 2 3 4.2 Dam Break 100k 4 3 4.3 Bunny Drop 80k 4 10 7.8 Table 2: Breakdown of a frame (percentages) for two examples. Constraint Solve includes collision handling with static objects, and Velocity Update includes vorticity confinement and viscosity calcu- lation. Step Armadillo Splash Dam Break Integrate 1 1 Create Hash Grid 8 6 Detect Neighbors 28 28 Constraint Solve 38 51 Velocity Update 25 14 10 Acknowledgments The authors would like to thank NVIDIA for supporting our re- search, especially Nuttapong Chentanez, Tae-Yong Kim and Simon Schirm for their valuable feedback, Simon Green for his rendering work, and Adam Moravansky and Christian Sigg for their support. We also wish to thank the anonymous reviewers for their generous comments and suggestions. The Bunny and Armadillo models are used courtesy of the Stanford Computer Graphics Laboratory. References AKI NC I , N. , I HMS E N, M. , AKI NC I , G. , S OL E NT HAL E R , B. , AND TE S C HNE R , M. 2012. Versatile rigid-fluid coupling for incom- pressible sph. ACM Trans. Graph. 31, 4 (July), 62:1Ð62:8. 4 To appear in ACM TOG 32(4). Figure 5: Dam break scenario at t=6.0, Left: vorticity confine- ment disabled. Right: vorticity confinement enabled, note the visi- bly higher splash. Figure 6: Convergence of our method over multiple iterations at t = 1.0 in the dam break scenario. AL DU ?AN, I . , AND OTADUY, M. A. 2011. Sph granular flow with friction and cohesion. In Proceedings of the 2011 ACM SIGGRAPH/Eurographics Symposium on Computer Animation, ACM, New York, NY, USA, SCA Õ11, 25Ð32. BE C KE R , M. , AND TE S C HNE R , M. 2007. Weakly compressible sph for free surface flows. In Proceedings of the 2007 ACM SIG- GRAPH/Eurographics symposium on Computer animation, Eu- rographics Association, Aire-la-Ville, Switzerland, Switzerland, SCA Õ07, 209Ð217. BE L L , N. , YU, Y. , AND MUC HA, P. J . 2005. Particle-based sim- ulation of granular materials. In Proceedings of the 2005 ACM SIGGRAPH/Eurographics symposium on Computer animation, ACM, New York, NY, USA, SCA Õ05, 77Ð86. BODI N, K. , LAC OUR S I E R E , C. , AND S E RVI N, M. 2012. Con- straint fluids. IEEE Transactions on Visualization and Computer Graphics 18, 3 (Mar.), 516Ð526. BR AC KB I L L , J . U. , AND RUP P E L , H. M. 1986. Flip: A method for adaptively zoned, particle-in-cell calculations of fluid flows in two dimensions. J. Comput. Phys. 65, 2 (Aug.), 314Ð343. BR I DS ON, R. , F E DKI W, R. , AND M ?UL L E R - FI S C HE R , M. 2006. Fluid simulation: Siggraph 2006 course notes fedkiw and muller-fischer presenation videos are available from the citation page. In ACM SIGGRAPH 2006 Courses, ACM, New York, NY, USA, SIGGRAPH Õ06, 1Ð87. CL AVE T, S. , BE AUDOI N, P. , AND POUL I N, P. 2005. Particle- based viscoelastic fluid simulation. In Proceedings of the 2005 ACM SIGGRAPH/Eurographics symposium on Computer ani- mation, ACM, New York, NY, USA, SCA Õ05, 219Ð228. F E DKI W, R. , S TAM, J . , AND J E NS E N, H. W. 2001. Visual sim- ulation of smoke. In Proceedings of the 28th annual conference on Computer graphics and interactive techniques, ACM, New York, NY, USA, SIGGRAPH Õ01, 15Ð22. GR E E N, S. 2008. Cuda particles. nVidia Whitepaper 2, 3.2, 1. HONG, J . - M. , LE E , H. - Y. , YOON, J . - C. , AND KI M, C. - H. 2008. Bubbles alive. In ACM SIGGRAPH 2008 papers, ACM, New York, NY, USA, SIGGRAPH Õ08, 48:1Ð48:4. LE NT I NE , M. , AANJANE YA, M. , AND FE DKI W, R. 2011. Mass and momentum conservation for fluid simulation. In Proceed- ings of the 2011 ACM SIGGRAPH/Eurographics Symposium on Computer Animation, ACM, New York, NY, USA, SCA Õ11, 91Ð 100. MONAGHAN, J . J . 1992. Smoothed particle hydrodynamics. An- nual Review of Astronomy and Astrophysics 30, 1, 543Ð574. MONAGHAN, J . J . 1994. Simulating free surface flows with sph. J. Comput. Phys. 110, 2 (Feb.), 399Ð406. MONAGHAN, J . J . 2000. Sph without a tensile instability. J. Comput. Phys. 159, 2 (Apr.), 290Ð311. M ?UL L E R , M. , CHARYPAR , D. , AND GROS S , M. 2003. Particle- based fluid simulation for interactive applications. In Proceed- ings of the 2003 ACM SIGGRAPH/Eurographics symposium on Computer animation, Eurographics Association, Aire-la-Ville, Switzerland, Switzerland, SCA Õ03, 154Ð159. M ?UL L E R , M. , HE I DE L B E R GE R , B. , HE NNI X, M. , AND RAT- C L I FF , J . 2007. Position based dynamics. J. Vis. Comun. Image Represent. 18, 2 (Apr.), 109Ð118. RAVE E NDR AN, K. , WOJ TAN, C. , AND TUR K, G. 2011. Hybrid smoothed particle hydrodynamics. In Proceedings of the 2011 ACM SIGGRAPH/Eurographics Symposium on Computer Ani- mation, ACM, New York, NY, USA, SCA Õ11, 33Ð42. SC HE C HT E R , H. , AND BR I DS ON, R. 2012. Ghost sph for animat- ing water. ACM Trans. Graph. 31, 4 (July), 61:1Ð61:8. SMI T H, R. 2006. Open dynamics engine v0.5 user guide. SOL E NT HAL E R , B. , AND GROS S , M. 2011. Two-scale particle simulation. ACM Trans. Graph. 30, 4 (July), 81:1Ð81:8. SOL E NT HAL E R , B. , AND PAJAROL A, R. 2009. Predictive- corrective incompressible sph. In ACM SIGGRAPH 2009 pa- pers, ACM, New York, NY, USA, SIGGRAPH Õ09, 40:1Ð40:6. VAN DE R LAAN, W. J . , GR E E N, S. , AND SAI NZ , M. 2009. Screen space fluid rendering with curvature flow. In Proceedings of the 2009 symposium on Interactive 3D graphics and games, ACM, New York, NY, USA, I3D Õ09, 91Ð98. YU, J . , AND TUR K, G. 2013. Reconstructing surfaces of particle- based fluids using anisotropic kernels. ACM Trans. Graph. 32, 1 (Feb.), 5:1Ð5:12. ZHU, Y. , AND BR I DS ON, R. 2005. Animating sand as a fluid. In ACM SIGGRAPH 2005 Papers, ACM, New York, NY, USA, SIGGRAPH Õ05, 965Ð972. 5