Home

29. Stam 1999 solver · 2026 illustration · Modern experiments

Dual tank

Demos are animating. Following your system setting.
Colors and knobs

These rewrite the live file. Copy HTML picks up whatever you set.

About

Split screen, two copies of the GDC Stam solver. Left is the usual dye tank, labeled bounded. Right runs the same step, then a display-only strain: dye is sampled through an inward spiral and an axial stretch, and the false-color speed scale goes as 1/(T*−t). Orange is the boosted speed, teal is slow. The code on the right is still stable. The legend is not. Stir either side. T* is the horizon slider.

Move or drag on either half to stir that tank. Right-hand color scale grows as t approaches T*.

Reference

J. Stam, "Stable Fluids", SIGGRAPH 1999, for the solver; the strained display on the right is a 2026 illustration

Browser APIs

  • Canvas 2D
  • ImageData
  • Pointer Events
  • requestAnimationFrame
  • matchMedia

If WebGL is missing, the demo draws a message on the canvas instead of a blank frame. By default a demo holds a still frame when your system asks for reduced motion; the Motion buttons override that, and a standalone copy animates with #motion=play on its URL.

Source

(function () {
  const PARAMS = {
    color: "#ff6a2a",
    accent: "#3ce0c8",
    horizon: 6
  };
  function hexToRgb(hex) {
    const n = parseInt(String(hex).replace("#", ""), 16);
    if (isNaN(n)) return [255, 106, 42];
    return [n >> 16 & 255, n >> 8 & 255, n & 255];
  }

  const canvas = document.getElementById("c");
  const fallback = document.getElementById("fallback");
  const ctx = canvas.getContext("2d");
  if (!ctx) {
    fallback.hidden = false;
    fallback.textContent = "Canvas 2D is not available in this browser.";
    return;
  }

  const N = 48, ITER = 4, dim = N + 2, size = dim * dim;
  function ix(i, j) { return i + dim * j; }

  function makeTank() {
    function field() { return new Float32Array(size); }
    const tnk = {
      u: field(), v: field(), u0: field(), v0: field(),
      dens: field(), dens0: field(), p: field(), div: field()
    };
    function setBnd(b, x) {
      for (let i = 1; i <= N; i++) {
        x[ix(0, i)] = b === 1 ? -x[ix(1, i)] : x[ix(1, i)];
        x[ix(N + 1, i)] = b === 1 ? -x[ix(N, i)] : x[ix(N, i)];
        x[ix(i, 0)] = b === 2 ? -x[ix(i, 1)] : x[ix(i, 1)];
        x[ix(i, N + 1)] = b === 2 ? -x[ix(i, N)] : x[ix(i, N)];
      }
      x[ix(0, 0)] = 0.5 * (x[ix(1, 0)] + x[ix(0, 1)]);
      x[ix(0, N + 1)] = 0.5 * (x[ix(1, N + 1)] + x[ix(0, N)]);
      x[ix(N + 1, 0)] = 0.5 * (x[ix(N, 0)] + x[ix(N + 1, 1)]);
      x[ix(N + 1, N + 1)] = 0.5 * (x[ix(N, N + 1)] + x[ix(N + 1, N)]);
    }
    function linSolve(b, x, x0, a, c) {
      const inv = 1 / c;
      for (let k = 0; k < ITER; k++) {
        for (let j = 1; j <= N; j++) {
          for (let i = 1; i <= N; i++) {
            x[ix(i, j)] = (x0[ix(i, j)] + a * (
              x[ix(i - 1, j)] + x[ix(i + 1, j)] + x[ix(i, j - 1)] + x[ix(i, j + 1)]
            )) * inv;
          }
        }
        setBnd(b, x);
      }
    }
    function diffuse(b, x, x0, diff, dt) {
      const a = dt * diff * N * N;
      linSolve(b, x, x0, a, 1 + 4 * a);
    }
    function advect(b, d, d0, uu, vv, dt) {
      const dt0 = dt * N;
      for (let j = 1; j <= N; j++) {
        for (let i = 1; i <= N; i++) {
          let x = i - dt0 * uu[ix(i, j)];
          let y = j - dt0 * vv[ix(i, j)];
          if (x < 0.5) x = 0.5;
          if (x > N + 0.5) x = N + 0.5;
          if (y < 0.5) y = 0.5;
          if (y > N + 0.5) y = N + 0.5;
          const i0 = x | 0, i1 = i0 + 1;
          const j0 = y | 0, j1 = j0 + 1;
          const s1 = x - i0, s0 = 1 - s1;
          const t1 = y - j0, t0 = 1 - t1;
          d[ix(i, j)] =
            s0 * (t0 * d0[ix(i0, j0)] + t1 * d0[ix(i0, j1)]) +
            s1 * (t0 * d0[ix(i1, j0)] + t1 * d0[ix(i1, j1)]);
        }
      }
      setBnd(b, d);
    }
    function project(uu, vv, pp, dv) {
      for (let j = 1; j <= N; j++) {
        for (let i = 1; i <= N; i++) {
          dv[ix(i, j)] = -0.5 * (uu[ix(i + 1, j)] - uu[ix(i - 1, j)] + vv[ix(i, j + 1)] - vv[ix(i, j - 1)]) / N;
          pp[ix(i, j)] = 0;
        }
      }
      setBnd(0, dv);
      setBnd(0, pp);
      linSolve(0, pp, dv, 1, 4);
      for (let j = 1; j <= N; j++) {
        for (let i = 1; i <= N; i++) {
          uu[ix(i, j)] -= 0.5 * N * (pp[ix(i + 1, j)] - pp[ix(i - 1, j)]);
          vv[ix(i, j)] -= 0.5 * N * (pp[ix(i, j + 1)] - pp[ix(i, j - 1)]);
        }
      }
      setBnd(1, uu);
      setBnd(2, vv);
    }
    tnk.step = function (dt) {
      const u = tnk.u, v = tnk.v, u0 = tnk.u0, v0 = tnk.v0;
      diffuse(1, u0, u, 0.00012, dt);
      diffuse(2, v0, v, 0.00012, dt);
      project(u0, v0, tnk.p, tnk.div);
      advect(1, u, u0, u0, v0, dt);
      advect(2, v, v0, u0, v0, dt);
      project(u, v, tnk.p, tnk.div);
      diffuse(0, tnk.dens0, tnk.dens, 0.00008, dt);
      advect(0, tnk.dens, tnk.dens0, u, v, dt);
      for (let i = 0; i < size; i++) tnk.dens[i] *= 0.994;
    };
    tnk.splat = function (gx, gy, dx, dy, amount) {
      const i = Math.max(1, Math.min(N, gx | 0));
      const j = Math.max(1, Math.min(N, gy | 0));
      for (let oj = -1; oj <= 1; oj++) {
        for (let oi = -1; oi <= 1; oi++) {
          const ii = i + oi, jj = j + oj;
          if (ii < 1 || jj < 1 || ii > N || jj > N) continue;
          const wgt = oi === 0 && oj === 0 ? 1 : 0.45;
          tnk.dens[ix(ii, jj)] += amount * wgt;
          tnk.u[ix(ii, jj)] += dx * wgt;
          tnk.v[ix(ii, jj)] += dy * wgt;
        }
      }
    };
    tnk.swirl = function (force, dye) {
      const c = N * 0.5 + 1;
      for (let j = 1; j <= N; j++) {
        for (let i = 1; i <= N; i++) {
          const dx = i - c, dy = j - c;
          const r2 = dx * dx + dy * dy;
          const g = Math.exp(-r2 * 0.012);
          tnk.u[ix(i, j)] += -dy * force * g;
          tnk.v[ix(i, j)] += dx * force * g;
          tnk.dens[ix(i, j)] += dye * g * (0.4 + 0.6 * ((i + j) & 1));
        }
      }
    };
    tnk.sample = function (arr, x, y) {
      if (x < 0.5) x = 0.5;
      if (x > N + 0.5) x = N + 0.5;
      if (y < 0.5) y = 0.5;
      if (y > N + 0.5) y = N + 0.5;
      const i0 = x | 0, i1 = i0 + 1;
      const j0 = y | 0, j1 = j0 + 1;
      const s1 = x - i0, s0 = 1 - s1;
      const t1 = y - j0, t0 = 1 - t1;
      return s0 * (t0 * arr[ix(i0, j0)] + t1 * arr[ix(i0, j1)]) +
             s1 * (t0 * arr[ix(i1, j0)] + t1 * arr[ix(i1, j1)]);
    };
    return tnk;
  }

  const left = makeTank();
  const right = makeTank();
  left.swirl(0.35, 8);
  right.swirl(0.35, 8);

  const pixelsL = new ImageData(N, N);
  const pixelsR = new ImageData(N, N);
  const offL = document.createElement("canvas");
  const offR = document.createElement("canvas");
  offL.width = offR.width = N;
  offL.height = offR.height = N;
  const ctxL = offL.getContext("2d");
  const ctxR = offR.getContext("2d");
  if (!ctxL || !ctxR) {
    fallback.hidden = false;
    fallback.textContent = "Canvas 2D is not available in this browser.";
    return;
  }

  const motionQuery = window.matchMedia("(prefers-reduced-motion: reduce)");
  let reduced = motionQuery.matches;
  const pointer = { x: 0, y: 0, px: 0, py: 0, down: false, side: -1 };
  let last = 0, tAbs = 0, raf = 0;
  let primed = false;

  function whichSide(e) {
    const rect = canvas.getBoundingClientRect();
    const nx = (e.clientX - rect.left) / rect.width;
    const ny = (e.clientY - rect.top) / rect.height;
    const side = nx < 0.5 ? 0 : 1;
    const localX = side === 0 ? nx / 0.5 : (nx - 0.5) / 0.5;
    return { side: side, x: localX * N + 1, y: ny * N + 1 };
  }
  canvas.addEventListener("pointerdown", function (e) {
    const g = whichSide(e);
    pointer.down = true; pointer.side = g.side;
    pointer.x = pointer.px = g.x; pointer.y = pointer.py = g.y;
    canvas.setPointerCapture(e.pointerId);
    wake();
  });
  canvas.addEventListener("pointerup", function () { pointer.down = false; wake(); });
  canvas.addEventListener("pointerleave", function () { pointer.side = -1; wake(); });
  canvas.addEventListener("pointermove", function (e) {
    const g = whichSide(e);
    pointer.px = pointer.x; pointer.py = pointer.y;
    pointer.x = g.x; pointer.y = g.y;
    if (pointer.down) pointer.side = g.side;
    else pointer.side = g.side;
    wake();
  });

  function strainAt(t) {
    const T = Math.max(2, PARAMS.horizon);
    const cycle = T + 0.9;
    const tm = t % cycle;
    if (tm >= T) return 1;
    return Math.min(22, 1 / Math.max(0.05, 1 - tm / T));
  }

  function paintLeft() {
    const data = pixelsL.data;
    const rgb = hexToRgb(PARAMS.accent);
    for (let j = 1; j <= N; j++) {
      for (let i = 1; i <= N; i++) {
        const d = left.dens[ix(i, j)];
        const spd = Math.hypot(left.u[ix(i, j)], left.v[ix(i, j)]);
        const glow = Math.min(1, d * 0.045 + spd * 0.08);
        const o = ((j - 1) * N + (i - 1)) * 4;
        data[o] = 12 + glow * rgb[0];
        data[o + 1] = 14 + glow * rgb[1];
        data[o + 2] = 16 + glow * rgb[2];
        data[o + 3] = 255;
      }
    }
    ctxL.putImageData(pixelsL, 0, 0);
  }

  function paintRight(strain) {
    const data = pixelsR.data;
    const hot = hexToRgb(PARAMS.color);
    const cool = hexToRgb(PARAMS.accent);
    const c = N * 0.5 + 1;
    for (let j = 1; j <= N; j++) {
      for (let i = 1; i <= N; i++) {
        const dx = i - c, dy = j - c;
        const r = Math.hypot(dx, dy) + 0.001;
        const th = Math.atan2(dy, dx);
        const r0 = r * (1 + (strain - 1) * 0.62);
        const th0 = th - Math.log(r + 1) * (strain - 1) * 0.48;
        const x0 = c + r0 * Math.cos(th0);
        const y0 = c + (r0 * Math.sin(th0)) / (1 + (strain - 1) * 1.15);
        const d = right.sample(right.dens, x0, y0);
        const uu = right.sample(right.u, x0, y0);
        const vv = right.sample(right.v, x0, y0);
        const spd = Math.hypot(uu, vv) * strain;
        const heat = Math.min(1, (strain - 1) * 0.08 + spd * 0.07);
        const k = Math.min(1, d * 0.05 + spd * 0.04);
        const o = ((j - 1) * N + (i - 1)) * 4;
        data[o] = 10 + k * (cool[0] + (hot[0] - cool[0]) * heat);
        data[o + 1] = 12 + k * (cool[1] + (hot[1] - cool[1]) * heat);
        data[o + 2] = 14 + k * (cool[2] + (hot[2] - cool[2]) * heat);
        data[o + 3] = 255;
      }
    }
    ctxR.putImageData(pixelsR, 0, 0);
  }

  function labels(w, h, strain, t) {
    const T = Math.max(2, PARAMS.horizon);
    const tm = t % (T + 0.9);
    ctx.save();
    ctx.font = Math.max(10, (h * 0.028) | 0) + "px ui-monospace, 'IBM Plex Mono', monospace";
    ctx.textBaseline = "top";
    ctx.fillStyle = "rgba(11,12,14,0.55)";
    ctx.fillRect(6, 6, w * 0.28, h * 0.13);
    ctx.fillRect(w * 0.5 + 6, 6, w * 0.42, h * 0.16);
    ctx.fillStyle = "#c8d2dc";
    ctx.fillText("BOUNDED", 12, 10);
    ctx.fillText("same Stam 48×48", 12, 10 + h * 0.045);
    ctx.fillStyle = "#ff6a2a";
    ctx.fillText("DISPLAY STRAIN", w * 0.5 + 12, 10);
    ctx.fillStyle = "#c8d2dc";
    ctx.fillText("t → T*   scale 1/(T*−t)", w * 0.5 + 12, 10 + h * 0.045);
    ctx.fillText("×" + strain.toFixed(1) + (tm >= T ? "  reset" : ""), w * 0.5 + 12, 10 + h * 0.09);
    ctx.restore();
  }

  function resize() {
    const dpr = Math.min(window.devicePixelRatio || 1, 2);
    const w = Math.max(1, Math.floor(canvas.clientWidth * dpr));
    const h = Math.max(1, Math.floor(canvas.clientHeight * dpr));
    if (canvas.width !== w || canvas.height !== h) { canvas.width = w; canvas.height = h; }
  }

  function frame(now) {
    const dt = Math.min(0.033, last ? (now - last) * 0.001 : 0.016);
    last = now;
    if (!reduced) tAbs += dt;
    const strain = strainAt(tAbs);

    if (!reduced) {
      left.swirl(0.018, 0.35);
      right.swirl(0.018, 0.35);
    }
    if (pointer.side === 0 || (pointer.down && pointer.side === 0)) {
      const dx = (pointer.x - pointer.px) * 12;
      const dy = (pointer.y - pointer.py) * 12;
      left.splat(pointer.x, pointer.y, dx, dy, pointer.down ? 90 : 22);
    }
    if (pointer.side === 1 || (pointer.down && pointer.side === 1)) {
      const dx = (pointer.x - pointer.px) * 12;
      const dy = (pointer.y - pointer.py) * 12;
      right.splat(pointer.x, pointer.y, dx, dy, pointer.down ? 90 : 22);
    }
    pointer.px = pointer.x; pointer.py = pointer.y;

    if (!reduced) {
      left.step(dt);
      right.step(dt);
    } else if (!primed) {
      left.step(0.016);
      right.step(0.016);
    }
    primed = true;

    paintLeft();
    paintRight(strain);
    resize();
    const w = canvas.width, h = canvas.height;
    const gap = Math.max(2, (w * 0.006) | 0);
    const hw = (w - gap) >> 1;
    ctx.imageSmoothingEnabled = false;
    ctx.fillStyle = "#0b0c0e";
    ctx.fillRect(0, 0, w, h);
    ctx.drawImage(offL, 0, 0, hw, h);
    ctx.drawImage(offR, hw + gap, 0, w - hw - gap, h);
    ctx.fillStyle = "#4a5560";
    ctx.fillRect(hw, 0, gap, h);
    labels(w, h, strain, tAbs);
    if (!reduced) raf = requestAnimationFrame(frame);
  }

  // Motion. The system's reduced-motion preference decides by default. A
  // standalone copy animates anyway with "#motion=play" in its URL (or holds a
  // frame with "#motion=still"); an embedding page can post {motion: "play"}.
  let motionMode = "auto";
  function applyMotion() {
    const next = motionMode === "still" || (motionMode === "auto" && motionQuery.matches);
    if (next === reduced) return;
    reduced = next;
    last = 0;
    if (reduced) { cancelAnimationFrame(raf); raf = 0; frame(performance.now()); }
    else if (!raf) raf = requestAnimationFrame(frame);
  }
  function setMotion(mode) {
    if (mode === "play" || mode === "still" || mode === "auto") { motionMode = mode; applyMotion(); }
  }
  function readMotionHash() {
    const m = /motion=(play|still|auto)/.exec(location.hash);
    let mode = m ? m[1] : "auto";
    // Embedded in the gallery, the page sets window.__demoMotion before any frame loads.
    try { if (!m && window.parent !== window && typeof window.parent.__demoMotion === "string") mode = window.parent.__demoMotion; } catch (e) { /* cross-origin parent */ }
    setMotion(mode);
  }
  function wake() { if (reduced) frame(performance.now()); }
  if (motionQuery.addEventListener) motionQuery.addEventListener("change", applyMotion);
  else if (motionQuery.addListener) motionQuery.addListener(applyMotion);
  window.addEventListener("hashchange", readMotionHash);
  window.addEventListener("message", function (e) { if (e.data && typeof e.data.motion === "string") setMotion(e.data.motion); });
  window.addEventListener("resize", wake);

  if (reduced) frame(performance.now());
  else raf = requestAnimationFrame(frame);
  readMotionHash();
})();