The user wants a real-time fluid simulation on a full-screen canvas with mouse interaction where dragging pushes colored dye. I'll build this as a single-file index.html with CSS and JS inline, using no external resources.

I'll implement a stable fluid simulation using the Navier-Stokes equations on a grid. This is a classic technique used in games and visual effects (like Jos Stam's "Real-Time Fluid Dynamics for Games"). I'll use a semi-Lagrangian advection approach.

Let me build this step by step:

1. **Grid-based fluid simulation**: Velocity field on a staggered grid (or collocated), with pressure projection
2. **Dye advection**: Color dye that gets advected by the velocity field
3. **Mouse interaction**: Dragging adds velocity and injects dye
4. **Full-screen canvas**: Responsive canvas

I'll implement Jos Stam's stable fluid solver:
- Add forces (mouse input)
- Advect velocity
- Diffuse velocity  
- Project (make it incompressible via pressure solve)
- Advect dye
- Diffuse dye

For the color, I'll use a swirling rainbow dye pattern that changes color based on direction/duration.

Let me write a complete implementation:

Key components:
- Velocity grid (u, v components)
- Density/dye grid (RGB channels)
- Jacobi iteration for diffusion and pressure
- Semi-Lagrangian advection
- Mouse drag tracking

I'll use a reasonable grid size for performance (maybe 100-150 cells on each axis) and scale it to full screen.

Let me write this out carefully.

Actually, for performance with JavaScript on a full-screen canvas, I should be careful about grid resolution. Let me use a fixed grid (like 128x128 or 100x100 internal simulation) and render it scaled up to the full canvas.

Let me think about color dye. The user wants "swirling colored dye." I'll make it so that dragging injects dye with colors that shift over time (like a rainbow/hue rotation), creating beautiful swirling patterns.

Here's my plan:

```javascript
// Grid size
const N = 128;
const size = N + 2; // with boundary cells

// Arrays
let u = new Float32Array(size * size);    // x-velocity
let v = new Float32Array(size * size);    // y-velocity
let u_prev = new Float32Array(size * size);
let v_prev = new Float32Array(size * size);
let dens_r = new Float32Array(size * size); // red dye
let dens_g = new Float32Array(size * size); // green dye
let dens_b = new Float32Array(size * size); // blue dye
let dens_r_prev, dens_g_prev, dens_b_prev;
```

Actually, let me simplify. I'll use one density array for grayscale or I can use 3 arrays for RGB. Let me use 3 arrays for RGB to get the colored dye effect.

For the mouse, I'll track position and add velocity + dye at the mouse location when dragging.

Let me write the complete code now.

Actually, for a clean and performant implementation, let me use a grid that's about 150 internal cells, which should be fast enough at 60fps. The canvas can be any resolution - we just render the internal grid scaled up.

Let me code this up:

```javascript
const SIM_SIZE = 120; // internal simulation size
const SIZE = SIM_SIZE + 2; // with boundary cells

// Index function
function IX(x, y) { return x + y * SIZE; }
```

Wait, I need to be careful - this might be slow with large arrays and many Jacobi iterations. Let me use a moderate grid like 80-100 cells for better performance.

Actually, modern JS engines are pretty fast. Let me use 100 cells and see. I can adjust.

Let me write the full implementation now. I'll use:

- Grid: 100x100 internal (102x102 with boundaries)
- Jacobi iterations: ~20 for diffusion, ~40 for pressure
- Viscosity and diffusion parameters
- Mouse interaction for adding velocity and dye
- Rainbow color dye that changes hue over time

Here's the full implementation:

The fluid solver has these steps each frame:
1. Add velocity from mouse
2. Diffuse velocity (viscosity)
3. Project velocity (pressure solve to maintain incompressibility)
4. Advect velocity  
5. Project velocity again
6. Add dye from mouse
7. Diffuse dye
8. Advect dye

Actually, the classic Stam solver order is:
1. Add forces
2. Advect
3. Diffuse
4. Project

But there are variations. Let me follow the standard approach.

Let me also think about the rendering. I'll render the dye density to the canvas using putImageData or fillRect per cell. Since the internal grid is smaller than the canvas, I'll scale it up.

For smooth rendering with color, I'll create an ImageData and set pixel colors based on the dye RGB values.

Let me write this all out now.

Actually, for even better visuals, I can use the velocity magnitude/curl to enhance the dye colors or add some turbulence-like effects.

Let me finalize the code. Here's what I'll do:

- Canvas: full screen, resizes with window
- Internal grid: 128x128
- Mouse drag adds velocity in the direction of movement AND injects colored dye
- Dye color shifts hue based on a global time variable, creating a rainbow swirl effect
- Render using ImageData for performance

Let me be careful with the math. The projection step (making the velocity field divergence-free) involves solving a Poisson equation for pressure using Jacobi iteration.

The discrete divergence at cell (i,j):
div = -0.5 * h * (u[i+1,j] - u[i-1,j] + v[i,j+1] - v[i,j-1])

Actually, the standard formulation:

For a cell-centered grid with cell size h:
- divergence at (i,j) = (u[i+1,j] - u[i-1,j])/(2h) + (v[i,j+1] - v[i,j-1])/(2h)
- gradient of pressure: grad_x = (p[i+1,j] - p[i-1,j])/(2h)

In Stam's implementation (which uses a unit cell size), things simplify. Let me follow his approach closely.

Actually, I'll implement this based on the well-known "Real-Time Fluid Dynamics for Games" paper by Jos Stam. The key insight is using a staggered grid approach implicitly through the linear solver.

Let me write a clean implementation:

```javascript
// Constants
const N = 100;  // grid size
const size = (N + 2) * (N + 2);
const dt = 0.1;
const diff = 0.0001;
const visc = 0.0001;

// Main simulation arrays
let u = new Float32Array(size);     // x velocity
let v = new Float32Array(size);     // y velocity  
let uPrev = new Float32Array(size);
let vPrev = new Float32Array(size);
let dens = new Float32Array(size);  // density (we'll expand for color)
let densPrev = new Float32Array(size);
```

Hmm, for colored dye I need 3 density arrays (R, G, B). Let me use:
- densR, densG, densB

Or I can use a single density but render with a color based on something else... No, let me use 3 density arrays for true colored dye.

Actually, to keep it simpler and focused on the visual effect, let me use:
- 3 dye arrays (r, g, b)
- Mouse adds dye with a hue that rotates over time

Let me finalize and write the code. I'll aim for something performant and visually stunning.

One more consideration: the Jacobi solver. I need to be careful about numerical stability. Using 20-40 iterations should be enough for visual quality.

Alright, let me write the complete index.html now. Here's the full plan:

**HTML structure:**
- Full-screen canvas
- Maybe a subtle instruction text

**JavaScript:**
- Fluid solver class
- Mouse tracking
- Animation loop
- Resize handler

**CSS:**
- Full screen, no margins, no scrollbars
- Canvas fills the viewport

Let me write it all out.

I need to be very precise in my implementation. Let me think through the key algorithms:

**Jacobi solver for diffusion:**
```
For each iteration:
  For each interior cell (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)])) / (1+4*a)
```
where a = diff * dt * N * N (for diffusion)

**Projection step:**
```
// Compute divergence
For each cell:
  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

// Solve pressure Poisson
For k iterations:
  For each cell:
    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

// Subtract gradient
For each cell:
  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)])
```

**Advection (semi-Lagrangian):**
```
For each cell:
  // Trace backward
  x = i - dt*N*u[IX(i,j)]
  y = j - dt*N*v[IX(i,j)]
  // Clamp to boundaries
  x = clamp(x, 0.5, N+0.5)
  y = clamp(y, 0.5, N+0.5)
  // Interpolate
  i0 = floor(x), j0 = floor(y)
  i1 = i0+1, j1 = j0+1
  s1 = x-i0, s0 = 1-s1
  t1 = y-j0, t0 = 1-t1
  result = s0*(t0*d[IX(i0,j0)]+t1*d[IX(i0,j1)]) + s1*(t0*d[IX(i1,j0)]+t1*d[IX(i1,j1)])
```

**Mouse force:**
When mouse is down and moving:
- dx = mouseX - prevMouseX (scaled to grid coords)
- dy = mouseY - prevMouseY
- Add velocity at mouse position proportional to (dx, dy)
- Add dye at mouse position

For the color, I'll inject RGB values based on the mouse movement direction and a global hue that rotates slowly.

Now, I also need to handle boundaries. The boundary conditions ensure no flow through walls. Stam's paper handles this by setting the velocity at boundaries to be the negative of the adjacent interior cell.

Let me simplify boundary handling: after each operation, I'll set boundaries.

For the canvas rendering:
- `putImageData` is faster than `fillRect` per cell
- I'll create an ImageData buffer matching the canvas size
- For each grid cell, I'll fill the corresponding pixel region with the dye color
- The background will be black/dark

Alright, I think I have enough to write a solid implementation. Let me code it up.

Actually wait - let me reconsider the grid size. At 100x100, each frame needs to process 10,000 cells multiple times (diffusion, advection, projection with ~20 Jacobi iterations each). That's several hundred thousand operations per frame. With Float32Array and modern JS, this should be fine.

But rendering 100x100 cells onto a full-screen canvas (say 1920x1080) means each cell covers about 19x11 pixels. That might look blocky. Let me use bilinear interpolation when rendering, or use a larger grid.

Actually, for full HD at 60fps with a large internal grid, JS might struggle. Let me use a 150-200 grid for better visual quality and rely on the fact that modern browsers are fast.

Hmm, let me compromise. I'll use N=128 for the internal simulation. That's 16,384 cells. Each operation needs ~16K operations. With 60 Jacobi iterations total (20 for diffusion + 40 for pressure), that's about 1M operations per frame, which should be fine.

Wait, actually the Jacobi iterations for pressure solve are the bottleneck. 40 iterations × 16K cells = 640K operations just for pressure. Plus diffusion (20 × 16K = 320K) and advection (16K each for u, v, and dye). Total around 1-2M operations per frame. On a modern machine, this should be okay for 60fps.

Let me use N=128 and implement it.

One more thing: I should also handle the window resize event to resize the canvas and re-render.

Alright, let me write the complete code now. I'll be very careful and thorough.

Actually, I realize there's an important optimization: I should use `requestAnimationFrame` and compute dt from the actual frame time rather than using a fixed dt.

And for the mouse, I should track the velocity of mouse movement (not just position) to create appropriate fluid forces.

Let me also add some nice features:
- Fade the dye slightly each frame so it doesn't accumulate forever
- Add some initial "ambient" dye to make it look alive
- The color of injected dye rotates through the color wheel

OK I'm going to write this now. Let me be complete and correct.

For the rendering, I'll use `createImageData` to get a pixel buffer, then fill it based on the dye values, then `putImageData` to the canvas. This is efficient.

Each grid cell maps to multiple canvas pixels. I'll compute the mapping and fill all pixels within each cell's rectangle.

Actually, for even smoother rendering, I can do bilinear interpolation of the dye values when rendering to canvas pixels. But that might be slow. Let me just use nearest-neighbor (blocky) rendering for performance, and use a larger N (like 150-200) to make the blocks less noticeable.

With N=150, each block is about 1920/150 ≈ 13 pixels wide on a full HD display. Not too bad.

Let me finalize: N=128, rendering with nearest-neighbor (just fill rects). I'll use ImageData for speed.

Actually, you know what, let me use `createImageData` and compute the pixel values directly. For each canvas pixel, find which grid cell it maps to, and use that dye color. This way each block gets filled smoothly.

Wait, that's just nearest-neighbor interpolation done in JavaScript. Let me just compute:

```
cellW = canvas.width / N
cellH = canvas.height / N
```

And fill each cell as a rectangle in the ImageData. This is faster than calling fillRect for each cell.

Actually, the fastest approach: for a canvas of size W×H and grid N×N:
- cellW = W / N, cellH = H / N
- For each grid cell (i,j), fill the rectangle from (i*cellW, j*cellH) to ((i+1)*cellW, (j+1)*cellH) in the ImageData

This involves setting pixel values in the ImageData's Uint8ClampedArray. Let me be efficient about it.

For each cell (i,j) with color (r,g,b):
```
for (let py = j*cellH; py < (j+1)*cellH; py++) {
  for (let px = i*cellW; px < (i+1)*cellW; px++) {
    let idx = (py * W + px) * 4;
    data[idx] = r;
    data[idx+1] = g;
    data[idx+2] = b;
    data[idx+3] = 255;
  }
}
```

This is a lot of pixel writes. For N=128 and canvas ~1920×1080, that's about 2M pixels being set. It might be slow to loop through all pixels in JS. Let me use `fillRect` approach instead - create a small offscreen canvas for each cell color and use `drawImage`, or just use `fillRect` on the main canvas.

Actually, `ctx.fillRect` is hardware accelerated and would be quite fast for 128×128 = 16K rectangles. Let me use that approach. For each cell, set fillStyle and fillRect.

But wait - 16K fillRect calls per frame might have overhead. Let me think...

Alternative: I can use a small offscreen canvas at the simulation resolution (N×N) and then draw it scaled up to the main canvas. The offscreen canvas would have one pixel per grid cell. Then `drawImage` with scaling would handle the rest via GPU.

That's clever! Let me do that:
1. Create a small canvas of size N×N
2. Set each pixel's color based on the dye value
3. Draw the small canvas scaled to full screen on the main canvas
4. The browser handles the upscaling efficiently

This way I only set N×N = 16K pixels (using ImageData on the small canvas) instead of 2M.

Yes, that's the way to go. Let me implement this.

Actually, for even better visual quality, I can use the CSS `image-rendering: auto` (default) which does bilinear filtering when scaling up. This will give a smooth look instead of blocky.

Alright, here's my final plan:

**JavaScript Architecture:**
- `IX(x, y)` - index function for 2D arrays
- `FluidSim` class with:
  - Constructor: allocates arrays, sets up
  - `addForce(x, y, dx, dy)` - adds velocity at position
  - `addDye(x, y, r, g, b, amount)` - adds dye at position
  - `step(dt)` - performs one simulation step
  - `diffuse(b, x, x0, diff, dt)` - diffuse a field
  - `advect(b, d, d0, u, v, dt)` - advect a field
  - `project(u, v, p, div)` - projection step
  - `setBounds(b, type)` - set boundary conditions
  - `render(ctx)` - render to canvas
- Mouse event handlers
- Animation loop

Let me write it now!

Actually, for the dye rendering, I want beautiful colors. Let me use HSV-to-RGB conversion:
- Hue: changes based on mouse direction and global time
- Saturation: 1.0 (full)
- Value: based on mouse speed

Or simpler: have a global hue that slowly rotates, and use it for all newly injected dye. This creates a rainbow swirl as the user drags.

I'll also add a subtle fade so old dye gradually dissipates (giving it a lifespan).

Alright, coding time. I need to be very careful with the math.

For the projection step, let me follow Stam's formulation exactly:

The projection solves ∇²p = ∇·w (where w is the intermediate velocity), then subtracts ∇p.

In discrete form:
- div[i,j] = (u[i+1,j] - u[i-1,j] + v[i,j+1] - v[i,j-1]) / (2 * h)
  where h = 1/N

Wait, in Stam's code, the cell size is implicitly 1/N. Let me check...

In Stam's GDC 2003 paper, the spatial step h = 1/N. So:
- div = (u[i+1,j] - u[i-1,j] + v[i,j+1] - v[i,j-1]) * N / 2
- Jacobi for pressure: p[i,j] = (div[i,j] + p[i-1,j] + p[i+1,j] + p[i,j-1] + p[i,j+1]) / 4
- Subtract gradient: u[i,j] -= (p[i+1,j] - p[i-1,j]) * N / 2

Wait, let me re-derive this. The standard projection is:
1. Compute w (divergent velocity field)
2. Solve ∇²p = ∇·w
3. u = w - ∇p

Discretized on a uniform grid with spacing h:
- ∇·w at (i,j) = (w_x[i+1,j] - w_x[i-1,j])/(2h) + (w_y[i,j+1] - w_y[i,j-1])/(2h)
- ∇²p at (i,j) = (p[i+1,j] + p[i-1,j] + p[i,j+1] + p[i,j-1] - 4p[i,j])/h²

So the linear system is:
(p[i+1,j] + p[i-1,j] + p[i,j+1] + p[i,j-1] - 4p[i,j]) = h * (w_x[i+1,j] - w_x[i-1,j] + w_y[i,j+1] - w_y[i,j-1]) / 2

Rearranging for Jacobi:
p[i,j] = (p[i+1,j] + p[i-1,j] + p[i,j+1] + p[i,j-1] - h * div_w[i,j]) / 4

where div_w = (w_x[i+1,j] - w_x[i-1,j] + w_y[i,j+1] - w_y[i,j-1]) / 2

With h = 1/N:
p[i,j] = (p[i+1,j] + p[i-1,j] + p[i,j+1] + p[i,j-1] - div_w[i,j] / N) / 4

Then gradient subtraction:
u[i,j] = w_x[i,j] - (p[i+1,j] - p[i-1,j]) / (2h) = w_x[i,j] - (p[i+1,j] - p[i-1,j]) * N / 2
v[i,j] = w_y[i,j] - (p[i,j+1] - p[i,j-1]) * N / 2

Hmm, I'm getting a bit confused with the signs. Let me look at a reference implementation approach.

Actually, let me simplify. Many JS fluid sim implementations use this approach:

**Projection step:**
```
// div = -(u[i+1]-u[i-1] + v[i,j+1]-v[i,j-1]) * 0.5 * h
// With h=1 for simplicity:
// div = -(u[i+1]-u[i-1] + v[i,j+1]-v[i,j-1]) * 0.5

// But with N scaling:
// div = -(u[i+1]-u[i-1] + v[i,j+1]-v[i,j-1]) * 0.5 * N  (if grid has N cells per unit)
```

Actually, you know what, let me just follow the well-tested implementation pattern. The key insight is that the scaling factors need to be consistent. Let me use the approach where the grid coordinates go from 0 to N+1 (with the interior being 1 to N), and cell size is implicitly 1.

In this formulation (which is what Stam uses):
- The grid has N interior cells
- Cell size is 1 (so the total domain is N×N cells of size 1 each)

The formulas become:
- div[i,j] = -0.5 * (u[i+1,j] - u[i-1,j] + v[i,j+1] - v[i,j-1])
- p[i,j] = (div[i,j] + p[i-1,j] + p[i+1,j] + p[i,j-1] + p[i,j+1]) / 4
- u[i,j] -= 0.5 * (p[i+1,j] - p[i-1,j])
- v[i,j] -= 0.5 * (p[i,j+1] - p[i,j-1])

This is clean and simple. Let me use this.

But wait - this means the velocities are in "cells per time step" units. I need to make sure the advection and forces use consistent units.

For advection, the semi-Lagrangian backtrace:
x = i - dt * u[i,j]
y = j - dt * v[i,j]

Where u and v are in cells per unit time, dt is the time step. If u=1, that means 1 cell per time unit. If dt=0.1, we trace back 0.1 cells. This makes sense.

For forces from mouse:
When the mouse moves by dx pixels from one frame to the next, I need to convert that to velocity in grid cells. If the canvas width is W pixels and the grid is N cells wide, then dx pixels = dx * N / W cells. The velocity (cells per time unit) would be (dx * N / W) / dt.

Let me set things up:
- Mouse dx in pixels → convert to grid cells: dx_grid = dx * N / canvasWidth
- Force = dx_grid / dt (since velocity = distance / time)
- But I'll just add this as a velocity impulse: u += force * dt = dx_grid

So effectively, I add dx_grid to the velocity at the mouse position. But since I'm adding force, not setting velocity, I'll multiply by some strength factor.

OK, I think I'm overcomplicating this. Let me just write reasonable code and tune the parameters. Many fluid sim demos just add mouse delta directly as velocity with some scaling, and it looks great.

Let me now write the actual code. I'll be very careful to get it right.

One more thing about the dye: I want RGB channels. So I'll have 3 density arrays and advect/diffuse each one separately. This triples the work for dye, but it's worth it for the color effect.

For the color injection, I'll use a hue that changes over time:
```javascript
let hue = (Date.now() / 1000) % 1; // cycles every second
let [r, g, b] = hsvToRgb(hue, 1, 1);
```

But I want different hues for different drags, creating a nice rainbow swirl. I can use the mouse position to influence the hue:
- hue = (mouseX / canvasWidth + globalTime * 0.1) % 1

This way, dragging left creates red dye, right creates blue, etc., and it slowly shifts.

OK, I'm going to write the complete code now. No more planning!

Let me also handle window resize properly and use devicePixelRatio for sharp rendering on Retina displays.

Hmm, using devicePixelRatio means the canvas pixel dimensions are larger than CSS dimensions, which complicates mouse coordinate mapping. Let me keep it simple: CSS full-screen canvas, no devicePixelRatio scaling. The internal simulation resolution is independent of the screen resolution.

Actually for the rendering, using the small offscreen canvas approach, devicePixelRatio doesn't matter much since we're upscaling from the simulation grid anyway. Let me keep it simple.

Here's the complete code I'll write:

```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>
/* full screen, no scroll */
</style>
</head>
<body>
<canvas id="canvas"></canvas>
<script>
// Fluid simulation code
</script>
</body>
</html>
```

Let me write the full implementation now. I'll focus on correctness and visual quality.

One key optimization: I'll pre-allocate all working arrays in the constructor to avoid GC pressure during animation.

Let me also think about the boundary conditions. I want "walls" on all 4 sides (no fluid escapes). The boundary condition for velocity is:
- At left wall (x=0): u[0,j] = 0 (no flow through wall; actually for the staggered grid it's different)
- For a collocated grid, the no-slip condition means: at boundaries, the velocity component normal to the wall is zero, and the tangential component is set to the negative of its neighbor (so the average vanishes at the wall)

For the dye, boundaries have zero gradient (no dye escapes).

OK let me just write the code. I'll follow the standard Stam implementation closely.

Here we go:

Let me think about one more thing. The `dt` in the simulation. For stability with semi-Lagrangian advection, dt shouldn't be too large. Typically dt = 0.1 to 0.5 works well. But since I'm using requestAnimationFrame, the actual frame time varies (typically ~16ms at 60fps). I should cap dt and use sub-stepping if needed.

I'll use a fixed dt of ~0.016 (matching 60fps) and run one simulation step per frame. If the frame time is larger, I might need to sub-step. But for simplicity, let me use a small fixed dt and accept that the simulation might run slightly slower on low-fps devices.

Alright, NOW I'm writing the code:

Actually, let me reconsider the array organization. I'll use typed arrays for performance:

- `u`, `v`: velocity components (size = (N+2)*(N+2))
- `uPrev`, `vPrev`: previous velocity (for advection)
- `densR`, `densG`, `densB`: dye density (size = (N+2)*(N+2))
- `densRPrev`, `densGPrev`, `densBPrev`: previous dye
- `p`: pressure (for projection)
- `div`: divergence (for projection)

All as Float32Array.

Now for the functions, I'll keep them as standalone functions operating on these arrays, passing N and size as needed. Or I can make a class.

Let me use a class for organization.

Here's the full implementation. I'm going to be very detailed and careful:

Actually, I realize I should double-check my Jacobi solver for the pressure projection. The standard formula:

For the Poisson equation ∇²p = ∇·w:

Discretized:
(p[i+1,j] + p[i-1,j] + p[i,j+1] + p[i,j-1] - 4p[i,j]) / h² = (w_x[i+1,j] - w_x[i-1,j] + w_y[i,j+1] - w_y[i,j-1]) / (2h)

With h=1:
p[i+1,j] + p[i-1,j] + p[i,j+1] + p[i,j-1] - 4p[i,j] = 0.5 * (w_x[i+1,j] - w_x[i-1,j] + w_y[i,j+1] - w_y[i,j-1])

Let me define:
div[i,j] = 0.5 * (u[i+1,j] - u[i-1,j] + v[i,j+1] - v[i,j-1])

Then the equation is:
4p[i,j] - p[i+1,j] - p[i-1,j] - p[i,j+1] - p[i,j-1] = -div[i,j]

Wait no:
p[i+1,j] + p[i-1,j] + p[i,j+1] + p[i,j-1] - 4p[i,j] = div[i,j]

For Jacobi iteration:
4p[i,j] = p[i+1,j] + p[i-1,j] + p[i,j+1] + p[i,j-1] - div[i,j]
p[i,j] = (p[i+1,j] + p[i-1,j] + p[i,j+1] + p[i,j-1] - div[i,j]) / 4

Then gradient correction:
u_new[i,j] = u[i,j] - 0.5 * (p[i+1,j] - p[i-1,j])
v_new[i,j] = v[i,j] - 0.5 * (p[i,j+1] - p[i,j-1])

This looks right. Now let me write the full code.

One more consideration: I should use "setBounds" after each operation to enforce boundary conditions. For velocity at walls:
- Left/right walls: u[0,j] = -u[1,j], ensuring zero average = no flow through wall
- Top/bottom walls: v[i,0] = -v[i,1]
- Same for the far walls: u[N+1,j] = -u[N,j], v[i,N+1] = -v[i,N]

For density at walls, we want zero gradient (no flux): just copy the neighboring interior value.

OK here we go. Writing the complete code now.

Actually, let me think about the layout of the arrays a bit more. I'll use a flat array of size (N+2)*(N+2). Index function:
```javascript
function IX(x, y) { return x + y * (N+2); }
```

For loops, interior cells go from x=1 to x=N, y=1 to y=N. The boundary cells are x=0, x=N+1, y=0, y=N+1.

Now let me write the complete code. I'll be very careful.

Let me write it step by step, then assemble it into the final file.

**addForce function:**
Adds velocity impulse at a specific grid position (x,y in grid coordinates):
```javascript
addForce(x, y, fx, fy) {
    if (x < 1 || x > N || y < 1 || y > N) return;
    let idx = IX(Math.round(x), Math.round(y));
    u[idx] += fx;
    v[idx] += fy;
}
```

But for smoother results, I should spread the force to neighboring cells using bilinear interpolation (the reverse of what we do for advection). That way the force is distributed smoothly.

```javascript
addForce(x, y, fx, fy) {
    x = Math.max(0.5, Math.min(N + 0.5, x));
    y = Math.max(0.5, Math.min(N + 0.5, y));
    let i0 = Math.floor(x), j0 = Math.floor(y);
    let i1 = i0 + 1, j1 = j0 + 1;
    let sx = x - i0, sy = y - j0;
    
    // Distribute to 4 surrounding cells
    u[IX(i0,j0)] += fx * (1-sx) * (1-sy);
    u[IX(i1,j0)] += fx * sx * (1-sy);
    u[IX(i0,j1)] += fx * (1-sx) * sy;
    u[IX(i1,j1)] += fx * sx * sy;
    
    // Same for v
    v[IX(i0,j0)] += fy * (1-sx) * (1-sy);
    v[IX(i1,j0)] += fy * sx * (1-sy);
    v[IX(i0,j1)] += fy * (1-sx) * sy;
    v[IX(i1,j1)] += fy * sx * sy;
}
```

Similarly for dye.

OK, I think I have everything planned. Let me write the complete HTML file now. I'll aim for about 400-500 lines of clean, well-commented code.

For the visual style:
- Dark background
- Dye colors pop against dark background
- Smooth, fluid-like appearance

Let me think about color injection. I want the dye to look beautiful:
- Hue rotates with global time
- Different mouse positions get different colors
- Over time, the colors mix and create beautiful patterns

For rendering, I'll map dye RGB values (which could be > 1.0 due to accumulation) to 0-255:

```javascript
let r = Math.min(255, Math.max(0, densR[idx] * 255));
let g = Math.min(255, Math.max(0, densG[idx] * 255));
let b = Math.min(255, Math.max(0, densB[idx] * 255));
```

I'll also apply a slight overall fade each frame to prevent dye from accumulating too much.

Now, one important thing: when the user drags, I need to inject dye at each intermediate position between the previous and current mouse position. I'll interpolate along the line.

Also, I want to inject the dye in a small radius around the mouse, not just at a single cell. This makes it look smoother.

Let me use a Gaussian-like spread: add dye to a 3x3 or 5x5 neighborhood around the mouse position with weights decreasing by distance.

OK, I'm going to write this now. Here's the final code.

Actually, for even better visuals, let me add some "smoke" lighting effect. Instead of just showing raw dye, I can use the density as alpha against a dark background, making it look more volumetric. Or I can use the velocity magnitude to create highlights.

Let me keep the rendering simple: just show the dye colors on a dark background. The dye colors will be vibrant and the fluid motion will do the rest.

I'll also add the Jacobi iteration count as a tunable parameter. More iterations = more accurate physics but slower. 20-30 iterations for the pressure solve is usually enough for visual quality.

Hmm actually, now I want to reconsider the grid size. 128 might be a bit large leading to slow performance. Let me go with N=100 for a nice round number and good performance. 100×100 = 10,000 cells, which is quite manageable.

Wait, the user said "full-screen canvas" and they'll be dragging the mouse. With 100 cells across a 1920px screen, each cell is 19px wide. That's a reasonable resolution for fluid dynamics. The upscaling will smooth things out.

Let me finalize with N=100. If it's too slow, it can be reduced, but I think 100 will work fine.

OK, writing now for real. Full code:

Let me also think about the Jacobi solver efficiency. I need to avoid creating new arrays or objects inside the animation loop. All working arrays should be pre-allocated.

For the diffusion and projection, I'll use the same temporary arrays (swap buffers). The diffusion function needs:
- x: current field
- x0: previous field  
- a: diffusion coefficient * dt * N*N

For the projection, I need:
- p: pressure (zero-initialized)
- div: divergence

I'll allocate these once and reuse.

Let me structure the code:

```javascript
class FluidSim {
    constructor(N) {
        this.N = N;
        this.size = (N + 2) * (N + 2);
        
        // Velocity
        this.u = new Float32Array(this.size);
        this.v = new Float32Array(this.size);
        this.uPrev = new Float32Array(this.size);
        this.vPrev = new Float32Array(this.size);
        
        // Dye (RGB)
        this.dr = new Float32Array(this.size);
        this.dg = new Float32Array(this.size);
        this.db = new Float32Array(this.size);
        this.drPrev = new Float32Array(this.size);
        this.dgPrev = new Float32Array(this.size);
        this.dbPrev = new Float32Array(this.size);
        
        // Working arrays for projection
        this.p = new Float32Array(this.size);
        this.div = new Float32Array(this.size);
        
        // Offscreen canvas for rendering
        this.renderCanvas = document.createElement('canvas');
        this.renderCanvas.width = N;
        this.renderCanvas.height = N;
        this.renderCtx = this.renderCanvas.getContext('2d');
        this.imageData = this.renderCtx.createImageData(N, N);
        
        // Parameters
        this.viscosity = 0.0000001;
        this.diffusion = 0.0000001;
        this.dt = 0.016;
        this.dyeFade = 0.9995; // slight fade each frame
    }
    
    // Methods...
}
```

Wait, I had `visc = 0.0001` earlier which is very large. Let me think about appropriate values.

For viscosity, the diffusion coefficient in the PDE is ν∇²u. The stability condition for the Jacobi solver is ν*dt/(h²) ≤ some value. With h = 1/N = 0.01 (for N=100), h² = 0.0001. So ν*dt/h² = ν*0.016/0.0001 = ν*160. For stability, this should be reasonably small.

Actually, in Stam's implementation, the diffusion parameter `a = diff * dt * N * N`. For N=100, dt=0.016, N²=10000, so a = diff * 160. If diff = 0.00001, a = 0.0016. The Jacobi solver handles this well for a < 1 or so.

For the fluid to look nice and swirly, I want fairly low viscosity. Let me use:
- viscosity diffusion: small, like 0.000001 → a ≈ 0.00016
- dye diffusion: also small, like 0.000001

Actually wait, for the dye, I want some diffusion so it spreads nicely. Let me use slightly higher diffusion for dye.

Let me set:
- visc = 0.000001
- dyeDiff = 0.00001

These give small `a` values which means few Jacobi iterations needed.

Actually, you know what, for visual quality, let me use:
- visc = 0.0 (no viscosity - more swirly)
- dyeDiff = 0.000001 (tiny diffusion for dye)

With almost no viscosity, the fluid will be very swirly/chaotic, which looks great!

For the projection (pressure) solve, I want enough iterations to maintain incompressibility. 30-40 iterations should be good.

Let me also add a "fade" step to the dye to give it a finite lifetime. I'll reduce all dye values by a small percentage each frame (like multiply by 0.9995). This prevents the screen from becoming a uniform color.

OK, truly writing the code now. No more deliberation.

One last thing: for the mouse interaction, when dragging, I want to add:
1. Velocity: proportional to mouse movement direction and speed
2. Dye: at the mouse position with a color

For the velocity, I'll track mouse position changes between frames. If the mouse moved by (dx, dy) pixels since last frame, I'll add:
- u_add = dx * forceScale (converted to grid coords)
- v_add = dy * forceScale

For the dye, I'll add it at the current mouse position. The dye color will be based on:
- Hue: based on mouse position and global time
- Inject in a small radius

For the hue generation, I'll use:
```javascript
hueToRgb(h) {
    let r, g, b;
    h = ((h % 1) + 1) % 1; // wrap to [0,1)
    let i = Math.floor(h * 6);
    let f = h * 6 - i;
    let q = 1 - f;
    switch(i % 6) {
        case 0: r = 1; g = f; b = 0; break;
        case 1: r = q; g = 1; b = 0; break;
        case 2: r = 0; g = 1; b = f; break;
        case 3: r = 0; g = q; b = 1; break;
        case 4: r = f; g = 0; b = 1; break;
        case 5: r = 1; g = 0; b = q; break;
    }
    return [r, g, b];
}
```

Now I'm truly going to write the complete HTML file.

I'll structure the code clearly:

1. HTML boilerplate
2. CSS (full screen, no overflow)
3. JavaScript:
   a. Utility functions (HSV to RGB, IX index)
   b. FluidSim class
   c. Mouse tracking
   d. Main loop
   e. Event handlers

Let me write it all out now. This is going to be a substantial file.

About the offscreen canvas rendering: I'll create an ImageData for the offscreen canvas (N×N pixels), fill it from the dye arrays, then putImageData. Then on the main canvas, I'll draw the offscreen canvas scaled to fill the screen. The main canvas will have `imageSmoothingEnabled = true` for bilinear filtering.

Actually, I need to be careful: `createImageData(N, N)` creates an N×N ImageData. Each pixel is 4 bytes (RGBA). I'll set RGB from the dye and alpha to 255.

Let me also add a subtle background color/gradient instead of pure black. Maybe a very dark blue-black.

For the background, I'll fill the main canvas with a dark color before drawing the scaled simulation on top. Or I can just use a dark CSS background and the simulation canvas will be transparent where there's no dye.

Hmm, actually the dye values are always non-negative. Where there's no dye, the color is (0,0,0) = black. So the background is black. That looks fine. I can add a very subtle dark blue tint.

OK, writing now. Final answer incoming.

Let me think about one more detail: the "addForce" function should spread the force to a neighborhood. Using just one cell creates sharp features. I'll spread to a 3x3 or 5x5 area using a Gaussian-like weight.

Actually, the semi-Lagrangian advection naturally smooths things out, so even point forces look OK. But for the dye, I definitely want to spread it to a small radius.

For dye injection, I'll add dye to all cells within a radius of ~2-3 cells from the mouse position, with the amount decreasing with distance.

For velocity, I'll add it to a 3x3 area around the mouse position.

OK here's my plan for the `addAt` function:

```javascript
addAt(x, y, fx, fy, r, g, b, radius = 2) {
    let i0 = Math.max(1, Math.floor(x - radius));
    let i1 = Math.min(N, Math.floor(x + radius));
    let j0 = Math.max(1, Math.floor(y - radius));
    let j1 = Math.min(N, Math.floor(y + radius));
    
    for (let j = j0; j <= j1; j++) {
        for (let i = i0; i <= i1; i++) {
            let dx = i - x;
            let dy = j - y;
            let dist = Math.sqrt(dx*dx + dy*dy);
            if (dist > radius) continue;
            let weight = Math.exp(-dist*dist / (radius*radius*0.5));
            let idx = IX(i, j);
            this.u[idx] += fx * weight;
            this.v[idx] += fy * weight;
            this.dr[idx] += r * weight;
            this.dg[idx] += g * weight;
            this.db[idx] += b * weight;
        }
    }
}
```

This is clean and effective.

Now let me truly write the complete code. I'm going to focus and produce the output.

One concern: I'm going to have a long swap section after advection where I swap prev and current arrays. This can be done efficiently by swapping references.

For the diffusion, I'll use Jacobi iteration:
- Read from `x0` (previous state)
- Write to `x` (current state)
- Swap and repeat

The standard pattern:
```
for (iter = 0; iter < ITER; iter++) {
    for (all cells) {
        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);
    }
    setBounds(x);
}
```

For advection:
```
for (all cells) {
    // backtrace
    x = i - dt*N*u[IX(i,j)];
    y = j - dt*N*v[IX(i,j)];
    // clamp
    // interpolate from d0 at (x,y)
    d[IX(i,j)] = interpolated value;
}
setBounds(d);
```

OK I really need to stop planning and write the code. Let me do it now.

Let me think about the overall structure of the step function:

```javascript
step(dt) {
    // 1. Velocity step
    // Diffuse velocity
    this.diffuse(1, this.uPrev, this.u, this.viscosity, dt); // diffuse u
    this.diffuse(2, this.vPrev, this.v, this.viscosity, dt); // diffuse v
    
    // Project
    this.project(this.uPrev, this.vPrev, this.u, this.v);
    
    // Advect velocity
    this.advect(1, this.u, this.uPrev, this.uPrev, this.vPrev, dt);
    this.advect(2, this.v, this.vPrev, this.uPrev, this.vPrev, dt);
    
    // Project again
    this.project(this.u, this.v, this.uPrev, this.vPrev);
    
    // 2. Dye step
    // Diffuse dye
    this.diffuse(0, this.drPrev, this.dr, this.dyeDiff, dt);
    this.diffuse(0, this.dgPrev, this.dg, this.dyeDiff, dt);
    this.diffuse(0, this.dbPrev, this.db, this.dyeDiff, dt);
    
    // Advect dye
    this.advect(0, this.dr, this.drPrev, this.u, this.v, dt);
    this.advect(0, this.dg, this.dgPrev, this.u, this.v, dt);
    this.advect(0, this.db, this.dbPrev, this.u, this.v, dt);
    
    // Fade dye
    for (let i = 0; i < this.size; i++) {
        this.dr[i] *= this.dyeFade;
        this.dg[i] *= this.dyeFade;
        this.db[i] *= this.dyeFade;
    }
}
```

Wait, there's an issue. After `diffuse(1, uPrev, u, ...)`:
- Input: current u array
- Output: uPrev array (diffused version)

Then `advect(1, u, uPrev, uPrev, vPrev, dt)`:
- Input: uPrev (diffused field to advect)
- Output: u (final advected field)
- Uses uPrev and vPrev as velocity

So the flow is:
1. u, v → diffuse → uPrev, vPrev
2. uPrev, vPrev → project → (uPrev, vPrev become divergence-free)
3. uPrev, vPrev → advect → u, v
4. u, v → project → (u, v become divergence-free)

Actually wait, the standard Stam approach is:
1. Add forces
2. Advect velocity (semi-Lagrangian)  
3. Diffuse velocity
4. Project

But I had the order backwards. Let me reconsider.

Standard steps from Stam's paper:
1. Add force: w₁ = w + Δt f
2. Advect: w₂ = advect(w₁)
3. Diffuse: w₃ = diffuse(w₂)
4. Project: w₄ = project(w₃)

Where w = (u, v) is the velocity field.

Actually, I think the order doesn't matter much for visual quality. Many implementations swap steps 2 and 3. Let me use: diffuse → advect → project. That's fine.

But for the actual arrays, I need to be careful about which is the source and which is the destination. Let me trace through:

Initial state: u, v have current velocity

Step 1 - Diffuse:
- Source: u (current velocity)
- Destination: uPrev (diffused velocity)
- After: uPrev = diffuse(u, visc)
- Similarly vPrev = diffuse(v, visc)

Step 2 - Project:
- Input: uPrev, vPrev (diffused velocity)
- We solve for pressure to make (uPrev, vPrev) incompressible
- After project: uPrev, vPrev are divergence-free

Step 3 - Advect:
- Source: uPrev, vPrev (projected velocity from step 2)
- Destination: u, v
- Uses uPrev, vPrev as advection velocity
- After: u = advect(uPrev, using uPrev/vPrev as vel)
- After: v = advect(vPrev, using uPrev/vPrev as vel)

Step 4 - Project:
- Input: u, v
- After: u, v are divergence-free
- These are now ready for the next frame

Wait, but in step 3, I'm using uPrev and vPrev as the velocity field for advection (the field that tells us how the fluid moves). After step 2, uPrev/vPrev are divergence-free, which is good. But I'm advecting uPrev using uPrev as the transport velocity. That means the field is being transported by itself, which is correct for the momentum equation (the nonlinear term u·∇u).

After advection, the result goes into u, v. Then we project to ensure incompressibility.

This is correct! Now let me also trace through the dye:

Step 5 - Diffuse dye:
- Source: dr, dg, db (current dye)
- Destination: drPrev, dgPrev, dbPrev
- After: drPrev = diffuse(dr, dyeDiff)

Step 6 - Advect dye:
- Source: drPrev (diffused dye)
- Destination: dr
- Uses u, v (final velocity from step 4) as transport velocity
- After: dr = advect(drPrev, using u/v)

This makes sense! The dye is advected by the final (incompressible) velocity field.

Now, for the `diffuse` function, it needs to know the type for boundary conditions:
- type 0: scalar (dye) - zero gradient at boundaries
- type 1: u velocity - no-slip at vertical walls
- type 2: v velocity - no-slip at horizontal walls

And it does Jacobi iterations to solve the diffusion equation.

For the `advect` function:
- type 0: scalar (dye)
- type 1: u velocity
- type 2: v velocity

OK, now I have a complete understanding. Let me write the code.

Let me also handle the case where dt is zero or very small (first frame).

Alright, here's the final complete code. I'm going to be very precise.

Oh wait, one more issue: in the diffuse function, when I do multiple Jacobi iterations, I need to swap buffers between iterations. The standard way:

```
for (k = 0; k < ITER; k++) {
    for (each cell)
        x_new[IX(i,j)] = (x0[IX(i,j)] + a*(...)) / (1+4*a);
    setBounds(x_new, b);
    // For next iteration, x_new becomes x_old
    // Need to swap or copy
}
```

This requires either swapping references or copying. I'll use two temporary buffers and swap.

Actually, let me simplify. In each diffuse call, I have a source array (x0) and a destination array (x). I'll iterate:
```
for (k = 0; k < ITER; k++) {
    for (each cell)
        x[IX(i,j)] = (x0[IX(i,j)] + a*(x[IX(i-1,j)]+...)) / (1+4*a);
    setBounds(x, b);
    // Swap x and x0 for next iteration
}
```

But swapping means I need to swap array references. I can do `[x, x0] = [x0, x]` at the end of each iteration. The final result might end up in either x or x0 depending on the number of iterations. Let me handle this carefully.

Actually, the standard Stam implementation does:
```
for (k = 0; k < ITER; k++) {
    for (each cell)
        x[IX(i,j)] = (x0[IX(i,j)] + a*(x[IX(i-1,j)]+...)) / (1+4*a);
    setBounds(x, b);
    // swap
    [x, x0] = [x0, x];
}
// After the loop, result is in x0 (since we swapped after the last iteration)
// Note: if ITER is even, result is in x0, but we swap one more time at end
```

Hmm, this is getting complicated. Let me handle the swapping differently.

Simpler approach: after each iteration, copy x into x0 (or just keep track and swap at the end if needed).

Actually, the simplest approach: always write to x, and at the end of each iteration, copy x to x0:
```
for (k = 0; k < ITER; k++) {
    for (each cell)
        x[IX(i,j)] = (x0[IX(i,j)] + a*(x[IX(i-1,j)]+...)) / (1+4*a);
    setBounds(x, b);
    x0.set(x); // copy x to x0
}
```

This is slightly slower due to the copy but very simple. With N=100, copying 10K floats 20 times is 200K copies - negligible.

But wait, in my implementation, the "x0" is the original unfiltered field, and "x" is the filtered result. The Jacobi method needs the original field in each iteration... no, it uses the current estimate from neighboring cells.

Actually, let me re-read Stam's paper. In the diffusion solver:

The equation is (I - a∇²)x = x0, where x0 is the original field and x is the diffused field.

The Jacobi iteration for this is:
x_new[i,j] = (x0[i,j] + a*(x[i-1,j] + x[i+1,j] + x[i,j-1] + x[i,j+1])) / (1+4a)

Where x0 is the ORIGINAL field (before diffusion), not the previous iteration's result. The x on the RHS is from the previous iteration.

So in each iteration:
- x0: constant, the original unfiltered field
- x: the current estimate, gets updated

So the pattern is:
```
// x0 is the original field, x starts as a copy of x0
x.set(x0); // initial guess
for (k = 0; k < ITER; k++) {
    for (each cell)
        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);
    setBounds(x, b);
}
```

This is correct! x0 stays constant, x is updated each iteration. No swapping needed.

Similarly for the pressure projection: div stays constant, p is updated each iteration.

Excellent, this simplifies things.

OK NOW I'm going to write the code. For real this time. Let me produce the final output.

I need to be careful about the indexing in the diffuse and advect loops. The interior is [1..N] and boundaries are [0] and [N+1].

In the loop:
```javascript
for (let j = 1; j <= N; j++) {
    for (let i = 1; i <= N; i++) {
        let idx = IX(i, j);
        // ...
    }
}
```

Now for rendering, the offscreen canvas is N×N pixels, mapping to cells 1..N (interior cells). So pixel (i-1, j-1) in the offscreen canvas corresponds to cell (i, j) in the simulation.

```javascript
render() {
    let data = this.imageData.data;
    for (let j = 1; j <= N; j++) {
        for (let i = 1; i <= N; i++) {
            let idx = IX(i, j);
            let px = (i - 1);
            let py = (j - 1);
            let pi = (py * N + px) * 4;
            
            let r = Math.min(255, this.dr[idx] * 255);
            let g = Math.min(255, this.dg[idx] * 255);
            let b = Math.min(255, this.db[idx] * 255);
            
            data[pi] = r;
            data[pi+1] = g;
            data[pi+2] = b;
            data[pi+3] = 255;
        }
    }
    this.renderCtx.putImageData(this.imageData, 0, 0);
}
```

Then on the main canvas:
```javascript
ctx.drawImage(fluid.renderCanvas, 0, 0, canvas.width, canvas.height);
```

This upscales with bilinear filtering (since imageSmoothingEnabled defaults to true).

Alright, I'm going to write everything now in one go. Here is the complete index.html:

Actually, I realize there might be an issue with `putImageData` - it bypasses all compositing and scaling. So I need to putImageData on the offscreen canvas first, then drawImage from the offscreen canvas to the main canvas.

That's what I described above. The offscreen canvas is N×N. I putImageData on it. Then I drawImage from it to the main canvas scaled up. This works.

Let me also add a slight background color. Before drawing the simulation, I'll fill the main canvas with a very dark color (like #0a0a15).

Now, for the mouse interaction, I need to convert mouse pixel coordinates to grid coordinates:
```javascript
gridX = (mouseX / canvas.width) * N + 0.5; // +0.5 because cell centers start at 1
gridY = (mouseY / canvas.height) * N + 0.5;
```

Wait: the grid interior cells are at x=1..N. The canvas pixel 0 maps to the left edge of cell 0 (which is the boundary). Canvas width maps to the right edge of cell N+1. So:
```javascript
gridX = (mouseX / canvas.width) * (N + 2); // maps to [0, N+2]
```

But I want to map to interior cells [1, N]. Hmm, let me think...

Canvas goes from 0 to canvasWidth. Grid goes from 0 to N+1 (cells). Cell centers are at 0.5, 1.5, ..., N+0.5.

Mapping: `gridX = (mouseX / canvasWidth) * (N + 2)`

This maps mouseX=0 to gridX=0 and mouseX=canvasWidth to gridX=N+2. The interior cells (1 to N) map to grid coordinates 1 to N. So:
- gridX = mouseX / canvasWidth * N + 1

Wait:
- mouseX = 0 → gridX = 1 (left edge of leftmost interior cell)
- mouseX = canvasWidth → gridX = N + 1 (right edge of rightmost interior cell)

Hmm, that gives range [1, N+1]. Let me use:
- gridX = mouseX / canvasWidth * (N) + 0.5

This gives:
- mouseX = 0 → gridX = 0.5 (left edge)
- mouseX = canvasWidth → gridX = N + 0.5 (right edge)

The interior cell centers are at 1, 2, ..., N. So this maps canvas properly to grid coordinates. Good.

Actually let's make it simpler:
```javascript
gridX = 1 + (mouseX / canvasWidth) * N;
gridY = 1 + (mouseY / canvasHeight) * N;
```

This maps [0, canvasWidth] → [1, N+1], which slightly overflows. That's fine because I clamp in the addAt function.

Let me just use a nice clean mapping and clamp:
```javascript
gridX = 1 + (mouseX / canvas.width) * N;
gridY = 1 + (mouseY / canvas.height) * N;
```

OK I'm done planning. Writing the actual code now.

Let me also include a nice feature: when the user isn't dragging, maybe add some ambient slow dye at random positions to keep things interesting. Actually no, let's keep it interactive-only. The user drags to create the dye.

But I could add some initial dye patterns so the canvas isn't completely black on load. Like a few colored blobs.

Actually, let me add a simple startup: when the page loads, inject some dye at a few random positions to get things started.

OR, I can create an initial velocity field with some swirls so the fluid is "alive" even before the user interacts. Let me add a few initial velocity vortices.

Actually, you know what, let me keep it simple. The page loads with a black canvas and a subtle hint text "Drag to create fluid". When the user drags, colorful dye appears and swirls through the fluid. The dye naturally diffuses and fades.

One thing I should be careful about: the dye values should be capped so they don't grow unboundedly. The fade factor helps, but if the user keeps dragging in the same spot, dye could accumulate. I'll cap each channel at some max value (e.g., 5.0, which when scaled to 0-255 gives 255).

OK, I truly need to write this now. Let me produce the complete code. I'll aim for clean, well-organized JavaScript.

Let me also add keyboard controls: press 'R' to reset, 'Space' to pause/unpause.

Actually no, let me keep it simple. Just the fluid simulation with mouse drag.

Here goes - the final implementation:

I realize I should also handle the case where the mouse is down but not moving (stationary). In this case, I should still inject some dye so the user sees something. Let me add a small amount of dye even when stationary, and velocity only when moving.

Alright. Code time. For real. No more changes.

Let me think about performance one more time. With N=100:
- diffuse (velocity): 20 iters × 2 components × 10K cells = 400K operations
- project: 30 iters × 10K cells = 300K operations  
- advect (velocity): 2 components × 10K cells = 20K operations
- project again: 300K
- diffuse (dye): 20 iters × 3 channels × 10K cells = 600K
- advect (dye): 3 channels × 10K cells = 30K
- fade: 3 × 10K = 30K
Total: ~1.7M operations per frame

For comparison, a simple operation in JS takes ~1ns, so 1.7M operations ≈ 1.7ms. Plus overhead. Easily within 16ms budget for 60fps. Should be very smooth.

Great, let me write the complete code now.

One more thing: for the rendering, I need to handle the case where the dye values might be very small (near zero) and should render as black. The Math.min(255, val*255) handles this correctly. But I should also ensure values aren't negative.

Also, maybe I should apply a small gamma correction or use sqrt to make the darker areas look better. Something like:
```javascript
let intensity = Math.sqrt(dr[idx]); // or some other curve
let r = Math.min(255, intensity * 255);
```

This makes low-density areas more visible. I'll use this.

OK HERE IS THE CODE. Writing now:

Actually wait, I want to add one more visual enhancement. The velocity field can be used to create a "glow" effect - areas with high velocity get a slight boost in brightness. This makes the fluid look more energetic. I can compute the velocity magnitude at each cell and add it to the dye intensity.

```javascript
let velMag = Math.sqrt(u[idx]*u[idx] + v[idx]*v[idx]);
let boost = 1 + velMag * 0.5; // up to 50% boost for high velocity
```

This is a nice touch. Let me include it.

OK NOW FOR REAL THE CODE:

I'll produce the complete file with all these features. Let me write it out carefully.

The workspace is empty. Let me create the index.html file with the fluid simulation.

Let me write the complete code now. I'll be very careful with the implementation.

I'll use a grid of N=100 interior cells for good performance and visual quality.

Let me write out the complete file:

The file has been written. Let me verify it's correct by reading it back and also double-check the logic.

Let me review the key parts of the simulation:

1. Grid: N=120 interior cells, SIZE = (122)*(122) = 14884
2. Jacobi iterations: 16 for diffusion, 40 for projection
3. The simulation step follows the standard Stam approach
4. Mouse interaction converts pixel coords to grid coords and injects both velocity and dye
5. Rendering uses an off-screen canvas (N×N) that gets upscaled

Let me verify the file was created properly.

The file is 437 lines and looks correct. Let me verify the simulation logic is sound by reviewing the key parts once more.

The file is complete. Let me verify there are no syntax issues or logic bugs by checking a few key sections. Let me also make sure the index.html is the only file and self-contained.

The JavaScript syntax check passed. The file is complete and valid. Let me provide a summary of what was built.