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. I asked it to 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.

After some finagling with parameters, it came together pretty quickly. It has issues. The FFT isn’t a good solver for non-periodic walls, so the boundaries leak fluid as a result of the projection step. Memory usage could be optimized. I’d have organized the shaders differently, probably consolidated some of the force computation. The FFT could be reworked to perform better. But overall it’s honestly not terrible, and that’s the unsettling part.

Reflecting on the results

The creation of other notebooks on this site involved varying degrees of toil. I worked quite hard on some of them, slowly building tools and techniques over the course of years. But this notebook I vaguely willed into existence while making dinner. I certainly feel conflicted about the result. On one hand, I want to believe that the result is judged by the result. On the other hand, I want to believe that meaning arises from the process.

In an ideal world, this would be moot because of the value this site has for other people. Truth is though, I’m not sure this site has much meaning to other people, and I’m definitely not sure it’s useful.

Most of the notebooks on this on this site link to original notebooks on observablehq.com which were definitely created by manual, human labor. I translated them to Observable Notebook Kit largely by way of AI. I improved many things along the way. The camera interactions are better. I translated a number of notebooks from WebGL to WebGPU. The presentation is cleaned up. I added explanatory figures which I definitely would not have taken the time to hand-craft with SVG.

Does this retroactively cheapen the hard work I put into these notebooks? I don’t know. Probably. Am I more pleased than I’ve ever been before about creating things with limited time and resources that match my vision for them? Yeah, definitely.

The future is going to get weird. The view that AI worthless nonsense because it doesn’t work is, I think, wildly outdated and only becoming more so. The view that it’s destructive is probably right. But it’s creative as well and enables ideas and approaches that could not have existed before. I’m not saying just vibe-code everything and pump out garbage. But there are things I can do now that I simply would not have done a year ago becuase I’m willing to devote entire weekends to messing around with low-level click handlers. Maybe constructing worlds that match your vision will become a new norm. This feels more old-school geocities-style to me than like an apocalyptic vibe-coded hellscape.

So I don’t know. I’m fearful but hopeful.

---

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:

1. Add forces: Apply external forces (mouse input) to the velocity field
2. Advect: Move the velocity field along itself using semi-Lagrangian advection
3. Diffuse: Apply viscous diffusion (done in frequency domain as multiplication)
4. Project: Make velocity divergence-free by subtracting the gradient of pressure

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:

1. Add forces (mouse interaction)
2. Advect velocity (semi-Lagrangian)
3. Advect dye
4. Project + diffuse (FFT-based)
5. Vorticity confinement
6. Buoyancy
7. Boundary enforcement

FFT-Based Projection and Diffusion

The FFT is used to efficiently solve the pressure projection and diffusion in frequency domain:

1. Split the interleaved velocity into separate and complex buffers
2. Forward FFT on both components
3. 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 ()
4. Inverse FFT to return to spatial domain
5. 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.