Fluid Simulation For Computer Graphics_ A Tutorial in Grid Based

Document Sample
Fluid Simulation For Computer Graphics_ A Tutorial in Grid Based Powered By Docstoc
					    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 fluid simulator for computer graph-      simulations.
ics applications. Many research papers on fluid 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 fluid as discrete blobs
lation for Computer Graphics.[Bridson 2009]” We base a large por-       of fluid. 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 benefit 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.       fixed points inside of the fluid. At each fixed point, we store quan-
Furthermore, Bridson’s text does not cover particle based methods,      tities such as the velocity of the fluid as it flows by, or the density of
like SPH, which are quickly becoming commonplace within the             the fluid 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 fixed 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 fluid may mean something dif-
ferent than you might usually think. In physics, fluids fall into two    Here we will describe the governing equations for fluid motion, and
categories incompressible and compressible flow. Incompressible          also describe some of the special notation used in fluid simulation
flow is a liquid, such as water or juice. Compressible flow, 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 final para-
flow 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 fluid. Note that there is no such thing as 100 percent       fluid simulation, this section should be read in full.
incompressible fluid. All fluids, 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 fluids 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 fluids. 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:
    †                                                                             ·u=0                               (2)
In these equations, u is fluid velocity. The time variable is t The
density of the fluid is represented by ρ. For water, ρ ≈ 1000 m3 .
The pressure inside the fluid (in force units per unit area) is repre-
sented by p. Note that p and ρ are different things(the first is the
Greek letter rho while the second is p). Body forces (usually just
gravity) are represented by F . Finally, ν is the fluid’s coefficient of
kinematic viscosity.
In our simulator, we currently don’t take fluid viscosity into ac-
count. For inviscid fluids 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 fluids 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 define what exactly we mean by ”grid” in grid based
Lastly, we will quickly discuss some unique notation used in fluid         simulation. Throughout the simulation, we must store many dif-
simulation. In simulation, we take small discrete time steps in or-       ferent quantities (velocity, pressure, fluid concentration, etc.) at
der simulate some phenomenon. Consider the velocity field 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 field. 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 fluid 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 fluid simulation, assuming one wants to simulate n frames of
    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)

                                                                           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 fluids, 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 field, 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 field.
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 field 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 sacrifice 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 fluid 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.
However, this staggered grid is not without drawbacks. In order                                     ∆t = kCF L                                (12)
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 modification
The above notation with half-indices is clearly useful since it            where umax is calculated with:
greatly simplifies 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
                                                                           account the effect that the body forces will have on the simulation’s
                       uy [i][j][k] = vi,j− 1 ,k                    (9)    current timestep.
4.3.2    Advection                                                         cell? This requires extrapolation for our MAC grid. In our expe-
                                                                           rience, simply clamping grid indices is fine 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 fluids. In this section, we will develop a numerical routine such
                                                ∂t                         that our fluid satisfies both the incompressibility condition:
In this section we will develop a computational function for
advect(Qn , ∆t, ∂Q ).
                                                                                                                     · 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
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 finally 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
                                                                                                                            1 pi+1,j − pi,j
                                                                                         un+1,j = un 1 ,j − ∆t
                                                                                          i+ 1
             Calcluate the spatial position of Qi,j,k , store it in X                             2
                                                                                                   i+            2          ρ     ∆h
             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-
                                                                           tain fluid. 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 first 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 fluid. Therefore, we must specify our boundary condi-
We will use this advection scheme throughout our simulator. One
common use is advecting fluid 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 fluid.
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 defined by βi,j . Through
for the pressure.                                                          our definition 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 fluid incompressible. This             for storing our matrix. We store a main linked list, sorted first 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),
                                                                           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 find 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 finite differ-         called the 7 Point Laplacian Matrix. Bridson recommends the Mod-
ences developer earlier:                                                    ified 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
               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]
                         ∆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 finally 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 fluid 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 fluid 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
                                                   
                                                                              • rbridson/mpcg/ Open Source Matlab
                                  .      p2   −D2                          Implementation of Specialized Form of Conjugate Gradient
    β2,1     −Ω2                 . 
                                  .  .          .   
                                        .  =    .
                                                     
    .                   .                  .       .
    .                 − ..
      .                         βn−1,n 
                                                
                                                                          4.3.4    Grid Update
       βn,1    ...    βn,n−1     −Ωn       pn      −Dn
                                                        (31)               4.4     Tracking the Water Surface In Grid Bases Simula-
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 field 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 fluid           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, first 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 first,
For the level set method, we define 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 define 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 satisfied:              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 define Φ(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 first be calculated at
or Catmull-Romm is a fine 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 define our signed dis-            tion, we have not discussed how to update the signed distance as
tance function D(X) as:                                                   the fluid’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 fluid’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 fluid, 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-specific techniques are developed frequently. Typically,
people classify these methods into two groups: PDE based ap-              In this section we describe an alternative approach to fluid 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 fluid. We no longer simulate our fluid 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 briefly 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 difficul-
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
nique was originally introduced for astrophysical simulations[?],          on particle i due to ”something.” Also, recall from basic physics
but has also found a lot of uses in computer graphics [?]. This            (F = M ) that Ai = Mii .
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 fluid dynam-           fluid 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:

                                                                                                                Mj WRij                    (36)

                                                                           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 influence
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 scientifically 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 final 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 fluid 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 fluids 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 fluid 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
                                                                         final production.
Additionally, this fixed grid requires us to change our kernel func-
tion. Oftentimes, implementors in graphics define their kernel func-      While there are benefits 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 fluid 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)

                                                                         Where h is a user specified 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 fluid 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 fluids. ACM Transactions
                                                                         on Graphics (Proc. of ACM SIGGRAPH 2004) 23, 3, 463–468.
Other active areas of fluid simulation research in the graphics com-
munity include, smoke simulation, fire simulation, simulation of       G UENDELMAN , E., S ELLE , A., L OSASSO , F., AND F EDKIW, R.
highly viscous fluids, 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 flow of fluid with a free
lations, and cloth simulations exist as well.
                                                                         surface. the physics of fluids 8.
                                                                      H ARLOW, F., AND W ELCH , J., 1965. Numerical calculation of
                                                                         time-dependent viscous incompressible flow of fluid with a free
                                                                         surface. the physics of fluids 8.
                                                                      H ARLOW, F., AND W ELCH , J., 1965. Numerical calculation of
                                                                         time-dependent viscous incompressible flow of fluid with a free
                                                                         surface. the physics of fluids 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 flowing.
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 fluid:
Ron Fedkiw’s Group                                                       Animating the interplay between rigid bodies and fluid. 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 fluids. 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-fluid 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.