International audienceRasterization pipelines are ubiquitous today. They can be found in most of our personal computers as well as in smaller, hand-held devices--like smart phones--with lower-end hardware. However, simulating particle-based liquids requires sorting the particles which is cumbersome when using a rasterization pipeline. In this chapter, we describe a method to simulate liquids without having to sort the particles. Our method was specifically designed for these architectures and low shader model specifications (starting from shader model 3 for 3D liquids). Instead of sorting the particles, we splat them onto a grid (i.e. a 3D or 2D texture) and solve the inter-particle dynamics directly on the grid. Splatting is simple to perform in a rasterization pipeline, but can also be costly. Thanks to the simplified pass on the grid, we only need to splat the particles once. The grid also provides additional benefits: we can easily add artificial obstacles for the particles to interact with, we can ray cast the grid directly to render the liquid surface, and we can even gain a speed up over sort-based liquid solvers--such as the optimized solver found in the DirectX 11 SDK

Guay, Martin

Hal - Université Grenoble Alpes

Simple, Rasterization-based LiquidsMartin GuayTo cite this version:Martin Guay. Simple, Rasterization-based Liquids. Carsten Dachsbacher. GPU Pro 5, A KPeters/CRC Press, pp.55-64, 2013, 978-1-4822-0864-1. <10.1201/b16721-6>. <hal-01015306>HAL Id: hal-01015306https://hal.inria.fr/hal-01015306Submitted on 26 Jun 2014HAL is a multi-disciplinary open accessarchive for the deposit and dissemination of sci-entific research documents, whether they are pub-lished or not. The documents may come fromteaching and research institutions in France orabroad, or from public or private research centers.L’archive ouverte pluridisciplinaire HAL, estdestine´e au de´poˆt et a` la diffusion de documentsscientifiques de niveau recherche, publie´s ou non,e´manant des e´tablissements d’enseignement et derecherche franc¸ais ou e´trangers, des laboratoirespublics ou prive´s.Simple, Rasterization-based LiquidsMartin GuayGPU Pro 5 Chapter-Author Preprint-1Rasterization pipelines are ubiquitous today. They can be found in mostof our personal computers as well as in smaller, hand-held devices—like smartphones—with lower-end hardware. However, simulating particle-based liq-uids requires sorting the particles which is cumbersome when using a raster-ization pipeline.In this chapter, we describe a method to simulate liquids without havingto sort the particles. Our method was specifically designed for these archi-tectures and low shader model specifications (starting from shader model 3for 3D liquids). Instead of sorting the particles, we splat them onto a grid(i.e. a 3D or 2D texture) and solve the inter-particle dynamics directly onthe grid. Splatting is simple to perform in a rasterization pipeline, but canalso be costly. Thanks to the simplified pass on the grid, we only need tosplat the particles once.The grid also provides additional benefits: we can easily add artificialobstacles for the particles to interact with, we can ray cast the grid directlyto render the liquid surface, and we can even gain a speed up over sort-basedliquid solvers—such as the optimized solver found in the DirectX 11 SDK.1 IntroductionSimulating liquids requires dealing with two phases of fluid—the liquid andthe air—which can be tricky to model as special care may be required forthe interface between the two phases depending on the fluid model. In com-puter graphics, there are mainly two popular formulations for fluids : stronglyincompressible and weakly incompressible.The strong formulation is usually more complex as it requires a hard con-straint (e.g. solving a Poisson eq.), but is more accurate and therefore morevisually pleasing. Because it is more complex, it is often used along simple,regular grid discretizations [Stam 99]. For liquids, several intermediate stepsare required for the surface to behave adequately [Enright et al. 02]. Imple-menting these steps using rasterization APIs is challenging. For instance,[Crane et al. 07] only partially implements them and the fluid behaves morelike a single phase fluid.Furthermore, the strong formulation requires a surface representationlike a level-set density field which requires its own set of specificities (re-initialisation). Again, in [Crane et al. 07] the level-set method is only par-tially implemented and had to be hacked into staying at a certain height;2preventing them from generating such scenarios as the water jet shown inFig. 3.The weak formulation on the other hand, requires only a simple soft con-straint to keep the fluid from compressing. It is much simpler, but also less ac-curate. It is often used along particle discretizations and mesh-free numericalschemes like Smooth Particles Hydrodynamics (SPH) [Desbrun and Cani 96].The advantage of the weak formulation along particles is really for liquids.This combination allowed reproducing the behavior of liquids without com-puting any surface boundary conditions similar to [Enright et al. 02].Additionally, the particles can be used directly to render the liquid surfaceand there is no need for a level-set. The drawback however, is that particlesrequire finding their neighbors, in order to compute forces ensuring they keepat a minimal distance. Typically buckets or other spacial sorting algorithmsare used to cluster the particles into groups [Amada et al. 04, Rozen et al. 08,Bayraktar et al. 08], which can be cumbersome to implement using rasteri-zation APIs.Instead of sorting the particles, our method makes use of rasterizationcapabilities. First we rasterize, or splat, the particle density onto a grid.Then we use simple finite difference to compute all the interacting forces—including the soft incompressibility constraint—on the grid, in a single pass.Some have considered splatting before [Kolb and Cuntz 05], but had tosplat for each force (Pressure and Viscosity), while we only need to splatonce—for all the forces. Finally, the particles are corrected and moved bysampling the grid—which in turn can also be used directly to render theliquid surface by ray casting the splatted density. Overall, our method allowstreating each particle independently while making sure they automaticallyrepulse one-another and avoid rigid obstacles in the domain.2 Simple Liquid ModelOur liquid model is best described as follows: particles are allowed to movefreely in the domain while their mass and kinetic energy remains conserved.For instance, if we add gravity it should translate into velocity u. To keepthe particles from interpenetrating, we add an incompressibility constraint Pderived from the density ρ of particles and the resulting force is the negativegradient of P :3Figure 1: Illustrating the scalar density field used to defined the pressureconstraint. The pressure force is proportional to the density gradient pushingthe particles towards minimum density. For simplicity, we show the idea in2D with density shown as a height field—which, can also be used as the liquidsurface in Shallow Water simulations (see Fig. 4).DρDt= 0, (1)DuDt= −∇P + g + fext, (2)where g gravity and fext accounts for external forces such as user inter-actions. The terms DρDtand DuDtaccount for the transport of density andvelocity. There is never any density added nor removed, it is only trans-ported and held by the particles leading to DρDt= 0. Energy added to thevelocity—like gravity—needs to be conserved as well, resulting in equation(2). Next we define the incompressibility constraint P , and how to computethe density of particles ρ.2.1 Pressure ConstraintTo keep the fluid from compressing and the particles from inter-penetrating,we penalize high density: P = kρ [Desbrun and Cani 96], where k is a stiff-ness parameter that makes the particles repulse one another more strongly4(but can also make the simulation unstable if too large). Hence, by minimiz-ing the density, the particles will move in the direction that reduces densitythe most and thereby avoiding each other. At the same time, gravity andboundary conditions like the walls, will act as to keep the particles fromsimply flying away.Keeping the particles from interpenetrating is crucial in particle-basedliquid simulations. To give strict response to the particles that are close tocolliding, we can make the density nonlinear leading to the pressure con-straint [Becker and Teschner 07]: P = k(ρρ0)γ, where γ is an integer (5 inour case). Dividing by the initial density ρ0 comes in handy for stability; itbecomes easier to control the magnitude of the force through the parameterk and thereby keeping the simulation stable.Setting the initial ρ0 can be tricky. It should be proportional an evalua-tion of ρ at an initial rest configuration; i.e. with a uniform distribution ofparticles. In practice however, we can set it manually. We approximate ρ byperforming a convolution, which we pre-compute on a texture by rasterizing,or splatting the kernels of each particle.3 SplattingTo evaluate the density of particles smoothly, we perform a convolution: thedensity is the weighted average of the surrounding discrete samples; in thiscase the particles. The weight is given by a kernel function which falls offexponentially with distance, as shown in Fig. 2. Instead of sampling thenearby particles, we rasterize the kernel function centered at each particle.The final result is a smooth density grid (texture)—like the one shown inFig. 1—that is equivalent to a convolution evaluation at each point on thetexture. We could say also that we now have a virtual particle on each gridcell.3.1 Rasterizing Kernel FunctionsTo update the velocity on the grid, we need to transfer both the density ofthe particles—to compute pressure—and their velocity. Hence, the first stepin our algorithm is to rasterize the smooth kernel function (red in Fig. 2)and the weighted velocity of each particle. We render the particles as pointsand create quad slices—spanning a cube—in the geometry shader. For each5Figure 2: Splatting consists in rasterizing the smooth kernel function, the 2Dcase shown here in red. In Fig. 1, we see the sum of all the kernel functions;a scalar field representing the density of particles.corner vertex xi, we write the distance d = ∥xi − xp∥ to the center of theparticle xp, and let the rasterizer perform the interpolation between vertices.Then, in a pixel shader, we render the smooth kernel value w(d, r) to thealpha channel, and the weighted velocity w(d, r)up to the other 3 channels—in an additive fashion. Finally, the density on the grid can be sampled bymultiplying the sum of kernel weights by the mass, and the velocity, bydividing the sum of weighted velocities by the sum of weights:ρ(xi) = mi∑pw (∥xi − xp∥, r) u(xi) =∑pw (∥xi − xp∥, r)up∑pw (∥xi − xp∥, r),where i denotes texture indices, p particle indices, and r the kernel radius.We used the following convolution kernel:w(d, r) =(1− d2r2)3.Next we update the velocity field on the grid as to make the particles movein a direction that keeps them from compressing; by computing a pressureforce from the density of particles and adding it to the velocity.64 Grid PassIn the grid pass, we update the splatted velocity field with the pressure force,gravity, and artificial pressure for obstacles (see Section 6). We compute thepressure gradient using a finite difference approximation and add forces tothe velocity field using forward Euler integration:un+1 = un −∆t(P (xi+1)− P (xi−1)∆x,P (xj+1)− P (xj−1)∆y), (3)where ∆t is the time step, ∆x the spatial resolution of grid, and n thetemporal state of the simulation.While we update the velocity, we set the velocity on the boundary cellsof the grid to a no-slip boundary condition by setting the component of thevelocity which is normal to the boundary to 0. This is a simple boundarytest; before writing the final velocity value, we check if the neighbor is aboundary and set the component of the velocity in that direction to 0.5 Particle UpdateWe update the position and velocity of particles following the Particle-In-Cell (PIC) and Fluid-In-Particle (FLIP) approaches that mix particles andgrids [Zhu and Bridson 05]. The main idea with these numerical schemes isthat instead of sampling the grid to assign new values (e.g. velocities) to theparticles, we can recover only the differences to their original values.5.1 Particle VelocityIn PIC, the particle’s velocity is taken directly from the grid, which tendsto be very dissipative, viscous and leads to damped flow. For more livelyfeatures and better energy conservation, FLIP assigns only the differencein velocities; the difference between the splatted velocity and the updatedsplatted velocity discussed in Section 4.By combining both PIC and FLIP, the liquid can be made very viscous likemelting wax, as it can also be made very energetic like water. A parameterr lets the user control the amount of each:7up = run+1(xp) + (1− r)(up −∆u),with ∆u = un(xp)− un+1(xp),where un and un+1 are grid velocities before and after the velocity updatein Section 4, and xp,up are the particle’s position and velocity. Using a bit ofPIC (r = 0.05) is useful in stabilizing the simulation performed with explicitEuler integration, which can become unstable if the time step is too large.5.2 Particle PositionWhile we update the particle velocity, we also update the particle positions.We integrate the particle position using two intermediate steps of Runge-Kutta 2 (RK2); each time sampling the velocity on the grid. The secondorder scheme is only approximate in our case because the velocity field onthe grid is kept constant during integration. At each intermediate step, wekeep the particles from leaving the domain by clamping their positions nearthe box boundaries:xn+1p = ∆tun+1(xn+ 12p ) with xn+ 12p = 0.5∆tun+1(xnp ).Note that we never modify the density value of the particles.6 Rigid ObstaclesWe can prevent the particles from penetrating rigid obstacles in the domainby adding an artificial pressure constraint where the objects are. This followsthe same logic as with the particle density; we define a smooth distance fieldto the surface of the object which can also be viewed as a density field. Wecan use analytical functions for primitives like spheres to approximate theshape of the objects. This avoids rasterizing volumes or voxelizing mesheson the GPU. Alternatively, one could approximate shapes using Metaballs[Blinn 82], and implicit surfaces which naturally provide a distance field.No matter which approach we choose, the gradient of these distance fieldsρObstacle can be computed analytically or numerically and plugged into thevelocity update formula covered in Section 4.8Figure 3: A 3D liquid simulation with obstacles in the domain implementedusing the rasterization pipeline. The simulation runs at 35 FPS using 125kparticles on a Quadro 2000M graphics card.When looking at the ShallowWater Equations (SWE), we find similaritieswith the equations outlined in this chapter. In fact, by looking at the densityfield as the height field of a liquid surface, we can imagine using our methoddirectly for height field simulations shown in Fig. 4. On the other hand, thereis usually a ground height term in the SWE. This in turn can be interpretedin 3D as a distance field for rigid objects as we mentioned above.7 ExamplesWe show a few simulation examples implemented with HLSL, compiled aslevel 4 shaders. They include a 3D liquid with rigid objects in the domain,a 2D shallow water height field simulation, and a 2D simulation comparingwith the optimized Direct Compute implementation of SPH available in theDirectX 11 SDK. Measures include both simulation and rendering. We splatparticles with a radius of 3 cells on 16-bit floating point textures without anysignificant loss in visual quality. In general we used a grid size close to d√NPtexels per axis, where d is the domain dimension (2 or 3) and NP is the totalnumber of particles.In Fig. 3, we rendered the fluid surface by raycasting the density directly.We perform a fixed number of steps and finish with an FXAA anti-aliasingpass (frame buffer size 1024× 768). We used 125k particles on a 963 textureand the simulation performs at 35 frames per second (FPS) on a Quadro2000M graphics card.9Figure 4: The shallow water equations describe the evolution of water heightover time. By looking at the height as the density of a 2D fluid, we seethe equations becoming similar. Hence, our method can be used directly tosimulate height fields. This figure shows a SWE simulation with 42k particlesusing our method.In Fig. 4, we see a shallow water simulation using 42k particles on a 2562grid. The simulation and rendering together run at 130 FPS using the samegraphics card.We compared our solver qualitatively with the SPH GPU implementationin the DirectX 11 SDK (Fig. 5). Their solver is implemented using the DirectCompute API, has shared memory optimizations and particle neighbor searchacceleration. Ours uses HLSL shaders only. In table 6, we compare theperformance of both methods with different particle quantities. We can seethat our method scales better with the number of particles involved for agiven grid size and splatting radius on the hardware we used.8 ConclusionIn GPU Pro 2, we described the $1 fluid solver: by combining the simplicity ofweakly incompressible fluids, with the simplicity of grids, we could simulatea single phase fluid (smoke or fire) in a single pass [Guay et al. 11]. In thischapter, we described the $1 liquid solver for rasterization APIs by combiningthe simplicity of the particles for dealing with liquids, with the simplicity ofthe grids to compute the forces. This is useful in getting a liquid solverrunning quickly on platforms that do not necessarily implement computeAPIs.10Figure 5: Comparison between an optimized SPH solver implemented withCompute Shaders on the left and our method implemented with rasterizationAPIs. Our method performs at 296 FPS using while the optimized SPH solverruns at 169 FPS.Particle Grid Size Our Method DX SDK 11 SpeedupAmount Dim(φ,u) (FPS) (FPS) Ratio64, 000 2562 296 169 1.7532, 000 2562 510 325 1.5516, 000 2562 830 567 1.45Figure 6: Comparison between our method and the optimized SPH solverfound in the DirectX 11 SDK.11References[Amada et al. 04] T. Amada, M. Imura, Y. Yasumoto, Y. Yamabe, andK. Chihara. “Particle-based Fluid simulation on gpu.” In ACM Work-shop on General-Purpose Computing on Graphics Processors, pp. 228–235. ACM, 2004.[Bayraktar et al. 08] S. Bayraktar, U. Gu¨du¨kbay, and B. O¨zgu¨c¸. “GPU-Based Neighbor-Search Algorithm for Particle Simulations.” journal ofgraphics, gpu, and game tools 14:1 (2008), 31–42.[Becker and Teschner 07] M. Becker and M. Teschner. “Weakly compress-ible SPH for free surface flows.” In Proceedings of Eurographics/ ACMSIGGRAPH Symposium on Computer Animation 2007, pp. 209–218,2007.[Blinn 82] J. F. Blinn. “A generalization of algebraic surface drawing.” ACMTransactions on Graphics (TOG) 1:3 (1982), 235–256.[Crane et al. 07] K. Crane, I. Llamas, and S. Tariq. Real-time simulationand rendering of 3d fluids, GPU Gems, 3, Chapter 30. Addison Wesley,2007.[Desbrun and Cani 96] M. Desbrun and M.P. Cani. “Smoothed particles: Anew paradigm for animating highly deformable bodies.” In Proceedingsof Eurographics Workshop on Animation and Simulation, pp. 61–76,1996.[Enright et al. 02] D. Enright, S. Marschner, and R. Fedkiw. “Animationand rendering of complex water surfaces.” In Proceedings of the 29thSIGGRAPH, pp. 736–744, 2002.[Guay et al. 11] M. Guay, F. Colin, and R. Egli. Simple and Fast Fluids,GPU Pro, 2, Chapter VII-3, pp. 433–444. A.K. Peters Ltd, 2011.[Kolb and Cuntz 05] A. Kolb and N. Cuntz. “Dynamic particle coupling forgpu-based fluid simulation.” In Proceedings of the 18th Symposium onSimulation Technique, pp. 722–727, 2005.12[Rozen et al. 08] T. Rozen, K. Boryczko, and W. Alda. “GPU Bucket SortAlgorithm with Applications to Nearest-neighbour Search.” In Jour-nal of the 16th Int. Conf. in Central Europe on Computer Graphics,Visualization and Computer Vision 16:1-3 (2008), 161–167.[Stam 99] J. Stam. “Stable Fluids.” In Proceedings of the 26th ACM SIG-GRAPH conference, pp. 121–128, 1999.[Zhu and Bridson 05] Y. Zhu and R. Bridson. “Animating Sand as a Fluid.”In Proceedings of ACM SIGGRAPH 2005, pp. 965–972, 2005.13

Simple, Rasterization-based Liquids

INRIA a CCSD electronic archive server

HAL Id: hal-01015306https://hal.inria.fr/hal-01015306Submitted on 26 Jun 2014HAL is a multi-disciplinary open accessarchive for the deposit and dissemination of sci-entific research documents, whether they are pub-lished or not. The documents may come fromteaching and research institutions in France orabroad, or from public or private research centers.L’archive ouverte pluridisciplinaire HAL, estdestinée au dépôt et à la diffusion de documentsscientifiques de niveau recherche, publiés ou non,émanant des établissements d’enseignement et derecherche français ou étrangers, des laboratoirespublics ou privés.Simple, Rasterization-based LiquidsMartin GuayTo cite this version:Martin Guay. Simple, Rasterization-based Liquids. Carsten Dachsbacher. GPU Pro 5, A K Pe-ters/CRC Press, pp.55-64, 2013, 978-1-4822-0864-1. ￿10.1201/b16721-6￿. ￿hal-01015306￿Simple, Rasterization-based LiquidsMartin GuayGPU Pro 5 Chapter-Author Preprint-1Rasterization pipelines are ubiquitous today. They can be found in mostof our personal computers as well as in smaller, hand-held devices—like smartphones—with lower-end hardware. However, simulating particle-based liq-uids requires sorting the particles which is cumbersome when using a raster-ization pipeline.In this chapter, we describe a method to simulate liquids without havingto sort the particles. Our method was specifically designed for these archi-tectures and low shader model specifications (starting from shader model 3for 3D liquids). Instead of sorting the particles, we splat them onto a grid(i.e. a 3D or 2D texture) and solve the inter-particle dynamics directly onthe grid. Splatting is simple to perform in a rasterization pipeline, but canalso be costly. Thanks to the simplified pass on the grid, we only need tosplat the particles once.The grid also provides additional benefits: we can easily add artificialobstacles for the particles to interact with, we can ray cast the grid directlyto render the liquid surface, and we can even gain a speed up over sort-basedliquid solvers—such as the optimized solver found in the DirectX 11 SDK.1 IntroductionSimulating liquids requires dealing with two phases of fluid—the liquid andthe air—which can be tricky to model as special care may be required forthe interface between the two phases depending on the fluid model. In com-puter graphics, there are mainly two popular formulations for fluids : stronglyincompressible and weakly incompressible.The strong formulation is usually more complex as it requires a hard con-straint (e.g. solving a Poisson eq.), but is more accurate and therefore morevisually pleasing. Because it is more complex, it is often used along simple,regular grid discretizations [Stam 99]. For liquids, several intermediate stepsare required for the surface to behave adequately [Enright et al. 02]. Imple-menting these steps using rasterization APIs is challenging. For instance,[Crane et al. 07] only partially implements them and the fluid behaves morelike a single phase fluid.Furthermore, the strong formulation requires a surface representationlike a level-set density field which requires its own set of specificities (re-initialisation). Again, in [Crane et al. 07] the level-set method is only par-tially implemented and had to be hacked into staying at a certain height;2preventing them from generating such scenarios as the water jet shown inFig. 3.The weak formulation on the other hand, requires only a simple soft con-straint to keep the fluid from compressing. It is much simpler, but also less ac-curate. It is often used along particle discretizations and mesh-free numericalschemes like Smooth Particles Hydrodynamics (SPH) [Desbrun and Cani 96].The advantage of the weak formulation along particles is really for liquids.This combination allowed reproducing the behavior of liquids without com-puting any surface boundary conditions similar to [Enright et al. 02].Additionally, the particles can be used directly to render the liquid surfaceand there is no need for a level-set. The drawback however, is that particlesrequire finding their neighbors, in order to compute forces ensuring they keepat a minimal distance. Typically buckets or other spacial sorting algorithmsare used to cluster the particles into groups [Amada et al. 04, Rozen et al. 08,Bayraktar et al. 08], which can be cumbersome to implement using rasteri-zation APIs.Instead of sorting the particles, our method makes use of rasterizationcapabilities. First we rasterize, or splat, the particle density onto a grid.Then we use simple finite difference to compute all the interacting forces—including the soft incompressibility constraint—on the grid, in a single pass.Some have considered splatting before [Kolb and Cuntz 05], but had tosplat for each force (Pressure and Viscosity), while we only need to splatonce—for all the forces. Finally, the particles are corrected and moved bysampling the grid—which in turn can also be used directly to render theliquid surface by ray casting the splatted density. Overall, our method allowstreating each particle independently while making sure they automaticallyrepulse one-another and avoid rigid obstacles in the domain.2 Simple Liquid ModelOur liquid model is best described as follows: particles are allowed to movefreely in the domain while their mass and kinetic energy remains conserved.For instance, if we add gravity it should translate into velocity u. To keepthe particles from interpenetrating, we add an incompressibility constraint Pderived from the density ρ of particles and the resulting force is the negativegradient of P :3Figure 1: Illustrating the scalar density field used to defined the pressureconstraint. The pressure force is proportional to the density gradient pushingthe particles towards minimum density. For simplicity, we show the idea in2D with density shown as a height field—which, can also be used as the liquidsurface in Shallow Water simulations (see Fig. 4).DρDt= 0, (1)DuDt= −∇P + g + fext, (2)where g gravity and fext accounts for external forces such as user inter-actions. The terms DρDtand DuDtaccount for the transport of density andvelocity. There is never any density added nor removed, it is only trans-ported and held by the particles leading to DρDt= 0. Energy added to thevelocity—like gravity—needs to be conserved as well, resulting in equation(2). Next we define the incompressibility constraint P , and how to computethe density of particles ρ.2.1 Pressure ConstraintTo keep the fluid from compressing and the particles from inter-penetrating,we penalize high density: P = kρ [Desbrun and Cani 96], where k is a stiff-ness parameter that makes the particles repulse one another more strongly4(but can also make the simulation unstable if too large). Hence, by minimiz-ing the density, the particles will move in the direction that reduces densitythe most and thereby avoiding each other. At the same time, gravity andboundary conditions like the walls, will act as to keep the particles fromsimply flying away.Keeping the particles from interpenetrating is crucial in particle-basedliquid simulations. To give strict response to the particles that are close tocolliding, we can make the density nonlinear leading to the pressure con-straint [Becker and Teschner 07]: P = k(ρρ0)γ, where γ is an integer (5 inour case). Dividing by the initial density ρ0 comes in handy for stability; itbecomes easier to control the magnitude of the force through the parameterk and thereby keeping the simulation stable.Setting the initial ρ0 can be tricky. It should be proportional an evalua-tion of ρ at an initial rest configuration; i.e. with a uniform distribution ofparticles. In practice however, we can set it manually. We approximate ρ byperforming a convolution, which we pre-compute on a texture by rasterizing,or splatting the kernels of each particle.3 SplattingTo evaluate the density of particles smoothly, we perform a convolution: thedensity is the weighted average of the surrounding discrete samples; in thiscase the particles. The weight is given by a kernel function which falls offexponentially with distance, as shown in Fig. 2. Instead of sampling thenearby particles, we rasterize the kernel function centered at each particle.The final result is a smooth density grid (texture)—like the one shown inFig. 1—that is equivalent to a convolution evaluation at each point on thetexture. We could say also that we now have a virtual particle on each gridcell.3.1 Rasterizing Kernel FunctionsTo update the velocity on the grid, we need to transfer both the density ofthe particles—to compute pressure—and their velocity. Hence, the first stepin our algorithm is to rasterize the smooth kernel function (red in Fig. 2)and the weighted velocity of each particle. We render the particles as pointsand create quad slices—spanning a cube—in the geometry shader. For each5Figure 2: Splatting consists in rasterizing the smooth kernel function, the 2Dcase shown here in red. In Fig. 1, we see the sum of all the kernel functions;a scalar field representing the density of particles.corner vertex xi, we write the distance d = ∥xi − xp∥ to the center of theparticle xp, and let the rasterizer perform the interpolation between vertices.Then, in a pixel shader, we render the smooth kernel value w(d, r) to thealpha channel, and the weighted velocity w(d, r)up to the other 3 channels—in an additive fashion. Finally, the density on the grid can be sampled bymultiplying the sum of kernel weights by the mass, and the velocity, bydividing the sum of weighted velocities by the sum of weights:ρ(xi) = mi∑pw (∥xi − xp∥, r) u(xi) =∑p w (∥xi − xp∥, r)up∑p w (∥xi − xp∥, r),where i denotes texture indices, p particle indices, and r the kernel radius.We used the following convolution kernel:w(d, r) =(1− d2r2)3.Next we update the velocity field on the grid as to make the particles movein a direction that keeps them from compressing; by computing a pressureforce from the density of particles and adding it to the velocity.64 Grid PassIn the grid pass, we update the splatted velocity field with the pressure force,gravity, and artificial pressure for obstacles (see Section 6). We compute thepressure gradient using a finite difference approximation and add forces tothe velocity field using forward Euler integration:un+1 = un −∆t(P (xi+1)− P (xi−1)∆x,P (xj+1)− P (xj−1)∆y), (3)where ∆t is the time step, ∆x the spatial resolution of grid, and n thetemporal state of the simulation.While we update the velocity, we set the velocity on the boundary cellsof the grid to a no-slip boundary condition by setting the component of thevelocity which is normal to the boundary to 0. This is a simple boundarytest; before writing the final velocity value, we check if the neighbor is aboundary and set the component of the velocity in that direction to 0.5 Particle UpdateWe update the position and velocity of particles following the Particle-In-Cell (PIC) and Fluid-In-Particle (FLIP) approaches that mix particles andgrids [Zhu and Bridson 05]. The main idea with these numerical schemes isthat instead of sampling the grid to assign new values (e.g. velocities) to theparticles, we can recover only the differences to their original values.5.1 Particle VelocityIn PIC, the particle’s velocity is taken directly from the grid, which tendsto be very dissipative, viscous and leads to damped flow. For more livelyfeatures and better energy conservation, FLIP assigns only the differencein velocities; the difference between the splatted velocity and the updatedsplatted velocity discussed in Section 4.By combining both PIC and FLIP, the liquid can be made very viscous likemelting wax, as it can also be made very energetic like water. A parameterr lets the user control the amount of each:7up = run+1(xp) + (1− r)(up −∆u),with ∆u = un(xp)− un+1(xp),where un and un+1 are grid velocities before and after the velocity updatein Section 4, and xp,up are the particle’s position and velocity. Using a bit ofPIC (r = 0.05) is useful in stabilizing the simulation performed with explicitEuler integration, which can become unstable if the time step is too large.5.2 Particle PositionWhile we update the particle velocity, we also update the particle positions.We integrate the particle position using two intermediate steps of Runge-Kutta 2 (RK2); each time sampling the velocity on the grid. The secondorder scheme is only approximate in our case because the velocity field onthe grid is kept constant during integration. At each intermediate step, wekeep the particles from leaving the domain by clamping their positions nearthe box boundaries:xn+1p = ∆tun+1(xn+ 12p ) with xn+ 12p = 0.5∆tun+1(xnp ).Note that we never modify the density value of the particles.6 Rigid ObstaclesWe can prevent the particles from penetrating rigid obstacles in the domainby adding an artificial pressure constraint where the objects are. This followsthe same logic as with the particle density; we define a smooth distance fieldto the surface of the object which can also be viewed as a density field. Wecan use analytical functions for primitives like spheres to approximate theshape of the objects. This avoids rasterizing volumes or voxelizing mesheson the GPU. Alternatively, one could approximate shapes using Metaballs[Blinn 82], and implicit surfaces which naturally provide a distance field.No matter which approach we choose, the gradient of these distance fieldsρObstacle can be computed analytically or numerically and plugged into thevelocity update formula covered in Section 4.8Figure 3: A 3D liquid simulation with obstacles in the domain implementedusing the rasterization pipeline. The simulation runs at 35 FPS using 125kparticles on a Quadro 2000M graphics card.When looking at the ShallowWater Equations (SWE), we find similaritieswith the equations outlined in this chapter. In fact, by looking at the densityfield as the height field of a liquid surface, we can imagine using our methoddirectly for height field simulations shown in Fig. 4. On the other hand, thereis usually a ground height term in the SWE. This in turn can be interpretedin 3D as a distance field for rigid objects as we mentioned above.7 ExamplesWe show a few simulation examples implemented with HLSL, compiled aslevel 4 shaders. They include a 3D liquid with rigid objects in the domain,a 2D shallow water height field simulation, and a 2D simulation comparingwith the optimized Direct Compute implementation of SPH available in theDirectX 11 SDK. Measures include both simulation and rendering. We splatparticles with a radius of 3 cells on 16-bit floating point textures without anysignificant loss in visual quality. In general we used a grid size close to d√NPtexels per axis, where d is the domain dimension (2 or 3) and NP is the totalnumber of particles.In Fig. 3, we rendered the fluid surface by raycasting the density directly.We perform a fixed number of steps and finish with an FXAA anti-aliasingpass (frame buffer size 1024× 768). We used 125k particles on a 963 textureand the simulation performs at 35 frames per second (FPS) on a Quadro2000M graphics card.9Figure 4: The shallow water equations describe the evolution of water heightover time. By looking at the height as the density of a 2D fluid, we seethe equations becoming similar. Hence, our method can be used directly tosimulate height fields. This figure shows a SWE simulation with 42k particlesusing our method.In Fig. 4, we see a shallow water simulation using 42k particles on a 2562grid. The simulation and rendering together run at 130 FPS using the samegraphics card.We compared our solver qualitatively with the SPH GPU implementationin the DirectX 11 SDK (Fig. 5). Their solver is implemented using the DirectCompute API, has shared memory optimizations and particle neighbor searchacceleration. Ours uses HLSL shaders only. In table 6, we compare theperformance of both methods with different particle quantities. We can seethat our method scales better with the number of particles involved for agiven grid size and splatting radius on the hardware we used.8 ConclusionIn GPU Pro 2, we described the $1 fluid solver: by combining the simplicity ofweakly incompressible fluids, with the simplicity of grids, we could simulatea single phase fluid (smoke or fire) in a single pass [Guay et al. 11]. In thischapter, we described the $1 liquid solver for rasterization APIs by combiningthe simplicity of the particles for dealing with liquids, with the simplicity ofthe grids to compute the forces. This is useful in getting a liquid solverrunning quickly on platforms that do not necessarily implement computeAPIs.10Figure 5: Comparison between an optimized SPH solver implemented withCompute Shaders on the left and our method implemented with rasterizationAPIs. Our method performs at 296 FPS using while the optimized SPH solverruns at 169 FPS.Particle Grid Size Our Method DX SDK 11 SpeedupAmount Dim(φ,u) (FPS) (FPS) Ratio64, 000 2562 296 169 1.7532, 000 2562 510 325 1.5516, 000 2562 830 567 1.45Figure 6: Comparison between our method and the optimized SPH solverfound in the DirectX 11 SDK.11References[Amada et al. 04] T. Amada, M. Imura, Y. Yasumoto, Y. Yamabe, andK. Chihara. “Particle-based Fluid simulation on gpu.” In ACM Work-shop on General-Purpose Computing on Graphics Processors, pp. 228–235. ACM, 2004.[Bayraktar et al. 08] S. Bayraktar, U. Güdükbay, and B. Özgüç. “GPU-Based Neighbor-Search Algorithm for Particle Simulations.” journal ofgraphics, gpu, and game tools 14:1 (2008), 31–42.[Becker and Teschner 07] M. Becker and M. Teschner. “Weakly compress-ible SPH for free surface flows.” In Proceedings of Eurographics/ ACMSIGGRAPH Symposium on Computer Animation 2007, pp. 209–218,2007.[Blinn 82] J. F. Blinn. “A generalization of algebraic surface drawing.” ACMTransactions on Graphics (TOG) 1:3 (1982), 235–256.[Crane et al. 07] K. Crane, I. Llamas, and S. Tariq. Real-time simulationand rendering of 3d fluids, GPU Gems, 3, Chapter 30. Addison Wesley,2007.[Desbrun and Cani 96] M. Desbrun and M.P. Cani. “Smoothed particles: Anew paradigm for animating highly deformable bodies.” In Proceedingsof Eurographics Workshop on Animation and Simulation, pp. 61–76,1996.[Enright et al. 02] D. Enright, S. Marschner, and R. Fedkiw. “Animationand rendering of complex water surfaces.” In Proceedings of the 29thSIGGRAPH, pp. 736–744, 2002.[Guay et al. 11] M. Guay, F. Colin, and R. Egli. Simple and Fast Fluids,GPU Pro, 2, Chapter VII-3, pp. 433–444. A.K. Peters Ltd, 2011.[Kolb and Cuntz 05] A. Kolb and N. Cuntz. “Dynamic particle coupling forgpu-based fluid simulation.” In Proceedings of the 18th Symposium onSimulation Technique, pp. 722–727, 2005.12[Rozen et al. 08] T. Rozen, K. Boryczko, and W. Alda. “GPU Bucket SortAlgorithm with Applications to Nearest-neighbour Search.” In Jour-nal of the 16th Int. Conf. in Central Europe on Computer Graphics,Visualization and Computer Vision 16:1-3 (2008), 161–167.[Stam 99] J. Stam. “Stable Fluids.” In Proceedings of the 26th ACM SIG-GRAPH conference, pp. 121–128, 1999.[Zhu and Bridson 05] Y. Zhu and R. Bridson. “Animating Sand as a Fluid.”In Proceedings of ACM SIGGRAPH 2005, pp. 965–972, 2005.13

https://hal.inria.fr/hal-01015306/file/AuthorVersion.pdf

Simple, Rasterization-based Liquids

Abstract

Similar works

Full text

Available Versions

Hal - Université Grenoble Alpes

INRIA a CCSD electronic archive server