Eulerian Fluids: Stable Fluids

The rolling smoke in a film, the fireball of an explosion, the churn of a river surface — a huge share of the fluid you see in modern computer graphics traces back to one remarkably short SIGGRAPH paper: Jos Stam's "Stable Fluids" (1999). Before it, fluid solvers in animation were fragile: push the timestep a little too hard and the simulation would explode into a screenful of NaNs. Stam's trick made the solver unconditionally stable — you could take big, movie-sized timesteps and it would never blow up. That single property is why grid-based fluids became practical for production, and the method still sits at the heart of tools like FumeFX, Phoenix and Houdini's smoke solver.

This page builds the solver from its idea outward: the Eulerian way of storing a fluid, the equations that govern it, and Stam's four-step recipe — add forces, advect, diffuse, project — with the two clever pieces (semi-Lagrangian advection and pressure projection) that make it work. It assumes you have met the Navier–Stokes equations and know your way around computer animation.

Two ways to watch a fluid: Eulerian vs Lagrangian

There are two fundamentally different vantage points for describing a moving fluid, and choosing one is the first design decision in any solver.

Lagrangian — follow the parcels. You track a swarm of moving particles, each carrying its own velocity and other quantities, and you ride along with them. This is the world of SPH particle fluids.
Eulerian — watch fixed windows. You lay down a stationary grid of cells and, at each cell, store the fluid quantities right now — the velocity \mathbf{u} passing through and the smoke density \rho sitting there. The fluid flows past the cells; the cells never move.

Stable Fluids is thoroughly Eulerian. Picture a fixed checkerboard laid over your smoke: every square remembers two things, a little velocity arrow (which way and how fast the air is moving there) and a greyness (how much smoke is in that square). Every timestep the solver rewrites those two fields. The great advantage of a grid is that spatial derivatives — the gradients and Laplacians all over Navier–Stokes — become simple differences between neighbouring cells, and pressure can be solved globally across the whole grid at once. The cost is that fine detail can only be as sharp as the grid is fine, a limitation we will meet head-on later.

The equations we are solving

A moving fluid of constant density obeys the incompressible Navier–Stokes equations. The first says how the velocity field \mathbf{u} changes moment to moment; the second is the incompressibility constraint that must hold at all times:

\frac{\partial \mathbf{u}}{\partial t} = -(\mathbf{u}\cdot\nabla)\mathbf{u} \;-\; \frac{1}{\rho}\nabla p \;+\; \nu\nabla^2\mathbf{u} \;+\; \mathbf{f}, \qquad \nabla\cdot\mathbf{u} = 0.

Read the top equation term by term — each term is one thing the fluid does, and Stam's solver handles each with its own dedicated step:

-(\mathbf{u}\cdot\nabla)\mathbf{u} — advection: the fluid carries itself (and everything in it) along its own flow.
-\tfrac{1}{\rho}\nabla p — the pressure gradient: pressure pushes fluid from crowded regions to empty ones.
\nu\nabla^2\mathbf{u} — diffusion (viscosity): neighbouring parcels drag on each other and smear velocity out; \nu is the viscosity.
\mathbf{f} — external forces: gravity, buoyancy, the fan you switch on, a stir of the mouse.

The constraint \nabla\cdot\mathbf{u} = 0 says the velocity field is divergence-free: no cell may be a net source or sink. Divergence at a cell would mean more fluid flowing out than in (or vice-versa) — fluid appearing from nowhere or vanishing into it. For an incompressible fluid that is forbidden, and enforcing it is the subtle heart of the algorithm.

Operator splitting: solve one term at a time

Trying to satisfy that whole coupled system in one shot is hard. Stam's move — operator splitting — is to march through the terms one after another, each a small, well-understood sub-problem, feeding the output of one into the next. Over a single timestep \Delta t the field passes through four stations:

1. Add forces. Nudge the velocity by \mathbf{f}\,\Delta t — buoyancy, gravity, user pushes.
2. Advect. Transport the field along the current velocity — this is the -(\mathbf{u}\cdot\nabla)\mathbf{u} term, done by the semi-Lagrangian trick below.
3. Diffuse. Apply viscosity \nu\nabla^2\mathbf{u} by an implicit solve (an inviscid smoke sim can skip or minimise this).
4. Project. Subtract the pressure gradient so that the result is divergence-free — this is where \nabla\cdot\mathbf{u}=0 is enforced.

After the four stations the velocity field is ready for the next frame, and the smoke density is advected along that same clean velocity so the smoke rides the flow. Splitting is an approximation — the terms genuinely interact — but with small enough \Delta t it is accurate, and it turns one intractable problem into four tractable ones.

The key trick: semi-Lagrangian advection

Advection is "move the field along the flow", and the naïve way to do it — step each quantity forward along its velocity — is exactly what used to blow up: a big timestep sends stuff flying past several cells and the scheme goes unstable. Stam's inspired reversal is to work backward.

To find the new value at grid cell \mathbf{x}, ask: where was the fluid that is now arriving at \mathbf{x}? Follow the velocity backward for one timestep to the departure point \mathbf{x} - \mathbf{u}(\mathbf{x})\,\Delta t, sample the old field there (interpolating between the surrounding cells, since it rarely lands exactly on a grid point), and copy that value to \mathbf{x}.

Why is this unconditionally stable? Because the new value is an interpolation of values that already existed in the old field. An average of existing numbers can never be larger than the largest of them, so no matter how big \Delta t is, the field can never grow without bound. The backtrace might overshoot wildly on a huge timestep, but the worst that happens is you sample the wrong cell — you get an inaccurate answer, never an exploding one. That trade — sacrifice some accuracy to buy guaranteed stability — is the whole reason Stable Fluids earned its name.

Projection: the Helmholtz–Hodge decomposition

After forces, advection and diffusion, the velocity field is generally not divergence-free — those steps don't respect the incompressibility constraint. The projection step fixes that, and its justification is a clean theorem from vector calculus.

Any sufficiently smooth vector field \mathbf{w} splits uniquely into two pieces: a divergence-free part \mathbf{u} and the gradient of a scalar field p, \mathbf{w} = \mathbf{u} + \nabla p with \nabla\cdot\mathbf{u} = 0.
The two parts are orthogonal, so the divergence-free part is recovered by simply subtracting the gradient: \mathbf{u} = \mathbf{w} - \nabla p.

The scalar p is exactly the pressure. To find it, take the divergence of both sides of \mathbf{w} = \mathbf{u} + \nabla p. Since \nabla\cdot\mathbf{u} = 0, the divergence-free part drops out and we are left with a Poisson equation for the pressure:

\nabla^2 p = \nabla\cdot\mathbf{w}.

On the grid this is a big sparse linear system (each cell's pressure couples to its neighbours), solved iteratively — Gauss–Seidel, Jacobi or a multigrid/conjugate-gradient method. Once p is known, subtract its gradient, \mathbf{u} = \mathbf{w} - \nabla p, and the velocity is divergence-free. Projection is what guarantees fluid neither appears nor vanishes — it is conservation of mass, enforced globally, every single frame.

Worked example: one step on a tiny grid

Take a small 3\times 3 patch of cells with a puff of smoke sitting in the bottom-middle cell and, for now, zero velocity everywhere. Walk it through one Stable-Fluids step.

Add forces. Warm smoke is buoyant, so apply an upward force \mathbf{f} = (0, +g) to the smoke-filled cell. After \mathbf{u} \mathrel{+}= \mathbf{f}\,\Delta t that cell now holds an upward velocity arrow; its neighbours are still at rest.

Advect. Move the smoke density along \mathbf{u}. For the cell above the puff, backtrace down its (borrowed) upward velocity to a departure point that lands inside the smoky cell, sample the density there, and copy it up. The grey blob has shifted one step upward — the smoke is rising, exactly as buoyancy should make it.

Project. That upward push has made the velocity field diverge: fluid is streaming out of the top of the puff cell with nothing filling in from below, so \nabla\cdot\mathbf{u} \ne 0 there — a phantom source. Solve \nabla^2 p = \nabla\cdot\mathbf{u} across the nine cells; the pressure comes out high where fluid is piling up and low where it is thinning. Subtracting \nabla p curls the straight-up arrows into a gentle circulation — fluid is drawn in from the sides to replace what rose. Now every cell has as much flowing in as out. The smoke rose, and not a single gram of fluid was created or destroyed. That is one frame; repeat a few hundred times and you have a rising plume of smoke.

Grid-based fluids in computational physics predate Stam by decades — but physicists cared about accuracy and were willing to take tiny, carefully-bounded timesteps (the Courant–Friedrichs–Lewy condition: a parcel must not cross more than one cell per step). Animation has the opposite priorities — it wants big timesteps, it wants to never crash mid-render, and it can trade a little physical accuracy for a lot of robustness because the goal is a plausible-looking fluid, not a certified-correct one. Stam's genius was recognising that the graphics community's different value function made the unconditionally-stable semi-Lagrangian scheme — too diffusive for a physicist — exactly the right tool. It is a lovely case of a field's constraints, not just its techniques, driving the breakthrough.

The very interpolation that makes advection stable is also its curse: every backtrace averages nearby cells, and averaging is smoothing. Step after step, that repeated blurring is numerical dissipation — it bleeds energy out of the fine-scale motion and rounds off exactly the small vortices and crisp filaments that make smoke look alive. A raw Stable-Fluids sim tends to look too smooth, too viscous, as if the air were syrup, even when you set the viscosity \nu to zero.

The fix is to put the lost detail back. Two standard remedies: vorticity confinement (measure where curl is being lost and inject a small force that pushes it back up, re-energising the swirls), and higher-order advection — BFECC and MacCormack schemes, which run the semi-Lagrangian step forward and backward to estimate its own error and cancel most of it, giving far sharper smoke at a modest extra cost while keeping the stability. If your smoke looks like fog, reach for one of these before you reach for a finer grid.