Stable Fluids is one of my favorite papers. Its blend of mathematics, numerical analysis, computer graphics, and physics has been a tremendous source of inspiration for me: a few elegant ideas become a simulation one can immediately see and interact with.
A moving cloud of smoke looks complicated. Its solver need not be. Jos Stam’s Stable Fluids [1] is built from a handful of operations with clear geometric meanings: carry a field along trajectories, smooth it, and enforce incompressibility by an orthogonal projection. Its achievement is an accessible synthesis for interactive graphics, built on a much older numerical-analysis tradition.
The trick is not to solve the entire fluid equation at once. Split it into manageable pieces, and choose each piece to avoid a familiar numerical instability. On a periodic domain, the pressure solve becomes particularly elegant: each Fourier coefficient is projected onto a line perpendicular to its wave vector. Stam’s later implementation notes [4] and Bridson’s book [5] are useful companions.
Velocity, pressure, and the constraint
Work on a periodic square \([0,L)^2\). The velocity \(u(t,x)\in\mathbb R^2\) solves the incompressible Navier–Stokes equations
\[\partial_tu+(u\cdot\nabla)u =-\nabla p+\nu\Delta u+f, \qquad \nabla\cdot u=0.\]The force \(f\) pushes the fluid and the viscosity \(\nu\ge0\) smooths velocity variations. Pressure is the multiplier enforcing incompressibility: a parcel may move or deform, but its volume cannot change. Dye is a separate passive scalar \(a\), evolving by
\[\partial_ta+u\cdot\nabla a=\kappa\Delta a+s,\]where \(\kappa\ge0\) is dye diffusivity and \(s\) a dye source. Since \(\nabla\cdot u=0\), the advection term can equally be written \(\nabla\cdot(au)\). Dye does not feed back into the velocity in the simple model used here.
First understand the projection
Suppose \(w\) is a provisional velocity that may compress the fluid. Remove a gradient:
\[\mathsf P w=w-\nabla\phi, \qquad \Delta\phi=\nabla\cdot w, \qquad \int\phi=0.\]The Poisson equation is compatible because a periodic divergence has integral zero. The field \(\mathsf Pw\) is divergence free by construction.
Proposition 1 (the pressure step is an orthogonal projection). In periodic \(L^2\), \(\mathsf Pw\) is the unique closest divergence-free field to \(w\). In particular,
\[\|w\|_{L^2}^2 =\|\mathsf Pw\|_{L^2}^2+\|\nabla\phi\|_{L^2}^2.\]Proof: gradients are perpendicular to incompressible motion
For any divergence-free \(z\), periodic integration by parts gives \(\int\nabla\phi\cdot z=-\int\phi\nabla\cdot z=0\). Hence the decomposition \(w=\mathsf Pw+\nabla\phi\) is orthogonal. More generally,
\[\|w-z\|_{L^2}^2 =\|\nabla\phi\|_{L^2}^2+\|\mathsf Pw-z\|_{L^2}^2,\]so the unique minimum occurs at \(z=\mathsf Pw\). \(\square\)
This is the Helmholtz–Hodge decomposition, or Leray projection. Harlow and Welch’s MAC method (1965) already used staggered velocities and a pressure Poisson equation to enforce incompressibility [6]. Chorin and Temam developed the projection/fractional-step approach [2], [3]; Chorin’s 1968 paper also uses implicit alternating-direction schemes, not merely explicit Euler. Stam explicitly follows Chorin and Marsden’s book for the decomposition [7].
One time step, four operations
Let \(h>0\) be the time step. A first-order splitting in the order of Stam’s 1999 presentation is
\[\begin{aligned} w_1&=u^n+hf^n &&\text{add forces},\\ w_2(x)&=\mathcal I[w_1](x-hw_1(x)) &&\text{advect},\\ (I-h\nu\Delta)w_3&=w_2 &&\text{diffuse},\\ u^{n+1}&=\mathsf Pw_3 &&\text{project}. \end{aligned}\]Here \(\mathcal I\) is periodic bilinear interpolation of the grid values. The pressure potential \(\phi\) in the projection has the units of \(hp\); one should not identify it directly with the continuous pressure without that time-step scaling.
Why trace backwards? To fill a grid point at the next time, ask where the arriving fluid came from. The first-order answer is \(x-hw_1(x)\). Interpolate the old field there. Forward tracing would leave the grid and create a separate problem of depositing particles back onto it.
Dye is advected by the resulting projected velocity, with the same backward-tracing operation, then optionally diffused. A mouse drag supplies a localized velocity impulse and a positive dye source. This is an externally forced simulation, not an unforced energy-decay experiment.
Earlier characteristic methods
Backward tracing also has a substantial pre-1999 history. Robert (1982) combined semi-Lagrangian transport with semi-implicit integration for an atmospheric shallow-water model [8]. Staniforth and Côté’s 1991 review documents the wider meteorological use of these methods [9]: taking large stable steps was already an important motivation, although the equations and stability analyses differ.
An especially close antecedent is Pironneau’s transport–diffusion algorithm (1982) [10, §3]. It transports the evolving velocity itself along backward characteristics, then performs an implicit incompressible Stokes solve. On our periodic domain with constant viscosity, that latter solve is equivalent, before spatial discretization, to implicit diffusion followed by pressure projection: the two operators commute. The distinction is in the realization—Pironneau’s variational finite-element formulation versus Stam’s grid interpolation and separate diffusion/projection steps. Pironneau proves an unconditional energy estimate for his analyzed schemes; it should not be transferred uncritically to the interpolated scheme below.
Graphics had its own predecessors: Kass and Miller (1990) used an implicit solver for a shallow-water approximation [11], while Foster and Metaxas (1997) supplied Stam’s explicit three-dimensional gas-simulation baseline [12]. Stam’s contribution is the simple, robust combination that makes volumetric fluid animation interactive, accepting numerical dissipation. He cites these graphics works and Harlow–Welch; Robert, Pironneau, and Staniforth–Côté are historical comparisons here, not sources cited in his paper.
Fourier projection: one small matrix per wave
For a Fourier mode \(e^{ik\cdot x}\), with \(k=(2\pi/L)(j_1,j_2)\), differentiation becomes multiplication by \(ik\). Thus, for \(k\ne0\),
\[\widehat\phi(k)=-\frac{i\,k\cdot\widehat w(k)}{|k|^2}, \qquad \widehat{\nabla\phi}(k) =\frac{k(k\cdot\widehat w(k))}{|k|^2}.\]
Proposition 2 (spectral pressure and diffusion). The periodic projection is
\[\widehat{\mathsf Pw}(k) =\left(I-\frac{kk^\top}{|k|^2}\right)\widehat w(k) \quad(k\ne0), \qquad \widehat{\mathsf Pw}(0)=\widehat w(0).\]Implicit velocity diffusion multiplies each coefficient by \((1+h\nu\lvert k\rvert^2)^{-1}\). The two operations commute and can be combined into a single Fourier pass.
Proof: solve the scalar Poisson equation mode by mode
The Fourier Poisson equation is \(-\lvert k\rvert^2\widehat\phi=ik\cdot\widehat w\). Solving it gives the formula above; subtracting its gradient gives the projector. The matrix \(I-kk^\top/\lvert k\rvert^2\) removes the component parallel to \(k\), so its output has zero Fourier divergence. The zero mode is the mean velocity, already divergence free, and remains unchanged. Finally the diffusion equation gives \((1+h\nu\lvert k\rvert^2)\widehat w_3=\widehat w_2\). Both multipliers act on each mode independently, and the diffusion multiplier is scalar. \(\square\)
In practice: transform the two velocity components, apply these multipliers, inverse transform. An FFT makes the cost \(O(M\log M)\) for \(M\) grid cells; no iterative pressure solver is hidden inside this step. The experiment uses a genuine mixed-radix FFT on a \(45\times45\) grid.
What “stable” really means
Proposition 3 (three separate stability facts). Bilinear semi-Lagrangian advection preserves the range of each transported scalar component. Implicit spectral diffusion is an \(L^2\) contraction for every \(h>0\). Pressure projection is also an \(L^2\) contraction.
Proof: convex weights, bounded multipliers, and orthogonality
A bilinear interpolated value is a convex combination of four old values, even if the backtrace crosses many grid cells. It therefore lies between the old minimum and maximum. Every diffusion multiplier lies in \((0,1]\), so Parseval’s identity gives its \(L^2\) contraction. Proposition 1, or the orthogonal matrix in Proposition 2 at each Fourier mode, gives the projection contraction. \(\square\)
These are statements in different norms. They do not combine into a universal monotone kinetic-energy theorem for the entire discrete split scheme. In particular, bilinear advection is not generally a grid-mass-preserving operator, nor automatically an \(\ell^2\) contraction. Increasing \(h\) can make a bounded animation very inaccurate: the backtrace is less precise, splitting error grows, and interpolation introduces numerical diffusion. Stability is not an accuracy certificate.
One more discretization detail matters for dye. Spectral implicit diffusion is \(L^2\)-stable but its truncated-grid kernel need not preserve positivity. We instead diffuse dye with the periodic five-point Laplacian \(\Delta_h\). Its implicit inverse is positive, and the same FFT diagonalizes it, with multiplier
\[\left[1+\frac{4h\kappa}{\Delta x^2} \left(\sin^2\frac{\pi j_1}{N} +\sin^2\frac{\pi j_2}{N}\right)\right]^{-1}.\]The maximum-principle argument is elementary: at a negative minimum of a hypothetical solution of \((I-h\kappa\Delta_h)a=b\ge0\), one has \(\Delta_h a\ge0\), so the left side is negative, a contradiction. This diffusion step also preserves the mean. The advection step can still change dye mass; the display reports it honestly rather than renormalizing it away.
Numerical experimentation
Paint and stir, adjust the brush radius and stirring strength, compare viscosity levels, and increase the time step. Zero stirring strength paints without pushing. Every stroke is a discrete impulse, and its actual grid-integrated dye increment is recorded: the reference mass is the initial mass plus all these increments. The signed difference from that reference exposes the advection error. Turn projection off to see what is lost, then turn it back on and take one step: the spectral divergence falls back to roundoff. This is an educational low-resolution periodic solver, not a turbulence benchmark.
There is a useful companion viewpoint in Moser’s interpolation of measures. There, one prescribes a density change and solves for a flux having the right divergence. Here, one starts with a velocity and removes its divergence. In both constructions, a Poisson equation turns conservation of mass into a practical operation.
Bibliography
- Jos Stam. Stable Fluids. SIGGRAPH, 121–128, 1999. Paper; PDF.
- Alexandre J. Chorin. Numerical Solution of the Navier–Stokes Equations. Mathematics of Computation 22(104), 745–762, 1968. Paper.
- Roger Temam. Sur l’approximation de la solution des équations de Navier–Stokes par la méthode des pas fractionnaires (I). Archive for Rational Mechanics and Analysis 32, 135–153, 1969. Paper.
- Jos Stam. Real-Time Fluid Dynamics for Games. Game Developers Conference, 2003. Author’s publications.
- Robert Bridson. Fluid Simulation for Computer Graphics. Second edition, A K Peters/CRC Press, 2015. Author’s materials.
- Francis H. Harlow and J. Eddie Welch. Numerical Calculation of Time-Dependent Viscous Incompressible Flow of Fluid with Free Surface. Physics of Fluids 8(12), 2182–2189, 1965. Paper.
- Alexandre J. Chorin and Jerrold E. Marsden. A Mathematical Introduction to Fluid Mechanics. Second edition, Springer, 1990. Book.
- André Robert. A Semi-Lagrangian and Semi-Implicit Numerical Integration Scheme for the Primitive Meteorological Equations. Journal of the Meteorological Society of Japan, Series II 60(1), 319–325, 1982. Paper.
- Andrew Staniforth and Jean Côté. Semi-Lagrangian Integration Schemes for Atmospheric Models—A Review. Monthly Weather Review 119(9), 2206–2223, 1991. Paper.
- Olivier Pironneau. On the Transport-Diffusion Algorithm and Its Applications to the Navier–Stokes Equations. Numerische Mathematik 38, 309–332, 1982. Paper.
- Michael Kass and Gavin Miller. Rapid, Stable Fluid Dynamics for Computer Graphics. Computer Graphics (SIGGRAPH ’90) 24(4), 49–57, 1990. Paper.
- Nick Foster and Dimitris Metaxas. Modeling the Motion of a Hot, Turbulent Gas. SIGGRAPH, 181–188, 1997. Paper.