[In this technical article originally printed in Game Developer magazine, Neversoft cofounder Mick West looks at how to efficiently implement fluid effects  from smoke to water and beyond  in video games, with example code.]
Fluid effects, such as rising smoke and turbulent water flow, are everywhere in nature but are seldom implemented convincingly in computer games. The simulation of fluids (which covers both liquids and gases) is computationally very expensive.
It's also mentally expensive, with even introductory papers on the subject relying on the reader to have math skills at least at the undergraduate calculus level.
In this twopart article, I will attempt to address both these problems from the perspective of a game programmer who's not necessarily conversant with vector calculus. I'll explain how certain fluid effects work without using advanced equations and without too much new terminology.
I'll also describe one way of implementing the simulation of fluids in an efficient manner without the expensive iterative diffusion and projection steps found in other implementations.
A working demonstration in ZIP form with source code accompanies this article; example output from this can be seen in Figure 1.
FIGURE 1 Smoke output is shown from the accompanying code.
There are several ways of simulating the motion of fluids, but they all generally divide into two common styles: grid methods and particle methods. In a grid method, the fluid is represented by dividing up the space a fluid might occupy into individual cells and storing how much of the fluid is in each cell.
In a particle method, the fluid is modeled as a large number of particles that move around and react to collisions with the environment, interacting with nearby particles. Let's focus first on simulating fluids with grids.
The simplest way to discuss the grid method is in respect to a regular twodimensional grid, although the techniques apply equally well in three dimensions. At the most basic level, to simulate fluid in the space covered by a grid you need two grids: one to store the density of liquid or gas at each point and another to store the velocity of the fluid.
Figure 2 shows a representation of this, with each point having a velocity vector and containing a density value (not shown). The actual implementation of these grids in C/C++ is most efficiently done as onedimensional arrays. The amount of fluid in each cell is represented as a float.
The velocity grid (also referred to as a velocity field, or vector field) could be represented as an array of 2D vectors, but for coding simplicity it's best represented as two separate arrays of floats, one for X and one for Y.
FIGURE 2 The fluid density moves over a field of velocities, with a density stored at each point.
In addition to these two grids, we can have any number of other matching grids that store various attributes. Again, each will be stored as a matching array of floats, which can store factors such as the temperature of the fluid at each point or the color of the fluid (whereby you can mix multiple fluids together).
You can also store more esoteric quantities such as humidity, for example, if you were simulating steam or cloud formation.
Anonymous 
26 Jun 2008 at 12:44 pm PST

Semilagrangian advection has been used since a long time in the field of numerical weather forecast. SemiLagrangian schemes presented in the literature have used either two or three time levels for time integration (Staniforth & Cot e 1991); the latter are more expensive but provide greater accuracy.
In graphic application one can use a faster (second order in time) trapezoidal rule for the integration : x_i  x_i* = (u(x_i, t+dt) + u(x_i*, t)) / 2 This is an implicit equation for x_i* (the position back in time) which could be solved by estimating u using the 2 first term of the Taylor expansion around x_i x_i* = x_i  dt ( 1 + dt/2 Grad u(x_i,t)^1 * ((u(x_i,t) + u(x_i, t+dt)/2) where all terms on the RHS are computed at grid point x_i. The biggest problem with Stam algoritm is the use of linear interpolation, which is not appropriated for rendering of fluid. One early drawback of semiLagrangian schemes was that, since the basis functions of common highorder interpolants such as cubic Lagrange or Hermite are not positive definite, a scalar that is initially bounded between 0 and 1 will not necessarily remain so. However, this problem can be easily fixed by a local clipping, as proposed by Bermejo and Staniforth (1992). This local clipping produces a quasimonotone scheme. The term quasimonotone means that the scalar field can develop wiggles, but it cannot exceed its range at the previous time step. This is different from traditional clipping, where undershootsand overshoots are reset to minimum and maximum background values chosen, based on physical considerations, at the beginning of the computation. As shown by Bermejo (2001) this clipping is equivalent to a linear combination of a higher order interpolant and a loworder (monotone) Lagrange interpolant. To guarantee (discrete) conservation of scalar, two approaches are proposed by Priestly (1993) and Bermejo and Conde (2002). Personnally I prefer the one from Priestly (1993), since it can be used with the quasimonotone interpolation scheme previously described. The complete reference is Priestly, A. 1993 A QuasiConservative version of the SemiLagrangian Advection Scheme. Mon. Wea. Rev., 121 


Stéphane Gaudreault 
According to google, patent number: 6266071 is issue since July 24, 2001.
http://www.google.com/patents?id=ySgIAAAAEBAJ&printsec=abstract&z oom=4 


Anonymous 
Normally, you should always avoid the case where two trajectories arrive at the same point back in time. The semilagrangian method is stable even for a CFL number larger than 1, but this does not mean that it is stable for any timestep. In fact, you are limited by your framerate, but there is also a theorical limit. When the advection velocity varies, the timeintegration interval should be bounded for both numerical accuracy and computational stability. A stability condition for semiLagrangian schemes in terms of the socalled Lipschitz number or deformational CFL (DeCFL) number has been discussed by Smolarkiewicz and Pudykiewicz (1992) and Lin and Rood (1996).
see http://ams.allenpress.com/perlserv/?request=getabstract&doi=10.1 175%2F15200493(1996)124%3C2046%3AMFFSLT%3E2.0.CO%3B2 Using this stability condition for the choice of the timestep, any two trajectories cannot cross each other. In this case, there is no need to keep a list of points of origin and mass fraction. 


Anonymous 
Could you post a sample that shows it as a liquid?



Lou Hayt 
During my education at FullSail I made a liquid dynamics project that was based on Stam's articlesource code. I used the densities to form a dynamic height map (this was mentioned in this article), added a floating boat (not polished though), waves, rain effects etc you can download a zip file with my source code together with Stam's article and source code here: http://louhayt.com/Documents/5.%20Fluid%20Dynamics%20Project%20(O
ne%20month%20class%20project).zip from my understanding this code is light enough to perform well on mobile devices, I was happy to read the possible improvements suggested in this article! :) 


Richard Mathie 
Hey why is the pressure term completely neglected in the naiver stokes equation, it should read
Du/Dt =du/dt + (u.grad)u= + nu grad^2 u it dosent matter in the end as the "vorticity confinement, and projection" steps are effectively taking care of it, but the physical reason they are needed is the pressure! 


Richard Mathie 
odd that didnt come out right try
Du/Dt=du/dt + (u.grad)u= ## 1/rho grad P## + nu grad^2 u 


Bryan Rayner 
Using this article as a guideline, I've made a simple fluid simulation for a physics assignment.
However, there are a few problems which arise later  particularly dealing with the decompression step. A special case exists wherein a thick wave gradually gains speed  the cells before it cause the wave to move faster towards them  since the wave is more than one cell thick, it isn't balanced off by the rear. As it moves to a faster and faster velocity, it begins to move across multiple cells at the same time  resulting in an ever increasing velocity which is not slowed in the slightest. This causes multiple glitches in the system  as soon as a thick wave starts, it will grow and fill the entire fluid (as the decompression algorithm will continue to spread it out in a given direction.) My quick fix  dissalow cells that have more than a given amount of fluid to decompress  the backward advection will pull them out, and they will eventually slow down enough to spread out to a normal state. 

