Vorticity Stir: Navier-Stokes

// Vorticity-streamfunction Navier-Stokes on a periodic grid.
// State: w (vorticity), psi (streamfunction), velocities derived as u=dpsi/dy, v=-dpsi/dx.
// Equations: dw/dt + (u.grad)w = nu * lap(w);  lap(psi) = -w.
let GW, GH, W, H, scale;
let w, w2, psi, u, v;
let img, pix;
let buf, bufCtx;
let cw, ch;
let lastMX = -1, lastMY = -1, lastDown = false;
let quietFrames = 0;

function alloc(width, height) {
  W = width; H = height; cw = width; ch = height;
  scale = Math.max(4, Math.floor(width / 100));
  GW = Math.max(80, Math.floor(width / scale));
  GH = Math.max(60, Math.floor(height / scale));
  const N = GW * GH;
  w = new Float32Array(N);
  w2 = new Float32Array(N);
  psi = new Float32Array(N);
  u = new Float32Array(N);
  v = new Float32Array(N);
}

function buildImg(ctx) {
  img = ctx.createImageData(GW, GH);
  pix = img.data;
  for (let i = 3; i < pix.length; i += 4) pix[i] = 255;
  buf = new OffscreenCanvas(GW, GH);
  bufCtx = buf.getContext('2d');
}

function seed() {
  // a couple of opposing vortices so the first frame is alive
  const cx = GW * 0.5, cy = GH * 0.5;
  for (let y = 0; y < GH; y++) {
    for (let x = 0; x < GW; x++) {
      const dx1 = (x - cx + 14), dy1 = (y - cy);
      const dx2 = (x - cx - 14), dy2 = (y - cy);
      const a = Math.exp(-(dx1 * dx1 + dy1 * dy1) / 60);
      const b = Math.exp(-(dx2 * dx2 + dy2 * dy2) / 60);
      w[y * GW + x] = 2.2 * a - 2.2 * b;
    }
  }
}

function init({ ctx, width, height }) {
  alloc(width, height);
  buildImg(ctx);
  seed();
}

// Inject vorticity along a drag segment: curl = dvx/dy - dvy/dx, approximated
// by depositing a dipole pair perpendicular to drag direction.
function inject(px, py, dx, dy) {
  const gx = px / scale, gy = py / scale;
  const len = Math.hypot(dx, dy);
  if (len < 0.5) return;
  const nx = -dy / len, ny = dx / len; // perpendicular
  const amp = Math.min(8, len * 0.6);
  const r = 5;
  for (let oy = -r; oy <= r; oy++) {
    for (let ox = -r; ox <= r; ox++) {
      const d2 = ox * ox + oy * oy;
      if (d2 > r * r) continue;
      const k = Math.exp(-d2 * 0.18);
      // dot with perpendicular gives sign of curl on each side of the drag axis
      const side = (ox * nx + oy * ny);
      const xi = ((gx + ox) | 0);
      const yi = ((gy + oy) | 0);
      if (xi < 1 || yi < 1 || xi >= GW - 1 || yi >= GH - 1) continue;
      w[yi * GW + xi] += amp * k * side * 0.18;
    }
  }
}

// Jacobi/Gauss-Seidel solve of lap(psi) = -w with zero boundaries (no-slip-ish).
function solvePsi(iters) {
  for (let k = 0; k < iters; k++) {
    for (let y = 1; y < GH - 1; y++) {
      const row = y * GW;
      for (let x = 1; x < GW - 1; x++) {
        const i = row + x;
        psi[i] = 0.25 * (psi[i - 1] + psi[i + 1] + psi[i - GW] + psi[i + GW] + w[i]);
      }
    }
  }
}

function velocities() {
  for (let y = 1; y < GH - 1; y++) {
    const row = y * GW;
    for (let x = 1; x < GW - 1; x++) {
      const i = row + x;
      u[i] = 0.5 * (psi[i + GW] - psi[i - GW]);
      v[i] = -0.5 * (psi[i + 1] - psi[i - 1]);
    }
  }
}

// Semi-Lagrangian advection + explicit diffusion of w.
function advectDiffuse(dt) {
  const nu = 0.06;
  for (let y = 1; y < GH - 1; y++) {
    const row = y * GW;
    for (let x = 1; x < GW - 1; x++) {
      const i = row + x;
      let sx = x - u[i] * dt;
      let sy = y - v[i] * dt;
      if (sx < 0.5) sx = 0.5; else if (sx > GW - 1.5) sx = GW - 1.5;
      if (sy < 0.5) sy = 0.5; else if (sy > GH - 1.5) sy = GH - 1.5;
      const x0 = sx | 0, y0 = sy | 0;
      const fx = sx - x0, fy = sy - y0;
      const j = y0 * GW + x0;
      const a = w[j], b = w[j + 1], c = w[j + GW], d = w[j + GW + 1];
      const advected = a * (1 - fx) * (1 - fy) + b * fx * (1 - fy) + c * (1 - fx) * fy + d * fx * fy;
      const lap = w[i - 1] + w[i + 1] + w[i - GW] + w[i + GW] - 4 * w[i];
      w2[i] = advected + nu * lap * dt * 0.25;
      // mild ambient decay so the field doesn't grow without bound
      w2[i] *= 0.9995;
    }
  }
  const t = w; w = w2; w2 = t;
}

function tick({ ctx, dt, frame, width, height, input }) {
  if (width !== cw || height !== ch) { alloc(width, height); buildImg(ctx); seed(); }

  // mouse drag injection
  if (input.mouseDown) {
    if (lastDown && lastMX >= 0) inject(input.mouseX, input.mouseY, input.mouseX - lastMX, input.mouseY - lastMY);
    lastMX = input.mouseX; lastMY = input.mouseY; lastDown = true;
  } else {
    lastDown = false; lastMX = -1; lastMY = -1;
  }

  const sub = 2;
  const stepDt = 1.0;
  for (let s = 0; s < sub; s++) {
    solvePsi(12);
    velocities();
    advectDiffuse(stepDt);
  }

  // diverging colormap: red positive, blue negative, near-white at zero
  let maxAbs = 0.6;
  for (let i = 0; i < w.length; i += 31) { const a = Math.abs(w[i]); if (a > maxAbs) maxAbs = a; }
  const inv = 1 / maxAbs;
  for (let i = 0, p = 0; i < w.length; i++, p += 4) {
    let t = w[i] * inv;
    if (t > 1) t = 1; else if (t < -1) t = -1;
    const m = Math.abs(t);
    // base latte color (stirred coffee crema)
    const baseR = 238, baseG = 222, baseB = 200;
    let r, g, b;
    if (t >= 0) { r = baseR - (baseR - 178) * m; g = baseG - (baseG - 52) * m; b = baseB - (baseB - 48) * m; }
    else { r = baseR - (baseR - 46) * m; g = baseG - (baseG - 88) * m; b = baseB - (baseB - 168) * m; }
    pix[p] = r; pix[p + 1] = g; pix[p + 2] = b;
  }

  // quiet-field detector: if vorticity stays low for ~10s @60fps, drop a swirl
  if (maxAbs < 0.65) quietFrames++; else quietFrames = 0;
  if (quietFrames > 600 && !input.mouseDown) {
    const cx = GW * (0.3 + 0.4 * Math.random());
    const cy = GH * (0.3 + 0.4 * Math.random());
    const sgn = Math.random() < 0.5 ? -1 : 1;
    const r = 7;
    for (let oy = -r; oy <= r; oy++) {
      for (let ox = -r; ox <= r; ox++) {
        const d2 = ox * ox + oy * oy;
        if (d2 > r * r) continue;
        const xi = (cx + ox) | 0, yi = (cy + oy) | 0;
        if (xi < 1 || yi < 1 || xi >= GW - 1 || yi >= GH - 1) continue;
        w[yi * GW + xi] += sgn * 1.6 * Math.exp(-d2 * 0.12);
      }
    }
    quietFrames = 0;
  }

  bufCtx.putImageData(img, 0, 0);
  ctx.imageSmoothingEnabled = true;
  ctx.drawImage(buf, 0, 0, GW, GH, 0, 0, width, height);

  ctx.fillStyle = "rgba(40, 28, 22, 0.78)";
  ctx.font = "12px system-ui, sans-serif";
  ctx.fillText("drag to stir — red +omega, blue -omega", 12, height - 12);
}