Stable Fluid
Dynamics Solver
A real-time 3D fluid simulator implementing Jos Stam’s Stable Fluids: a semi-Lagrangian Navier-Stokes solver that is unconditionally stable at any time step.
Overview
Fluid motion is governed by the incompressible Navier-Stokes equations, which describe how a velocity field evolves under advection (the fluid carrying itself along), viscous diffusion, external forces, and a pressure field that keeps the fluid from compressing. The catch: a straightforward explicit solver is only stable for tiny time steps. Push the time step or the forces too far and the simulation literally explodes.
This project implements Jos Stam's Stable Fluids, whose key contribution is making every stage of the solver unconditionally stable: the simulation never blows up regardless of time step, which is what makes it practical for real-time use.
The first equation says velocity changes due to advection, viscosity, and forces; the second says the fluid never compresses, so the velocity field must stay divergence-free.
The solver, one step at a time
Each time step splits the equations into four simpler sub-steps applied in sequence (operator splitting):
Add external forces
Gravity, buoyancy, and user-driven impulses are added directly to the velocity field: \(\mathbf{u}\leftarrow\mathbf{u}+\Delta t\,\mathbf{f}\).
Semi-Lagrangian backtrace
Instead of pushing quantities forward (unstable), trace each cell backwards through the velocity field to find where its fluid came from, and interpolate the old field there.
Implicit viscosity solve
Viscosity is applied by solving a sparse linear system \((I-\nu\Delta t\nabla^2)\mathbf{u}'=\mathbf{u}\) rather than stepping explicitly, which stays stable for any viscosity or time step.
Pressure projection
Solve a Poisson equation for pressure and subtract its gradient, projecting the velocity back onto the space of divergence-free flows. This is where the swirling comes from.
By the Helmholtz-Hodge decomposition, any velocity field u* splits into a divergence-free part plus the gradient of a scalar field. Solving this sparse system for \(p\) and subtracting \(\nabla p\) removes exactly the compressible part.
The staggered MAC grid
Velocities live on a staggered Marker-and-Cell grid: pressure is stored at cell centers, while each velocity component is stored on the cell faces perpendicular to it. The finite-difference derivatives needed for divergence and pressure gradients then line up exactly where they're evaluated, avoiding the checkerboard artifacts of collocated grids and keeping the pressure solve well-conditioned.