Physics Fluid Dynamics WebGL GLSL

Solving Navier–Stokes at 60 FPS: Stable Fluids & Helmholtz–Hodge Projection in WebGL

October 31, 2026 • ☕ 11 min read

Live GPU Navier–Stokes solver. Drag across the canvas to inject momentum and dye, or open ⚙ Controls to switch the field view between Dye ($q$), Vorticity ($\omega = \nabla \times \vec u$), Pressure ($p$), and Velocity ($\vec u$). Click Kármán Street to drop a solid cylinder into a wind tunnel and watch an alternating von Kármán vortex street shed in real time.

In our previous posts on Schwarzschild Gargantua, Morris–Thorne Wormholes, Kerr Spacetime, and Binary Black Hole Mergers, the GPU was stateless: every pixel integrated its own independent ordinary differential equation from the camera to infinity and forgot the result on the next frame.

Simulating a turbulent fluid flips that paradigm upside down. Instead of tracing independent rays through a fixed metric, the GPU must evolve a coupled system of nonlinear partial differential equations (PDEs) where every grid cell talks to its neighbors across time. Using ping-pong floating-point Framebuffer Objects (FBOs), Jos Stam's Stable Fluids (1999) method, and Fedkiw, Stam & Jensen's Vorticity Confinement (2001), a single HTML file can execute ~30 full-grid shader passes per frame at 60 fps—turning your browser into a real-time computational fluid dynamics (CFD) wind tunnel.

0. Anatomy of a Real-Time GPU Fluid Solver

The simulation above solves the 2D incompressible Euler / Navier–Stokes equations coupled to a passive scalar dye concentration $q(\vec x, t)$:

Momentum PDE: $\partial_t \vec u + (\vec u \cdot \nabla)\vec u = -\frac{1}{\rho}\nabla p - \nu_d \vec u + \vec f_{\text{ext}} + \vec f_{\text{vort}}$, evolving the 2D velocity field $\vec u = (u_x, u_y)$.
Incompressibility: $\nabla \cdot \vec u = \partial_x u_x + \partial_y u_y = 0$ everywhere—fluid volume is neither created nor compressed.
Ping-Pong FBOs: Because WebGL cannot read from and write to the same texture simultaneously, velocityFBO, dyeFBO, and pressureFBO each hold a read and write RGBA16F texture pair swapped after every pass.
Advection: Semi-Lagrangian backward-characteristic lookup $\vec x_{\text{back}} = \vec x - \Delta t\,\vec u(\vec x)$, unconditionally stable for any Courant number $\text{CFL} = |\vec u|\Delta t/\Delta x$.
Pressure Projection: 26 Jacobi relaxation passes per frame solving the Poisson equation $\nabla^2 p = \frac{\rho}{\Delta t}\nabla \cdot \vec u^*$, followed by subtracting $\frac{\Delta t}{\rho}\nabla p$.

1. Operator Splitting the Navier–Stokes Equations

Written out in full, an incompressible fluid with velocity $\vec u(\vec x, t)$, pressure $p(\vec x, t)$, and passive dye concentration $q(\vec x, t)$ obeys:

Incompressible Navier–Stokes + Passive Scalar Transport $$\frac{\partial \vec u}{\partial t} + (\vec u \cdot \nabla)\vec u = -\frac{1}{\rho}\nabla p + \nu \nabla^2 \vec u + \vec f, \qquad \nabla \cdot \vec u = 0, \qquad \frac{\partial q}{\partial t} + (\vec u \cdot \nabla)q = -\gamma_q q + S_q$$ Notice that there is no time-derivative equation $\partial_t p$ for pressure! In an incompressible flow, the speed of sound is infinite: pressure acts as an instantaneous global Lagrange multiplier that adjusts itself everywhere at once to keep $\nabla \cdot \vec u = 0$.

Instead of integrating all terms simultaneously, Jos Stam's operator-splitting pipeline advances the state across each time step $\Delta t$ through a sequence of simpler fractional steps, where the output texture of each shader pass becomes the input texture of the next:

Fractional-Step Operator Splitting Pipeline per Frame $\Delta t$ $$\vec u^{(0)} \;\xrightarrow{\text{1. Splat + Vorticity}}\; \vec u^{(1)} \;\xrightarrow{\text{2. Self-Advect}}\; \vec u^* \;\xrightarrow{\text{3. Poisson } \nabla^2\tilde p = \nabla\cdot\vec u^*}\; \vec u^{(n+1)} = \vec u^* - \nabla\tilde p \;\xrightarrow{\text{4. Dye Advect}}\; q^{(n+1)}$$ Crucially, the Helmholtz–Hodge pressure projection (Step 3) runs after velocity self-advection (Step 2) and before passive dye advection (Step 4), guaranteeing that the dye is transported by a freshly projected, strictly divergence-free velocity field.

2. Semi-Lagrangian Advection: Why Backward Rays Beat the CFL Limit

If you discretize the nonlinear advection term $(\vec u \cdot \nabla)\vec u$ with standard forward finite differences, explicit Euler blows up whenever a fluid parcel crosses more than one grid cell per time step ($\text{CFL} = |\vec u|\Delta t/\Delta x > 1$). Fast mouse swipes or high-velocity jets would instantly NaN the simulation.

Stam's key insight was to use the method of characteristics backward in time. The advection equation $\frac{Dq}{Dt} = \partial_t q + (\vec u \cdot \nabla)q = 0$ states that $q$ is constant along particle trajectories $\frac{d\vec x}{dt} = \vec u$. To find the value arriving at grid cell center $\vec x$ at time $t + \Delta t$, we ask: where was this fluid parcel at time $t$?

Unconditionally Stable Backward-Characteristic Advection (with Implicit Damping) $$\vec x_{\text{back}} = \vec x - \Delta t\,\vec u(\vec x, t), \qquad q(\vec x, t + \Delta t) = \frac{\mathrm{bilerp}\big(q(\cdot, t),\; \vec x_{\text{back}}\big)}{1 + \gamma\,\Delta t}$$ Because bilinear interpolation is a convex combination of four grid neighbors ($w_k \ge 0$, $\sum_{k=1}^4 w_k = 1$), the new value is strictly bounded by $\min q(\cdot, t) \le q(\vec x, t+\Delta t) \le \max q(\cdot, t)$ for any $\Delta t > 0$. Dividing by $1 + \gamma\Delta t$ applies an unconditionally stable backward-Euler decay.

💻 GLSL: Semi-Lagrangian Backward-Characteristic Advection PassGLSL

void main() {
  // Zero velocity and dye inside the solid Kármán cylinder obstacle
  if (uObstacle > 0.5) {
    vec2 dp = vUv - uObsPos;
    dp.x *= uAspect;
    if (length(dp) < uObsRad) {
      gl_FragColor = vec4(0.0, 0.0, 0.0, 1.0);
      return;
    }
  }

  // Absorb dye near walls in closed mode, or fade near outlet in wind-tunnel mode
  float wallFade = smoothstep(0.0, 0.025, vUv.x)
                 * (uObstacle > 0.5 ? smoothstep(1.0, 0.94, vUv.x) : smoothstep(1.0, 0.975, vUv.x))
                 * smoothstep(0.0, 0.025, vUv.y) * smoothstep(1.0, 0.975, vUv.y);

  // Trace characteristic backward in time: x_back = x - dt * u(x)
  vec2 vel   = sampleBilerp(uVelocity, vUv, uVelTexel).xy;
  vec2 coord = clamp(vUv - uDt * vel * uVelTexel, uSourceTexel * 0.5, 1.0 - uSourceTexel * 0.5);

  // Hardware linear texture fetch (with 4-tap manual bilerp fallback) + nonlinear decay
  vec4 result = sampleBilerp(uSource, coord, uSourceTexel);
  float peak  = max(result.r, max(result.g, result.b));
  float decay = 1.0 + (uDissipation + 0.55 * peak) * uDt;

  gl_FragColor = vec4((result.rgb / decay) * mix(1.0, wallFade, 0.15), 1.0);
}

3. Vorticity Confinement: Defeating Numerical Viscosity

Unconditional stability comes with a catch: bilinear interpolation averages four neighboring texels on every single frame. Taylor-expanding that interpolation reveals an artificial numerical viscosity $\nu_{\text{num}} \sim \frac{\Delta x^2}{\Delta t}$ that damps out small-scale rotational eddies. Drag the Vorticity Confinement · $\varepsilon_{\text{curl}}$ slider in the controls down to 0 and watch how quickly the fine Kelvin–Helmholtz billows smooth out into laminar syrup.

To restore crisp sub-grid turbulence without giving up an Eulerian grid, Fedkiw, Stam & Jensen (2001) introduced vorticity confinement. In 2D, the curl of the velocity field is a scalar signed vorticity $\omega = (\nabla \times \vec u)_z = \partial_x u_y - \partial_y u_x$. The shader locates vortex cores by computing the unit gradient of vorticity magnitude $\hat N = \nabla|\omega| / |\nabla|\omega||$ (which points radially inward toward each vortex center) and injects a tangential paddle-wheel force $\vec f_{\text{vort}} \propto \hat N \times \vec\omega$:

2D Scalar Vorticity $\omega$ and Fedkiw–Stam–Jensen Vorticity Confinement Force $$\omega = \frac{\partial u_y}{\partial x} - \frac{\partial u_x}{\partial y}, \qquad \hat N = \frac{\nabla |\omega|}{\big|\nabla |\omega|\big| + \epsilon}, \qquad \vec f_{\text{vort}} = \varepsilon_{\text{curl}}\,\Delta x\,\big(\hat N \times \vec\omega\big) = \varepsilon_{\text{curl}}\,\Delta x\,\begin{pmatrix} \hat N_y\,\omega \\ -\hat N_x\,\omega \end{pmatrix}$$ Switch the field visualization pill to Vorticity ($\omega$) in the controls above: counter-clockwise eddies ($\omega > 0$) glow fiery orange-white, while clockwise eddies ($\omega < 0$) glow electric cyan.

💻 GLSL: Scalar Curl Pass & Vorticity Confinement ForceGLSL

/* --- Pass A: Scalar Vorticity omega = du_y/dx - du_x/dy --- */
float L = texture2D(uVelocity, vUv - vec2(uTexel.x, 0.0)).y;
float R = texture2D(uVelocity, vUv + vec2(uTexel.x, 0.0)).y;
float B = texture2D(uVelocity, vUv - vec2(0.0, uTexel.y)).x;
float T = texture2D(uVelocity, vUv + vec2(0.0, uTexel.y)).x;
float vorticity = 0.5 * ((R - L) - (T - B));

/* --- Pass B: Paddle-Wheel Confinement Force f = eps * (N x omega) --- */
float wL = abs(texture2D(uCurl, vUv - vec2(uTexel.x, 0.0)).x);
float wR = abs(texture2D(uCurl, vUv + vec2(uTexel.x, 0.0)).x);
float wB = abs(texture2D(uCurl, vUv - vec2(0.0, uTexel.y)).x);
float wT = abs(texture2D(uCurl, vUv + vec2(0.0, uTexel.y)).x);
float wC = texture2D(uCurl, vUv).x;

vec2 gradMag = 0.5 * vec2(wR - wL, wT - wB);
vec2 N       = gradMag / (length(gradMag) + 1e-4);
vec2 force   = uCurlStrength * vec2(N.y * wC, -N.x * wC);
vec2 vel     = texture2D(uVelocity, vUv).xy + force * uDt;

4. Helmholtz–Hodge Decomposition & the Poisson Pressure Solver

After force injection and nonlinear self-advection, the intermediate velocity field $\vec u^*$ has non-zero divergence ($\nabla \cdot \vec u^* \neq 0$): some cells have net inflow, others net outflow. How do we strip out the compressible part while keeping every bit of the swirl?

The answer is the Helmholtz–Hodge Decomposition Theorem: any smooth vector field $\vec u^*$ decomposes uniquely into the sum of a divergence-free (solenoidal) vector field $\vec u$ ($\nabla \cdot \vec u = 0$) and a curl-free scalar gradient $\nabla \tilde p$ ($\nabla \times \nabla \tilde p = \vec 0$):

Helmholtz–Hodge Decomposition & Discrete Jacobi Poisson Iteration $$\vec u^* = \vec u + \nabla \tilde p \quad\xrightarrow{\;\nabla\cdot(\dots)\;}\quad \nabla^2 \tilde p = \nabla \cdot \vec u^*, \qquad \vec u = \vec u^* - \nabla \tilde p$$ $$\tilde p^{(k+1)}_{i,j} = \frac{\tilde p^{(k)}_{i+1,j} + \tilde p^{(k)}_{i-1,j} + \tilde p^{(k)}_{i,j+1} + \tilde p^{(k)}_{i,j-1} - D_{i,j}}{4}, \qquad D_{i,j} = \frac{u^*_{x,i+1,j} - u^*_{x,i-1,j} + u^*_{y,i,j+1} - u^*_{y,i,j-1}}{2}$$ Here $\tilde p \equiv \frac{\Delta t}{\rho h}p$ absorbs the grid spacing $h$ and time step $\Delta t$. Warm-starting $\tilde p^{(0)}$ from the previous frame's pressure field allows 26 Jacobi ping-pong passes per frame to converge smoothly across the entire grid.

Seeing pressure and vortex streets in action: click the Kármán Street preset, then switch the view pill to Pressure ($p$). You will see a bright golden high-pressure stagnation zone piled up against the front face of the cylinder ($\partial p/\partial n = 0$) and a trail of blue low-pressure pockets marking the cores of the alternating vortices shedding behind it!

💻 GLSL: Jacobi Poisson Iteration & Gradient Projection (with Neumann Obstacles)GLSL

/* --- Jacobi Iteration Pass (run N = 26 times per frame via ping-pong FBO swap) --- */
float pC = texture2D(uPressure, vUv).x;
// Pure Neumann boundary condition (dp/dn = 0): solid obstacle neighbors mirror center pressure pC
float L  = isSolid(uvL) ? pC : texture2D(uPressure, uvL).x;
float R  = (uObstacle > 0.5 && uvR.x > 1.0) ? 0.0 : (isSolid(uvR) ? pC : texture2D(uPressure, uvR).x);
float B  = isSolid(uvB) ? pC : texture2D(uPressure, uvB).x;
float T  = isSolid(uvT) ? pC : texture2D(uPressure, uvT).x;
float div = texture2D(uDivergence, vUv).x;

float pNew = (L + R + B + T - div) * 0.25;

/* --- Helmholtz-Hodge Gradient Subtraction Pass: u = u* - grad(p) --- */
vec2 vel = texture2D(uVelocity, vUv).xy - 0.5 * vec2(R - L, T - B);

5. The Multi-Pass Ping-Pong FBO Pipeline

Unlike a single-pass ray tracer, each animation frame orchestrates 31 sequential WebGL draw calls across five offscreen framebuffers before compositing to the screen:

Figure 1: Each frame executes 5 advection/calculus passes plus $N_{\text{iter}} = 26$ Jacobi relaxation passes on RGBA16F ping-pong FBOs before rendering the final field.

6. Two Paradigms of GPU Physics: Ray-Tracing vs. Ping-Pong PDEs

Stepping back to compare the relativistic ray tracers from Parts 1–4 with this Navier–Stokes solver shows the two fundamental ways to use a fragment shader as a physics engine:

Dimension	Relativistic Ray Tracers (Black Holes / Wormhole)	Eulerian PDE Solver (Stable Fluids)
Mathematical problem	Initial-value ODEs ($d^2 x^\mu / d\lambda^2 + \Gamma^\mu_{\alpha\beta} \dot x^\alpha \dot x^\beta = 0$)	Coupled nonlinear PDEs + elliptic boundary-value problem ($\nabla^2 p = \nabla\cdot\vec u^*$)
Frame-to-frame state	Stateless (pure function of time uniform `uTime` and camera pose)	Stateful (`RGBA16F` ping-pong FBO textures persist and evolve across frames)
Draw calls per frame	1 fullscreen geodesic pass + 1 bloom pass	5 operator passes + $N_{\text{iter}}$ (26–48) Jacobi Poisson passes
Pixel coupling	Embarrassingly parallel (pixels never communicate)	Tightly coupled via 5-point finite-difference stencils ($L, R, B, T, C$)
Key stability trick	Symplectic Velocity-Verlet + Mino time / Fermat unit tangent	Backward-characteristic advection + Helmholtz–Hodge projection

Try dragging the Poisson Jacobi Iterations · $N_{\text{iter}}$ slider down to 2 while in Kármán Street mode, then slide it back up to 40—watching the flow transition from a compressible, leaky gas back into a tight, incompressible vortex street shows why Helmholtz–Hodge projection is the heartbeat of computational fluid dynamics.