This page came about after an online discussion of what Claude Code was capable of, or perhaps not capable of. I feared it might actually do an alright job of implementing Stable Fluids, so I pretty blindly threw Claude Code at the task, giving it the pdftotext-extracted text of the Stable Fluids paper. I connected my mcp-observable-notebookkit-debug tool so that it would be able to interact with the simulation and analyze the <canvas> content. After some finagling over parameter calibration to get things in the ballpark, I asked that it implement some initial conditions for fluid instabilities and sharpen things up a little with the monotonic cubic interpolation scheme from Fedkiw, Stam, and Jensen’s Visual Simulation of Smoke. To instruct it about the method, I took a screenshot of the Appendix of the paper and dropped it into the terminal window.
I want to reject AI, but I politely decline to be religious about it, either supporting it unconditionally or rejecting it out of hat. I view it as a tool, nothing more, nothing less. You can produce good things with it, and you can produce bad things with it. Honestly it seems to have done a decent job here. The FFT isn’t an appropriate fit given our application of walls, so the boundaries leak fluid as a result of the projection step. Memory usage could be optimized. I think I’d have organized the shaders a little differently, maybe consolidating some force computation or something. I think the FFT could be reworked to perform better. But gosh overall it’s just not terrible.
A WebGPU implementation of Jos Stam’s classic Stable Fluids algorithm (SIGGRAPH 1999). The method solves the incompressible Navier-Stokes equations using an unconditionally stable semi-Lagrangian advection scheme and FFT-based pressure projection.
Click and drag to add velocity to the fluid, or select an initial condition to observe classic fluid instabilities:
- Interactive: Start with a blank canvas and add velocity/dye with mouse interaction
- Rayleigh-Taylor: Dense fluid (dye) on top falls into lighter fluid below, forming characteristic mushroom patterns
- Kelvin-Helmholtz: Two layers moving in opposite directions create shear at the interface, forming rolling wave patterns
Algorithm
The Stable Fluids algorithm solves the incompressible Navier-Stokes equations:
subject to the incompressibility constraint .
Each timestep consists of four stages:
- Add forces: Apply external forces (mouse input) to the velocity field
- Advect: Move the velocity field along itself using semi-Lagrangian advection
- Diffuse: Apply viscous diffusion (done in frequency domain as multiplication)
- Project: Make velocity divergence-free by subtracting the gradient of pressure
The key insight is that projection can be done efficiently in frequency domain. The velocity is decomposed into components parallel and perpendicular to the wavenumber . The divergence-free projection removes the parallel component:
This makes the algorithm unconditionally stable regardless of timestep size.
Implementation Details
Grid Layout
Velocities are cell-centered (collocated), not staggered. Cell stores both velocity components at physical position . Velocity is stored as vec4<f32> with layout (u_re, u_im, v_re, v_im) to accommodate FFT operations. In the spatial domain, the imaginary parts are zero.
Simulation Loop
Each timestep executes these stages in order:
- Add forces (mouse interaction)
- Advect velocity (semi-Lagrangian)
- Advect dye
- Project + diffuse (FFT-based)
- Vorticity confinement
- Buoyancy
- Boundary enforcement
FFT-Based Projection and Diffusion
The FFT is used to efficiently solve the pressure projection and diffusion in frequency domain:
- Split the interleaved velocity into separate and complex buffers
- Forward FFT on both components
- Project and diffuse in frequency domain:
- Compute wavenumber with frequencies scaled by
- Remove component parallel to :
- Apply implicit diffusion by multiplying by
- Zero the DC component ()
- Inverse FFT to return to spatial domain
- Merge back into interleaved velocity buffer
Forcing Model
The force uses a velocity-matching model rather than direct impulse injection:
where is computed from mouse velocity, and the falloff is Gaussian: .
Dye injection uses a hard circle (step function) rather than Gaussian for sharper smoke edges.
Semi-Lagrangian Advection
For each grid cell, trace backward through the velocity field to find the departure point:
Then sample the source field at using monotone cubic interpolation (Fedkiw et al., “Visual Simulation of Smoke”, SIGGRAPH 2001). The factor of converts velocity from normalized to grid units.
The monotone cubic interpolation uses Hermite interpolation with slope limiting to prevent overshoot:
where the slopes are set to zero when they differ in sign from . This preserves monotonicity and eliminates spurious oscillations that can arise with standard cubic interpolation.
Vorticity Confinement
Two-pass process to counteract numerical dissipation of rotational motion:
Pass 1: Compute vorticity (curl) using central differences:
Pass 2: Apply confinement force in the direction perpendicular to the vorticity gradient:
In 2D, the cross product gives . This force points toward vortex centers, amplifying rotational structures.