Gamasutra - Practical Fluid Dynamics: Part 1

[In this technical article originally printed in Game Developer magazine, Neversoft co-founder 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 two-part 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.

Grids and Particles

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 two-dimensional 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 one-dimensional 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.

Comments

Anonymous

26 Jun 2008 at 12:44 pm PST

Semi-lagrangian advection has been used since a long time in the field of numerical weather forecast. Semi-Lagrangian 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 semi-Lagrangian schemes was that, since the basis functions of common high-order 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 quasi-monotone 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 under-shootsand over-shoots 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 low-order (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 quasi-monotone interpolation scheme previously described. The complete reference is

Priestly, A. 1993 A Quasi-Conservative version of the Semi-Lagrangian Advection Scheme. Mon. Wea. Rev., 121

Stéphane Gaudreault

26 Jun 2008 at 12:55 pm PST

According to google, patent number: 6266071 is issue since July 24, 2001.

http://www.google.com/patents?id=ySgIAAAAEBAJ&printsec=abstract&zoom=4

Anonymous

26 Jun 2008 at 2:11 pm PST

Normally, you should always avoid the case where two trajectories arrive at the same point back in time. The semi-lagrangian 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 time-integration interval should be bounded for both numerical accuracy and computational stability. A stability condition for semi-Lagrangian schemes in terms of the so-called 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=get-abstract&doi=10.1175%2F1520-0493(1996)124%3C204
6%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

26 Jun 2008 at 3:47 pm PST

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

Lou Hayt

29 Jun 2008 at 8:45 am PST

During my education at FullSail I made a liquid dynamics project that was based on Stam's article\source 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(One%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! :)