Advection

The fundamental operation in grid-based fluid dynamics is advection. Advection is basically the moving of things on the grid, but more specifically, it's moving the quantities stored in one array by the movement vectors stored in the velocity arrays.

It's quite simple to understand what's happening if you think of each point on the grid as an individual particle, with some attribute (density) and a velocity. Then you are familiar with the process of moving a particle by adding the velocity vector to the position vector. On the grid, however, the possible positions are fixed, so all we can do is move (advect) the quantity (density) from one grid point to another.

In addition to advecting the density value, we also need to advect all the other quantities associated with the point. This would include additional attributes such as temperature and color, but also the velocity of the point itself. The process of moving a velocity field over itself is referred to as self-advection.

The grid does not represent a series of discrete quantities, density or otherwise; it actually represents (inaccurately) a smooth surface, with the grid points just being sample points on that surface.

Think of the points as being X,Y vertices of a 3D mesh, with the density field being the Z height. Thus, you can pick any X and Y position on the mesh, and find the Z value at that point by interpolating between the closest four points. Similarly while advecting a value across the grid, the destination point will not fall directly on a grid point, and you'll have to interpolate the value into the four grid points closest to the target position.

FIGURE 3 Forward advection: The value in P moves forward to A, B, C, and D. This dissipates it when moving diagonally.

In Figure 3, point P has a velocity V, which, after a time step of t, will put it in position P´=P+t*V. This point falls between the points A, B, C, and D, and so a bit of P has to go into each of them. Generally, t*V will be significantly smaller than the width of a cell, so one of the points A, B, C, or D will be P itself. There are various inaccuracies when advecting an entire grid like this, particularly that quantities dissipate when moving in a direction that is not axis-aligned. But this inaccuracy can be turned to our advantage.

Stam's Advection

Programmers looking into grid-based fluid dynamics for the first time will most often come across the work of Jos Stam and Ron Fedkiw, particularly Stam's paper "Real-Time Fluid Dynamics for Games," which he presented at the 2003 Game Developers Conference. He describes a very short procedure for making a grid-based fluid simulator.

In particular, he shows how to implement the advection step using what he calls a "linear backtrace," which simply means that instead of moving the point forward in space, he inverts the velocity and finds the source point in the opposite direction, essentially back in time. He then takes the interpolated density value from that source (which, again, will lay between four actual grid points), and moves the value into the point P. See Figure 4 for an example.

Stam's approach produces visually pleasing results, yet suffers from a number of problems. First, the specific collection of techniques discussed may be covered by an existing patent (U.S. #6,266,071), although as Stam notes, backtracing dates back to 1952. Check with a lawyer if this is a concern to you.

On a more practical note, the advection alone as Stam describes it simply does not work accurately unless the velocity field is smooth in a way termed mass conserving, or incompressible.

Consider the case of a vector field in which all the velocities are zero except for one. The velocity cannot move (advect) forward through the field, since there's nothing ahead of it to "pull" it forward. Instead, the velocity simply bleeds backward. The resultant velocity field will terminate at the original point, and any quantities moving through this field will end up there.

We can solve this particular problem by adding a step to the algorithm termed projection. Projection essentially smooths out the velocity by making it incompressible, allowing the backtracing advection to work perfectly and making the paths formed by the velocity "swirly," as if it were real water. The problem with this approach is that projection is quite expensive, requiring 20 iterations over the velocity field in order to "relax" it to a usable state.

Another performance problem with Stam's approach is that he uses a diffusion step, which also involves 20 iterations over a field. It's necessary to allow the gas to spread out from areas of high density to areas of low density. If the diffusion step were missing, solid blocks of the fluid would remain solid as they moved over the velocity field. Diffusion is an important cosmetic step.