Fluid Simulation For Computer Graphics: A Tutorial in Grid Based and Particle Based Methods Colin Braley∗ Adrian Sandu† Virginia Tech Virginia Tech Figure 1: Fluid Simulation Examples Abstract are typically highly accurate, although relatively slow compared to particle based solutions. Particle based simulations are usually In this paper we present a tutorial on the implementation of both a much faster, but they typically do not look as good as grid based grid based and a particle based ﬂuid simulator for computer graph- simulations. ics applications. Many research papers on ﬂuid simulation are read- ily available, but these papers often assume a very sophisticated Some readers may not know the difference between grid based and mathematical background not held by many undergraduates. Fur- particle based simulations. The best description of this can be found thermore, these papers tend to gloss over the implementation de- in page 6 of [Bridson 2009]. I will attempt to paraphrase this de- tails, which are very important to people trying to implement a scription here. working system. Fluid can be simulated from 2 viewpoints, Lagrangian or Eulerian. Recently, Robert Bridson release the wonderful book, ”Fluid Simu- In the Lagrangian viewpoint, we simulate the ﬂuid as discrete blobs lation for Computer Graphics.[Bridson 2009]” We base a large por- of ﬂuid. Each particle has various properties, such as mass, veloc- tion of our own grid-based simulator off of this text. However, this ity, etc. The beneﬁt of this approach is that conservation of mass text is very dense and theory intensive, and this document serves as comes easily. The Eulerian viewpoint, on the other hand tracks easy version for those who want to implement a simulator quickly. ﬁxed points inside of the ﬂuid. At each ﬁxed point, we store quan- Furthermore, Bridson’s text does not cover particle based methods, tities such as the velocity of the ﬂuid as it ﬂows by, or the density of like SPH, which are quickly becoming commonplace within the the ﬂuid as it passes by. The Eulerian approach corresponds to grid graphics community. This work provides an introduction to SPH based techniques. Grid based techniques have the advantage of hav- as well. ing higher numerical accuracy, since it is easier to work with spatial derivatives on a ﬁxed grid, as opposed to an unstructured cloud of Keywords: Fluids, Physically Based Animation, SPH, Grid Based particles. However, grid based techniques often suffer from mass loss, and are often slower than particle based simulations. Finally, grid based simulations often do much better tracking smooth water 1 Introduction surfaces, whereas particle based approaches often have issues with these smooth surfaces. 2 Introduction 3 Governing Equations In the context of physics, the word ﬂuid may mean something dif- ferent than you might usually think. In physics, ﬂuids fall into two Here we will describe the governing equations for ﬂuid motion, and categories incompressible and compressible ﬂow. Incompressible also describe some of the special notation used in ﬂuid simulation ﬂow is a liquid, such as water or juice. Compressible ﬂow, on the literature. For someone only interested in a basic implementation, other hand, corresponds to gas such as air or steam. Compressible this section can be skimmed with the exception of the ﬁnal para- ﬂow is called compressible, because you can easily change the vol- graph on notation. However, for a full understanding of the why in ume of this ﬂuid. Note that there is no such thing as 100 percent ﬂuid simulation, this section should be read in full. incompressible ﬂuid. All ﬂuids, even water can change volume to some degree. If they could not, there would be no way to audi- First, we will describe the general incompressible Navier Stokes bly yell under water. However, we simply choose to ignore com- equations. We will attempt to intuitively describe the vector cal- pressibility in ﬂuids like water that are nearly incompressible, and culus operators involved, and those without a knowledge of vector instead we refer to them simply as incompressible. calculus may wish to consult the appendix of [Bridson 2009] for review. Here are the incompressible Navier Stokes equations: There are many, many ways to simulate ﬂuids. In graphics, the most common two techniques are grid based simulations, and par- ticle based simulations. (Very recently, new techniques such as the ∂u 1 Lattice-Boltzmann method have been introduced to graphics, but +u+ p=F +ν · u (1) ∂t ρ they are beyond the scope of this paper.) Grid based simulations ∗ e-mail: firstname.lastname@example.org † e-mail:email@example.com ·u=0 (2) In these equations, u is ﬂuid velocity. The time variable is t The kg density of the ﬂuid is represented by ρ. For water, ρ ≈ 1000 m3 . The pressure inside the ﬂuid (in force units per unit area) is repre- sented by p. Note that p and ρ are different things(the ﬁrst is the Greek letter rho while the second is p). Body forces (usually just gravity) are represented by F . Finally, ν is the ﬂuid’s coefﬁcient of kinematic viscosity. In our simulator, we currently don’t take ﬂuid viscosity into ac- count. For inviscid ﬂuids like water, viscosity does usually not play a large role in the look of an animation. However, if you wish to include viscosity, see chapter 8 of [Bridson 2009]. The fundamen- tal work on viscous ﬂuids in graphics is [Goktekin et al. 2004], and can be taken as a starting point for implementation. When the viscosity term is dropped from the incompressible Navier Stokes equations, we get the following set of equations: Figure 2: Our 2D Eulerian Solver ∂u 1 + p=F (3) ∂t ρ ·u=0 (4) 4.2 Data Structures These equations are simpler, and they are the equations we will con- While we have discussed Lagrangian vs. Eulerian viewpoints, we sider for the rest of this paper. These are called the Euler Equations. have yet to deﬁne what exactly we mean by ”grid” in grid based Lastly, we will quickly discuss some unique notation used in ﬂuid simulation. Throughout the simulation, we must store many dif- simulation. In simulation, we take small discrete time steps in or- ferent quantities (velocity, pressure, ﬂuid concentration, etc.) at der simulate some phenomenon. Consider the velocity ﬁeld u be- various points in space. Clearly, we will lay them out in some form ing simulated. In a grid based simulation, we store u as a discreetly of regular grid. However, not just any grid will do. It turns out that sampled vector ﬁeld. However, we need some type of notation to some grids work much better than others. describe which grid element we are referring to. To do this, we use Most peoples intuition is to go with the simplest approach: store ev- subscripts like the following ua,b,c to refer to the vector at a, b, c ery quantity on the same grid. However, for reasons which we will (Note that our indexing scheme is actually more complicated, but soon discuss, this is not a good approach. Back in the 1950’s the this is discussed in the beginning of the grid based simulation sec- seminal paper [Harlow and Welch 1965a] by Harlow and Welch de- tion.) However, we also need a way to indicate which timestep we veloped the innovate MAC Grid technique. (Note that MAC stands are referring to. To do this, we use superscripts, such as uk . a,b,c for Marker-and-Cell.) Among many things, this paper developed a The previous equation would indicate the velocity at index a, b, c new way to track liquid movement through the grid, called marker at timestep k. While some would consider this an egregious abuse particles, and a new type of grid, a staggered grid. Marker particles of notation, this syntax is extremely convenient in ﬂuid simulation. are still used in some simulators for their simplicity, but they are no In order to maintain clarity, we will explicitly state when we are longer state of the art. However, the staggered grid developed by raising a quantity to a power (as opposed to indicating a timestep). Harlow and Welch is still used in many, many simulators. Furthermore, timesteps are written in bold (ab ), whereas exponents are written in a regular script (ab ). This grid is called staggered because it stores different quantities at different locations. In two dimensions, a single cell in a MAC Grid might look as follows: 4 Grid Based Simulation 4.1 Overview Here we will present a high level version of the algorithm for grid based ﬂuid simulation, assuming one wants to simulate n frames of animation. 1. Initialize Grid with some Fluid 2. for( i from 1 to n ) Let t = 0.0 While t < tf rame Calculate ∆t Figure 3: Two Dimensional MAC Cell Advect Fluid Pressure Projection (Pressure Solve) In three dimensions, a MAC cell would look like this: Advect Free Surface Note that in these images p represents pressure, and u is velocity. t = t + ∆t We see in these images that pressure is stored in the center of every grid cell, while velocity is stored on the faces of the cells. Note that Write frame i to disk these velocity samples are the normal component of the velocity at uz [i][j][k] = wi,j,k− 1 (10) 2 Therefore, for a grid of nx , ny , nz cells, we store the pressure in a nx , ny , nz array, the x component of the velocity in a nx +1, ny , nz array, the y component of the velocity in a nx , ny + 1, nz array, and the z component of the velocity in a nx , ny , nz + 1 array. 4.3 Algorithm 4.3.1 Choosing a Timestep Figure 4: Three Dimensional MAC Cell When simulating ﬂuids, we want to simulate as fast as possible without losing numerical accuracy. Therefore, we want to chose a timestep that is as large as possible, but not large enough to desta- each cell face. We will now describe why the grid is arranged in bilize our simulation. The CFL condition helps us do this. The CFL this manner. says to chose a value of ∆t small enough so that when any quantity Consider a quantity w sampled at discrete locations is moved from the center of some cell through the velocity ﬁeld, it w0 , w1 , . . . wi−1 , wi , wi+1 , . . . wn−1 , wn along the real line. will only move ∆h distance. This makes sense intuitively, seeing Imagine we want to estimate the derivative at some sample point i. as if a particle was allowed to move in any larger amounts than this We must use central differences of some form. Immediately, one you would effectively be ignoring some parts of the velocity ﬁeld. can see that: Therefore, our equation for ∆t is as follows: ∂w wi+1 − wi−1 ∆h ≈ (5) ∆t = (11) ∂x i 2∆x umax Recall from numerical analysis that this is O(∆x2 ) accurate. How- As you can see, this requires us to know the maximum velocity in ever, there is clearly a huge problem with this equation. It ignores the velocity ﬁeld at any given time. There are 2 ways to get this the actual value of w at wi ! We clearly need a way to estimate this value, either by doing a linear search through all of the velocities, derivative without ignoring the actual value at the point which we or by keeping track of the maximum velocity throughout the sim- are trying to estimate. We could forward or backward differences, ulation. These details are discussed in the Implementation section. but these are biased and only O(∆x) accurate. Instead, we choose In computer graphics, we are often willing to sacriﬁce strict nu- to stagger our grid to make these accurate central differences work. merical accuracy for the sake of increased computational speed. In Our staggered central difference looks like this: many situations, a practitioner is not worried about if the ﬂuid being simulated is one-hundred percent accurate. Instead, we want plau- wi+ 1 − wi− 1 sible looking results. Therefore, in many applications you can get ∂w 2 2 away with using a timestep larger than that prescribed by the CFL ≈ (6) ∂x i 2∆x condition. For instance, in [Foster and Fedkiw 2001], the authors were able to use a timestep 5 times bigger than that dictated by the This formula is still O(x2 ) accurate. CFL condition. Either way, it is good practice to let the user be able to scale the CFL based timestep by a factor of their choice, kCF L . It turns out that, later on in our Pressure Projection stage, this stag- In this case, our equation becomes: gered grid is very useful, as it allows us to accurately estimate cer- tain derivatives that we require. ∆h However, this staggered grid is not without drawbacks. In order ∆t = kCF L (12) umax to evaluate a pressure value at an arbitrary point which is not an exact grid point, trilinear interpolation (or bilinear in the 2D case) is required. If we want to evaluate velocity anywhere in the grid, a separate trilinear interpolation is required for each component of the κmax − κmin velocity! Therefore, in 3D we need to to 3 trilinear interpolations, f (v, vmax , vmin ) = κmin (v − vmin ) (13) vmax − vmin and in the 2D case we need to do 2 bilinear interpolations! Clearly, this is slightly unwieldy. However, the added accuracy makes up for the complications in interpolation. In [Bridson 2009], a slightly more robust treatment of the CFL condition is presented. In this text, Bridson suggests a modiﬁcation The above notation with half-indices is clearly useful since it where umax is calculated with: greatly simpliﬁes our formulae. However, it is obvious that half- indices can’t be used in an actual implementation. [Bridson 2009] recommends the following formulae for converting these half- umax = max (|u|) + ∆h |F | (14) indices to array indices for a real system: p[i][j][k] = pi,j,k (7) where F is whatever body forces are to be applied (usually just gravity), and max u is simply the largest velocity value currently ux [i][j][k] = ui− 1 ,j,k (8) on the grid. This solution is slightly more robust in that it takes into 2 account the effect that the body forces will have on the simulation’s uy [i][j][k] = vi,j− 1 ,k (9) current timestep. 2 4.3.2 Advection cell? This requires extrapolation for our MAC grid. In our expe- rience, simply clamping grid indices is ﬁne in the advection code, Central to any grid based method is our ability to advect both scalar but more advanced techniques do exist. However, we have found and vector quantities through our simulation grid. Advection can that in practice these advanced extrapolation techniques do little to be informally described as follows: ”Given some quantity Q on our visually augment the simulation. simulation grid, how will Q change ∆t later?” More formally, we can describe advection as: 4.3.3 Pressure Solve ∂Q n So far, we have done nothing to deal with the incompressbility of Qn+1 = advect(Qn , ∆t, ) (15) our ﬂuids. In this section, we will develop a numerical routine such ∂t that our ﬂuid satisﬁes both the incompressibility condition: In this section we will develop a computational function for n advect(Qn , ∆t, ∂Q ). ∂t · un+1 = 0 (20) Consider a P on our simulation grid. Using our central differencing as well as our boundary conditions: schemes described previously, we can trivially calculate ∂Q . Using ∂t this derivative, along with our grid information, we can develop un+1 · n = usolid · n ˆ ˆ (21) a technique to advect quantities through the grid. This technique sometimes called a backwards particle trace. Since we are using a particle, this is also commonly referred to as Semi-Lagrangian Additionally, this section ﬁnally allows us to show the reason for Advection. It is important to note that no particle is ever created, our staggered MAC grid discussed in our previous section. First, and the particle is purely conceptual. This is what leads to the Semi we will work out the individual equations to make a single grid in Semi-Lagrangian Advection. cell satisfy our two conditions above. Then, we will show how this information can be used to make the entire grid incompressible. Note that our algorithm can not be done in place, and requires an extra copy of the pertinent data in our simulation grid. Our algo- Consider a 2D MAC cell at location i, j. Per the Euler equations, rithm works as follows: on every step we must update our cells velocity by the following equations. First, we present them in 2D where our velocity is rep- 1. For each grid cell with index i, j, k resented by u =< u, v >. Calculate − ∂Q ∂t 1 pi+1,j − pi,j un+1,j = un 1 ,j − ∆t i+ 1 (22) Calcluate the spatial position of Qi,j,k , store it in X 2 i+ 2 ρ ∆h ∂Q Calculate Xprev = X − ∂t ∗ ∆t n+1 n 1 pi,j+1 − pi,j vi,j+ 1 = vi,j+ 1 − ∆t (23) Set the gridpoint for Qn+1 that is nearest to Xprev equal 2 2 ρ ∆h to Qi,j,k Here are the equivalent equations in 3D for u =< u, v, w >. 2. Set Q = Qn+1 This algorithm is very simple, and fairly accurate. However, if you 1 pi+1,j,k − pi,j,k are familiar with numerical analysis, you will recognize that this un+1,j,k = un 1 ,j,k − ∆t i+ 1 i+ (24) 2 2 ρ ∆h algorithm uses the simple time integrator Forward Euler. This inte- grator is not very accurate. We recommend at least using an integra- tor such as RK2 or better (Runge Kutta Order-2). In our simulator, n+1 n 1 pi,j+1,k − pi,j,k vi,j+ 1 ,k = vi,j+ 1 ,k − ∆t (25) we tested out 5 different integrators, and found the following O(h3 ) 2 2 ρ ∆h accurate scheme to work the best: n+1 n 1 pi,j,k+1 − pi,j,k wi,j,k+ 1 = wi,j,k+ 1 − ∆t (26) κ1 = f (Q )n (16) 2 2 ρ ∆h Just in case these don’t seem complicated enough already, there is 1 another issue that must be attended to. These equations are only ap- κ2 = f (Qn + ∆tκ1 ) (17) plied to components of the velocity that border a grid cell that con- 2 tain ﬂuid. Getting these conditions correct was one of the hardest things in the actual programming of our simulation, and we recom- 3 mend that an implementor try to ﬁrst program this in 2D for easier κ3 = f (Qn + ∆tκ2 ) (18) 4 debugging. A careful reader will also notice that these equations may require n+1 n2 3 4 the pressures of grid cells that lie either outside of the grid, or out- Q = Q + ∆tκ1 + ∆tκ2 + ∆tκ3 (19) 9 9 9 side of the ﬂuid. Therefore, we must specify our boundary condi- tions. We will use this advection scheme throughout our simulator. One common use is advecting ﬂuid velocity itself. Another use is ad- There are two primary types of boundary conditions in grid based vecting temperatures or material properties in advanced simulators. simulation, Dirichlet and Neumann. We will use Dirichlet condi- tions for free surface boundaries, indicating that we will specify the A careful reader might have noticed one issue with the advection value of the quantity at and boundary case. Therefore, we simply psuedocode. How would we perform an advection for a boundary assume that pressure is 0 in any region of air outside of the ﬂuid. The more complicated boundary is with solid walls. Here we will we should store is as a sparse matrix. Each row has at most 4 non- use a Neumann boundary condition. Using the above pressure up- zero entries in the 2D case, and 6 non-zero entries in the 3D case. date equations, we substitute in the solids velocity (0 for simula- tions without moving solids), and then we arrive at a single linear Furthermore, is is clear that A is symmetric. Every entry at index equation for our pressure. Rearranging our terms allows us to solve i, j that is not on the main diagonal is deﬁned by βi,j . Through for the pressure. our deﬁnition of β, it is clear that βi,j = βj,i . Therefore, we only need to store half of the entries in A. We use the following scheme Now we will work out how to make our ﬂuid incompressible. This for storing our matrix. We store a main linked list, sorted ﬁrst by means that, for every velocity component on the grid, we want to column, then by row. We store in each linked list node the following satisfy: < i, j, pij >. The column index is stored as i, the row index is stored by j, and the corresponding pressure value is stored as p. This storage scheme has many pros and cons. We are storing the ·u=0 (27) minimum amount of data, as we are storing only half of the matrix’s non-zero entries. However, this memory saving comes at the cost Note that this divergence operator can be expanded to: of speed. Accessing an arbitrary matrix entry is O( 1 n) = O(n), 2 which can be slow. For small simulations where memory usage is ∂u ∂v ∂w not a concern, we recommend including the option to store entries ·u= + + (28) ∂x ∂y ∂z in a dense matrix. Therefore, it is clear that we simply want to ﬁnd a way so that each Finally, a reader highly experienced in numerical analysis will real- component of our spatial derivatives equals zero. Recalling our cen- ize that A has a form that is common to many other matrices. In 2D, tral differences used earlier, we can approximate these divergences A is often called the 5 point Laplacian Matrix, whereas in 3D it is using the following numerical routines, which use our ﬁnite differ- called the 7 Point Laplacian Matrix. Bridson recommends the Mod- ences developer earlier: iﬁed Incomplete Cholesky Conjugate Gradient Level 0 algorithm [Bridson 2009]. Essentially, this is simply the conjugate gradient algorithm with a special preconditioner designed for this particular matrix. ui+ 1 ,j,k − ui− 1 ,j,k vi,j+ 1 ,k − vi,j− 1 ,k wi,j,k+ 1 − wi,j,k− 1 If the reader is interested in implementing their own con- ( ·u)i,j,k ≈ 2 2 + 2 2 + 2 jugate gradient solver, we recommend the paper [Shewchuk 2007] 2 ∆h ∆h ∆h as a good starting point. However, in our implementation we are (29) not planning on dealing with enormous bodies of water so we im- Finally, we have all the quantities necessary for our pressure update. plemented both the regular Conjugate Gradient algorithm, as well We have developed equations for how the pressure affects the veloc- as Parallel Successive Over-Relaxation (Parallel SOR), and the Ja- ity, and we also have numerical equations to estimate the pressure cobi Method. We found the parallel SOR to be faster than the con- gradient. Using this information, we can create a linear equation jugate gradient implementation, but this is probably because we for the new pressure in every grid cell. We can then combine these are not using a preconditioner. However, conjugate gradient was equations together into a system of simultaneous linear equations slightly faster than the Jacobi method in our tests. Note that we which we can solve for the whole grid, and ﬁnally complete our used OpenMP for parallelizing our SOR implementation. We hope pressure update. to add a preconditioner to our conjugate gradient implementation in the near future. Eventually, we hope to end up with a system of equations of the form: While we have described the characteristics of the linear system to Ax = b (30) solve, and what kind of solvers to use, we realize that most readers will not want to spend their time writing a highly optimized im- Every row of A corresponds to one equation for one ﬂuid cell. In plementation of a specialized conjugate gradient solver. Therefore, this formulation, we will setup our matrix such that b is simply our we will quickly direct the reader towards a few good linear algebra negative divergences for every ﬂuid cell. When written out, our packages that have routines that suit our purposes: linear system takes the following form: • Boost µBLAS • SparseKit −Ω1 β1,2 ... β1,n p1 −D1 • http://people.cs.ubc.ca/ rbridson/mpcg/ Open Source Matlab . p2 −D2 Implementation of Specialized Form of Conjugate Gradient β2,1 −Ω2 . . . . . = . . . . . . − .. . βn−1,n pn−1 −Dn−1 4.3.4 Grid Update βn,1 ... βn,n−1 −Ωn pn −Dn (31) 4.4 Tracking the Water Surface In Grid Bases Simula- tion In this equation, Di is the divergence through cell i, Ωi is the num- ber of non-solid neighbors of cell i, and βij takes values based on While we have discussed the basic mechanisms for a grid based the below equation: simulation, we have not discussed how to unify all these ideas into a full working simulator that can output data that encapsulates the 1 if cell i is a neighbor of cell j position of a moving water surface. βij = (32) 0 otherwise. There are many ways to approach the problem of tracking the move- ment of water through a simulation grid. The simplest way, in- Our matrix A has many unique properties we can exploit both in our troduced all the way back in Harlow and Welch’s seminal paper choice of linear solver and in our storage of A itself. It is immedi- [Harlow and Welch 1965b]. This approach is relatively simple, ately clear that A is sparse(most of its entires are zero), indicating and still quite useful. Here we store a collection of many discrete marker particles in our simulation, each representing a water parti- 2. for each grid point P directly at the free-surface, set the signed cle. Every timestep, we advect them through the velocity ﬁeld by distance to 0 ∆t. Also, we store an enumeration value inside of each cell indi- cating whether the cell contains liquid, air, or solid. Once a ﬂuid 3. Loop over each grid-point Gi,j,k at which signed distance is marker particle moves into a cell, we mark it as liquid. This is nec- unknown essary for our pressure solve. After each timestep, we can output Loop over each grid-point P that neighbors Gi,j,k as these particles to disk. However, the question remains as to how long as the signed distance at P is known to render these particles. One approach is to use a implicit surface function to generate a water surface from these particles. We dis- Find the distance from Gi,j,k to the surface points. cuss this approach later, in section the section on surfacing SPH If this distance is closer than P ’s signed distance, mark P as simulations. Unfortunately, this approach can lead to blobby, ugly unknown once again surfaces. Therefore, we turn to a level-set based approach, ﬁrst in- troduced in [?]. Take the minimum value of the distances of the neigh- bors, and determine if Gi,j,k is inside or outside, and set the Level set methods are currently the best way to achieve smooth sign of the distance based on this high quality free surfaces in liquid simulation. However, they are far more computationally expensive than the above implicit surface There are two main techniques for implementing such an algorithm. approach. Here we will present a brief introduction to these tech- These techniques are the fast marching method, and the fast sweep- niques. We recommend that an interested reader refer to [Fedkiw ing method. and Sethian 2002] for more detailed information. The fast marching method loops over the closest grid points ﬁrst, For the level set method, we deﬁne a new value, φi,j,k , at the center and then those that are farther away. This technique works rapidly of all of our simulation cells. We deﬁne our liquid free surface to by storing the unknown grid points in a priority queue data struc- exist at locations where the following equation is satisﬁed: ture. This algorithm runs in O(nlog(n)) when the priority queue is implemented with a heap. A detailed description can be found in [Sethian 1999]. φ(X) = 0 (33) The other technique is the fast sweeping method. The fast sweeping Where X is a position vector. Note that we can deﬁne Φ(X) at non- method takes the opposite approach from the fast marching method. grid cell locations through any type of interpolation, either trilinear Here, we allow the signed distance function to ﬁrst be calculated at or Catmull-Romm is a ﬁne choice. our farthest away points. We then have this information propagate back towards the surface. Fast marching is great because it is O(n), Furthermore, we say that locations that satisfy φ(X) < 0 to be and is very simple. Furthermore, it works well with narrow band inside of the water, and φ(X) < 0 to be outside of the water. To methods, discussed in [Bridson 2009] and [Fedkiw et al. 2001a]. represent φ, we use a function called the signed distance function. While we have discussed how to compute a signed distance func- Given an arbitrary set S of S points, we deﬁne our signed dis- tion, we have not discussed how to update the signed distance as tance function D(X) as: the ﬂuid’s free surface moves. While this may seem complicated, this step is quite trivial. Since our φ values are stored in the cen- ter of our grid cells, we can simply advect these values according DS (X) = minp∈S X − p (34) to the ﬂuid’s velocity. However, it turns out that advection does not perfectly preserve signed distance. Therefore, we periodically recalculate our signed distance every few timesteps. Bridson rec- Clearly, for some arbitrary point X, the magnitude of the signed comends that we recalculate our signed distance once per frame distance is the distance to the nearest point in the set S. Signed (note that typically many timesteps of ∆t are required per frame) distance is also useful because of another property. If we want to [Bridson 2009]. check whether a grid cell is inside or outside of the ﬂuid, all we must do is examine the sign of the signed distance. TODO: Finish this section. I am still working on it since my level set implementation is not complete. Thus far we have ignored an important problem with signed dis- tance: how to compute it. At the beginning of a simulation, we can assume our signed distance function is already computed on 5 Particle Based Simulation the grid. 5.1 Overview There are many ways to calculate signed distance, and new problem-speciﬁc techniques are developed frequently. Typically, people classify these methods into two groups: PDE based ap- In this section we describe an alternative approach to ﬂuid simula- proaches and Geometric Approaches. PDE Based techniques tion. Here we discuss Lagrangian techniques. In this approach, we approximate something called the Eikonal Equation, φ = have a set of discrete particles that move through space to represent 1. These techniques are mathematically and computationally in- our ﬂuid. We no longer simulate our ﬂuid on a grid structure. volved, and are often overkill for graphics work. Instead, we will This approach has many pros and cons compared to grid based tech- discuss brieﬂy the geometric approaches. Our discussion will not niques. In general, particle based approaches are less accurate than go into much depth; for a more in depth treatment of geometric their grid based counterparts. This is primarily due to the difﬁcul- algorithms for computing signed distance see [?]. ties in dealing with spatial derivatives on an unstructured particle Our algorithms for computing signed distance take the following cloud. However, particle based simulations are typically much eas- general form: ier to program and understand. Furthermore, particle based tech- niques are much faster, and can be used in real time applications 1. Set the signed distance of each grid-point to ”unknown” such as video games. We will describe Smoothed Particle Hydrodynamics. This tech- Note that in this equation, Asomething refers to the acceleration i nique was originally introduced for astrophysical simulations[?], on particle i due to ”something.” Also, recall from basic physics F but has also found a lot of uses in computer graphics [?]. This (F = M ) that Ai = Mii . A description is especially valuable, since SPH is not discussed in Bridson’s text[Bridson 2009], and crucial implementation details In order for our simulation to progress, we need a way to calculate are scattered through both astophyics, computational ﬂuid dynam- ﬂuid density at some arbitrary point. ics, and computer graphics literature. Certain particle properties, such as mass, are given initial values at the beginning of the simulation and are not expected to change. However, other properties must be recalculated every step. Con- sider the pressure property. Here is how to determine the new pres- sure every time step: However, we have a discrete cloud off particles, so we must use a discrete summation to approximate this integral. This leads us to the equation: n Mj WRij (36) j=i Here, Rij is equal to the Euclidean distance between particle i and particle j. Figure 5: Our Simple 2D SPH Implementation This function W (d) is known as a kernel function. This function takes a single scalar parameter, which is a distance between two 5.2 Data Structures particles, and returns a scalar ∈ [0, 1]. Typically, a kernel func- tion maps particles that are farther away to values closer to 0. This makes sense, since particles far away will not have a large inﬂuence First, we must outline what information we need to store for our on a particle. simulation. Clearly, we need a data structure to list all of our parti- cles. Since particles must frequently be added to the simulation, we Once particles are far enough away from the source, the kernel will choose a linked list. However, the question remains as to what function drops to 0, and therefore these particles no longer have information is stored in each particle. to be considered. We will exploit this fact in a later section when developing acceleration structures for SPH simulations. Clearly, we need to store position, velocity, mass, density, and pres- sure. It turns out to be useful to store color and a force vector as The question of what kernel function is best is still a very open well, so we will store these. We will refer to these quantities with research question. However, since our simulations are targeted at the following variables throughout our discussion: begin visually pleasing, rather than scientiﬁcally accurate, we do not care much about this. The following kernel function has been • X Position used extensively in research and practical applications. This is the Gaussian Kernel. • V Velocity • M Mass 1 r2 W (d) = 3 exp( ) (37) π h3 2 h2 • d Density Here r is a the distance between two particles, h is our smoothing • ρ Pressure width. Once particles are greater than distance 2h away, they will no longer affect the particles in question. Clearly, larger values of • C =< Cred , Cgreen , Cblue > Color h will make for a more realistic simulation, albeit at the expense of computational speed. • F Force Finally, we present full psueo-code for an SPH simulation: Note that for our color, each component of the color is ∈ [0, 1]. Each particle can be trivially implemented as a C structure. Nota- • Initialize all particles tionally, we will refer to particles with the variable P , and individ- • Set t = 0 ual particles using subscript notation. For instance, the i-th particle would be Pi . Finally, we will refer to particle quantities in a similar • Choose a ∆t manner. For example, the mass of the 12th particle would be M12 . • for i from 0 to n 5.3 Algorithm for j from 1 to numparticles Get list Lj of neighbors for Pj Our ﬁnal goal is to satisfy the following condition: Calculate Densityj for Pj using Lj For all particles Pi : Calculate P ressurej for Pj using Lj ∂V Calculate acceleration Aj for Pj using Densityj and = Apressure + Aviscosity + Agravity + Aexternal i i i i (35) P ressurej ∂t i Move Pj using Aj and ∆t using Euler step t = t + ∆t • Cleanup all data structures • Exit 5.4 Acceleration Structures As described above, there is a clear computational bottleneck in our application. We must calculate interaction forces between each and Figure 7: Raytraced GPU Based 3D SPH from University of Tokyo every particle. This is an O(n2 ) process, which is not computa- tionally viable. By using spatial data structures, we can reduce our computation time to O(n) in the typical case. In the worst case, 5.5 Surface Tracking where all of the particles are in one cell, our algorithm still runs O(n2 ) Advanced techniques using quadtrees in 2D, or octrees in At our current stage, all SPH results in is an unorganized point 3D, can eliminate this possibility. cloud of ﬂuid particles. This is unacceptable for most applications. Usually the compuiter graphics practitioner desires a way to render In addition to storing all of our particles in a linked list, we also these ﬂuids using a off-the-shelf 3D renderer. In this section, we store them in a spatial grid data structure. Our grid cells extend will outline a few techniques for transforming the result of our SPH by a distance of R in each dimension. Therefore, in order to cal- simulations into a renderable form. culate the forces on a particular particle, one must only examine 9 grid cells in the 2D case, or 27 grid cells in the 3D case. This T and probably the easiest technique, is to sample our SPH results is because, for any grid cells far enough away, our kernel function onto a uniform grid. Here we can step through the uniform grid will evaluate to 0 and their contributions will not be included on the points, and sample the ﬂuid density on these points. Then, we can current particle. In our experience, including this spatial grid can use this uniform grid as input to an application that performs iso- decrease simulation time by over an order of magnitude for large surface volume rendering. Here we have a choice between direct enough simulations. volume rendering, such as that presented in [Colin Braley 2009], and marching cubes[Lorenson and Cline 1982]. Direct volume ren- However, the addition of this grid data structure requires additional dering, typically done through volume raycasting, has the advan- book-keeping during the simulation process. Whenever a particle tage of speed. However, there is no easy way to integrate this result- is moved, one must remove it from its current grid cell, and add it ing image with other generated images, and therefore this technique to the grid cell it belongs in. Unfortunately, there is no way to do is only suitable for creating previews. Marching cubes, on the other this simulation in place, and one must maintain two copies of the hand, produces triangle meshes. These meshes are suitable for use simulation grid. in a 3D animation program, and this is therefore a viable option for ﬁnal production. Additionally, this ﬁxed grid requires us to change our kernel func- tion. Oftentimes, implementors in graphics deﬁne their kernel func- While there are beneﬁts to sampling our SPH onto a grid, there are tions using a piecewise function. If the particles are less than some other techniques that often produce better results. Usually, these distance h apart, they evaluate some type of spline. If the particles techniques use a special function for each particle that, whe com- are farther away, the kernel instead evalutes to 0. However, more bined with the other particles, produces a ﬂuid surface. advanced kernels exist for these types of simulations, and they have recently been used with success in graphics. One of the most com- The function usually used for this purpose was introduced by [Blinn mon advanced kernels is the cubic spline kernel. 1982]. This technique is often referred to as meta-balls or blobbies in the graphics community. X − xi F (X) = Σn k( i ) (38) h Where h is a user speciﬁed parameter representing the smooth- ness of the surface, and k is a kernel function like the one repre- sented above. Incremenetal improvements have been made to the the above function throughout the years, the most imporortant of which is presetned in [Williams 2008]. Figure 6: Our 2D SPH Implementation with Lookup Grid Finally, SPH can be further optimized in another way. SPH is clearly very data parallel. Therefore, each particle can be simu- lated in a separate thread with relative ease. Because of this, many Figure 8: Metaball Based Surface Reconstruction high performance SPH implementations are done on the GPU. 5.6 Extensions F EDKIW, R., S TAM , J., AND J ENSEN , H. W., 2001. Visual simu- lation of smoke. This document has only scratched the surface in discussing the cur- F OSTER , N., AND F EDKIW, R. 2001. Practical animation of liq- rent state of the art ﬂuid simulation techniques. For a thorough uids. In SIGGRAPH ’01: Proceedings of the 28th annual con- explanation of grid based methods, see [Bridson 2009]. No com- ference on Computer graphics and interactive techniques, ACM parable resource for SPH and particle based technique exists, but Press, New York, NY, USA, 23–30. the author recommends reading the papers of Nils Theurey, and the SIGGRAPH 2006 Fluid Simulation Course Notes as a starting G OKTEKIN , T. G., BARGTEIL , A. W., AND O’B RIEN , J. F. 2004. point. A method for animating viscoelastic ﬂuids. ACM Transactions on Graphics (Proc. of ACM SIGGRAPH 2004) 23, 3, 463–468. Other active areas of ﬂuid simulation research in the graphics com- munity include, smoke simulation, ﬁre simulation, simulation of G UENDELMAN , E., S ELLE , A., L OSASSO , F., AND F EDKIW, R. highly viscous ﬂuids, and coupled simulations. Coupled simula- 2005. Coupling water and smoke to thin deformable and rigid tions combine two or more simulations and get them to interact shells. In SIGGRAPH ’05: ACM SIGGRAPH 2005 Papers, plausibly. Recently, Ron Fedkiw’s group achieved 2-way coupled ACM, New York, NY, USA, 973–981. SPH and grid based simulations in [?]. Furthermore, examples of H ARLOW, F., AND W ELCH , J., 1965. Numerical calculation of couplings with thin shells, rigid body simulations, soft body simu- time-dependent viscous incompressible ﬂow of ﬂuid with a free lations, and cloth simulations exist as well. surface. the physics of ﬂuids 8. H ARLOW, F., AND W ELCH , J., 1965. Numerical calculation of time-dependent viscous incompressible ﬂow of ﬂuid with a free surface. the physics of ﬂuids 8. H ARLOW, F., AND W ELCH , J., 1965. Numerical calculation of time-dependent viscous incompressible ﬂow of ﬂuid with a free surface. the physics of ﬂuids 8. L ORENSON , AND C LINE. 1982. Marching cubes. ACM Trans. Graph. 1, 3, 235–256. P ETER , M. C., M UCHA , P. J., B ROOKS , R., I II , V. H., AND T URK , G., 2002. Melting and ﬂowing. Figure 9: Two Way Coupled SPH and Grid Based Simulation by P ETER , M. C., M UCHA , P. J., AND T URK , G. 2004. Rigid ﬂuid: Ron Fedkiw’s Group Animating the interplay between rigid bodies and ﬂuid. In ACM Trans. Graph, 377–384. S ETHIAN , J. A. 1999. Level Set Methods and Fast Marching Meth- Acknowledgments ods: Evolving Interfaces in Computational Geometry, Fluid Me- chanics, Computer Vision, and Materials Science (Cambridge ... Thanks to Robert Hagan for editing this work. on Applied and Computational Mathematics), 2 ed. Cambridge University Press, June. References S HEWCHUK , J. R. 2007. Conjugate gradient without the agonizing pain. Tech. rep., Carnegie Mellon. BATTY, C., AND B RIDSON , R. 2008. Accurate viscous free sur- faces for buckling, coiling, and rotating liquids. In Proceedings S TAM , J. 1999. Stable ﬂuids. 121–128. of the 2008 ACM/Eurographics Symposium on Computer Ani- W ILLIAMS , B. W. 2008. Fluid Surface Reconstruction from Par- mation, 219–228. ticles. Master’s thesis, University of British Columbia. BATTY, C., B ERTAILS , F., AND B RIDSON , R. 2007. A fast varia- tional framework for accurate solid-ﬂuid coupling. ACM Trans. Graph. 26, 3, 100. B LINN , J. F. 1982. A generalization of algebraic surface drawing. ACM Trans. Graph. 1, 3, 235–256. B RIDSON , R. 2009. Fluid Simulation For Computer Graphics. A.K Peters. C OLIN B RALEY, ROBERT H AGAN , Y. C. D. G. 2009. Gpu based isosurface volume rendering using depth based coherence. Sig- graph Asia Technical Sketches. F EDKIW, R., AND S ETHIAN , J. 2002. Level Set Methods for Dy- namic and Implicit Surfaces. F EDKIW, R., S TAM , J., AND J ENSEN , H. W. 2001. Visual sim- ulation of smoke. In Proceedings of SIGGRAPH 2001, ACM Press / ACM SIGGRAPH, E. Fiume, Ed., Computer Graphics Proceedings, Annual Conference Series, ACM, 15–22.