physics

FLUID DYNAMICS

Real-time 2D Navier-Stokes solver. Stir the fluid, visualize vortices, and explore incompressible flow dynamics.

0.0001
0.0001
0.999
1.0×

How it works

Stable Fluids (Jos Stam)
Simulates incompressible 2D flow using the Navier-Stokes equations. Each frame: (1) add velocity from mouse, (2) diffuse velocity (viscosity), (3) project (enforce incompressibility), (4) advect velocity and density along flow field.

Reynolds Number: Re = ρvL/μ. Low Re = smooth laminar flow; high Re = chaotic turbulent vortices.

Controls: Drag mouse to stir. Viscosity controls resistance (honey-like vs inviscid). Obstacles create vortex shedding. Wind tunnel mode shows Karman vortex streets.
Developer Reference

Core Algorithm & Standalone Script

Standalone, zero-dependency JavaScript implementation powering this tool. Free to inspect, copy, and build upon.

const canvas = document.getElementById('fluidCanvas');
    const ctx = canvas.getContext('2d');

    // Grid resolution
    const GRID_X = 128;
    const GRID_Y = 128;
    const ITERATIONS = 20;

    // Physics state
    let u = new Float32Array(GRID_X * GRID_Y);
    let v = new Float32Array(GRID_X * GRID_Y);
    let u_prev = new Float32Array(GRID_X * GRID_Y);
    let v_prev = new Float32Array(GRID_X * GRID_Y);
    let density = new Float32Array(GRID_X * GRID_Y);
    let density_prev = new Float32Array(GRID_X * GRID_Y);

    // Obstacle masks
    let obstacles = new Uint8Array(GRID_X * GRID_Y);

    // Parameters
    let params = {
      viscosity: 0.0001,
      diffusion: 0.0001,
      decay: 0.999,
      colorScheme: 'smoke',
      flowMode: 'interactive',
      showCylinder: false,
      showWing: false,
      paused: false,
      simSpeed: 1
    };

    // Mouse tracking
    let mouseX = 0, mouseY = 0;
    let lastMouseX = 0, lastMouseY = 0;
    let mousePressed = false;

    // Resize canvas to device pixel ratio for crisp rendering
    function resizeCanvas() {
      const dpr = window.devicePixelRatio || 1;
      const rect = canvas.parentElement.getBoundingClientRect();
      const size = Math.min(rect.width, rect.height);

      canvas.width = size * dpr;
      canvas.height = size * dpr;
      ctx.scale(dpr, dpr);
      canvas.style.width = size + 'px';
      canvas.style.height = size + 'px';
    }

    resizeCanvas();
    window.addEventListener('resize', resizeCanvas);

    // Get array index
    const IX = (x, y) => Math.floor(x) + Math.floor(y) * GRID_X;

    // Boundary conditions
    function setBounds(b, x) {
      for (let i = 1; i < GRID_X - 1; i++) {
        for (let j = 1; j < GRID_Y - 1; j++) {
          if (obstacles[IX(i, j)]) continue;
          const idx = IX(i, j);

          // Enforce no-slip at obstacles
          if (obstacles[IX(i-1, j)] || obstacles[IX(i+1, j)] ||
              obstacles[IX(i, j-1)] || obstacles[IX(i, j+1)]) {
            if (b === 1) x[idx] = 0; // u
            if (b === 2) x[idx] = 0; // v
            if (b === 0) x[idx] = 0; // density
          }
        }
      }
    }

    // Diffusion solver (Gauss-Seidel)
    function diffuse(b, x, x0, diff, dt) {
      const a = dt * diff * (GRID_X - 2) * (GRID_Y - 2);

      for (let iter = 0; iter < ITERATIONS; iter++) {
        for (let i = 1; i < GRID_X - 1; i++) {
          for (let j = 1; j < GRID_Y - 1; j++) {
            const idx = IX(i, j);
            if (obstacles[idx]) continue;

            x[idx] = (x0[idx] + a * (
              x[IX(i+1, j)] + x[IX(i-1, j)] +
              x[IX(i, j+1)] + x[IX(i, j-1)]
            )) / (1 + 4 * a);
          }
        }
        setBounds(b, x);
      }
    }

    // Project (enforce divergence-free condition)
    function project(u_field, v_field, p, div) {
      const h = 1.0 / (GRID_X - 2);

      // Compute divergence
      for (let i = 1; i < GRID_X - 1; i++) {
        for (let j = 1; j < GRID_Y - 1; j++) {
          const idx = IX(i, j);
          div[idx] = -0.5 * h * (
            u_field[IX(i+1, j)] - u_field[IX(i-1, j)] +
            v_field[IX(i, j+1)] - v_field[IX(i, j-1)]
          );
          p[idx] = 0;
        }
      }

      setBounds(0, div);
      setBounds(0, p);

      // Solve Poisson equation
      for (let iter = 0; iter < ITERATIONS; iter++) {
        for (let i = 1; i < GRID_X - 1; i++) {
          for (let j = 1; j < GRID_Y - 1; j++) {
            const idx = IX(i, j);
            if (obstacles[idx]) continue;

            p[idx] = (div[idx] + (
              p[IX(i+1, j)] + p[IX(i-1, j)] +
              p[IX(i, j+1)] + p[IX(i, j-1)]
            )) / 4;
          }
        }
        setBounds(0, p);
      }

      // Subtract gradient
      for (let i = 1; i < GRID_X - 1; i++) {
        for (let j = 1; j < GRID_Y - 1; j++) {
          const idx = IX(i, j);
          u_field[idx] -= 0.5 * (p[IX(i+1, j)] - p[IX(i-1, j)]) / h;
          v_field[idx] -= 0.5 * (p[IX(i, j+1)] - p[IX(i, j-1)]) / h;
        }
      }

      setBounds(1, u_field);
      setBounds(2, v_field);
    }

    // Semi-Lagrangian advection
    function advect(b, d, d0, u_field, v_field, dt) {
      const dt0 = dt * (GRID_X - 2);

      for (let i = 1; i < GRID_X - 1; i++) {
        for (let j = 1; j < GRID_Y - 1; j++) {
          const idx = IX(i, j);
          if (obstacles[idx]) continue;

          // Trace backwards
          let x = i - dt0 * u_field[idx];
          let y = j - dt0 * v_field[idx];

          // Clamp
          x = Math.max(0.5, Math.min(GRID_X - 1.5, x));
          y = Math.max(0.5, Math.min(GRID_Y - 1.5, y));

          const i0 = Math.floor(x);
          const i1 = i0 + 1;
          const j0 = Math.floor(y);
          const j1 = j0 + 1;

          const sx = x - i0;
          const sy = y - j0;

          // Bilinear interpolation
          d[idx] = (1 - sx) * (1 - sy) * d0[IX(i0, j0)] +
                   sx * (1 - sy) * d0[IX(i1, j0)] +
                   (1 - sx) * sy * d0[IX(i0, j1)] +
                   sx * sy * d0[IX(i1, j1)];
        }
      }
      setBounds(b, d);
    }

    // Obstacle helpers
    function drawCircle(cx, cy, radius) {
      for (let i = 0; i < GRID_X; i++) {
        for (let j = 0; j < GRID_Y; j++) {
          const dx = (i - cx) / GRID_X;
          const dy = (j - cy) / GRID_Y;
          const dist = Math.sqrt(dx * dx + dy * dy);
          if (dist < radius) {
            obstacles[IX(i, j)] = 1;
          }
        }
      }
    }

    function drawWing(cx, cy) {
      for (let i = 0; i < GRID_X; i++) {
        for (let j = 0; j < GRID_Y; j++) {
          const dx = i - cx;
          const dy = j - cy;

          // Airfoil-like shape
          const nx = dx / GRID_X;
          const ny = dy / GRID_Y;
          const angle = Math.atan2(ny, nx);
          const r = Math.sqrt(nx * nx + ny * ny);

          if (r < 0.08 && Math.abs(ny) < 0.04 + 0.02 * (0.08 - r) / 0.08) {
            obstacles[IX(i, j)] = 1;
          }
        }
      }
    }

    function clearObstacles() {
      obstacles.fill(0);
    }

    // Add force/density from mouse
    function addForce(x, y, amountU, amountV, amountDensity) {
      const xi = Math.min(GRID_X - 2, Math.max(1, Math.floor(x * GRID_X)));
      const yi = Math.min(GRID_Y - 2, Math.max(1, Math.floor(y * GRID_Y)));

      const radius = 4;
      for (let i = Math.max(0, xi - radius); i <= Math.min(GRID_X - 1, xi + radius); i++) {
        for (let j = Math.max(0, yi - radius); j <= Math.min(GRID_Y - 1, yi + radius); j++) {
          const dist = Math.sqrt((i - xi) ** 2 + (j - yi) ** 2);
          const falloff = Math.max(0, 1 - dist / radius);
          const idx = IX(i, j);
          u[idx] += amountU * falloff;
          v[idx] += amountV * falloff;
          density[idx] = Math.min(1, density[idx] + amountDensity * falloff);
        }
      }
    }

    // Simulation step
    let p = new Float32Array(GRID_X * GRID_Y);
    let div = new Float32Array(GRID_X * GRID_Y);

    function step(dt = 0.016) {
      if (params.paused) return;

      // Add forces based on flow mode
      if (params.flowMode === 'windtunnel') {
        for (let j = 0; j < GRID_Y; j++) {
          u[IX(5, j)] += 0.5;
          density[IX(5, j)] = Math.max(density[IX(5, j)], 0.5);
        }
      } else if (params.flowMode === 'vortexpair') {
        // Two counter-rotating vortices
        const t = Date.now() * 0.001;
        const cx1 = GRID_X * 0.35;
        const cy1 = GRID_Y * 0.5;
        const cx2 = GRID_X * 0.65;
        const cy2 = GRID_Y * 0.5;

        for (let i = 1; i < GRID_X - 1; i++) {
          for (let j = 1; j < GRID_Y - 1; j++) {
            const idx = IX(i, j);
            const dx1 = (i - cx1) * 0.1;
            const dy1 = (j - cy1) * 0.1;
            const r1 = Math.sqrt(dx1 * dx1 + dy1 * dy1) + 0.1;

            const dx2 = (i - cx2) * 0.1;
            const dy2 = (j - cy2) * 0.1;
            const r2 = Math.sqrt(dx2 * dx2 + dy2 * dy2) + 0.1;

            u[idx] += (dy1 / r1 - dy2 / r2) * 0.01;
            v[idx] += (-dx1 / r1 + dx2 / r2) * 0.01;
            density[idx] *= 0.999;
          }
        }
      } else if (params.flowMode === 'rising') {
        // Hot air rising
        const xi = Math.floor(0.5 * GRID_X);
        for (let i = Math.max(0, xi - 5); i <= Math.min(GRID_X - 1, xi + 5); i++) {
          for (let j = Math.max(GRID_Y - 20, 0); j < GRID_Y; j++) {
            const idx = IX(i, j);
            v[idx] += 0.3;
            density[idx] = Math.min(1, density[idx] + 0.5);
          }
        }
      }

      // Velocity step
      diffuse(1, u_prev, u, params.viscosity, dt);
      diffuse(2, v_prev, v, params.viscosity, dt);

      u.set(u_prev);
      v.set(v_prev);

      project(u, v, p, div);

      u_prev.set(u);
      v_prev.set(v);

      advect(1, u, u_prev, u_prev, v_prev, dt);
      advect(2, v, v_prev, u_prev, v_prev, dt);

      project(u, v, p, div);

      // Density step
      diffuse(0, density_prev, density, params.diffusion, dt);
      advect(0, density, density_prev, u, v, dt);

      // Decay
      for (let i = 0; i < density.length; i++) {
        density[i] *= params.decay;
      }
    }

    // Color mapping
    function getColor(val) {
      val = Math.max(0, Math.min(1, val));

      switch (params.colorScheme) {
        case 'smoke':
          const gray = Math.floor(val * 255);
          return `rgb(${gray},${gray},${gray})`;

        case 'fire':
          if (val < 0.25) {
            const r = Math.floor(val * 4 * 255);
            return `rgb(${r},0,0)`;
          } else if (val < 0.5) {
            const r = 255;
            const g = Math.floor((val - 0.25) * 4 * 255);
            return `rgb(${r},${g},0)`;
          } else if (val < 0.75) {
            const r = 255;
            const g = 255;
            const b = Math.floor((val - 0.5) * 4 * 255);
            return `rgb(${r},${g},${b})`;
          } else {
            const r = Math.floor(255 + (val - 0.75) * 4 * 255);
            const g = 255;
            const b = 255;
            return `rgb(${r},${g},${b})`;
          }

        case 'ocean':
          if (val < 0.3) {
            const b = Math.floor(val / 0.3 * 139);
            return `rgb(0,0,${b})`;
          } else if (val < 0.7) {
            const g = Math.floor((val - 0.3) / 0.4 * 255);
            return `rgb(0,${g},139)`;
          } else {
            const r = Math.floor((val - 0.7) / 0.3 * 255);
            const g = 255;
            const b = 255;
            return `rgb(${r},${g},${b})`;
          }

        case 'rainbow':
          const hue = val * 360;
          return `hsl(${hue},100%,50%)`;

        case 'neon':
          if (val < 0.33) {
            const r = 255;
            const g = Math.floor(val / 0.33 * 136);
            return `rgb(${r},${g},0)`;
          } else if (val < 0.66) {
            const r = Math.floor(255 - (val - 0.33) / 0.33 * 255);
            const g = Math.floor(136 + (val - 0.33) / 0.33 * 119);
            return `rgb(${r},${g},0)`;
          } else {
            const r = 0;
            const g = 255;
            const b = Math.floor((val - 0.66) / 0.34 * 255);
            return `rgb(${r},${g},${b})`;
          }
      }
    }

    // Render
    function render() {
      const imageData = ctx.createImageData(canvas.width / (window.devicePixelRatio || 1),
                                          canvas.height / (window.devicePixelRatio || 1));
      const data = imageData.data;

      for (let i = 0; i < GRID_X; i++) {
        for (let j = 0; j < GRID_Y; j++) {
          const idx = IX(i, j);
          const val = density[idx];
          const pixelIdx = (j * GRID_X + i) * 4;

          // Simple color mapping to RGB
          const color = getColor(val);
          const matches = color.match(/\d+/g);
          data[pixelIdx] = parseInt(matches[0]);
          data[pixelIdx + 1] = parseInt(matches[1]);
          data[pixelIdx + 2] = parseInt(matches[2]);
          data[pixelIdx + 3] = 255;
        }
      }

      ctx.putImageData(imageData, 0, 0);

      // Draw obstacles
      if (params.showCylinder || params.showWing) {
        ctx.strokeStyle = 'rgba(255, 34, 0, 0.4)';
        ctx.lineWidth = 2;

        if (params.showCylinder) {
          const cx = canvas.width / (window.devicePixelRatio || 1) * 0.5;
          const cy = canvas.height / (window.devicePixelRatio || 1) * 0.5;
          const r = (canvas.width / (window.devicePixelRatio || 1)) * 0.08;
          ctx.beginPath();
          ctx.arc(cx, cy, r, 0, Math.PI * 2);
          ctx.stroke();
        }

        if (params.showWing) {
          const cx = canvas.width / (window.devicePixelRatio || 1) * 0.5;
          const cy = canvas.height / (window.devicePixelRatio || 1) * 0.5;
          const w = (canvas.width / (window.devicePixelRatio || 1)) * 0.12;
          const h = (canvas.height / (window.devicePixelRatio || 1)) * 0.05;
          ctx.beginPath();
          ctx.ellipse(cx, cy, w, h, 0, 0, Math.PI * 2);
          ctx.stroke();
        }
      }
    }

    // Animation loop
    function animate() {
      step(0.016 * params.simSpeed);
      render();
      requestAnimationFrame(animate);
    }

    animate();

    // Mouse events
    canvas.addEventListener('mousedown', (e) => {
      mousePressed = true;
      const rect = canvas.getBoundingClientRect();
      const dpr = window.devicePixelRatio || 1;
      mouseX = (e.clientX - rect.left) * dpr / rect.width;
      mouseY = (e.clientY - rect.top) * dpr / rect.height;
      lastMouseX = mouseX;
      lastMouseY = mouseY;
    });

    canvas.addEventListener('mousemove', (e) => {
      const rect = canvas.getBoundingClientRect();
      const dpr = window.devicePixelRatio || 1;
      mouseX = (e.clientX - rect.left) * dpr / rect.width;
      mouseY = (e.clientY - rect.top) * dpr / rect.height;

      if (mousePressed && params.flowMode === 'interactive') {
        const amountU = (mouseX - lastMouseX) * 300;
        const amountV = (mouseY - lastMouseY) * 300;
        addForce(mouseX, mouseY, amountU, amountV, 1);
      }

      lastMouseX = mouseX;
      lastMouseY = mouseY;
    });

    canvas.addEventListener('mouseup', () => {
      mousePressed = false;
    });

    canvas.addEventListener('mouseleave', () => {
      mousePressed = false;
    });

    // Touch support
    canvas.addEventListener('touchstart', (e) => {
      mousePressed = true;
      const rect = canvas.getBoundingClientRect();
      const touch = e.touches[0];
      const dpr = window.devicePixelRatio || 1;
      mouseX = (touch.clientX - rect.left) * dpr / rect.width;
      mouseY = (touch.clientY - rect.top) * dpr / rect.height;
      lastMouseX = mouseX;
      lastMouseY = mouseY;
    });

    canvas.addEventListener('touchmove', (e) => {
      e.preventDefault();
      const rect = canvas.getBoundingClientRect();
      const touch = e.touches[0];
      const dpr = window.devicePixelRatio || 1;
      mouseX = (touch.clientX - rect.left) * dpr / rect.width;
      mouseY = (touch.clientY - rect.top) * dpr / rect.height;

      if (mousePressed && params.flowMode === 'interactive') {
        const amountU = (mouseX - lastMouseX) * 300;
        const amountV = (mouseY - lastMouseY) * 300;
        addForce(mouseX, mouseY, amountU, amountV, 1);
      }

      lastMouseX = mouseX;
      lastMouseY = mouseY;
    });

    canvas.addEventListener('touchend', () => {
      mousePressed = false;
    });

    // Control event listeners
    document.getElementById('colorScheme').addEventListener('change', (e) => {
      params.colorScheme = e.target.value;
    });

    document.getElementById('flowMode').addEventListener('change', (e) => {
      params.flowMode = e.target.value;
      clearObstacles();

      if (params.flowMode === 'windtunnel' && params.showCylinder) {
        drawCircle(GRID_X * 0.5, GRID_Y * 0.5, 0.12);
      }
      if (params.flowMode === 'windtunnel' && params.showWing) {
        drawWing(GRID_X * 0.5, GRID_Y * 0.5);
      }
    });

    document.getElementById('viscosity').addEventListener('input', (e) => {
      params.viscosity = parseFloat(e.target.value);
      document.getElementById('viscosityValue').textContent = params.viscosity.toFixed(4);
    });

    document.getElementById('diffusion').addEventListener('input', (e) => {
      params.diffusion = parseFloat(e.target.value);
      document.getElementById('diffusionValue').textContent = params.diffusion.toFixed(4);
    });

    document.getElementById('decay').addEventListener('input', (e) => {
      params.decay = parseFloat(e.target.value);
      document.getElementById('decayValue').textContent = params.decay.toFixed(3);
    });

    document.getElementById('cylinderToggle').addEventListener('click', (e) => {
      params.showCylinder = !params.showCylinder;
      e.target.classList.toggle('active');
      clearObstacles();

      if (params.showCylinder) {
        drawCircle(GRID_X * 0.5, GRID_Y * 0.5, 0.12);
      }
      if (params.showWing) {
        drawWing(GRID_X * 0.5, GRID_Y * 0.5);
      }
    });

    document.getElementById('wingToggle').addEventListener('click', (e) => {
      params.showWing = !params.showWing;
      e.target.classList.toggle('active');
      clearObstacles();

      if (params.showCylinder) {
        drawCircle(GRID_X * 0.5, GRID_Y * 0.5, 0.12);
      }
      if (params.showWing) {
        drawWing(GRID_X * 0.5, GRID_Y * 0.5);
      }
    });

    document.getElementById('resetBtn').addEventListener('click', () => {
      u.fill(0);
      v.fill(0);
      u_prev.fill(0);
      v_prev.fill(0);
      density.fill(0);
      density_prev.fill(0);
      clearObstacles();
      document.getElementById('cylinderToggle').classList.remove('active');
      document.getElementById('wingToggle').classList.remove('active');
      params.showCylinder = false;
      params.showWing = false;
    });

    document.getElementById('simSpeedSlider').addEventListener('input', (e) => {
      params.simSpeed = parseFloat(e.target.value);
      document.getElementById('simSpeedValue').textContent = params.simSpeed.toFixed(1);
    });

    document.getElementById('pauseBtn').addEventListener('click', (e) => {
      params.paused = !params.paused;
      e.target.textContent = params.paused ? 'Resume' : 'Pause';
      e.target.classList.toggle('active');
    });