The user wants a real-time fluid simulation on a full-screen canvas. This is a classic cfd (computational fluid dynamics) simulation using the Navier-Stokes equations. Let me think about the implementation:

1. **Fluid simulation approach**: I'll implement a real-time fluid simulation using the method described by Jos Stam (Real-Time Fluid Dynamics for Games). This uses:
   - Velocity field (u, v)
   - Pressure field
   - Density/dye field (the colored dye we see)
   - Diffusion and advection steps
   - Projection step to enforce incompressibility

2. **Features**:
   - Full-screen canvas
   - Mouse dragging adds velocity and dye
   - Colorful swirling dye patterns
   - Multiple colors or gradient colors based on velocity/movement

Let me implement this carefully. I'll use a grid-based simulation.

Key components:
- `Nx × Ny` grid (e.g., 128×96 or similar, scaled to canvas)
- Velocity fields `u`, `v` (staggered grid or collocated)
- Dye density field `dens`
- Operations: advect, diffuse, project (Hodge decomposition)

Actually, let me use a simpler but visually impressive approach: a vorticity-based or simpler fluid solver. The classic Stam solver is the most reliable.

Let me implement:

1. **Grid**: `N` cells in each dimension
2. **Fields**: `u`, `v` (velocity), `dens` (density/dye), and temporary fields
3. **Solver steps**:
   - `addSource`: add velocity/dye from mouse
   - `diffuse`: Gauss-Seidel relaxation
   - `advect`: semi-Lagrangian backtrace
   - `project`: enforce divergence-free velocity
4. **Rendering**: draw density field as colored pixels, with color based on velocity magnitude or position

For colors, I'll make the dye have a color that depends on where it was introduced or use velocity direction to determine hue. Let me use a technique where different mouse positions/velocities produce different colors — maybe assign a hue based on initial position and let it swirl.

Actually, a nice approach: use multiple dye channels (R, G, B) each with their own advection, colored differently. Or use a single dye field and map velocity to color.

Let me think of a visually stunning approach:
- Use the dye field with a color that depends on velocity magnitude/direction
- Or use multiple dye components (RGB) each advected separately, with the mouse injecting different colors based on movement direction or position

I'll go with: the mouse injects dye with a color that varies based on mouse position (mapped to hue) and the dye advects and diffuses, creating beautiful swirling patterns. I'll also add a slight rainbow effect based on velocity.

Actually, let me keep it relatively simple but visually impressive. I'll use:
- A single dye density field
- Color computed per-pixel from: base hue from mouse position + variation from velocity
- Or: three dye fields (R, G, B) each with its own source pattern

Wait, the cleanest approach for a beautiful result: use velocity-to-color mapping. The dye shows the flow patterns beautifully. Let me use a velocity magnitude + direction to compute color via HSV.

Let me implement the full Stam solver.

Grid size: I'll use maybe 150x100 or so for performance, with the canvas scaled up.

Actually, let me think about performance. For a full-screen canvas with real-time requirements, I should:
- Use a modest grid (e.g., 128×96)
- Render to an offscreen canvas at grid resolution, then scale up
- Or use ImageData directly

Let me use a canvas that matches the grid size, then CSS-scale it to full screen with image-rendering: pixelated or smooth.

Let me code this up:

```javascript
// Fluid solver constants
const N = 128; // grid resolution (cells)
const dt = 0.1; // time step
const diff = 0.0001; // diffusion rate
const visc = 0.0001; // viscosity

// Fields (size (N+2)*(N+2) for boundary cells)
const size = (N+2)*(N+2);
const u = new Float32Array(size);
const v = new Float32Array(size);
const uPrev = new Float32Array(size);
const vPrev = new Float32Array(size);
const dens = new Float32Array(size);
const densPrev = new Float32Array(size);
```

Wait, for the dye with color, I need RGB dye fields. Let me have 3 dye fields: r, g, b.

Actually, let me reconsider the visual approach for maximum beauty:

**Approach**: Use a single velocity field. Add dye in RGB channels. Each mouse drag injects dye with colors determined by:
- The direction of movement (angle → hue)
- Or the position (mapped to hue)

Then advect each dye channel separately. The velocity field advects all dye. This creates colorful swirling patterns.

Let me implement:

```javascript
const N = 100; // grid
const cellSize = ... // canvas size / N

// Fields
const u, v: velocity
const r, g, b: dye densities
const tmp fields for swapping
```

For the solver, I need:
- `setBnd` (boundary conditions)
- `linSolve` (diffusion solver)
- `advect` (semi-Lagrangian)
- `project` (pressure solve + subtract gradient)

Let me write the complete implementation.

For mouse interaction:
- On mousedown/mousemove: record mouse position in grid coordinates
- Add velocity in the direction of mouse movement (velocity = delta pos / dt)
- Add dye at the mouse position with color based on movement direction or position

Let me also add a nice touch: when not dragging, maybe add some subtle ambient flow, or just let the existing dye fade slowly.

For rendering:
- Each cell maps to a pixel on an offscreen canvas
- Color = (r, g, b) with some saturation/gamma
- Scale the offscreen canvas to full screen via CSS

Let me also consider: to make it more colorful, I'll have the dye color be based on the velocity direction (angle) when the dye is injected. This gives rainbow swirls.

Actually, an even nicer approach: inject dye with a color that depends on the mouse velocity direction, so as you drag in circles you get rainbow patterns. Let me also add a slight drift/ambient to keep things interesting.

Let me also consider using velocity magnitude to affect brightness.

OK let me just write the code. I'll aim for a beautiful, performant implementation.

Let me think about the grid indexing:
- Index: `IX(i, j) = i + (N+2)*j` where i, j range from 0 to N+1
- Interior cells: i, j from 1 to N
- Boundary: i=0, i=N+1, j=0, j=N+1

Stam solver functions:

```javascript
function addSource(x, s) {
  for (let i = 0; i < size; i++) x[i] += s[i];
}

function diffuse(b, x, x0, diff, dt) {
  // Gauss-Seidel
  const a = dt * diff * N * N;
  for (let k = 0; k < 20; k++) {
    for (let i = 1; i <= N; i++) {
      for (let j = 1; j <= N; 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)])) / (1 + 4*a);
      }
    }
    setBnd(b, x);
  }
}

function advect(b, d, d0, u, v, dt) {
  // Semi-Lagrangian backtrace
  for (let i = 1; i <= N; i++) {
    for (let j = 1; j <= N; j++) {
      const x = i - dt * N * u[IX(i,j)];
      const y = j - dt * N * v[IX(i,j)];
      // Clamp to [0.5, N+0.5]
      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) {
  // Compute divergence
  for (let i = 1; i <= N; i++) {
    for (let j = 1; j <= N; j++) {
      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);
  // Solve pressure Poisson
  for (let k = 0; k < 20; k++) {
    for (let i = 1; i <= N; i++) {
      for (let j = 1; j <= N; j++) {
        p[IX(i,j)] = (div[IX(i,j)] + p[IX(i-1,j)] + p[IX(i+1,j)] + p[IX(i,j-1)] + p[IX(i,j+1)]) / 4;
      }
    }
    setBnd(0, p);
  }
  // Subtract gradient
  for (let i = 1; i <= N; i++) {
    for (let j = 1; j <= N; j++) {
      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);
}
```

Boundary conditions:
- b=1: u (horizontal velocity) - left/right: u=0? No...
- For velocity: u is stored at cell centers. Boundary: u on left wall (i=0) = -u on i=1 (reflect). v on left wall = v on i=1. etc.

Stam's setBnd:
```javascript
function setBnd(b, x) {
  // Vertical boundaries (i=0, i=N+1)
  for (let j = 1; j <= N; j++) {
    x[IX(0,j)] = b == 1 ? -x[IX(1,j)] : x[IX(1,j)];
    x[IX(N+1,j)] = b == 1 ? -x[IX(N,j)] : x[IX(N,j)];
  }
  // Horizontal boundaries (j=0, j=N+1)
  for (let i = 1; i <= 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)];
  }
  // Corners
  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)]);
}
```

Now for the dye with color. I'll have separate R, G, B dye fields, each with their own advection and diffusion. The mouse injects dye into all three, with the color determined by... let me think.

Approach: Each mouse drag, determine a hue based on the drag direction (angle of velocity). Inject dye into R, G, B channels accordingly. This creates rainbow swirls as you drag in different directions.

Alternatively, map position to hue. But direction-based is more interactive and fun.

Let me use: hue = atan2(dy, dx) → map to RGB.

Actually, let me also consider: to make it really beautiful, I could have the dye color be a function of the velocity field itself (velocity magnitude → brightness, velocity direction → hue). But that's harder to compute since velocity changes.

I'll go with: each injection point gets a color based on the injection velocity direction. I'll maintain 3 dye fields (r, g, b). The injection adds to each channel based on the color.

Let me also add: a slight diffusion of dye so it fades over time (decay). Or rather, the std Stam solver has the dye just advect and diffuse but not decay. I'll add a small decay factor so the dye eventually fades (multiplied by 0.99 each frame or similar).

Actually, for visual clarity, I should make the dye persist long enough to see patterns but fade eventually. Let me add a decay of like 0.995 per frame, or make diffusion handle it.

Hmm, in Stam's solver, the dye field doesn't naturally decay. I should add an explicit decay: `dens *= 0.99` each frame, or include it in the diffuse step.

Let me think about the overall loop:

```
Each frame:
1. Add velocity from mouse movement (if dragging)
2. Add dye from mouse (if dragging) - with color
3. Diffuse velocity (viscous)
4. Advect velocity (self-advection)
5. Project velocity (make divergence-free)
6. Diffuse dye (R, G, B)
7. Advect dye (R, G, B)
8. Apply decay to dye
9. Render
```

Wait, the order in Stam's solver is:
1. Add sources (velocity)
2. Diffuse velocity
3. Project velocity
4. Advect velocity
5. Project velocity (again)

For dye:
1. Add sources (dye)
2. Diffuse dye
3. Advect dye

Hmm, actually the standard Stam loop is:

For velocity:
- Add sources
- Diffuse
- Project
- Advect
- Project

For density (dye):
- Add sources
- Diffuse
- Advect

Let me structure it properly.

Actually, looking at Stam's paper more carefully:
- The steps are: add source, diffuse, advect, project
- For velocity: the projection is done after advection to enforce incompressibility
- The diffusion uses the previous velocity field

Let me use a clean implementation. I'll follow the standard approach.

For performance with 3 dye channels, I could:
- Use 3 separate arrays and process each through the same functions (but with different source arrays)
- Or process them in the same loops

Given that N=100 or 128, the grid has ~10,000-16,000 cells. Three dye fields = manageable. Let me use N=100 for good performance on all devices.

Actually, let me use N = 120 or 150 for nicer visuals, and optimize where possible. Let me go with N=128.

Hmm, but the advect step for 3 dye fields + velocity = 4 advects per frame. With N=128, that's 128*128*4 = 65,536 cells per advect, 4 advects = 262K operations. Plus diffusion (20 iterations * 4 fields * 16K = 1.3M). Plus projection (20 * 2 * 16K = 640K). Total ~2M operations per frame. Should be fine.

Let me go with N=120 for a good balance. Actually let me try N=100 first for safety, can increase.

Wait, I want this to look really good. Let me use N=150 and see. The canvas can be rendered at 150×150 and scaled up. Let me compute: 150*150 = 22,500 cells. With 4 fields (u, v, 3 dye but I'll process dye separately)... 

Actually let me think about this differently. I'll use:
- N = 120 (grid cells per dimension)
- Velocity: u, v (2 fields)
- Dye: r, g, b (3 fields)
- Temporaries: uPrev, vPrev, rPrev, gPrev, bPrev (5 fields)
- Plus: div, p (2 fields for projection)

Total: 12 float arrays of size (N+2)^2 ≈ 12 * 122*122 ≈ 12 * 14,884 ≈ 178,608 floats. That's fine.

For the solver iterations, I'll use fewer iterations for the dye (since it's just visual) and maybe 10-15 for velocity.

Let me also think about rendering. I'll create an ImageData on a canvas of size N×N, fill it with colors from the dye fields, then draw that canvas scaled to full screen.

For color: the dye fields r, g, b have values in [0, 1] ideally. I'll clamp and use as RGB. To make it look nicer, I might apply some gamma or curve.

Actually, for a more vivid look, let me use the dye values directly as RGB but maybe apply a slight boost.

Let me also consider: what if I use velocity to modulate color? Like, show velocity magnitude as brightness and velocity direction as hue? That's a common visualization. But the user asked for "colored dye" so the dye approach is more appropriate.

OK let me also add: when the mouse drags, the color of dye injected depends on the mouse velocity direction. This gives rainbow swirls.

To compute the injection color: 
- Get mouse velocity (dx/dt, dy/dt)
- angle = atan2(dy, dx)
- Map angle to RGB via HSV: hue = (angle + π) / (2π), saturation = 1, value = 1
- Convert HSV to RGB, inject that color

This means dragging up gives one color, dragging right gives another, etc. Beautiful rainbow swirls.

Let me implement HSV to RGB.

Alright, let me write the complete code now.

One more consideration: performance. I should use requestAnimationFrame and measure dt. The simulation dt should be fixed (or clamped) for stability. I'll use a fixed dt of 0.1 or so, and run the simulation at a fixed rate.

Actually, for stability of the advection step, the CFL condition requires dt * N * max_velocity < 1. If max velocity is ~5 cell units per frame and N=120, dt should be < 1/(120*5) ≈ 0.0017. Hmm, that's very small.

Wait, in Stam's formulation, the velocities are in cell units per time step. If I want the fluid to move at a visible rate, I need velocities of maybe 1-5 cells per time unit. With N=120 and dt=0.016 (60fps), CFL would require velocity < 1/(120*0.016) ≈ 0.5 cells per time unit. That's slow.

To get around this, I can:
1. Use a smaller N (e.g., 60)
2. Use a smaller dt but run multiple sub-steps per frame
3. Use the "linear backtrace" with clamping which handles some CFL violation

Actually, in the Stam solver as commonly implemented for games, the velocities are scaled such that dt * N is the effective step. Let me look at how it's typically done.

In many game implementations, they use:
- `dt` as a parameter
- The advect step uses `dt * N` as the scaling factor
- Velocities are in the range of maybe 0-10 or so
- N is the grid size

The CFL condition in this formulation: `dt * N * |v| < 1` where v is velocity in cell units. If dt=0.016, N=120, then |v| < 0.52 for stability. 

So if I want visible motion, I should either:
- Use a lower N (e.g., 60-80)
- Or use a smaller effective dt (sub-stepping)
- Or clamp the backtrace coordinates (which the Stam advect already does by clamping to [0.5, N+0.5])

The clamping helps with stability but can cause issues if velocities are too high.

Let me use N=80 for good visual quality and performance. With N=80:
- CFL: |v| < 1/(80*0.016) ≈ 0.78. So velocities should be < 0.78 cell units per frame.
- I'll normalize mouse input so that the added velocity is in a reasonable range.

Actually, you know what, let me just use a moderately sized grid and tune the parameters. Let me use N=96.

Or actually, let me think about this more carefully. In many real-time fluid demos, they use N around 80-120 and it works fine. The key is to scale the mouse input velocity appropriately.

Let me set up the parameters so that:
- Mouse velocity of ~100 pixels/second maps to a fluid velocity of ~0.5 cell units per frame (stable)
- The fluid advects dye nicely

I'll compute: if mouse moves 100 px/s, and canvas is ~1000px wide, that's 0.1 canvas-widths per second. In grid cells (N=96), that's 9.6 cells per second, or 0.16 cells per frame at 60fps. That's well within CFL.

So I'll scale mouse movement to grid coordinates and add velocity = (grid_delta_x / dt_in_seconds). With dt=0.016 and grid_delta=1 cell, velocity = 62.5 cells/frame. That's way too high!

I need to be careful. Let me think:
- Mouse movement in grid cells per second: if I want subtle but visible fluid motion, maybe 10-50 cells/sec is good.
- At 60fps, that's 0.17-0.83 cells per frame.
- Velocity in Stam's formulation = cells per frame (since dt is per frame effectively).

Wait, actually in the Stam solver as typically implemented for games, `dt` is set to 1.0 or some value, and the velocity is in "cells per time unit". Let me re-examine.

In the original Stam paper and many implementations:
- The domain is [0, 1]×[0, 1] or [0, N]×[0, N]
- Velocity is in units of "cells per time step" or "domain per time step"
- dt is the time step

In the implementation I outlined above, I used `dt * N` as the scaling in advect. This means:
- Velocity u is in "cells per unit time"
- The backtrace: x = i - dt * N * u[IX(i,j)]
- For stability: dt * N * |u| < 1

If I set dt = 1.0 (per frame), then |u| < 1/N for stability. That's very restrictive.

If I set dt = 0.016 (real seconds per frame), then |u| < 1/(N*0.016) = 3.125/N for N=200... wait, this doesn't seem right either.

Let me reconsider. The standard game implementation (like the one often attributed to Mike Ash or various GitHub projects) uses:

```javascript
function advect(b, d, d0, u, v, dt) {
  for (let i = 1; i <= N; i++) {
    for (let j = 1; j <= N; j++) {
      let x = i - dt * u[IX(i,j)];
      let y = j - dt * v[IX(i,j)];
      ...
    }
  }
}
```

Where velocities are in cell units and dt is in frames (or a scaled value). In this case, CFL is dt * |v| < 1. With dt=1 (per frame), |v| < 1 cell per frame. That's slow but stable.

Many implementations use dt as a parameter they tune, often around 0.1-0.5, and velocities around 0-5. 

Let me use a different scaling: I'll keep dt as a parameter (say 0.1-0.5) and scale velocities accordingly. Mouse input will be converted to velocity in the appropriate range.

Actually, I think the cleanest approach is:

1. Use grid coordinates where the domain is [0, N] × [0, N]
2. Velocities are in cells per second
3. dt is in seconds (real frame time, clamped to e.g. 1/30)
4. In advect: `x = i - dt * u` (velocity is already in cells/sec)
5. CFL: dt * |u| < 1 → |u| < 1/dt ≈ 30-60 cells/sec

So if I set mouse velocity to add up to ~20 cells/sec, that's stable and visible.

For mouse input: I track mouse position in grid coords. On each frame, I compute delta in grid coords and divide by dt to get velocity in cells/sec. I add this to the velocity field at the mouse position.

But I need to be careful about the magnitude. If the user moves the mouse fast (say 1000 px/s on a 1000px canvas = 1 canvas width/sec = N cells/sec = 96 cells/sec for N=96), that's too fast for CFL.

I should clamp the input velocity to some maximum (e.g., 30 cells/sec) and also smooth it.

Let me implement this:

```javascript
// Mouse handling
let mouseX = -1, mouseY = -1; // in grid coords (0 to N)
let prevMouseX = -1, prevMouseY = -1;
let isDragging = false;

// On mousemove while dragging:
//   mouseX, mouseY = current position in grid coords
//   dx = mouseX - prevMouseX (in grid cells)
//   dy = mouseY - prevMouseY
//   vx = dx / dt (cells per second)
//   vy = dy / dt
//   clamp vx, vy to maxVelocity
//   Add velocity to u, v at mouse grid position (with some radius)
//   Add dye at mouse position with color from velocity angle
//   prevMouseX = mouseX, prevMouseY = mouseY
```

For adding to the field, I need to add to the discrete grid cells around the mouse position. I'll use bilinear interpolation or just add to the nearest cell / a small radius.

Let me use a simple approach: add to the cell nearest the mouse, and also to neighboring cells with a Gaussian falloff.

Actually, for simplicity and good visual results, I'll add to a small radius (like 2-3 cells) with a Gaussian weight.

OK, I think I've planned enough. Let me write the code.

For the visual rendering, I want it to look gorgeous. Let me think about the color mapping:

Option A: Velocity-direction-based injection color (rainbow swirls)
Option B: Position-based injection color (spatial color map)
Option C: Velocity magnitude → brightness, with a base hue

I'll go with Option A primarily, but also add a slight option B for variety. Actually, let me just do A with a nice HSV mapping.

For rendering the dye field: r, g, b values → pixel color. I'll apply a slight gamma (gamma < 1 for brighter midtones) and maybe a contrast boost.

Let me also add: the dye slowly fades (multiply by 0.998 per frame or so) so the screen doesn't get too full.

Alright, here's my plan for the complete index.html:

```html
<!DOCTYPE html>
<html>
<head>
<style>
  body { margin: 0; overflow: hidden; background: #000; }
  canvas { display: block; width: 100vw; height: 100vh; }
  #info { ... }
</style>
</head>
<body>
<canvas id="canvas"></canvas>
<script>
  // Fluid simulation code
</script>
</body>
</html>
```

Let me write the full implementation now. I'll be thorough.

Let me think about the grid size once more. For a full-screen experience on various devices:
- N=80: 80*80=6400 cells, fast on all devices, decent visual
- N=120: 14400 cells, good on desktop, might be slow on mobile
- N=100: 10000 cells, good balance

I'll use N=100 and allow the simulation to run at 60fps. If performance is an issue, I can always lower it.

Actually, let me use N=120 for nicer visuals and add a performance fallback... no, let me keep it simple. N=100 is good.

Hmm, actually let me use a slightly different approach for better visuals: N=96 (divisible by many things) or just N=100. Let me go with N=100.

Wait, I realize I should think about the canvas sizing. The simulation grid is N×N. I'll render to an offscreen canvas of size N×N, then draw it scaled to the full-screen canvas. The CSS will make the canvas fill the viewport.

For the offscreen canvas: I'll use ImageData to set pixels directly, which is fast.

For scaling: I'll use `drawImage` with the offscreen canvas to the main canvas, scaled to fill. I'll use `imageSmoothingEnabled = false` for a crisp pixel look, or `true` for a blurred/smooth look. I think smooth (true) looks better for fluid simulations.

Actually, for a fluid sim, smoothing (bilinear interpolation when scaling) looks better because it gives a fluid, organic look rather than pixelated blocks.

Let me finalize the implementation:

```javascript
const N = 100;
const cellSize = 1; // grid is N×N, rendered to N×N offscreen then scaled

// Allocate fields
const size = (N+2)*(N+2);
const u = new Float32Array(size);
const v = new Float32Array(size);
const uPrev = new Float32Array(size);
const vPrev = new Float32Array(size);
const r = new Float32Array(size);
const g = new Float32Array(size);
const b = new Float32Array(size);
const rPrev = new Float32Array(size);
const gPrev = new Float32Array(size);
const bPrev = new Float32Array(size);
const div = new Float32Array(size);
const p = new Float32Array(size);

// Solver functions
function IX(i, j) { return i + (N+2)*j; }
function setBnd(b, x) { ... }
function diffuse(b, x, x0, diff, dt) { ... }
function advect(b, d, d0, u, v, dt) { ... }
function project(u, v, p, div) { ... }

// Simulation step
function step(dt) {
  // Velocity
  addVelocitySource(dt);
  diffuse(1, u, uPrev, visc, dt);
  diffuse(2, v, vPrev, visc, dt);
  project(u, v, p, div);
  swap(u, uPrev);
  swap(v, vPrev);
  advect(1, u, uPrev, uPrev, vPrev, dt);
  advect(2, v, vPrev, uPrev, vPrev, dt);
  project(u, v, p, div);
  
  // Dye (R, G, B)
  addDyeSource(dt);
  diffuse(0, r, rPrev, diff, dt);
  diffuse(0, g, gPrev, diff, dt);
  diffuse(0, b, bPrev, diff, dt);
  swap(r, rPrev); swap(g, gPrev); swap(b, bPrev);
  advect(0, r, rPrev, u, v, dt);
  advect(0, g, gPrev, u, v, dt);
  advect(0, b, bPrev, u, v, dt);
  // Decay
  for (let i = 0; i < size; i++) {
    r[i] *= 0.997; g[i] *= 0.997; b[i] *= 0.997;
  }
}

// Rendering
function render() {
  const imageData = offscreenCtx.createImageData(N, N);
  const data = imageData.data;
  for (let i = 1; i <= N; i++) {
    for (let j = 1; j <= N; j++) {
      const idx = (i-1 + (N)*(j-1)) * 4;
      const ri = r[IX(i,j)];
      const gi = g[IX(i,j)];
      const bi = b[IX(i,j)];
      data[idx] = Math.min(255, ri * 255);
      data[idx+1] = Math.min(255, gi * 255);
      data[idx+2] = Math.min(255, bi * 255);
      data[idx+3] = 255;
    }
  }
  offscreenCtx.putImageData(imageData, 0, 0);
  mainCtx.imageSmoothingEnabled = true;
  mainCtx.drawImage(offscreenCanvas, 0, 0, canvasWidth, canvasHeight);
}
```

Wait, I need to be more careful with the advect function. In Stam's solver, the advect step uses the velocity field to backtrace. The velocity field used should be the one from the previous step (before advection). Let me re-check the standard implementation.

Standard Stam velocity update:
```
1. Add sources to velocity (u0, v0) ← (u, v) + sources
2. Diffuse: u ← diffuse(u0), v ← diffuse(v0)  [actually: u = diffuse(u0), etc.]
3. Project: project(u, v)
4. Advect: u1 = advect(u, u0, u, v), v1 = advect(v, v0, u, v)  [using pre-advect velocity]
5. Project: project(u1, v1)
6. Swap: u = u1, v = v1
```

Wait, the advect in Stam's code is:
```c
advect(1, u, u0, u0, v0, dt);  // b=1 for u
advect(2, v, v0, u0, v0, dt);  // b=2 for v
```

So it advects u using u0, v0 (the pre-advect velocity) and stores in u. Then projects.

For the dye:
```c
advect(0, dens, dens0, u, v, dt);  // using current velocity
```

OK so the dye advection uses the current (post-projection) velocity field. Good.

Let me restructure my step function:

```
// 1. Save current velocity to prev
uPrev = u; vPrev = v;  // (via swap or copy)
// 2. Add sources to u, v (now u, v have sources added)
addVelocitySource();  // adds to u, v (using prev as base? or just adds to current?)
```

Hmm, I need to be careful. Let me follow Stam's structure exactly:

```c
// In Stam's code:
void FluidSim::step(float dt) {
  // Add velocity sources
  addSource(u, uSource);  // u += uSource
  addSource(v, vSource);  // v += vSource
  
  // Diffuse
  diffuse(1, u, uPrev, visc, dt);  // u = diffuse(uPrev + sources, ...)
  // Wait, this isn't right either.
}
```

Let me look at the actual Stam implementation more carefully.

From Stam's paper "Real-Time Fluid Dynamics for Games":

```
void Fluid::step(float dt) {
  // 1. Add forces (sources)
  addSource(u, uForce);
  addSource(v, vForce);
  
  // 2. Diffuse
  diffuse(1, u, u0, visc, dt);
  diffuse(2, v, v0, visc, dt);
  
  // 3. Project
  project(u, v, p, div);
  
  // 4. Advect
  advect(1, u, u0, u0, v0, dt);
  advect(2, v, v0, u0, v0, dt);
  
  // 5. Project
  project(u, v, p, div);
}
```

Where `u0` and `v0` are the "previous" velocity fields. The pattern is:
- `addSource(u, force)` adds to `u` (which is the current field)
- `diffuse(b, x, x0, diff, dt)` solves `x = x0 + diff*dt*laplacian(x)`, where x0 is the source. So x0 should be the field before diffusion.
- `advect(b, d, d0, u, v, dt)` advects d0 using u,v and stores in d.

So the flow is:
1. `u` starts as the current velocity
2. Add sources to `u` (now `u` has sources)
3. `diffuse(1, u, u0, visc, dt)` — but wait, what's u0 here? It should be the velocity BEFORE adding sources and BEFORE diffusion. So u0 should be a copy of u from the start of the step.

I think the standard approach is:
- At the start of step(), u and v are the current velocity
- Save them to u0, v0 (the "previous" fields)
- Add sources to u, v (so u, v now have sources added to last frame's velocity)
- Diffuse: solve for u using u0+sources as the source term... 

Actually, I think the confusion is that Stam uses u0 as a temporary that gets overwritten. Let me look at the diffuse function:

```c
void diffuse(int b, float *x, float *x0, float diff, float dt) {
  float a = dt * diff * N * N;
  // Gauss-Seidel iteration: x = (x0 + a*neighbors) / (1+4a)
  // x0 is the source (what we're diffusing from)
  // x is the result
}
```

So `diffuse(b, x, x0, diff, dt)` computes `x` as the diffusion of `x0`. After calling, `x` contains the diffused field.

So the flow is:
1. `u0 = u` (copy current velocity to u0) — but actually, we need u0 to be the velocity FIELD before sources are added.
2. Add sources to `u` (u = u + sources)
3. `diffuse(1, u, u0, visc, dt)` — diffuses u0 (pre-source velocity) and stores in u. But this loses the sources!

Hmm, that doesn't seem right. Let me look at actual implementations.

In many implementations I've seen:
```javascript
// Step:
addSource(u, uSource);  // u += source
addSource(v, vSource);  // v += source
swap(u, uPrev);  // uPrev = old u (before sources)
swap(v, vPrev);  // vPrev = old v (before sources)
diffuse(1, u, uPrev, visc, dt);  // u = diffuse(uPrev) — but uPrev is pre-source!
```

Wait, this diffuses the PRE-source velocity, which means the sources are lost. That can't be right.

Oh, I see. In some implementations, the addSource happens AFTER the swap, and the source is added to uPrev (which becomes the source term for diffusion). Or the source is added to u after diffusion.

Let me look at a well-known implementation. Here's the structure from a typical JS implementation:

```javascript
function step(dt) {
  // Add sources
  addSource(u, uSrc);
  addSource(v, vSrc);
  addSource(dens, densSrc);
  
  // Swap
  swap(u, uPrev);
  swap(v, vPrev);
  swap(dens, densPrev);
  
  // Diffuse
  diffuse(1, u, uPrev, visc, dt);
  diffuse(2, v, vPrev, visc, dt);
  diffuse(0, dens, densPrev, diff, dt);
  
  // Project
  project(u, v, p, div);
  
  // Advect
  advect(1, u, uPrev, uPrev, vPrev, dt);
  advect(2, v, vPrev, uPrev, vPrev, dt);
  advect(0, dens, densPrev, u, v, dt);
  
  // Project again
  project(u, v, p, div);
}
```

In this structure:
- `addSource` adds to `u` (current field)
- `swap(u, uPrev)` saves the current (with sources) to uPrev, and u becomes the old uPrev
- `diffuse(1, u, uPrev, visc, dt)` — uPrev has the sources, and u gets the diffused result. 

Wait, that means uPrev (which is the old u + sources) is the source term for diffusion, and u gets the result. So the sources ARE included in the diffusion. 

But then `swap(u, uPrev)` at the start of the next step... let me trace through:

Frame 1:
- Initial: u = field_A, uPrev = field_0 (zero or whatever)
- addSource(u, src) → u = field_A + src
- swap(u, uPrev) → u = field_0, uPrev = field_A + src
- diffuse(1, u, uPrev, visc, dt) → u = diffuse(field_A + src) = field_A diffused + src diffused
- (continue with project, advect, etc.)

Frame 2:
- u = result from frame 1
- uPrev = field_A + src (from frame 1 swap)
- addSource(u, src2) → u = result1 + src2
- swap(u, uPrev) → u = field_A + src (old uPrev), uPrev = result1 + src2
- diffuse(1, u, uPrev, visc, dt) → u = diffuse(result1 + src2)

Hmm, this seems to work. The uPrev from the previous frame (which had sources added) becomes the source for diffusion in the current frame. And the result of diffusion goes into u.

But wait, there's an issue: in frame 1, uPrev = field_0 (initial, maybe zero). After addSource and swap, uPrev = field_A + src, u = field_0. Then diffuse(1, u, uPrev, ...) computes u = diffuse(field_A + src). This is correct — we're diffusing the current velocity (with sources) from the previous state.

Hmm, actually I think there might be a subtle issue. Let me think again...

In Stam's original code (C), the pattern is:

```c
// Step function
void step() {
  // Add sources
  addSource(u, uSrc);
  addSource(v, vSrc);
  
  // Swap u <-> u0, v <-> v0
  // Now u0 = old u + sources, u = old u0 (previous step's u0)
  
  // Diffuse: u = diffuse(u0), v = diffuse(v0)
  diffuse(1, u, u0, visc, dt);
  diffuse(2, v, v0, visc, dt);
  
  // Project
  project(u, v, p, div);
  
  // Advect: using u, v as velocity, u0, v0 as source
  // u = advect(u0, u, v), v = advect(v0, u, v)
  advect(1, u, u0, u, v, dt);
  advect(2, v, v0, u, v, dt);
  
  // Project
  project(u, v, p, div);
}
```

Wait, in the advect call, the velocity used is `u` and `v` (the current fields after diffusion and projection). And the source being advected is `u0` and `v0` (which are the pre-diffusion, post-source fields). The result goes into `u` and `v`.

Hmm, but that means after advection, u contains the advected version of u0 (pre-diffusion velocity). That loses the diffusion effect!

I think the actual Stam implementation is different. Let me look at the actual code from the paper.

OK I found the actual Stam pseudocode. Here it is:

```
void FluidSim::step(float dt) {
  // Add forces
  addSource(u, uSrc);
  addSource(v, vSrc);
  
  // Diffuse
  diffuse(1, u, uPrev, visc, dt);
  diffuse(2, v, vPrev, visc, dt);
  
  // Project
  project(u, v, p, div);
  
  // Advect
  advect(1, u, uPrev, uPrev, vPrev, dt);
  advect(2, v, vPrev, uPrev, vPrev, dt);
  
  // Project
  project(u, v, p, div);
  
  // Swap u <-> uPrev, v <-> vPrev
  swap(u, uPrev);
  swap(v, vPrev);
}
```

Wait, this swaps at the END. So:

Frame 1:
- u = initial velocity, uPrev = some initial state (maybe zero)
- addSource: u += src (u now has sources)
- diffuse(1, u, uPrev, visc, dt): u = diffuse(uPrev) — but uPrev is the initial state, not u+sources!

This doesn't include the sources in the diffusion. That seems wrong.

Hmm, unless... the addSource happens AFTER the swap in the previous step. Let me re-trace:

Let's say at the start of frame 1, u = init_u, uPrev = zero.

If the swap happens at the END of frame 0 (or initialization):
- After swap: u = zero, uPrev = init_u

Frame 1:
- addSource(u, src1): u = zero + src1 = src1
- diffuse(1, u, uPrev, visc, dt): u = diffuse(init_u) — diffuse the initial velocity, result in u
- project(u, v, ...)
- advect(1, u, uPrev, uPrev, vPrev, dt): u = advect(init_u using u,v velocity)
- project(u, v, ...)
- swap(u, uPrev): u = (advected init_u), uPrev = src1

Frame 2:
- addSource(u, src2): u = advected_init_u + src2
- diffuse(1, u, uPrev, visc, dt): u = diffuse(src1) — diffuse the source from frame 1!
- ...

This is wrong! The velocity from frame 1 (advected init_u) is not being diffused; instead, src1 (a tiny source) is being diffused.

I think I'm confusing myself. Let me look at actual working implementations.

Let me look at the implementation from https://github.com/h-schellander/Fluid-Simulation or similar.

Actually, I recall now. The issue is that in Stam's implementation, the "source" added via addSource is the FORCE (like wind, mouse input), and it gets added to the velocity field. The diffusion then smooths the COMBINED field (old velocity + force). The advect then moves this combined field.

The correct flow (which I've seen in working implementations) is:

```javascript
function step(dt) {
  // 1. Add force sources to current velocity
  addSource(u, uSrc);
  addSource(v, vSrc);
  
  // 2. Diffuse: solve for new velocity = diffusion of (old_velocity + sources)
  //    We need old_velocity + sources as the source term
  //    Since we added sources to u, and swapped so uPrev = old_u:
  swap(u, uPrev);  // uPrev = old_u + sources, u = old_uPrev (zero or old)
  diffuse(1, u, uPrev, visc, dt);  // u = diffuse(old_u + sources)
  
  // 3. Project
  project(u, v, p, div);
  
  // 4. Advect: backtrace using current velocity (u,v), advect uPrev (which = old_u + sources)
  swap(u, uPrev);  // uPrev = current u (diffused+projected), u = old_u + sources
  advect(1, u, uPrev, u, v, dt);  // u = advect(old_u + sources using current velocity)
  // Wait, this advects the old velocity field, not the current one.
  
  // Hmm, this is confusing.
}
```

OK I think I need to just look at a concrete working implementation. Let me recall the one from the popular "Fluid Simulation" by Pavel Doulin or the one on GitHub.

Actually, let me just go with a well-tested implementation pattern. Here's one that I know works (from various sources):

```javascript
function step(dt) {
  const N2 = N;
  
  // Add sources (forces) to velocity
  addSource(u, uSrc);
  addSource(v, vSrc);
  
  // Swap so that uPrev holds the current velocity (with sources)
  // and u holds the previous uPrev (which will be overwritten)
  swapField(u, uPrev);
  swapField(v, vPrev);
  
  // Diffuse velocity
  diffuse(1, u, uPrev, visc, dt);
  diffuse(2, v, vPrev, visc, dt);
  
  // Project (make divergence-free)
  project(u, v, p, div);
  
  // Swap again: uPrev = velocity for advection source, u = result
  swapField(u, uPrev);
  swapField(v, vPrev);
  
  // Advect velocity
  advect(1, u, uPrev, uPrev, vPrev, dt);
  advect(2, v, vPrev, uPrev, vPrev, dt);
  
  // Project again
  project(u, v, p, div);
}
```

Let me trace this:

Initial state (frame 0 done, or initialization):
- u = current velocity field
- uPrev = some buffer (don't care)

Frame 1:
1. addSource(u, src1): u = init_u + src1
2. swap(u, uPrev): uPrev = init_u + src1, u = buffer
3. diffuse(1, u, uPrev, visc, dt): u = diffuse(init_u + src1) ✓
4. project(u, v, p, div): u, v are now divergence-free ✓
5. swap(u, uPrev): uPrev = diffuse(init_u + src1) (divergence-free), u = init_u + src1
6. advect(1, u, uPrev, uPrev, vPrev, dt): u = advect(diffused_velocity, using diffused_velocity as the advection field)
   - Wait, this advects the PREVIOUS velocity field (uPrev after step 5 = diffused velocity) using the CURRENT velocity (uPrev after step 5 = diffused velocity, passed as the velocity field parameter).
   - Actually, uPrev after step 5 is the diffused+projected velocity. And in the advect call, uPrev is passed as both the source field AND the velocity field (for both u and v components). So we're advecting the velocity field using itself. That's self-advection. ✓
7. project(u, v, p, div): final projection ✓
8. End of frame: u = advected, divergence-free velocity. uPrev = init_u + src1 (old buffer).

Frame 2:
1. addSource(u, src2): u = advected_velocity_frame1 + src2
2. swap(u, uPrev): uPrev = advected_velocity_frame1 + src2, u = init_u + src1 (old)
3. diffuse(1, u, uPrev, visc, dt): u = diffuse(advected_velocity_frame1 + src2) ✓
4. ... and so on.

This looks correct! The velocity field is properly diffused, projected, advected, and projected each frame. The sources are included in the diffusion step.

Now for the dye, the pattern is simpler:
```javascript
// Dye
addSource(r, rSrc);
addSource(g, gSrc);
addSource(b, bSrc);
swapField(r, rPrev);
swapField(g, gPrev);
swapField(b, bPrev);
diffuse(0, r, rPrev, diff, dt);
diffuse(0, g, gPrev, diff, dt);
diffuse(0, b, bPrev, diff, dt);
// No projection for dye
swapField(r, rPrev);
swapField(g, gPrev);
swapField(b, bPrev);
advect(0, r, rPrev, u, v, dt);
advect(0, g, gPrev, u, v, dt);
advect(0, b, bPrev, u, v, dt);
// Decay
for (...) { r[i] *= decay; g[i] *= decay; b[i] *= decay; }
```

Wait, but there's a question: should the dye advection use the velocity BEFORE or AFTER the velocity step? In Stam's implementation, the dye is advected using the current velocity field (after the velocity step). So I should do the velocity step first, then the dye step.

Actually, looking at Stam's code, he processes velocity and density in the same step function, with velocity processed first. The dye advection uses the final velocity field (after projection).

So the order is:
1. Velocity: add source → diffuse → project → advect → project
2. Dye: add source → diffuse → advect → decay

Let me now also think about the dye source. When the mouse drags, I add dye at the mouse position. The dye color should be based on the mouse velocity direction (for rainbow swirls).

For the mouse velocity:
- Track mouse position in grid coords: (mx, my) where mx, my ∈ [0, N]
- Previous position: (pmx, pmy)
- Delta: (dx, dy) = (mx - pmx, my - pmy)
- dt is the frame time in seconds
- Velocity: (vx, vy) = (dx/dt, dy/dt) in cells per second
- Clamp to maxVel (e.g., 100 cells/sec)
- Color: hue = atan2(vy, vx) mapped to [0, 1], convert HSV to RGB

For adding dye to the grid:
- Add to cell at (mx, my) and neighbors with a Gaussian radius
- The amount added depends on the mouse speed and a source strength parameter

For adding velocity:
- Add (vx, vy) to the velocity field at (mx, my) and neighbors
- This creates the push effect

Let me also consider: to make the simulation more interesting, I could add a small ambient velocity or noise. But let's keep it simple — the user's mouse creates all the motion.

One more thing: I should handle the case where the mouse is dragged quickly. The velocity should be clamped to avoid CFL violations. And the dye injection should be proportional to the drag speed (faster drag = more dye and stronger push).

Let me now write the full code.

Oh wait, one more consideration about the `advect` function signature. In my implementation:

```javascript
function advect(b, d, d0, u, v, dt) {
  // d: output field
  // d0: source field (to be advected)
  // u, v: velocity fields
  // dt: time step
  // b: boundary type (0 for scalar, 1 for u, 2 for v)
}
```

The backtrace: `x = i - dt * N * u[IX(i,j)]` — wait, should it be `dt * N` or just `dt`?

In Stam's original code, the domain is [0, N] × [0, N], and the velocity is in "cells per unit time". The backtrace is:
```
x = i - dt * u[IX(i,j)]
```
where dt is in time units and u is in cells per time unit. So if dt=0.016 (1/60) and u=10 cells/sec, then x = i - 0.16, which is fine.

But in many implementations, they use `dt * N` as the scaling. This is when the velocity is normalized to [0, 1] range (fraction of domain per time unit). Then `dt * N * u` gives the backtrace distance in cells.

I'll use the former approach: velocity in cells per second, dt in seconds, backtrace = `i - dt * u[IX(i,j)]`. This is cleaner.

For the CFL condition: `dt * |u| < 1` (in cells). With dt=0.016, |u| < 62.5 cells/sec. I'll clamp mouse velocity to maybe 50 cells/sec to be safe.

Wait, but I also need to consider that the velocity field gets amplified by the advect step (self-advection can create high velocities). The projection step helps, but I should also clamp velocities or use a small enough dt.

For safety, I'll clamp velocities to a maximum (e.g., 50 cells/sec) after each step, and use dt clamped to a maximum (e.g., 1/30).

Let me also consider the diffusion coefficient. For the velocity, visc (viscosity) should be small (e.g., 0.0001) so the fluid is mostly inviscid (like water/air). For the dye, diff should be small too (e.g., 0.00001) so the dye doesn't blur too much.

Actually, for nice visual results, a small amount of dye diffusion creates a pleasing blur. And a small amount of viscosity smooths the velocity field. Let me use:
- visc = 0.00005 (velocity diffusion)
- diff = 0.00001 (dye diffusion)
- decay = 0.998 (dye decay per frame)

These are tunable. Let me start with these and adjust.

For the Gauss-Seidel iterations, I'll use 10-20 iterations for velocity and 10 for dye. More iterations = better accuracy but slower. 20 is a good balance.

For the pressure solver (in project), I'll use 20 iterations.

Alright, let me write the complete code now. I'll make it clean and well-organized.

One more thing about rendering: I want the background to be dark (black or very dark) so the colored dye pops. I'll render the dye as RGB on a black background. Pixels with zero dye will be black.

For the color mapping: r, g, b values are in [0, 1] (or higher if multiple sources overlap). I'll clamp to [0, 1] and multiply by 255. To make it look nicer, I might apply a slight gamma (e.g., sqrt or pow 0.8) to brighten the midtones.

Actually, for the dye, values can exceed 1.0 if multiple injections overlap. I should clamp to [0, 1] for rendering. The decay will bring high values down over time.

Let me also think about the initial state: empty (all zeros) is fine. The user starts dragging and the dye appears.

For a nice touch, I could pre-seed some dye or velocity, but let's keep it clean — the user creates everything by dragging.

Alright, coding time. Let me write the full index.html.

Actually, let me also think about touch support for mobile. I'll add touch event handlers too.

Let me also add a subtle UI hint: "Drag to create fluid" text that fades out after a few seconds, or a small instruction.

OK, let me write the code now. I'll be comprehensive.

Let me think about one more thing: the `setBnd` function. For the dye (scalar field, b=0), the boundary condition is reflexive (Neumann): the value at the boundary equals the value at the adjacent interior cell. For velocity (b=1 for u, b=2 for v), the boundary condition is:
- u (horizontal velocity): on vertical walls (i=0, i=N+1), u = -u (reflect, meaning no flow through the wall). On horizontal walls (j=0, j=N+1), u = u (free slip).
- v (vertical velocity): on horizontal walls, v = -v. On vertical walls, v = v.

This is what Stam's setBnd does with the `b` parameter.

For the dye, b=0 means all boundaries are reflective (copy from interior).

OK, I'm confident in my understanding now. Let me write the code.

For the rendering, I'll create an offscreen canvas of size N×N, put the image data, then draw it to the main canvas scaled to fullscreen. I'll use `imageSmoothingEnabled = true` for a smooth, fluid look.

Let me also think about the canvas sizing. I'll make the main canvas fill the window. On resize, I'll update the canvas dimensions. The offscreen canvas stays at N×N.

For high-DPI displays, I could use devicePixelRatio, but for simplicity and performance, I'll keep it at 1:1 (CSS pixels = canvas pixels). The fluid sim resolution is independent of display resolution anyway.

Actually, for crisp rendering on high-DPI screens, I should set the canvas width/height to window.devicePixelRatio * window.innerWidth/Height, and scale the context. But this adds complexity. For a fluid sim, it's fine to render at CSS pixel resolution — the fluid will look smooth either way.

Let me use `canvas.width = window.innerWidth` and `canvas.height = window.innerHeight`, and CSS `width: 100vw; height: 100vh`.

Hmm, but for high-DPI, this will look blurry. Let me handle DPR:

```javascript
function resize() {
  const dpr = window.devicePixelRatio || 1;
  canvas.width = window.innerWidth * dpr;
  canvas.height = window.innerHeight * dpr;
  canvas.style.width = window.innerWidth + 'px';
  canvas.style.height = window.innerHeight + 'px';
  ctx.scale(dpr, dpr);
}
```

Actually, for this fluid sim, the rendering is just drawing a scaled image. The DPR scaling doesn't affect the fluid resolution. So I can keep it simple:

```javascript
canvas.width = window.innerWidth;
canvas.height = window.innerHeight;
```

And the offscreen canvas is N×N. The drawImage scales it. On high-DPI, the canvas is physically larger (more pixels) but the CSS size is the same, so it might look slightly blurry. But for a fluid sim, that's actually fine — it gives a smooth look.

Hmm, actually I want it to look good. Let me use DPR scaling:

```javascript
const dpr = window.devicePixelRatio || 1;
canvas.width = innerWidth * dpr;
canvas.height = innerHeight * dpr;
canvas.style.width = innerWidth + 'px';
canvas.style.height = innerHeight + 'px';
ctx.scale(dpr, dpr);
// Now drawing at (0,0, innerWidth, innerHeight) fills the canvas
```

But then the offscreen canvas is N×N, and I draw it to (0,0, innerWidth, innerHeight) in CSS pixels, which gets scaled by DPR. This should look fine.

Wait, actually, if I scale the context by DPR, then all drawing operations are in CSS pixels. So `drawImage(offscreen, 0, 0, innerWidth, innerHeight)` draws the offscreen canvas scaled to fill the viewport in CSS pixels, and the DPR scaling handles the rest.

But I need to be careful: after `ctx.scale(dpr, dpr)`, the coordinate system is in CSS pixels. So the canvas's logical size (for drawing) is innerWidth × innerHeight. The physical size is innerWidth*dpr × innerHeight*dpr.

OK this is getting complicated. Let me simplify: I'll just set the canvas to window dimensions and not worry about DPR. The fluid sim will look fine.

Actually, you know what, let me just set the canvas width/height to the window's innerWidth/innerHeight (in CSS pixels), and not use DPR. The image will be drawn at this resolution. On a 1000px wide screen, the canvas is 1000px wide. The offscreen canvas (N×N = 100×100) is scaled up to 1000×1000. With imageSmoothingEnabled = true, it'll look smooth.

On a Retina display (2x DPR), the canvas element is 1000×1000 physical pixels (if I set width=1000), but the screen is 2000×2000 physical pixels. So it'll be slightly blurry. But honestly, for a fluid sim, this is fine and even desirable (smooth look).

Let me go with the simple approach: canvas.width = innerWidth, canvas.height = innerHeight.

OK, NOW let me actually write the code. I've been planning for too long.

Let me write a clean, complete implementation:

```html
<!DOCTYPE html>
<html lang="en">
<head>
<meta charset="UTF-8">
<meta name="viewport" content="width=device-width, initial-scale=1.0">
<title>Fluid Simulation</title>
<style>
  * { margin: 0; padding: 0; box-sizing: border-box; }
  html, body { width: 100%; height: 100%; overflow: hidden; background: #0a0a0f; }
  canvas { display: block; width: 100vw; height: 100vh; cursor: crosshair; }
  #hint {
    position: fixed;
    bottom: 30px;
    left: 50%;
    transform: translateX(-50%);
    color: rgba(255,255,255,0.4);
    font-family: system-ui, sans-serif;
    font-size: 14px;
    pointer-events: none;
    transition: opacity 2s;
    text-align: center;
    letter-spacing: 0.5px;
  }
</style>
</head>
<body>
<canvas id="c"></canvas>
<div id="hint">Drag to push colorful dye through the fluid</div>
<script>
// ... full implementation
</script>
</body>
</html>
```

Now for the JS:

```javascript
const canvas = document.getElementById('c');
const ctx = canvas.getContext('2d');
const offscreen = document.createElement('canvas');
const offCtx = offscreen.getContext('2d');

// Grid
const N = 120;
const cellSize = 1;
const size = (N + 2) * (N + 2);

// Fields
const u = new Float32Array(size);
const v = new Float32Array(size);
const uPrev = new Float32Array(size);
const vPrev = new Float32Array(size);
const r = new Float32Array(size);
const g = new Float32Array(size);
const b = new Float32Array(size);
const rPrev = new Float32Array(size);
const gPrev = new Float32Array(size);
const bPrev = new Float32Array(size);
const div = new Float32Array(size);
const p = new Float32Array(size);

// Parameters
const dt = 0.016; // ~60fps time step (will be clamped)
const visc = 0.0001;
const diff = 0.00002;
const decay = 0.997;
const iter = 15;
const maxVel = 80; // max velocity in cells/sec

// Mouse state
let mouseX = -1, mouseY = -1;
let pmouseX = -1, pmouseY = -1;
let isDown = false;
let lastTime = 0;

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

function setBnd(b, x) {
  for (let j = 1; j <= N; j++) {
    x[IX(0, j)] = b === 1 ? -x[IX(1, j)] : x[IX(1, j)];
    x[IX(N + 1, j)] = b === 1 ? -x[IX(N, j)] : x[IX(N, j)];
  }
  for (let i = 1; i <= 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 diffuse(b, x, x0, diff, dt) {
  const a = dt * diff * N * N;
  for (let k = 0; k < iter; k++) {
    for (let i = 1; i <= N; i++) {
      for (let j = 1; j <= N; 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)]
        )) / (1 + 4 * a);
      }
    }
    setBnd(b, x);
  }
}

function advect(b, d, d0, u, v, dt) {
  const dtN = dt * N;
  for (let i = 1; i <= N; i++) {
    for (let j = 1; j <= N; j++) {
      let x = i - dtN * u[IX(i, j)];
      let y = j - dtN * v[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 = Math.floor(x);
      const i1 = i0 + 1;
      const j0 = Math.floor(y);
      const j1 = j0 + 1;
      const s1 = x - i0;
      const s0 = 1 - s1;
      const t1 = y - j0;
      const 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 i = 1; i <= N; i++) {
    for (let j = 1; j <= N; j++) {
      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);
  for (let k = 0; k < iter; k++) {
    for (let i = 1; i

The file wrote successfully. Let me verify it's structurally sound by quickly checking it exists and has the basic pieces, and start a local server so the user can preview it.

Let me verify the file and also test it quickly to make sure it has no obvious issues.

The file is written. Let me verify the structure and make sure the JS is clean. Let me also double-check that there's no issue with the inject velocity vs per-frame injection — the inject multiplies by dt twice which would make it very weak. Let me review and fix the injection to be a proper impulse.

Actually, looking more carefully: in `inject`, I pass `dx, dy` as cells/sec, then compute `amount = strength * (norm/30) * dt`, and add `dx * dt * w * 0.6` to u. That's velocity impulse = dt * velocity. That's wrong — I should add velocity impulse directly (not multiplied by dt) since the velocity field is already per-frame in the solver (the solver uses dt inside advection). Actually, the classic Stam solver treats u,v as the field, and you add forces directly to u,v each step (force = dt * accel). So adding `dx * dt` is actually correct if dx is acceleration-like. But I named it as velocity cells/sec, so it should just be `dx * w * 0.6` added to u (force as velocity change per step), then the solver handles it.

Let me reconsider. Standard Stam: `u += force * dt` where force is acceleration. If I think of the mouse velocity as the desired velocity, then the force needed is (target - current) and the integration adds force*dt. For simplicity and responsiveness, adding `dx * dt * w` where dx is desired velocity (cells/sec) gives a proportional response — effectively adding impulse proportional to velocity. That's reasonable and feels snappy. But I'm also multiplying by dt once more in `amount` for dye. Let me keep `dx*dt` for velocity (force-based) — that's fine and stable.

Actually wait — the issue: I compute `vx = dmx/dt` giving cells/sec. Then in inject I add `vx * dt * w * 0.6` to u. That gives `dmx * w * 0.6` — independent of dt. Good, that's actually a position-delta-based impulse, which is stable and dt-independent. 

And dye amount: `amount = strength * (norm/30) * dt` where norm is mouse speed in cells/sec. So dye scales with speed * dt = distance moved. That's also dt-independent for a steady drag. Good.

Everything is consistent. The simulation should work well.

Let me do a quick sanity check by opening the file to see if there were any syntax issues I may have missed. Let me re-read the relevant injection function region.

The simulation is built but the `inject()` function computes a color from `Math.atan2(dy, dx)` and then uses the resulting hue to derive RGB. Good.

One thing — when the user drags, the injection happens at both the current and previous mouse positions. This creates a stronger trail effect.

Let me double-check a potential issue: in `onMove`, when `down` is false we skip. But `pmx`/`pmy` are set once on mousedown. Then `cmx`/`cmy` are updated on move. When down is true we inject using `cmx - pmx` divided by `dt`. But `dt` here is the frame delta. If the user moves the mouse fast between frames, that's fine.

But there's a subtle issue: we call `inject()` with `dt` being the full frame dt. However the injection should happen continuously. The current approach is reasonable.

Let me also verify: in the render function, I have two render functions (render and render2) and render2 is called. That's a bit messy. Let me check the full render section.

Actually I wrote two versions and keep render2. Let me clean that up since it's just wasting memory. Actually it's fine — the unused `render` function is never called, so it's dead code but doesn't hurt. However I can remove it for cleanliness.

Actually, let me leave it. The user just wants the simulation to work. Let me test it quickly by looking at whether it runs. I can't actually run a browser, but I can check for JS errors by reading the code carefully.

Let me trace through the key functions:
1. `inject` — uses `Math.atan2(dy, dx)`, converts to HSV, injects dye + velocity
2. `stepVelocity` — diffuse, project, advect, project, clamp
3. `stepDye` — diffuse, advect, decay

The `inject` function computes hue from atan2 and maps to RGB. Good.

One issue: in `onMove`, I only update `cmx`/`cmy` if `down` is true. But what if the user mousedowns, then moves — `down` is true so it works. On mouseup, `down` becomes false and `cmx`/`cmy` are reset.

I think this is solid. Let me finalize.

Actually wait — I notice I have a small issue. The `inject` call uses `cmx - pmx` and `cmy - pmy` but these are in grid coordinates, not pixels. That's correct since we converted to grid coords. The velocity `vmx = (cmx - pmx) / dt` is in grid cells per second. Good.

Let me also verify the initial setup: `pmx = -1, pmy = -1`. The `onStart` sets `pmx = cmx; pmy = cmy`. 

I'm satisfied the code works. Let me clean up the dead `render` function though — it's wasteful and confusing.

The exact text includes the comment `// ---------- render ----------` with trailing blank line. Let me read the area around line ~95-130.

The output is telling me the `render` and `render2` functions are actually the same, with `render2` replacing `render`. Let me just delete the `render` function since `render2` is what's actually used - but the edit says "Could not find the exact text". Let me re-read the file to see the actual content around the render functions.

The `render` function removal failed because whitespace doesn't match. Let me check the actual render section in the file.

I see the exact text now. Let me remove the unused `render()` function.

OK, now let me check what's left and whether there are other unused parts. Let me read the end of the file to make sure it's clean.

Now, find the leftover `render` function body (with its closing brace) so I can delete it.

The user's request was simply to build a single-page fluid simulation. I've created it and fixed the duplicate render function. The file now reads cleanly. Let me verify the final file is complete and coherent, then report back.

The dead render function is now gone. Let me do a final sanity check on the whole file:

The edit left a stray `ctx.imageSmoothingEnabled = true;` and `ctx.imageSmoothingQuality = 'high';` before the `// redraw render...` comment. Let me clean that up so the render function and main loop are tidy.

The file now has a clean render2. Let me do a final sanity check that the file is well-formed and doesn't have JS errors:

**Typing in Node REPL and checking parse errors**

I'm in the Node REPL, working with a match object for cleaning up HTML. I'm slicing out script content and trying to parse it as a Function to catch syntax errors before evaluating it. I'll log any errors and exit, keeping it quick and focused on scripting details.

The error could be from the render2 function syntax. Let me check the render2 function and the surrounding lines more carefully.

Let me find where the syntax error is - probably the stray blank line issue or a brace problem.

Let me not rely on node to parse it; instead find the parse error manually. The error is "Unexpected token '}'". This likely means somewhere there is a mismatch. Let me read the full JS section and look for mismatched braces.

Actually, easier: let me just read the full file and look for obvious issues.

The JS code ends with `})();` — that's an IIFE wrapping the arrow function IIFE. But the outer one is `( () => {...} )()` — wait, let me look at the very top of the script:

The file looks correct. The Node parse error is likely a Node issue with the IIFE ending — actually `})();` is fine. Let me try again with a raw module parse instead of Function constructor.

I want to find the exact location of the parse error. Let me extract the script content to a temp file and feed it to node with line numbers on errors:

Let me find the offending `}`.

I see the issue now — my earlier edits left the file in a broken state: lines 308–318 contain orphaned code fragments from the old `render()` function (variable definitions and comments that were part of the original function body) but the `function render()` declaration and its opening lines were removed. So there are dangling statements outside any function. Let me look at what I need to do: from line 308 to 318, I need to replace the orphaned code with clean content.

Let me look at the exact text from line 308 through line 318.

208: ctx.imageSmoothingEnabled = true;
209: ctx.imageSmoothingQuality = 'high';
210: // ---------- render ----------
211:         const idx = ((i - 1) + (j - 1) * N) * 4;
212:         data[idx]     = r8;
213:         data[idx + 1] = g8;
214:         data[idx + 2] = b8;
215:         data[idx + 3] = 255;
216:       }
217:     }
218:     offCtx.putImageData(img, 0, 0);
219:     // smooth scaling for fluid look
220: // redraw render using CSS-pixel dims
221:   function render2() {

So my plan: replace lines 210–220 (the orphaned render code and stray comment) with nothing — just delete those lines — and keep the `function render2()` declaration that starts at line 221.

Wait, I need to be precise about what's on each line. Let me re-read with awk to see exact line numbers and content.

Actually the awk output showed:
308: // ---------- render ----------
309:         const idx = ((i - 1) + (j - 1) * N) * 4;
310:         data[idx]     = r8;
311:         data[idx + 1] = g8;
312:         data[idx + 2] = b8;
313:         data[idx + 3] = 255;
314:       }
315:     }
316:     offCtx.putImageData(img, 0, 0);
317:     // smooth scaling for fluid look
318: // redraw render using CSS-pixel dims
319:   function render2() {

Hmm, but wait — line 309 has `const idx = ((i - 1) + (j - 1) * N) * 4;` which is indented as if inside a function. And lines 314, 315 have closing braces. And line 318 has `// redraw render using CSS-pixel dims` which is a comment I put there. And then line 319 is `function render2() {`.

So the structure is: at line 308 there's a comment `// ---------- render ----------` followed by orphaned code that was part of the old `render()` function (but the function declaration and opening lines were removed), ending with a stray comment at line 318, then `function render2()` at line 319.

The issue: lines 308–318 are orphaned code outside any function. They reference variables `idx`, `r8`, `g8`, `b8`, `data`, `offCtx`, `N`, `img` that are scoped inside a function, but these lines are now at module/top level, causing syntax errors.

I need to delete lines 308–318 entirely (the orphaned code and the stray comment), leaving only `function render2() {` at line 319 (which would become line 308 after deletion).

Let me do this carefully. I'll construct an edit that replaces the orphaned block from `// ---------- render ----------` (line 308) through `// redraw render using CSS-pixel dims` (line 318) with nothing, preserving the `  function render2() {` line.

But wait — I need to double-check: the `  function render2() {` line is at 319, and it's indented with 2 spaces. After I delete 308–318, that line becomes line 308 and should work fine as a function declaration.

Actually, let me reconsider — maybe a cleaner approach is to replace lines 308–318 with the proper opening of `render2()`. But `render2()` already starts at line 319. So I just need to delete 308–318.

Let me construct the edit. The orphaned block starts at line 308: `// ---------- render ----------` and ends at line 318: `// redraw render using CSS-pixel dims`. I need to match this exact text including newlines.

From the awk output:
- Line 308: `  // ---------- render ----------` (with 2-space indent? Let me check)

Actually the awk output shows `308: // ---------- render ----------` without leading spaces shown, but that's because awk's printf "%s" strips them... no wait, printf "%s" preserves spaces. Let me look again: `308: // ---------- render ----------` — there are no leading spaces shown. But the original write had `  // ---------- render ----------` with 2 leading spaces (since it was inside the IIFE). Hmm, but the read output earlier showed:

From read offset=320 earlier:
```
  // ---------- render ----------
  function render() {
```

Wait no, that was before edits. After my edits, the structure changed. Let me just look at the raw lines with cat -A or similar to see exact whitespace.

Actually, I realize I should just look at the actual file content around those lines very carefully. Let me use read with offset to see lines 305-322 with exact content.

Actually, I already have a good picture. The issue is that lines 308-318 contain orphaned code. Let me just replace that block. The safest way: match from `// ---------- render ----------` (the comment at line 308) through `// redraw render using CSS-pixel dims` (line 318, which is a comment I introduced), and replace with empty string.

But I need the exact text. Let me read lines 305-322.

From the read output at offset=320 (which was before my last edit), I saw:
```
  // ---------- render ----------
  function render() {
    const img = offCtx.createImageData(N, N);
    ...
```

But that was the OLD state. After my edit that removed `function render() { ... }` up to the `const r8 = ...` line, the state changed. Then my second edit removed more. Let me just read the current file state around lines 305-325 to get exact text.

I'll read with offset=305 and limit=20.

The previous edit removed the `function render() { ...` wrapper lines but left the "stub" comment line `// ---------- render ----------` and leftover function body lines that were inside `render()`. So now there's orphaned code: lines 309-317 are a fragment that was inside `render()` but the function declaration is gone.

I need to remove the orphaned fragment (lines 308-318 area) entirely. Let me look at a bit more context.

The orphaned fragment is:

```
// ---------- render ----------
        const idx = ((i - 1) + (j - 1) * N) * 4;
        data[idx]     = r8;
        data[idx + 1] = g8;
        data[idx + 2] = b8;
        data[idx + 3] = 255;
      }
    }
    offCtx.putImageData(img, 0, 0);
    // smooth scaling for fluid look
// redraw render using CSS-pixel dims
  function render2() {
```

This is garbage. I'll replace from `// ---------- render ----------` through to (but not including) `  function render2() {`.Actually, it's cleaner to just replace that whole broken region with a clean `// ---------- render ----------` comment and then `function render2() {`.

Let me do the edit.

I see a syntax error in the file. The render2 function is missing its opening brace. The edit showed the closing brace of the script, but the actual function body starts on the line after. Let me check the actual file content around that area to see the full context and identify where the missing brace is.

Looking at the edit result, it shows:
```javascript
// ---------- render ----------
// redraw render using CSS-pixel dims
function render2() {
  const img = offCtx.createImageData(N, N);
  const data = img.data;
  ...
```

But the original file had:
```javascript
// ---------- render ----------
        const idx = ((i - 1) + (j - 1) * N) * 4;
        data[idx]     = r8;
        ...
      }
    }
```

So the render function body got partially removed, but a few lines of the old render remained (like the `data[idx]` lines). Then after those, `// redraw render using CSS-pixel dims` and `function render2() {` appeared. So the file now has garbage code between the comment and `render2()`. Let me read more of the file to see the full extent of the corruption.

Let me do this cleanly — extract just the script to a temp .js file and check it.

The JS parses cleanly now. Let me do a final review of the full file.

The user asked me to create a real-time fluid simulation. I already built it and delivered the file. My last response was just a read of the file for review — that's not what I should output to the user. Let me give a short summary of what was built and confirm the file is ready.