Navier-Stokes Equations in WebGL — real-time fluid simulation

1. The incompressible Navier-Stokes equations

The incompressible Navier-Stokes equations describe the motion of a constant-density Newtonian fluid. In vector form they are two coupled PDEs:

The continuity equation ∇ · u = 0 is the incompressibility constraint: fluid cannot compress, so divergence of the velocity field must be zero everywhere. This constraint drives the pressure solver.

The momentum equation has four terms on the right:

• Advection — −(u · ∇)u: fluid carries its own velocity (non-linear, the hard part)
• Diffusion — ν ∇²u: viscous smoothing of velocity differences
• Pressure gradient — −(1/ρ)∇p: pressure differences drive flow
• Body forces — f: external inputs (mouse splat, gravity)

2. Operator splitting — the Stable Fluids approach

Stam's method advances the solution by one timestep Δt using four sequential sub-steps, each applied to a field grid stored as a floating-point texture:

3. Semi-Lagrangian advection

The idea: instead of asking "where does the fluid at grid cell x go?", ask "where did the fluid at x come from?" Trace the velocity field backwards by Δt and sample the previous field there.

In GLSL (one channel of a floating-point texture = velocity component):

// Advect shader — gl_FragCoord.xy maps to UV in [0,1]²
uniform sampler2D uVelocity;   // current velocity field (vec2 encoded as RG)
uniform sampler2D uSource;     // field to advect (velocity or dye)
uniform vec2      uTexelSize;  // 1/resolution
uniform float      uDt;

void main() {
  vec2 uv   = gl_FragCoord.xy * uTexelSize;
  vec2 vel  = texture(uVelocity, uv).rg;
  vec2 prev = uv - uDt * vel * uTexelSize; // backtrace in texture space
  gl_FragColor = texture(uSource, prev);  // bilinear sample (FREE in GPU)
}

4. Diffusion — Jacobi iteration

Solving (I − ν Δt ∇²) u_new = u_old implicitly for u_new gives a sparse linear system. On a regular grid, the Laplacian discretises as:

// Jacobi shader (one iteration; ping-pong between two FBOs)
uniform sampler2D uX;        // current iterate x_k
uniform sampler2D uB;        // right-hand side b  (u_old)
uniform float      uAlpha;    // Δx² / (ν·Δt)
uniform float      uBeta;     // 1 / (4 + α)
uniform vec2      uTexelSize;

void main() {
  vec2 uv = gl_FragCoord.xy * uTexelSize;
  vec4 xL = texture(uX, uv - vec2(uTexelSize.x, 0));
  vec4 xR = texture(uX, uv + vec2(uTexelSize.x, 0));
  vec4 xB = texture(uX, uv - vec2(0, uTexelSize.y));
  vec4 xT = texture(uX, uv + vec2(0, uTexelSize.y));
  vec4 b  = texture(uB, uv);
  gl_FragColor = (xL + xR + xB + xT + uAlpha * b) * uBeta;
}

For low-viscosity fluids you can skip explicit diffusion entirely (ν = 0) — the numerical diffusion from advection alone provides sufficient smoothing and saves 20–40 Jacobi passes per frame.

5. Pressure projection — Helmholtz decomposition

After advection, the velocity field u̅ is generally divergent (∇ · u̅ ≠ 0) — incompatible with real fluid behaviour. The Helmholtz decomposition guarantees any vector field can be split into a divergence-free part and a gradient:

// Divergence shader
void main() {
  vec2 uv = gl_FragCoord.xy * uTexelSize;
  float vL = texture(uVelocity, uv - vec2(uTexelSize.x,0)).x;
  float vR = texture(uVelocity, uv + vec2(uTexelSize.x,0)).x;
  float vB = texture(uVelocity, uv - vec2(0,uTexelSize.y)).y;
  float vT = texture(uVelocity, uv + vec2(0,uTexelSize.y)).y;
  gl_FragColor.r = 0.5 * (vR - vL + vT - vB);  // central differences
}

// Gradient subtraction shader
void main() {
  vec2 uv = gl_FragCoord.xy * uTexelSize;
  float pL = texture(uPressure, uv - vec2(uTexelSize.x,0)).r;
  float pR = texture(uPressure, uv + vec2(uTexelSize.x,0)).r;
  float pB = texture(uPressure, uv - vec2(0,uTexelSize.y)).r;
  float pT = texture(uPressure, uv + vec2(0,uTexelSize.y)).r;
  vec2  vel = texture(uVelocity, uv).xy;
  vel -= 0.5 * vec2(pR - pL, pT - pB);
  gl_FragColor = vec4(vel, 0, 1);
}

The pressure Poisson solve uses the same Jacobi shader as diffusion — parameter α = −Δx², β = 1/4. 20–50 iterations are typically enough for smooth flow; a Gauss-Seidel or multigrid solver would converge faster.

6. WebGL implementation — ping-pong textures

A WebGL fluid solver maintains several floating-point render targets (FBO pairs for ping-pong). A minimal setup requires:

• velocity — 2-component (RG) float texture
• pressure — 1-component (R) float texture
• divergence — 1-component (R) float texture
• dye — 4-component (RGBA) float texture

// Minimal WebGL2 FBO helper
function createDoubleFBO(gl, w, h, format) {
  const make = () => {
    const tex = gl.createTexture();
    gl.bindTexture(gl.TEXTURE_2D, tex);
    gl.texImage2D(gl.TEXTURE_2D, 0, format.internal,
                   w, h, 0, format.format, gl.FLOAT, null);
    gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MIN_FILTER, gl.LINEAR);
    gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_WRAP_S, gl.CLAMP_TO_EDGE);
    gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_WRAP_T, gl.CLAMP_TO_EDGE);
    const fbo = gl.createFramebuffer();
    gl.bindFramebuffer(gl.FRAMEBUFFER, fbo);
    gl.framebufferTexture2D(gl.FRAMEBUFFER, gl.COLOR_ATTACHMENT0, gl.TEXTURE_2D, tex, 0);
    return { tex, fbo };
  };
  let read = make(), write = make();
  return {
    get read()  { return read;  },
    get write() { return write; },
    swap() { [read, write] = [write, read]; }
  };
}

// Per-frame update loop
function update(dt) {
  addForces(velocity, mouseForce);          // step 1
  for (let i = 0; i < 20; i++)
    jacobi(velocity, velocity, alpha, beta); // step 2: diffuse
  advect(velocity, velocity, dt);           // step 3: advect vel
  computeDivergence(velocity, divergence);  // step 4a
  pressure.read.clearToZero();               // step 4b: reset p
  for (let i = 0; i < 40; i++)
    jacobi(pressure, divergence, -dx*dx, 0.25); // 4c: pressure solve
  subtractGradient(velocity, pressure);     // step 4d: project
  advect(dye, velocity, dt);                // advect dye
  renderDye(dye);
}

7. Dye advection and visualisation

The velocity field itself is invisible. Two popular visualisation techniques:

• Dye advection: inject colour at the mouse position and advect it each frame using the velocity field. The dye accumulates and reveals vortices, turbulence and mixing patterns.
• Curl / vorticity colouring: compute ω = ∂v/∂x − ∂u/∂y (the 2D vorticity scalar) and map positive/negative values to two colours. Vortex cores show up as bright spots.

// Curl shader — vorticity confinement colour pass
void main() {
  vec2 uv = gl_FragCoord.xy * uTexelSize;
  float vL = texture(uVelocity, uv - vec2(uTexelSize.x, 0)).y;
  float vR = texture(uVelocity, uv + vec2(uTexelSize.x, 0)).y;
  float uB = texture(uVelocity, uv - vec2(0, uTexelSize.y)).x;
  float uT = texture(uVelocity, uv + vec2(0, uTexelSize.y)).x;
  float curl = 0.5 * ((vR - vL) - (uT - uB));  // ∂v/∂x − ∂u/∂y
  gl_FragColor = vec4(curl, 0.0, 0.0, 1.0);
}

// Vorticity confinement force — re-energises vortices to counter dissipation
void main() {
  vec2  uv  = gl_FragCoord.xy * uTexelSize;
  float wL  = abs(texture(uCurl, uv - vec2(uTexelSize.x, 0)).r);
  float wR  = abs(texture(uCurl, uv + vec2(uTexelSize.x, 0)).r);
  float wB  = abs(texture(uCurl, uv - vec2(0, uTexelSize.y)).r);
  float wT  = abs(texture(uCurl, uv + vec2(0, uTexelSize.y)).r);
  float w   = texture(uCurl, uv).r;
  vec2  eta = normalize(vec2(wT - wB, wR - wL) + 1e-5);
  vec2  force = eta * w * uConfinement * uTexelSize;
  vec2  vel = texture(uVelocity, uv).xy + force * uDt;
  gl_FragColor = vec4(vel, 0, 1);
}

8. Extensions and further reading

• 3D Navier-Stokes: extend each texture to a 3D texture (WebGL 2 TEXTURE_3D). Memory usage grows as N³; a 64³ grid is practical in real time.
• Better Poisson solvers: Jacobi needs many iterations. Gauss-Seidel (red-black ordering) converges 2× faster; multigrid methods converge in O(N) work independent of grid size.
• Vortex particle methods: represent vorticity as discrete point vortices rather than a grid — naturally adaptive and free of numerical diffusion.
• SPH fluid (particle-based): rather than an Eulerian grid, Smoothed Particle Hydrodynamics uses Lagrangian particles — see the SPH article for details.
• LBM comparison: the Lattice Boltzmann Method achieves similar results with a completely different mathematical foundation — mesoscopic kinetic theory rather than macroscopic continuum mechanics.