GPU Fluid Simulation with WebGL

The incompressible Navier-Stokes momentum equation is:

∂u/∂t = −(u·∇)u − (1/ρ)∇p + ν∇²u + f

where u = velocity field, p = pressure, ν = viscosity, f = external forces. With the incompressibility constraint ∇·u = 0.

Each frame the GPU executes these passes in order:

A shader cannot read and write the same texture simultaneously. Ping-pong uses two FBOs — one is read while the other is written, then they swap:

import * as THREE from 'https://cdn.jsdelivr.net/npm/three@0.168/build/three.module.js';

function createFBO(width, height) {
  return new THREE.WebGLRenderTarget(width, height, {
    type: THREE.FloatType,           // 32-bit float per channel
    format: THREE.RGBAFormat,
    minFilter: THREE.LinearFilter,
    magFilter: THREE.LinearFilter,
    wrapS: THREE.ClampToEdgeWrapping,
    wrapT: THREE.ClampToEdgeWrapping,
  });
}

class DoubleFBO {
  constructor(w, h) {
    this.read  = createFBO(w, h);
    this.write = createFBO(w, h);
  }
  swap() { [this.read, this.write] = [this.write, this.read]; }
}

const SIM_RES = 256;
const velocity = new DoubleFBO(SIM_RES, SIM_RES);
const pressure = new DoubleFBO(SIM_RES, SIM_RES);
const dye      = new DoubleFBO(SIM_RES, SIM_RES);

Semi-Lagrangian advection: "where did this particle come from one timestep ago?" — follow the velocity field backward and sample there:

// advection fragment shader
const advectFrag = /* glsl */`
  uniform sampler2D uVelocity; // velocity field to advect WITH
  uniform sampler2D uSource;   // field to advect (velocity or dye)
  uniform float uDt;
  uniform vec2  uTexelSize;    // 1/SIM_RES for each axis

  varying vec2 vUv;

  void main() {
    // Read velocity at current cell
    vec2 velocity = texture2D(uVelocity, vUv).xy;

    // Trace backward in time: prev_pos = current - velocity * dt
    vec2 prevUv = vUv - velocity * uDt * uTexelSize;

    // Sample the source field at the previous position
    gl_FragColor = texture2D(uSource, prevUv);
  }
`;

To run a shader pass: render a full-screen quad with the shader, writing to the target FBO:

function runPass(material, targetFBO) {
  quad.material = material;
  renderer.setRenderTarget(targetFBO);
  renderer.render(quadScene, quadCamera);
  renderer.setRenderTarget(null);
}

advectMat.uniforms.uVelocity.value = velocity.read.texture;
advectMat.uniforms.uSource.value   = velocity.read.texture;
advectMat.uniforms.uDt.value       = dt;
runPass(advectMat, velocity.write);
velocity.swap();

Compute divergence (how much the velocity field expands), then iteratively solve for pressure using Jacobi iterations:

// Divergence shader — reads velocity, writes scalar div into .r channel
const divergenceFrag = /* glsl */`
  uniform sampler2D uVelocity;
  uniform vec2 uTexelSize;
  varying vec2 vUv;
  void main() {
    vec2 ts = uTexelSize;
    float L = texture2D(uVelocity, vUv - vec2(ts.x, 0)).x;
    float R = texture2D(uVelocity, vUv + vec2(ts.x, 0)).x;
    float B = texture2D(uVelocity, vUv - vec2(0, ts.y)).y;
    float T = texture2D(uVelocity, vUv + vec2(0, ts.y)).y;
    float div = 0.5 * (R - L + T - B);
    gl_FragColor = vec4(div, 0, 0, 1);
  }
`;

// Pressure Jacobi iteration shader — run 30–50 times
const pressureFrag = /* glsl */`
  uniform sampler2D uPressure;
  uniform sampler2D uDivergence;
  uniform vec2 uTexelSize;
  varying vec2 vUv;
  void main() {
    vec2 ts = uTexelSize;
    float L = texture2D(uPressure, vUv - vec2(ts.x, 0)).r;
    float R = texture2D(uPressure, vUv + vec2(ts.x, 0)).r;
    float B = texture2D(uPressure, vUv - vec2(0, ts.y)).r;
    float T = texture2D(uPressure, vUv + vec2(0, ts.y)).r;
    float div = texture2D(uDivergence, vUv).r;
    // Jacobi: p_new = (L+R+B+T - div) / 4
    float pNew = (L + R + B + T - div) * 0.25;
    gl_FragColor = vec4(pNew, 0, 0, 1);
  }
`;

// Run 40 Jacobi iterations
for (let i = 0; i < 40; i++) {
  pressureMat.uniforms.uPressure.value = pressure.read.texture;
  runPass(pressureMat, pressure.write);
  pressure.swap();
}

Subtract the pressure gradient from velocity to make the field divergence-free:

const projectFrag = /* glsl */`
  uniform sampler2D uVelocity;
  uniform sampler2D uPressure;
  uniform vec2 uTexelSize;
  varying vec2 vUv;
  void main() {
    vec2 ts = uTexelSize;
    float pL = texture2D(uPressure, vUv - vec2(ts.x, 0)).r;
    float pR = texture2D(uPressure, vUv + vec2(ts.x, 0)).r;
    float pB = texture2D(uPressure, vUv - vec2(0, ts.y)).r;
    float pT = texture2D(uPressure, vUv + vec2(0, ts.y)).r;
    vec2 gradP = vec2(pR - pL, pT - pB) * 0.5;
    vec2 vel = texture2D(uVelocity, vUv).xy;
    gl_FragColor = vec4(vel - gradP, 0, 1);
  }
`;

Splat forces and dye onto the simulation when the user moves the mouse:

const splatFrag = /* glsl */`
  uniform sampler2D uTarget;
  uniform vec2  uPoint;    // normalised mouse position [0,1]
  uniform vec3  uColor;    // force vector (xy) or dye colour (xyz)
  uniform float uRadius;
  varying vec2 vUv;
  void main() {
    float d = distance(vUv, uPoint);
    float splat = exp(-d * d / uRadius);
    vec4 current = texture2D(uTarget, vUv);
    gl_FragColor = current + vec4(uColor * splat, 0);
  }
`;

// Call on mousemove to inject velocity + dye
function onMouseMove(e) {
  const x = e.clientX / innerWidth;
  const y = 1.0 - e.clientY / innerHeight; // flip y
  const dx = (e.clientX - prevMouseX) / innerWidth  * 15.0; // scale to sim speed
  const dy = (e.clientY - prevMouseY) / innerHeight * 15.0;

  splatMat.uniforms.uPoint.value.set(x, y);
  splatMat.uniforms.uRadius.value = 0.0015;

  // Velocity splat
  splatMat.uniforms.uColor.value.set(dx, -dy, 0);
  splatMat.uniforms.uTarget.value = velocity.read.texture;
  runPass(splatMat, velocity.write); velocity.swap();

  // Dye splat
  splatMat.uniforms.uColor.value.set(0.3, 0.7, 1.0); // cyan dye
  splatMat.uniforms.uTarget.value = dye.read.texture;
  runPass(splatMat, dye.write); dye.swap();
}

Advect the dye texture with the divergence-free velocity field, then display it:

// Each frame:
// 1. Advect velocity
advectMat.uniforms.uSource.value = velocity.read.texture;
runPass(advectMat, velocity.write); velocity.swap();

// 2. Add forces (if mouse is down)

// 3. Compute divergence
runPass(divergenceMat, divFBO);

// 4. Solve pressure (40 Jacobi iterations)
for (let k = 0; k < 40; k++) { /* ... */ }

// 5. Project
projectMat.uniforms.uVelocity.value = velocity.read.texture;
runPass(projectMat, velocity.write); velocity.swap();

// 6. Advect dye
advectMat.uniforms.uVelocity.value = velocity.read.texture;
advectMat.uniforms.uSource.value   = dye.read.texture;
runPass(advectMat, dye.write); dye.swap();

// 7. Display dye to screen
displayMat.uniforms.uTexture.value = dye.read.texture;
renderer.setRenderTarget(null);
renderer.render(quadScene, quadCamera);