Fluid Simulation: How to Create Mesmerizing Liquid Art With Code

Pour cream into coffee and watch it spiral. Blow smoke across a beam of light and see vortices form. Run your hand through water and watch turbulence cascade into smaller and smaller eddies. Fluids are everywhere, and they produce some of the most hypnotic, organic motion in nature. They also produce some of the most stunning generative art you can create with code.

This guide covers eight working fluid simulations you can build in your browser. Every example uses vanilla JavaScript and HTML Canvas. No physics libraries, no WebGL (though we'll mention where it helps) — just the math of fluid dynamics, simplified for real-time art. You'll go from simple diffusion to a full Navier-Stokes grid solver, and everything in between.

The physics: what makes fluids flow

Fluid simulation boils down to three operations, applied repeatedly to a grid of velocity and density values:

1. Advection — the fluid carries things (including itself) along its velocity field
2. Diffusion — properties spread out to neighboring cells (viscosity, heat dissipation)
3. Projection — enforce incompressibility (mass conservation) by removing divergence from the velocity field

Jos Stam's 1999 paper "Stable Fluids" showed how to do all three in a way that never explodes numerically, no matter how large the time step. That's the method most real-time fluid art uses — including everything in this guide. The core insight: instead of pushing fluid forward (which is unstable), trace backward to find where each cell's value came from.

1. Diffusion — the simplest fluid behavior

Before simulating full fluid flow, let's start with diffusion alone. Drop some "ink" on a grid and watch it spread. This is the foundation of smoke and dye effects.

const canvas = document.createElement('canvas');
canvas.width = 200; canvas.height = 200;
document.body.appendChild(canvas);
const ctx = canvas.getContext('2d');

const N = 200;
let density = new Float32Array(N * N);
let densityPrev = new Float32Array(N * N);
const diffusion = 0.0001;
const dt = 1.0;

function diffuse(x, x0, diff, dt) {
  const a = dt * diff * N * N;
  for (let k = 0; k < 20; k++) {
    for (let i = 1; i < N - 1; i++) {
      for (let j = 1; j < N - 1; j++) {
        const idx = i + j * N;
        x[idx] = (x0[idx] + a * (
          x[idx - 1] + x[idx + 1] +
          x[idx - N] + x[idx + N]
        )) / (1 + 4 * a);
      }
    }
  }
}

canvas.addEventListener('click', (e) => {
  const rect = canvas.getBoundingClientRect();
  const cx = Math.floor((e.clientX - rect.left));
  const cy = Math.floor((e.clientY - rect.top));
  for (let di = -3; di <= 3; di++)
    for (let dj = -3; dj <= 3; dj++) {
      const idx = (cx + di) + (cy + dj) * N;
      if (idx >= 0 && idx < N * N) density[idx] = 1.0;
    }
});

function render() {
  densityPrev.set(density);
  diffuse(density, densityPrev, diffusion, dt);
  const img = ctx.createImageData(N, N);
  for (let i = 0; i < N * N; i++) {
    const v = Math.min(255, density[i] * 255);
    img.data[i * 4] = v * 0.3;
    img.data[i * 4 + 1] = v * 0.6;
    img.data[i * 4 + 2] = v;
    img.data[i * 4 + 3] = 255;
  }
  ctx.putImageData(img, 0, 0);
  requestAnimationFrame(render);
}
render();

Click on the canvas to drop ink. The Gauss-Seidel solver (the nested loop with 20 iterations) is the workhorse — it relaxes the diffusion equation iteratively. More iterations = more accurate spreading, but 20 is plenty for visual art.

2. Advection — carrying values along a velocity field

Diffusion spreads things out. Advection moves things along. This is what makes fluid simulation look like fluid — values swirl and streak instead of just blurring. The trick from Stam's method: for each cell, trace backward through the velocity field to find where the value came from, then interpolate.

const canvas = document.createElement('canvas');
canvas.width = 300; canvas.height = 300;
document.body.appendChild(canvas);
const ctx = canvas.getContext('2d');

const N = 150;
let density = new Float32Array(N * N);
let vx = new Float32Array(N * N);
let vy = new Float32Array(N * N);

// Create a swirling velocity field
for (let i = 0; i < N; i++) {
  for (let j = 0; j < N; j++) {
    const x = (i / N - 0.5) * 2;
    const y = (j / N - 0.5) * 2;
    const r = Math.sqrt(x * x + y * y);
    const strength = Math.exp(-r * 2) * 3;
    vx[i + j * N] = -y * strength;
    vy[i + j * N] = x * strength;
  }
}

function advect(d, d0, u, v, dt) {
  const dt0 = dt * N;
  for (let i = 1; i < N - 1; i++) {
    for (let j = 1; j < N - 1; j++) {
      const idx = i + j * N;
      let x = i - dt0 * u[idx];
      let y = j - dt0 * v[idx];
      x = Math.max(0.5, Math.min(N - 1.5, x));
      y = Math.max(0.5, Math.min(N - 1.5, y));
      const i0 = Math.floor(x), i1 = i0 + 1;
      const j0 = Math.floor(y), j1 = j0 + 1;
      const s1 = x - i0, s0 = 1 - s1;
      const t1 = y - j0, t0 = 1 - t1;
      d[idx] = s0 * (t0 * d0[i0 + j0 * N] + t1 * d0[i0 + j1 * N])
             + s1 * (t0 * d0[i1 + j0 * N] + t1 * d0[i1 + j1 * N]);
    }
  }
}

// Drop initial ink ring
for (let i = 0; i < N; i++) {
  for (let j = 0; j < N; j++) {
    const x = (i / N - 0.5) * 2;
    const y = (j / N - 0.5) * 2;
    const r = Math.sqrt(x * x + y * y);
    if (r > 0.3 && r < 0.5) density[i + j * N] = 1.0;
  }
}

function render() {
  const prev = new Float32Array(density);
  advect(density, prev, vx, vy, 0.02);
  const scale = canvas.width / N;
  ctx.fillStyle = '#0a0a0a';
  ctx.fillRect(0, 0, canvas.width, canvas.height);
  const img = ctx.createImageData(N, N);
  for (let i = 0; i < N * N; i++) {
    const v = Math.min(1, density[i]);
    const r = v * 60, g = v * 180, b = v * 255;
    img.data[i * 4] = r;
    img.data[i * 4 + 1] = g;
    img.data[i * 4 + 2] = b;
    img.data[i * 4 + 3] = 255;
  }
  ctx.putImageData(img, 0, 0);
  requestAnimationFrame(render);
}
render();

The ring of ink spirals inward along the vortex. The bilinear interpolation (the s0, s1, t0, t1 calculation) is essential — without it you'd get blocky artifacts from sampling grid cells. This backward-trace advection is unconditionally stable, meaning you can use any time step without the simulation blowing up.

3. Full Navier-Stokes grid solver — real fluid dynamics

Now we combine diffusion, advection, and pressure projection into a complete 2D fluid solver. This is Jos Stam's "Stable Fluids" algorithm — the same method used in most real-time fluid art, from screensavers to music visualizers.

const canvas = document.createElement('canvas');
canvas.width = 400; canvas.height = 400;
document.body.appendChild(canvas);
const ctx = canvas.getContext('2d');

const N = 128;
const size = (N + 2) * (N + 2);
let u = new Float32Array(size);
let v = new Float32Array(size);
let u0 = new Float32Array(size);
let v0 = new Float32Array(size);
let dens = new Float32Array(size);
let dens0 = new Float32Array(size);
const visc = 0.0;
const diff = 0.00001;
const dt = 0.1;

function IX(i, j) { return i + (N + 2) * j; }

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 cRecip = 1.0 / c;
  for (let k = 0; k < 4; 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)])) * cRecip;
    setBnd(b, x);
  }
}

function diffuseStep(b, x, x0, diff, dt) {
  const a = dt * diff * N * N;
  linSolve(b, x, x0, a, 1 + 4 * a);
}

function advectStep(b, d, d0, u, v, dt) {
  const dt0 = dt * N;
  for (let j = 1; j <= N; j++)
    for (let i = 1; i <= N; i++) {
      let x = i - dt0 * u[IX(i,j)];
      let y = j - dt0 * v[IX(i,j)];
      x = Math.max(0.5, Math.min(N + 0.5, x));
      y = Math.max(0.5, Math.min(N + 0.5, y));
      const i0 = Math.floor(x), i1 = i0 + 1;
      const j0 = Math.floor(y), 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(u, v, p, div) {
  for (let j = 1; j <= N; j++)
    for (let i = 1; i <= N; i++) {
      div[IX(i,j)] = -0.5 * (u[IX(i+1,j)] - u[IX(i-1,j)] + v[IX(i,j+1)] - v[IX(i,j-1)]) / N;
      p[IX(i,j)] = 0;
    }
  setBnd(0, div); setBnd(0, p);
  linSolve(0, p, div, 1, 4);
  for (let j = 1; j <= N; j++)
    for (let i = 1; i <= N; i++) {
      u[IX(i,j)] -= 0.5 * N * (p[IX(i+1,j)] - p[IX(i-1,j)]);
      v[IX(i,j)] -= 0.5 * N * (p[IX(i,j+1)] - p[IX(i,j-1)]);
    }
  setBnd(1, u); setBnd(2, v);
}

let mouseX = 0, mouseY = 0, pmouseX = 0, pmouseY = 0, mouseDown = false;
canvas.addEventListener('mousemove', (e) => {
  const rect = canvas.getBoundingClientRect();
  pmouseX = mouseX; pmouseY = mouseY;
  mouseX = (e.clientX - rect.left) / canvas.width * N;
  mouseY = (e.clientY - rect.top) / canvas.height * N;
});
canvas.addEventListener('mousedown', () => mouseDown = true);
canvas.addEventListener('mouseup', () => mouseDown = false);

function step() {
  if (mouseDown) {
    const i = Math.floor(mouseX), j = Math.floor(mouseY);
    if (i > 0 && i <= N && j > 0 && j <= N) {
      dens[IX(i,j)] += 8;
      u[IX(i,j)] += (mouseX - pmouseX) * 20;
      v[IX(i,j)] += (mouseY - pmouseY) * 20;
    }
  }
  u0.set(u); v0.set(v);
  diffuseStep(1, u, u0, visc, dt);
  diffuseStep(2, v, v0, visc, dt);
  project(u, v, u0, v0);
  u0.set(u); v0.set(v);
  advectStep(1, u, u0, u0, v0, dt);
  advectStep(2, v, v0, u0, v0, dt);
  project(u, v, u0, v0);
  dens0.set(dens);
  diffuseStep(0, dens, dens0, diff, dt);
  dens0.set(dens);
  advectStep(0, dens, dens0, u, v, dt);
}

function render() {
  step();
  const img = ctx.createImageData(canvas.width, canvas.height);
  const scale = canvas.width / (N + 2);
  for (let y = 0; y < canvas.height; y++)
    for (let x = 0; x < canvas.width; x++) {
      const i = Math.floor(x / scale), j = Math.floor(y / scale);
      const d = Math.min(1, dens[IX(i,j)]);
      const idx = (y * canvas.width + x) * 4;
      img.data[idx] = d * 40;
      img.data[idx+1] = d * 120 + d * d * 80;
      img.data[idx+2] = d * 200 + d * d * 55;
      img.data[idx+3] = 255;
    }
  ctx.putImageData(img, 0, 0);
  requestAnimationFrame(render);
}
render();

Click and drag to inject dye and push the fluid. This is the full algorithm: diffuse velocities, project to remove divergence, advect velocities, project again, then diffuse and advect the density field. The boundary conditions (setBnd) reflect velocity at walls so the fluid bounces instead of leaking through edges.

4. Smoke simulation with buoyancy

Real smoke rises because it's hot. We can add temperature as a second field (alongside density) and apply an upward buoyancy force proportional to temperature. This produces realistic rising plumes, mushroom clouds, and convection rolls.

const canvas = document.createElement('canvas');
canvas.width = 400; canvas.height = 500;
document.body.appendChild(canvas);
const ctx = canvas.getContext('2d');

const N = 100;
const size = (N + 2) * (N + 2);
const IX = (i, j) => i + (N + 2) * j;

let u = new Float32Array(size), v = new Float32Array(size);
let u0 = new Float32Array(size), v0 = new Float32Array(size);
let dens = new Float32Array(size), dens0 = new Float32Array(size);
let temp = new Float32Array(size), temp0 = new Float32Array(size);
const buoyancy = 0.15;
const ambientTemp = 0.0;

function addSource() {
  // Smoke source at bottom center
  for (let i = N/2 - 5; i < N/2 + 5; i++) {
    const idx = IX(i, N - 2);
    dens[idx] += 0.5;
    temp[idx] += 2.0;
    v[idx] -= 3.0;
    // Add slight random horizontal push for turbulence
    u[idx] += (Math.random() - 0.5) * 0.5;
  }
}

function applyBuoyancy() {
  for (let j = 1; j <= N; j++)
    for (let i = 1; i <= N; i++) {
      const idx = IX(i, j);
      v[idx] -= buoyancy * (temp[idx] - ambientTemp);
    }
}

// Reuse diffuse/advect/project from example 3
// (Simplified inline for this example)
function linSolve(x, x0, a, c) {
  for (let k = 0; k < 4; 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)])) / c;
}

function project_step() {
  const div = new Float32Array(size), p = new Float32Array(size);
  for (let j = 1; j <= N; j++)
    for (let i = 1; i <= N; i++)
      div[IX(i,j)] = -0.5 * (u[IX(i+1,j)]-u[IX(i-1,j)]+v[IX(i,j+1)]-v[IX(i,j-1)]) / N;
  linSolve(p, div, 1, 4);
  for (let j = 1; j <= N; j++)
    for (let i = 1; i <= N; i++) {
      u[IX(i,j)] -= 0.5*N*(p[IX(i+1,j)]-p[IX(i-1,j)]);
      v[IX(i,j)] -= 0.5*N*(p[IX(i,j+1)]-p[IX(i,j-1)]);
    }
}

function advectField(d, d0, dt) {
  const dt0 = dt * N;
  for (let j = 1; j <= N; j++)
    for (let i = 1; i <= N; i++) {
      let x = i - dt0 * u[IX(i,j)];
      let y = j - dt0 * v[IX(i,j)];
      x = Math.max(0.5, Math.min(N + 0.5, x));
      y = Math.max(0.5, Math.min(N + 0.5, y));
      const i0 = Math.floor(x), j0 = Math.floor(y);
      const s = x - i0, t = y - j0;
      d[IX(i,j)] = (1-s)*((1-t)*d0[IX(i0,j0)]+t*d0[IX(i0,j0+1)])
                  + s*((1-t)*d0[IX(i0+1,j0)]+t*d0[IX(i0+1,j0+1)]);
    }
}

function step() {
  addSource();
  applyBuoyancy();
  // Velocity step
  const a = 0.1 * 0.0 * N * N;
  u0.set(u); v0.set(v);
  linSolve(u, u0, a, 1 + 4*a);
  linSolve(v, v0, a, 1 + 4*a);
  project_step();
  u0.set(u); v0.set(v);
  advectField(u, u0, 0.1);
  advectField(v, v0, 0.1);
  project_step();
  // Density + temperature
  dens0.set(dens); temp0.set(temp);
  advectField(dens, dens0, 0.1);
  advectField(temp, temp0, 0.1);
  // Decay
  for (let i = 0; i < size; i++) {
    dens[i] *= 0.995;
    temp[i] *= 0.99;
  }
}

function render() {
  step();
  ctx.fillStyle = '#000';
  ctx.fillRect(0, 0, canvas.width, canvas.height);
  const img = ctx.createImageData(canvas.width, canvas.height);
  const sx = canvas.width / (N + 2), sy = canvas.height / (N + 2);
  for (let y = 0; y < canvas.height; y++)
    for (let x = 0; x < canvas.width; x++) {
      const i = Math.floor(x / sx), j = Math.floor(y / sy);
      const d = Math.min(1, dens[IX(i,j)]);
      const t = Math.min(1, temp[IX(i,j)] * 0.5);
      const idx = (y * canvas.width + x) * 4;
      img.data[idx] = d * 80 + t * 175;
      img.data[idx+1] = d * 80 + t * 40;
      img.data[idx+2] = d * 120;
      img.data[idx+3] = 255;
    }
  ctx.putImageData(img, 0, 0);
  requestAnimationFrame(render);
}
render();

The smoke rises from the bottom, curls at the top, and slowly dissipates. The temperature field drives the buoyancy force, which creates the characteristic mushroom-cap shape. The slight random horizontal push at the source breaks symmetry and creates natural-looking turbulence. Warm areas glow orange-red; cooler smoke appears blue-gray.

5. Interactive ink drops — mouse-driven fluid

The most satisfying fluid interaction: click to drop colored ink, drag to push it around. Each click uses a different hue, creating marble-like patterns as the colors swirl together. This is the classic "fluid art" demo that launched a thousand screensavers.

const canvas = document.createElement('canvas');
canvas.width = 500; canvas.height = 500;
document.body.appendChild(canvas);
const ctx = canvas.getContext('2d');

const N = 128;
const s = N + 2;
const IX = (i, j) => i + s * j;
const total = s * s;

let u = new Float32Array(total), v = new Float32Array(total);
let r = new Float32Array(total), g = new Float32Array(total), b = new Float32Array(total);
let u0 = new Float32Array(total), v0 = new Float32Array(total);
let r0 = new Float32Array(total), g0 = new Float32Array(total), b0 = new Float32Array(total);
let hue = 0;

function solve(x, x0, a, c) {
  for (let k = 0; k < 4; 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)])) / c;
}

function proj() {
  const div = new Float32Array(total), p = new Float32Array(total);
  for (let j=1;j<=N;j++) for (let i=1;i<=N;i++)
    div[IX(i,j)] = -0.5*(u[IX(i+1,j)]-u[IX(i-1,j)]+v[IX(i,j+1)]-v[IX(i,j-1)])/N;
  solve(p, div, 1, 4);
  for (let j=1;j<=N;j++) for (let i=1;i<=N;i++) {
    u[IX(i,j)] -= 0.5*N*(p[IX(i+1,j)]-p[IX(i-1,j)]);
    v[IX(i,j)] -= 0.5*N*(p[IX(i,j+1)]-p[IX(i,j-1)]);
  }
}

function adv(d, d0, dt) {
  const dt0 = dt*N;
  for (let j=1;j<=N;j++) for (let i=1;i<=N;i++) {
    let x = i - dt0*u[IX(i,j)], y = j - dt0*v[IX(i,j)];
    x = Math.max(0.5, Math.min(N+0.5, x));
    y = Math.max(0.5, Math.min(N+0.5, y));
    const i0=Math.floor(x), j0=Math.floor(y), s1=x-i0, t1=y-j0;
    d[IX(i,j)] = (1-s1)*((1-t1)*d0[IX(i0,j0)]+t1*d0[IX(i0,j0+1)])
                + s1*((1-t1)*d0[IX(i0+1,j0)]+t1*d0[IX(i0+1,j0+1)]);
  }
}

function hslToRgb(h, s, l) {
  const c = (1 - Math.abs(2*l - 1)) * s;
  const x = c * (1 - Math.abs((h/60)%2 - 1));
  const m = l - c/2;
  let rr, gg, bb;
  if (h < 60) { rr=c; gg=x; bb=0; }
  else if (h < 120) { rr=x; gg=c; bb=0; }
  else if (h < 180) { rr=0; gg=c; bb=x; }
  else if (h < 240) { rr=0; gg=x; bb=c; }
  else if (h < 300) { rr=x; gg=0; bb=c; }
  else { rr=c; gg=0; bb=x; }
  return [(rr+m), (gg+m), (bb+m)];
}

let mx=0, my=0, pmx=0, pmy=0, down=false;
canvas.onmousemove = (e) => {
  const rc = canvas.getBoundingClientRect();
  pmx=mx; pmy=my;
  mx=(e.clientX-rc.left)/canvas.width*N;
  my=(e.clientY-rc.top)/canvas.height*N;
};
canvas.onmousedown = () => { down=true; hue=(hue+47)%360; };
canvas.onmouseup = () => down=false;

function step() {
  if (down) {
    const ci=Math.floor(mx), cj=Math.floor(my);
    const [cr,cg,cb] = hslToRgb(hue, 1, 0.5);
    for (let di=-2;di<=2;di++) for (let dj=-2;dj<=2;dj++) {
      const idx = IX(ci+di, cj+dj);
      if (ci+di>0 && ci+di<=N && cj+dj>0 && cj+dj<=N) {
        r[idx]+=cr*2; g[idx]+=cg*2; b[idx]+=cb*2;
        u[idx]+=(mx-pmx)*15; v[idx]+=(my-pmy)*15;
      }
    }
  }
  u0.set(u); v0.set(v);
  solve(u, u0, 0, 1); solve(v, v0, 0, 1);
  proj();
  u0.set(u); v0.set(v);
  adv(u, u0, 0.1); adv(v, v0, 0.1);
  proj();
  for (const ch of [[r,r0],[g,g0],[b,b0]]) {
    ch[1].set(ch[0]); adv(ch[0], ch[1], 0.1);
  }
  for (let i=0;i

const canvas = document.createElement('canvas');
canvas.width = 600; canvas.height = 300;
document.body.appendChild(canvas);
const ctx = canvas.getContext('2d');

const NX = 200, NY = 100;
const IX = (i, j) => i + (NX+2) * j;
const total = (NX+2) * (NY+2);

let u = new Float32Array(total), v = new Float32Array(total);
let u0 = new Float32Array(total), v0 = new Float32Array(total);
let dye = new Float32Array(total);
let dye0 = new Float32Array(total);

// Obstacle: circle at 1/4 width, center height
const obsX = NX * 0.25, obsY = NY * 0.5, obsR = NY * 0.08;

function isObstacle(i, j) {
  const dx = i - obsX, dy = j - obsY;
  return dx*dx + dy*dy < obsR*obsR;
}

function solve(x, x0, a, c) {
  for (let k=0;k<4;k++)
    for (let j=1;j<=NY;j++) for (let i=1;i<=NX;i++)
      if (!isObstacle(i,j))
        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)]))/c;
}

function adv(d, d0, uu, vv, dt) {
  const dt0x = dt*NX, dt0y = dt*NY;
  for (let j=1;j<=NY;j++) for (let i=1;i<=NX;i++) {
    if (isObstacle(i,j)) { d[IX(i,j)]=0; continue; }
    let x = i - dt0x*uu[IX(i,j)], y = j - dt0y*vv[IX(i,j)];
    x = Math.max(0.5, Math.min(NX+0.5, x));
    y = Math.max(0.5, Math.min(NY+0.5, y));
    const i0=Math.floor(x), j0=Math.floor(y), s=x-i0, t=y-j0;
    d[IX(i,j)] = (1-s)*((1-t)*d0[IX(i0,j0)]+t*d0[IX(i0,j0+1)])
                + s*((1-t)*d0[IX(i0+1,j0)]+t*d0[IX(i0+1,j0+1)]);
  }
}

function proj() {
  const div=new Float32Array(total), p=new Float32Array(total);
  for (let j=1;j<=NY;j++) for (let i=1;i<=NX;i++) {
    if (isObstacle(i,j)) continue;
    div[IX(i,j)] = -0.5*(u[IX(i+1,j)]-u[IX(i-1,j)]+v[IX(i,j+1)]-v[IX(i,j-1)])/NX;
  }
  solve(p, div, 1, 4);
  for (let j=1;j<=NY;j++) for (let i=1;i<=NX;i++) {
    if (isObstacle(i,j)) { u[IX(i,j)]=0; v[IX(i,j)]=0; continue; }
    u[IX(i,j)] -= 0.5*NX*(p[IX(i+1,j)]-p[IX(i-1,j)]);
    v[IX(i,j)] -= 0.5*NY*(p[IX(i,j+1)]-p[IX(i,j-1)]);
  }
}

function step() {
  // Inflow from left
  for (let j=1;j<=NY;j++) {
    u[IX(1,j)] = 4.0;
    if (j % 6 < 3) dye[IX(1,j)] = 0.8;
  }
  // Enforce obstacle
  for (let j=0;j 0 ? vc*180 : 0);
      img.data[pi+1] = d*100;
      img.data[pi+2] = d*60 + (curl < 0 ? vc*180 : 0);
      img.data[pi+3] = 255;
    }
  }
  ctx.putImageData(img, 0, 0);
  requestAnimationFrame(render);
}
render();

const canvas = document.createElement('canvas');
canvas.width = 500; canvas.height = 400;
document.body.appendChild(canvas);
const ctx = canvas.getContext('2d');

const numParticles = 300;
const h = 16;          // smoothing radius
const restDensity = 1;
const stiffness = 200;
const viscosity = 10;
const gravity = 0.5;
const dt = 0.016;

const particles = [];
for (let i = 0; i < numParticles; i++) {
  particles.push({
    x: 100 + Math.random() * 150,
    y: 50 + Math.random() * 150,
    vx: 0, vy: 0,
    density: 0, pressure: 0,
    fx: 0, fy: 0
  });
}

function poly6(r2) {
  if (r2 >= h * h) return 0;
  const diff = h * h - r2;
  return 315 / (64 * Math.PI * Math.pow(h, 9)) * diff * diff * diff;
}

function spikyGrad(r) {
  if (r >= h || r < 0.001) return 0;
  const diff = h - r;
  return -45 / (Math.PI * Math.pow(h, 6)) * diff * diff;
}

function viscLap(r) {
  if (r >= h) return 0;
  return 45 / (Math.PI * Math.pow(h, 6)) * (h - r);
}

function step() {
  // Compute density and pressure
  for (const p of particles) {
    p.density = 0;
    for (const q of particles) {
      const dx = q.x - p.x, dy = q.y - p.y;
      p.density += poly6(dx*dx + dy*dy);
    }
    p.pressure = stiffness * (p.density - restDensity);
  }

  // Compute forces
  for (const p of particles) {
    p.fx = 0; p.fy = gravity * p.density;
    for (const q of particles) {
      if (p === q) continue;
      const dx = q.x - p.x, dy = q.y - p.y;
      const r = Math.sqrt(dx*dx + dy*dy);
      if (r >= h || r < 0.001) continue;
      const nx = dx/r, ny = dy/r;
      // Pressure force
      const pf = -0.5 * (p.pressure + q.pressure) / q.density * spikyGrad(r);
      p.fx += pf * nx;
      p.fy += pf * ny;
      // Viscosity force
      const vf = viscosity / q.density * viscLap(r);
      p.fx += vf * (q.vx - p.vx);
      p.fy += vf * (q.vy - p.vy);
    }
  }

  // Integrate
  for (const p of particles) {
    const ax = p.fx / p.density;
    const ay = p.fy / p.density;
    p.vx += ax * dt;
    p.vy += ay * dt;
    p.x += p.vx * dt * 60;
    p.y += p.vy * dt * 60;

    // Boundary collision
    if (p.x < 5) { p.x = 5; p.vx *= -0.3; }
    if (p.x > canvas.width-5) { p.x = canvas.width-5; p.vx *= -0.3; }
    if (p.y < 5) { p.y = 5; p.vy *= -0.3; }
    if (p.y > canvas.height-5) { p.y = canvas.height-5; p.vy *= -0.3; }
  }
}

canvas.addEventListener('click', (e) => {
  const rect = canvas.getBoundingClientRect();
  const cx = e.clientX - rect.left, cy = e.clientY - rect.top;
  for (const p of particles) {
    const dx = p.x - cx, dy = p.y - cy;
    const d = Math.sqrt(dx*dx + dy*dy);
    if (d < 100) {
      const force = (100 - d) * 0.15;
      p.vx += (dx/d) * force;
      p.vy += (dy/d) * force;
    }
  }
});

function render() {
  step();
  ctx.fillStyle = 'rgba(5, 10, 20, 0.3)';
  ctx.fillRect(0, 0, canvas.width, canvas.height);

  for (const p of particles) {
    const speed = Math.sqrt(p.vx*p.vx + p.vy*p.vy);
    const t = Math.min(1, speed * 0.3);
    const r = 30 + t * 200;
    const g = 100 + t * 100;
    const b = 220 - t * 100;
    ctx.beginPath();
    ctx.arc(p.x, p.y, 4, 0, Math.PI * 2);
    ctx.fillStyle = `rgb(${r},${g},${b})`;
    ctx.fill();
  }
  requestAnimationFrame(render);
}
render();

const canvas = document.createElement('canvas');
canvas.width = 500; canvas.height = 500;
document.body.appendChild(canvas);
const ctx = canvas.getContext('2d');

const N = 256;
let current = new Float32Array(N * N);
let previous = new Float32Array(N * N);
const damping = 0.99;

function step() {
  const next = new Float32Array(N * N);
  for (let j = 1; j < N - 1; j++) {
    for (let i = 1; i < N - 1; i++) {
      const idx = i + j * N;
      next[idx] = (
        current[idx - 1] + current[idx + 1] +
        current[idx - N] + current[idx + N]
      ) / 2 - previous[idx];
      next[idx] *= damping;
    }
  }
  previous = current;
  current = next;
}

canvas.addEventListener('click', (e) => {
  const rect = canvas.getBoundingClientRect();
  const cx = Math.floor((e.clientX - rect.left) / canvas.width * N);
  const cy = Math.floor((e.clientY - rect.top) / canvas.height * N);
  for (let di = -4; di <= 4; di++)
    for (let dj = -4; dj <= 4; dj++) {
      const r = Math.sqrt(di*di + dj*dj);
      if (r < 5) {
        const idx = (cx+di) + (cy+dj) * N;
        if (idx >= 0 && idx < N*N)
          current[idx] = Math.cos(r * Math.PI / 5) * 2;
      }
    }
});

// Auto-drop ripples
setInterval(() => {
  const x = 30 + Math.random() * (N - 60);
  const y = 30 + Math.random() * (N - 60);
  for (let di = -3; di <= 3; di++)
    for (let dj = -3; dj <= 3; dj++) {
      const idx = Math.floor(x+di) + Math.floor(y+dj) * N;
      if (idx >= 0 && idx < N*N) current[idx] = 1.5;
    }
}, 2000);

function render() {
  step(); step(); // 2 substeps for faster propagation
  const img = ctx.createImageData(canvas.width, canvas.height);
  const scale = canvas.width / N;

  for (let y = 0; y < canvas.height; y++) {
    for (let x = 0; x < canvas.width; x++) {
      const gi = Math.floor(x / scale);
      const gj = Math.floor(y / scale);
      const idx = gi + gj * N;
      const val = current[idx];

      // Compute fake normal for lighting
      const dx = gi > 0 && gi < N-1 ? current[idx+1] - current[idx-1] : 0;
      const dy = gj > 0 && gj < N-1 ? current[idx+N] - current[idx-N] : 0;

      // Base water color + specular highlight
      const light = Math.max(0, -dx * 3 + dy * 2);
      const spec = Math.pow(Math.max(0, light), 8) * 200;
      const depth = 0.5 + val * 0.15;

      const pi = (y * canvas.width + x) * 4;
      img.data[pi]     = Math.min(255, 10 + spec + Math.max(0, val * 30));
      img.data[pi + 1] = Math.min(255, 40 * depth + spec * 0.8);
      img.data[pi + 2] = Math.min(255, 120 * depth + spec * 0.5 + Math.abs(val) * 60);
      img.data[pi + 3] = 255;
    }
  }
  ctx.putImageData(img, 0, 0);
  requestAnimationFrame(render);
}
render();