The user wants a real-time fluid simulation on a full-screen canvas. Dragging the mouse pushes swirling colored dye through the fluid. It needs to be a single index.html file with inline CSS and JS, no external resources.

This is a classic Jos Stam-style stable fluid solver task. We need:
1. Full-screen canvas
2. Fluid velocity field simulation
3. Dye/color advection
4. Mouse interaction that adds velocity and dye
5. Real-time rendering

Key aspects of a stable fluid solver:
- Grid of cells with velocity (u, v) and density/color
- Stable solver: diffuse + project + advect + project
- Jacobi iteration for diffusion and pressure projection
- Boundary conditions
- Dye advection using velocity field

We need to be careful about performance. A JavaScript implementation at reasonable resolution should work.

Since no external resources, we write all WebGL/canvas code in JS. We can use HTML5 Canvas 2D for simplicity and to avoid WebGL complexity, but it might be slower for high res. WebGL shaders would be much faster but more complex. Either way, single file.

Given the requirement: "real-time" and "single-page web app". We can use Canvas 2D with a moderate grid size (e.g. 128x128 or 256x256) and render scaled up. This yields a nice simulation.

But to make it look nicer and faster, we could implement with WebGL and framebuffers. However, the prompt says "no external resources, CDNs or network requests". WebGL with inline shaders is fine: no external resources. But Canvas 2D is simpler and self-contained. Which should we choose?

A JS-based solver on CPU is a known implementation. For 256x256 grid, Jacobi iterations (10-20 per step) might be borderline but okay on modern devices. Let's think about typical performance:
- 256x256 = 65536 cells
- Each step: 2 diffuse passes (u, v) each 20 iterations = 40 * 65536 = 2.6M operations
- 2 advect passes (u, v) = 131K
- projection: divergence + 20 Jacobi + gradient subtract = ~40 * 65536 = 2.6M
- dye: 3 color channels diffuse + advect = ~3*(1.3M+131K) ~ 4.3M
- Total ~9M JS ops per frame at 256x256. At 60fps: 540M ops/sec. JS might manage this in ~5-10ms? On modern engines, maybe. But a low-end device may struggle.

We might reduce resolution to 128x128 (16384 cells) for CPU stability. But visual less detailed.

Alternative: WebGL solver. It uses GPU and can handle 512x512 or even higher real-time. It requires writing vertex/fragment shaders and managing textures/framebuffers. The single file can include all shader code as strings. More complex but definitely "expert coding" quality. The user didn't specify implementation. We can use WebGL if we can keep it robust.

Given the request is somewhat simple and expects a single HTML file, Canvas 2D CPU solver is a common correct answer. Many implementations exist. Let's create a Jos Stam stable fluids CPU implementation.

But to improve visuals and performance, we can implement a grid-based simulation at e.g. 128x128 simulation, render interpolated to full screen via `drawImage(canvas, 0, 0, width, height)` or simply set canvas size small and CSS scaling. We can use a low-res simulation canvas and display via a second canvas scaled, or just use one canvas smaller internal resolution scaled via CSS. The canvas element width/height set to grid resolution, CSS width/height 100%, so browser scales with smoothing; but pixelated or smooth? We can set `image-rendering: pixelated` for retro; smooth interpolation via default. For dye swirls, smooth scaling may look nicer. Let's set canvas internal resolution maybe 1/4 of screen or fixed 512? We can determine based on screen size to keep cell count reasonable.

Better: set simulation grid to fixed e.g. 256 x (256 * aspect). Use lower resolution for CPU. Let's make grid height 256, width = Math.round(256 * aspect), max maybe 384. For 16:9 aspect, width ~455, cells ~116k. That's bigger. Let's set grid to 128 height, width = 128*aspect; cells ~29k. That is safe.

Actually we can make grid resolution dynamic: choose grid size such that total cells ~50k, e.g. base 200, width base*aspect. Maybe fine. For performance, try 128 base and see. But initial simulation may use a grid of 256 columns? Let's compute: if gridWidth=256, gridHeight=144 (for 16:9), cells=36864. Very good. Use gridHeight 180 for width 320 => 57600 cells. Should be okay. The simulation resolution can be 320x180, display scaled to full screen. That gives enough detail for swirls. I think 320x180 is safe for CPU.

Now, implementation details.

Fluid solver CPU in JS:
We'll create a `Fluid` class with:
- Nx, Ny: grid dimensions (not including boundary cells? We can include boundaries in arrays but solver loops interior; boundary sets edges).
- size = Nx*Ny
- arrays: u0, v0, u1, v1, r0, g0, b0, r1, g1, b1 or one Uint32 image.
- density/dye arrays for each color channel.

Functions:
- step(dt): diffuse velocities, project, advect velocities, project; diffuse dye, advect dye.
- addDensity(x, y, r, g, b, amount)
- addVelocity(x, y, dx, dy)
- render(ctx, imageData): compose dye into imageData and putImageData.

Mouse interaction:
- track mouse positions and deltas
- on mousemove/ pointermove while down: 
  - map client coords to grid coords.
  - add velocity scaled by delta.
  - add dye based on cursor movement amount/velocity; color cycling hue.
- Mouse may move fast; we can make a small spray, but maps to grid cell.

Boundary conditions:
- For velocity, boundary should reflect (set edge value to negative of adjacent interior) with sign depending on axis.
- For dye, boundary clamps to zero or same? Typically set_bound(b=0) for density, and for velocity: x component negates on vertical edges, y component negates on horizontal edges. Implementation: loop edges, set arr[boundary] = (axis == 1 ? -arr[interior] : arr[interior])? Actually for u (x velocity), set u[edge] = -u[interior] for vertical edges? Let's recall standard: 
```
for i in 1..N:
  x[0,i] = b==1 ? -x[1,i] : x[1,i]
  x[N+1,i] = b==1 ? -x[N,i] : x[N,i]
for j in 1..N:
  x[j,0] = b==2 ? -x[j,1] : x[j,1]
  x[j,N+1] = b==2 ? -x[j,N] : x[j,N]
```
But b=1 applies to u when setting horizontal boundaries? Wait, in Jos Stam code, for velocity array u (horizontal component), they call set_bnd(N, 1, u); for v, set_bnd(N, 2, v); for density, set_bnd(N, 0, x). In that code, `b==1` negates on left/right boundaries; `b==2` negates on top/bottom boundaries. So for u, b=1 means vertical walls reverse x velocity (normal flow). For v, b=2 means horizontal walls reverse y velocity. That is physically correct.

We need to handle grid dimensions possibly non-square. The N in code for square; we need Nx, Ny. The boundary loops need to adapt.

Implementation of linear_solver:
```
function lin_solve(b, x, x0, a, c, iter, Nx, Ny) {
  const invC = 1/c;
  for (let k=0; k<iter; k++) {
    for (let j=1; j<=Ny; j++) {
      for (let i=1; i<=Nx; i++) {
        const idx = IX(i,j);
        x[idx] = (x0[idx] + a*(x[IX(i+1,j)] + x[IX(i-1,j)] + x[IX(i,j+1)] + x[IX(i,j-1)])) * invC;
      }
    }
    set_bnd(b, x, Nx, Ny);
  }
}
```
Where IX(i,j) returns i + j*(Nx+2) perhaps.

Diffusion:
```
function diffuse(b, x, x0, diff, dt, iter, Nx, Ny) {
  const a = dt * diff * Nx * Ny; // for square N. For rectangle, maybe use sqrt? Stable if a not too high.
  lin_solve(b, x, x0, a, 1+4*a, iter, Nx, Ny);
}
```
For non-square, `a = dx^2 / (nu dt)?` in standard. But for visual, we can use `const a = dt * diff * Nx * Ny;` This is approximated. We could compute based on min(Nx,Ny). Fine.

Projection:
```
function project(u, v, p, div, iter, Nx, Ny) {
  for interior:
    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)]) / h; // h = 1?? Depending grid. In original with N, div = -0.5*h*(...), then p solve, then u -= 0.5*(p[i+1]-p[i-1])/h. If h=1/N?? Original uses N as cell count and grid [0,N], h=1/N. So div = -0.5*N*(...). Then after pressure solve, u -= 0.5*(p[i+1]-p[i-1])*N. If we use h=1 / min(Nx,Ny), or simply use 1.0? For the incompressibility effect, scale can be adjusted: using N factor ensures convergence but divergence field magnitude scales with N. We can use `scale = 1.0` to avoid instabilities at high resolution? Actually the standard code uses `h = 1/N` in continuous space, but if we use normalized grid 0..1. Let N = min(Nx,Ny). Then use `scale = 0.5 * N`. That works reasonably across resolutions. Then lin_solve for pressure with a=1, c=4? Actually poisson: `p = (div + p_neighbors)/4`. In original, `lin_solve(0, p, div, 1, 4, iter, Nx, Ny)`. Then `u[i,j] -= 0.5*(p[i+1]-p[i-1])*N`. That gives divergence-free field.

Advection:
```
function advect(b, d, d0, u, v, dt, Nx, Ny) {
  dt0 = dt * N; // if N used as grid size. Actually original: dt0 = dt*N. For rectangle, maybe dt0x = dt*Nx, dt0y = dt*Ny.
  for j=1..Ny, i=1..Nx:
    x = i - dt0x * u[idx]; y = j - dt0y * v[idx];
    clamp x,y in [0.5, Nx+0.5] etc.
    i0 = floor(x), i1=i0+1; s1=x-i0, s0=1-s1; etc.
    interpolate from d0 at four corners.
}
```
The original uses `dt0 = dt * N` because grid coordinates range 1..N and velocity is in grid cells per second? We can simulate time step such that dt0 ~ small. We'll determine `dt = 0.1` maybe, and `dt0x = dt*(Nx)`, dt0y = dt*(Ny). That's large if dt=0.1 and Nx=320: dt0=32, huge. Need careful scaling. Better to interpret velocities in cells per frame, and dt=1. Then dt0x = dt (maybe 1) for square? Actually if velocities are in grid cells per unit time, advection trace back: position = current - velocity * dt. If grid cell index spacing =1, then dt0 = dt. The original code used dt0 = dt*N because velocity field is scaled to simulation [0,1]. But if velocities are cell units, dt0 = dt. This is simpler: dt=0.2, dt0=0.2.

Velocity mouse input: add `force = (mouse dx, dy) * factor`. Dye add at cell.

We need to avoid instabilities from high velocity; we can clamp velocity. Add diffusion to velocity and dye (viscosity, dye diffusion). Keep dt small or use sub-steps.

Time stepping: on each animation frame, compute `dt` from elapsed; maybe clamp to max 0.05, and run multiple simulation substeps to keep stable (e.g. 4 substeps). We can use fixed dt per substep.

Let's design:

Grid setup:
```
let Nx = Math.floor(window.innerWidth / 4); // maybe /3
let Ny = Math.floor(window.innerHeight / 4);
```
Total grid cells ~ width*height/16. For 1920x1080 => 480x270 = 129600 cells. For CPU heavy but maybe okay if iterations low? Hmm 129k cells * solver per frame. Could be heavy. Use larger divisor e.g. 6 => 320x180=57.6k. Use divisor based on `Math.sqrt(targetCells/(aspect))`. Let's target ~40000-60000 cells. Compute base = Math.sqrt(target / aspect); Nx = Math.round(base*aspect); Ny=Math.round(base). For 16:9 target 50000: base = sqrt(50000/(16/9))=sqrt(28125)=167.7; Nx=298; Ny=168. Good. For 4:3, Nx=258,Ny=194. Use target = 50000. That's safe.

We can recompute grid on resize. But re-init arrays on every resize may reset simulation. Could be fine; or we can only resize display canvas and keep sim grid?? Simpler: reinitialize on resize.

Rendering:
- Create ImageData(Nx, Ny). For each interior cell, set RGB from dye arrays, clamp 0-255. Put image data to a canvas of size Nx,Ny, then drawImage to display canvas? We can use a single canvas with CSS sizing but `putImageData` uses canvas pixel dimensions. If canvas CSS is full-screen but internal width=window.innerWidth/scale, then putImageData writes to the scaled? Actually canvas internal resolution is independent of CSS. If we set canvas.width=320, canvas.height=180 but CSS width/height 100%, then the image appears scaled. `putImageData` writes at actual pixels. So we can have one canvas full-screen in CSS with low internal resolution. When mouse events, use `getBoundingClientRect` scale factor to map to grid coords. That's efficient.

We need to be careful: when CSS scales canvas, mouse coordinates must be converted accounting for css-to-pixel ratio. `const rect = canvas.getBoundingClientRect(); const gx = Math.floor((e.clientX - rect.left) / rect.width * Nx);` etc. It handles non-square aspect and CSS full-screen.

Use `image-rendering: auto` for smoothing (default). Might produce bilinear scaling -> smooth swirls. Good. Optionally `pixelated` maybe a more "simulation" look. Smooth is good for dye swirls. But bilinear stretch blurs details; acceptable.

Dye colors:
- Cycle hue based on time or based on mouse movement angle.
- Add dye at mouse position with gaussian falloff over a few cells. Map hue to RGB using HSL function. Let's implement conversion.
- Also perhaps dye color depends on velocity magnitude or direction.

Mouse interactions:
- Add velocity when moving. For a nicer swirl, maybe add a force with rotational component around a point? Like if moving fast, add angular vorticity. This can create swirls. A simple way: compute velocity vector (dx, dy) over cursor path, add that to u/v at cursor and in small radius. Additionally, add a tangential component perpendicular to drag direction to induce swirl. Or use "vorticity confinement" not needed. Actually simply adding velocity while fluids advects itself creates swirls indirectly due to incompressibility. So plain dragging is fine.
- Dye amount proportional to speed.
- Trailed injection: when moving fast, inject at interpolated positions between old and new mouse grid coordinates so no gaps.

Implementation plan:
- index.html with inline style and script.
- Style: html, body { margin:0; width:100%; height:100%; overflow:hidden; background:#000; } canvas { display:block; width:100%; height:100%; }
- Script: 
  - canvas = document.getElementById('c')
  - ctx = canvas.getContext('2d', { alpha: false }) maybe.
  - Grid setup function `resizeGrid()`.

Need to handle high-DPI? Canvas CSS full-screen, internal resolution smaller. OK.

Data structures: use Float32Array for velocity/dye arrays. Since values can exceed 1, fine. Dye density can be >255 clamped.

Fluid class functions:
```
function Fluid(nx, ny) {
  this.nx = nx; this.ny = ny;
  const size = (nx+2)*(ny+2);
  this.u0 = new Float32Array(size);
  this.v0 = new Float32Array(size);
  this.u1 = new Float32Array(size);
  this.v1 = new Float32Array(size);
  this.r0 = new Float32Array(size);
  this.g0 = new Float32Array(size);
  this.b0 = new Float32Array(size);
  this.r1 = new Float32Array(size);
  this.g1 = new Float32Array(size);
  this.b1 = new Float32Array(size);
}
```
Default values zero. Boundary zero.

`IX(i,j)`: `return i + j*(this.nx+2)`, with 0 <= i <= nx+1.

`set_bnd(b, x)`: edges and corners.

Implementation of boundaries for rectangle:
```
for (let i=1; i<=nx; i++) {
  x[IX(i,0)] = b==2 ? -x[IX(i,1)] : x[IX(i,1)];
  x[IX(i,ny+1)] = b==2 ? -x[IX(i,ny)] : x[IX(i,ny)];
}
for (let j=1; j<=ny; j++) {
  x[IX(0,j)] = b==1 ? -x[IX(1,j)] : x[IX(1,j)];
  x[IX(nx+1,j)] = b==1 ? -x[IX(nx,j)] : x[IX(nx,j)];
}
x[IX(0,0)] = 0.5*(x[IX(1,0)] + x[IX(0,1)]);
x[IX(0,ny+1)] = 0.5*(x[IX(1,ny+1)] + x[IX(0,ny)]);
x[IX(nx+1,0)] = 0.5*(x[IX(nx,0)] + x[IX(nx+1,1)]);
x[IX(nx+1,ny+1)] = 0.5*(x[IX(nx,ny+1)] + x[IX(nx+1,ny)]);
```
For velocity arrays, b=1 for u, b=2 for v, zero for dye.

`lin_solve(b, x, x0, a, c, iter)`:
```
const cRecip = 1/c;
for (let k=0; k<iter; k++) {
  for (let j=1; j<=ny; j++) {
    const row = j*(nx+2);
    const rowp = row + (nx+2);
    const rowm = row - (nx+2);
    for (let i=1; i<=nx; i++) {
      x[row+i] = (x0[row+i] + a*(x[row+i+1] + x[row+i-1] + x[rowp+i] + x[rowm+i])) * cRecip;
    }
  }
  set_bnd(b, x);
}
```
Manual indexing for speed, no IX function in inner loop.

`diffuse(b, x, x0, diff, dt, iter)`:
```
const a = dt * diff * nx * ny; // As original with N^2
```
For rectangular, use grid area; higher resolution -> higher a; lin_solve coefficient c = 1+4*a. Use 4 neighbors. This ensures that to diffuse at same physical rate independent of resolution, a scaled by Nx*Ny. Good. We can set diff = 0.0001? Actually if a = dt*diff*Nx*Ny. For N~300, dt ~ .02 per substep, a = 0.02*diff*50000 = diff*1000. If diff=0.0001, a=0.1. So diff can be very small. We can instead not square N; use a = dt * diff * (nx*ny)? This is fine. Let's use `a = dt * diff * (nx * ny)`. Actually original uses N*N. For dye, we set diff to e.g. 0.00001. With nx*ny=50000, dt=0.2, a=0.1. Good.

`project(u, v, p, div, iter)`:
```
// div
const h = 1/Math.max(nx,ny)? Actually scale factor N. Use scale = Math.max(nx,ny); but original used N as number of cells.
```
Let's choose `scale = Math.max(nx,ny)`. Wait original with an NxN grid uses h = 1/N and scale = 0.5*N in div and gradient. For rectangle, use `N = Math.max(nx, ny)` as reference. So:
```
for interior:
  div[idx] = -0.5 * N * (u[idx+1] - u[idx-1] + v[idx+(nx+2)] - v[idx-(nx+2)]);
  p[idx] = 0;
set_bnd(0, div); set_bnd(0, p);
lin_solve(0, p, div, 1, 4, iter);
for interior:
  u[idx] -= 0.5 * N * (p[idx+1] - p[idx-1]);
  v[idx] -= 0.5 * N * (p[idx+(nx+2)] - p[idx-(nx+2)]);
set_bnd(1,u); set_bnd(2,v);
```
But if N large (300), scale=150, velocities might be huge from projection unless velocity is normalized small. Mouse velocities are maybe 0.1-1. Div could be ~150*(...). Might be too much. Original mouse velocities were scaled by N factor? They used small dt=0.1 and N maybe 64. Mouse force divided by something. It might still work because velocity arrays after projection are corrected. If scale too large relative to velocity, perturbations get magnified? Let's test thought:
- Suppose mouse adds u=5 in cell. Then div ~ -0.5*N*(5-0) = -2.5*N. Pressure solve yields p ~ (div + neighbors)/4; p magnitude ~ N. Gradient subtract: u -= 0.5*N*(p_neighbor - p)/h? Wait gradient uses scale N too, so u -= 0.5*N*(p diff). p diff ~ div/... maybe O(N). So u correction O(N^2), huge and unstable. That's because our velocity and pressure scaling inconsistent with N. In original, velocities are in simulation units, i.e., value 5 means moves 5/N of the domain per time. With dt=0.1 and advection dt0=dt*N=1 means trace back by 5 cells? Hmm.

Let's re-derive. In original code, grid coordinate j indexes 1..N, and velocity in cells per time? `dt0 = dt * N` means velocity multiplied by N to convert cells? Actually if coordinate in grid cells is 1..N, then velocity there is in cells per time, so dt0 = dt. But original uses dt0 = dt*N, suggesting velocities are in domain units [0,1] per time; multiply by N to get grid cell displacement. So mouse input maybe added in domain units. That means to move one cell, u must be 1/N. So velocities are small (e.g., 0.01). Then projection scale N * u ~ O(0.01*N) maybe a few. That is stable. So we should treat velocity in domain units and advect with dt0=dt*N.

For dye arrays, likewise domain density. Dye add small amount maybe 0.5.

Thus:
- Mouse delta per frame in client pixels; map to grid cells; convert to domain velocity: delta_grid_x / Nx (if total domain width is 1). Add to u and v with scaling factor, e.g. `du = (dx / rect.width) / dt * strength`? Actually moving mouse across the full screen in one second at 60fps should produce domain velocity ~1. If we map mouse delta in pixels to grid delta: `gdx = e.movementX / rect.width * Nx` (grid cells per frame). Domain velocity = gdx / Nx / dt = movementX / rect.width / dt. So `u = movementX / rect.width / dt * strength`. Since dt between frames maybe 1/60, u ~ movementX*60/width. If movementX ~ 10 px at width 1920, u ~0.3. That's plausible. Add scaling factor 5 -> u ~1.5 domain velocity. OK.

- Dye add: at cursor position, amount proportional to speed in domain units clamped.

- Advection: `dt0x = dt * Nx; dt0y = dt * Ny`. For dt=0.1, Nx=300 => dt0x=30; with u=1.5, trace back 45 cells? That's too large, unstable. Wait if u domain 1.5 means moves 150% of domain per time, very fast. dt0=dt*Nx huge. To keep Courant stable, we need u*dt*Nx <= 1 cell roughly: u <= 1/(dt*Nx). With dt=0.1, u<=0.033. That's very slow. But original used dt=0.1, N=64 => u<=0.156, and mouse force might be small. But we want interactive swirling, velocities can be a few cells per frame; we can use smaller dt and/or many substeps, or use advection that allows larger displacement? Stable semi-Lagrangian advection is stable for large displacements (trace back far), but projection with dt0 large maybe okay if input small.

Hmm. Let's reconsider scaling. The simulation doesn't need to correspond physical units. We can define velocities in grid cell units per time. Then advection uses `dt0=dt`. Projection scale h=1 (since grid spacing=1). Then projection div uses `0.5*(u[i+1]-u[i-1]+v[j+1]-v[j-1])`. Pressure solve p; gradient subtract `0.5*(p[i+1]-p[i-1])`. This is the original code if N=1 domain. Then mouse velocities are actual grid cells per frame. If mouse moves across screen, gdx per frame maybe a few cells; add scaled to maybe 0.5-2. Good. dt perhaps 0.2 per substep. This is much simpler and stable.

But diffusion coefficient using nx*ny may be too large. For diffusion in grid units, use `a = dt * diff * (nx*ny)` same as original? Original used N*N because h=1/N; diffusion coefficient in code `a = dt * diff * N * N`. If we use grid units h=1, then a = dt*diff. For rectangle, not N*N. So we can use a = dt*diff. Good.

Projection in grid units:
```
div[idx] = -0.5 * (u[idx+1] - u[idx-1] + v[idx+stride] - v[idx-stride]);
p solve = lin_solve(0,p,div,1,4,iter);
u[idx] -= 0.5*(p[idx+1]-p[idx-1]);
v[idx] -= 0.5*(p[idx+stride]-p[idx-stride]);
```
No N factor. This is physically consistent if pressure gradient computed per grid cell. This prevents large projection blow-up.

Advection: `x = i - dt*u[idx]; y = j - dt*v[idx];` dt maybe 0.5. This is stable since SL. Use bilinear interpolation.

Diffusion: `a = dt*diff; c = 1+4*a;` With dt=0.5, diff=0.0001 => a=5e-5, diffusion minimal. But we want some dye spread. We can set diff to 0.001. Use maybe 10-20 iterations.

The entire solver in grid units is simpler and often used in CPU implementations.

Mouse injection: For each segment between previous and current grid positions, inject velocity and dye along line with spacing <1 cell; use linear interpolation. This ensures continuous trails even with fast mouse.

Let's implement solver in grid units.

Need to ensure incompressibility projection is correct for zero boundary velocities. Dye boundary clamp/continuity? For dye, set_bnd 0 gives black edges because dye clamps zero outside? Actually `set_bnd(0, x)` sets edge equal to interior, meaning no-flux (Neumann). Then dye won't leak but smoothing near edges maybe odd. For density, original uses `b=0`, so boundary = interior, okay. We'll use same.

Velocity boundary: horizontal walls reverse v (b=2 for v), vertical walls reverse u (b=1 for u). This creates sealed box.

Implementation details:
- Use temporary arrays p and div of size. We'll allocate in Fluid class.

Class methods:
- `step(dt, viscosity, diff, iterations, substeps)`:
  For substep:
    let sdt = dt / substeps;
    // velocity
    diffuse(1, this.u1, this.u0, viscosity, sdt, iter); diffuse(2, this.v1, this.v0, viscosity, sdt, iter);
    project(this.u1, this.v1, this.u0, this.v0, iter); // Note after diffuse, u0/v0 used as temp p/div
    // swap u0<->u1? Then advect: advect(1, this.u0, this.u1, ...)
    // Actually original: diffuse x into x array? Let's design with explicit arrays.
```

velocity step:
  diffuse u0 -> uTmp (using u1 as target? Let's maintain uPrev and uCurr and swap)
  swap, so u0 is updated after diffuse.
  diffuse v0 -> vTmp and swap.
  project(u0,v0, p, div)
  advect u0 -> uTmp with velocity (u0,v0); swap.
  advect v0 -> vTmp with velocity (u0,v0) (after first advect we have updated u in u0 and old v? Need careful: original advects both using the before-advect velocity). Standard:
```
diffuse(u, u0); diffuse(v, v0);
project(u,v, p, div); // now u,v are diffused and divergence-free
advect(u, u0, u0, v0); advect(v, v0, u0, v0);
project(u, v, p, div);
```
Where u0,v0 are original before step, u,v after diffuse/project, advect uses u0,v0 (midpoint), result overwrites u,v? Actually original:
```
 diffuse(N, 1, u, u0, visc, dt );
 diffuse(N, 2, v, v0, visc, dt );
 project(N, u, v, u0, v0 );
 advect(N, 1, u, u0, u0, v0, dt );
 advect(N, 2, v, v0, u0, v0, dt );
 project(N, u, v, u0, v0 );
```
Here u0/v0 are used as prev arrays; after diffuse target is u/v. Then project using temporary u0/v0. Then advect reads u0/v0 (which now are pressure/div temps? Wait after project, u0/v0 overwritten with pressure/div, not original. Hmm standard implementation uses u0/v0 as previous before step, but project clobbers them. Actually in Stam's implementation, project returns u0/v0 as temp arrays (p and div). The advect step reads u0/v0 which are now undefined? Let's not overthink. Many implementations maintain separate arrays: u, v as current, uPrev, vPrev as prev, use temp arrays for project. The sequence:
```
diffuseVelocity(visc); // writes new velocity into prev arrays? 
projectVelocity(); // makes current divergence-free
advectVelocity(); // uses prev (projected) velocity to move current

```
Implementation from六百行：```
fluid.step() {
  diffuse(1, this.u, this.u0, visc); // u0 = diffused u
  diffuse(2, this.v, this.v0, visc);
  project(this.u0, this.v0, this.u, this.v); // result in u0/v0 divergence-free; u/v used as p/div tmp
  // now velocity divergence-free is in u0/v0; original current is in u/v (old)
  advect(1, this.u, this.u0, this.u0, this.v0); // advect u0 (now divergence-free) using u0/v0 and store in u
  advect(2, this.v, this.v0, this.u0, this.v0);
  project(this.u, this.v, this.u0, this.v0);
}
```
This is one canonical approach. Need separate arrays.

Let's define:
- `u`, `v`: current velocity
- `uPrev`, `vPrev`: previous/temp
- `r`,`g`,`b`: current dye
- `rPrev`,`gPrev`,`bPrev`: previous dye

Methods:
`step(dt)`:
```
diffuse(1, uPrev, u, viscosity, dt, iter);
diffuse(2, vPrev, v, viscosity, dt, iter);
project(uPrev, vPrev, u, v, iter); // result divergence-free in uPrev,vPrev; p,div in u,v
advect(1, u, uPrev, uPrev, vPrev, dt, iter); // u = advect(uPrev) using uPrev,vPrev
advect(2, v, vPrev, uPrev, vPrev, dt, iter);
project(u, v, uPrev, vPrev, iter); // result divergence-free in u,v; p,div in uPrev,vPrev
for each color:
  diffuse(0, rPrev, r, dyeDiff, dt, iter);
  advect(0, r, rPrev, u, v, dt, iter);
```
At the end, current dye in r/g/b. Note after diffuse dye, rPrev is diffused, then advect target r reads rPrev and u,v. Good.

Need to be careful that during velocity advect, u is overwritten for x but uses old vPrev for y; since both advect use uPrev,vPrev as velocity, no problem. After advect, u and v are updated; then final project.

Now `diffuse(b, x, x0, diff, dt, iter)`:
- `a = dt * diff * (this.nx * this.ny)`? Since we're using grid units (h=1), physically diffusion coefficient in grid units should be a = dt*diff. But to get similar spread across resolutions maybe a = dt*diff * (this.nx*this.ny)/scale? If diff is resolution-independent, original a = dt*diff*N*N because domain physical size is 1. If we want diff to correspond to domain, use `a = dt * diff * Math.max(this.nx,this.ny)**2 / 10000?` Hmm.

The user's resolution varies. It may be okay if dye spreads more on higher resolution, because sim grid cells smaller. We can make diff slightly larger. Let's use a = dt * diff * (this.nx * this.ny) * 0.0001? This is messy. Simpler: define a "diffuse factor" f = dt * diff, where diff is in "grid cell squared per frame". So diff=0.001 means small spread in grid cells. We'll tune diff to ~0.0005 maybe. For 50k grid, a small. Good enough.

I think use `a = dt * diffusion * (this.nx * this.ny) / 50000`. That makes diff resolution independent relative to domain [0,1], similar to original. If N=300, factor ~1.8. So a = dt * diffusion * 1.8. If diffusion=0.001 and dt=0.2, a=0.00036. Very small (nearly no diffusion), maybe not enough? Actually diffusion coefficient of 0.001 in domain units is small. We may want more visible spread. Set diffusion=0.005 -> a=0.0018. After 20 iterations, some spread. Good.

Actually we can tune later.

For viscosity similarly.

Advection with dt: Use `dt0x = dt` if u is grid cells per time. For dt maybe 0.2, trace back 0.2 cells. If mouse velocity u is 1 cell/frame, movement fine. Use multiple substeps? If u can be 5 cells/frame (fast drag), dt=0.2 => trace back 1 cell if substeps=5? We can set dt total per frame ~0.5 and substeps=1 maybe enough. SL handles large tracebacks but interpolates. Large u*dt might cross many cells; still stable if clamp. We'll clamp trace coordinate to interior. Use dt maybe 0.2 per step; mouse velocities clamped maybe 5; traceback 1 cell. Good.

However, time step per frame should be based on elapsed. We can run with fixed `dt = 0.2` each animation frame, independent of frame rate? That would make simulation speed depend on frame rate. Better compute `dt` from timestamp and run N substeps of fixed `sdt=0.02` or accumulate. But fixed at 60fps assumed. We can compute `elapsedDelta` and cap; run one step with dt = elapsedDelta * timeScale. To ensure stability, we can limit dt max (e.g., 0.05) and substeps = Math.max(1, Math.ceil(dt / maxStep)); sdt = dt/substeps. Then each substep uses small dt. This robust.

Set maxStep = 0.05 (20 substeps per second). If dt=0.016 (60fps), substeps=1, sdt=0.016. If frame drops dt=0.1 -> substeps=2, sdt=0.05. Good.

Now rendering:
- We need to map dye density to color. Dye density floats will accumulate and can saturate. We can clamp to 1 and multiply by 255. But to get vivid colors, maybe normalize? Simpler clamp. If user drags a lot, dye saturates white. Could add "decay" to dye each frame to prevent blow-up/whiteout. We can multiply dye by 0.995 each frame or after step. That keeps scene dynamic. Also velocity can decay to avoid perpetual motion? Maybe not, but small damping helps.

Dye decay: after step, multiply current dye arrays by 0.99 (or 1 - dt*decayRate). This prevents saturation. Also velocity damping maybe 0.99. For incompressible fluid, momentum may persist; damping helps stable.

Use a background base color maybe dark black. Dye emits color. We'll add color where mouse is and let it advect. Decay means trails fade.

Color generation:
- Use HSL to RGB. Hue rotates over time. We can use `hue = (hue + 1) % 360`; on dye injection add color based on hue. Also make color somewhat pastel? For dyes maybe bright saturated; use lightness ~0.55, saturation 1. Convert HSL to RGB.
- Dye amount = clamp(speed*factor, 0.1, 2.0). Multiply rgb by amount. If dye accumulates heavily, clamp per-channel 0-255 in rendering. But also we can maintain dye density >1? Maybe cap at 2 or 3 to avoid blow-up.

Mouse velocity injection:
- The simulation grid coordinates from mouse: `const mx = (e.clientX - rect.left) / rect.width * nx; const my = (e.clientY - rect.top) / rect.height * ny;` (0..nx, 0..ny). We can treat as float. Previous grid pos `prevMx, prevMy`.
- On move, compute `dx = mx - prevMx; dy = my - prevMy;` This is grid cells per raw frame. To convert to velocity per simulation time, divide by dt? Since dt is frame time. But if we inject each animation frame, we can inject `forceX = dx * force` and `forceY = dy * force` as cells per simulation step? If simulation dt is 0.016 and we inject along segment, the velocity should be ~dx/dt maybe. If mouse moved across half the grid between frames (fast), dx~160 cells, dt=0.016 => velocity ~10000, too huge. Mouse shouldn't move that fast; but if it does, we need clamp. It's okay to limit max injection velocity and spread along line.

A better approach: Add velocity proportional to `e.movementX / rect.width * nx` per frame. That's grid delta per frame. If we divide by dt, too large; but our simulation adds velocity each frame, which is used over dt. If we add a velocity equal to grid delta per frame, then in one frame, advect moves by u * dt (with dt ~ frame time). That would be negligible movement relative to the cursor (u*dt = dx*dt vs dx). To make cursor actually push fluid at the speed of the cursor, we should add velocity = dx / dt. But then if dt matches the frame, the injected particles move dx in one frame. So yes, `force = dx / dt * forceScale`. For 60fps, dt ~0.017, forceScale maybe 0.5, velocity ~30*dx. If dx=3 cells/frame -> u=45 cells/sec -> in one dt moves 0.75 cells. Hmm still less than the cursor. Actually speed = dx/dt cells per second; in one dt of simulation (which equals frame time), displacement = speed*dt = dx. So to match, force = dx/dt. We can scale down to e.g. 0.5 so fluid lags. Then displacement per frame = 0.5*dx, trail behind cursor. That's fine.

Injecting velocity per substep? We inject once per frame at current maybe. Then solver handles substeps.

Clamp injection velocity to max ~20*Nx? Actually 20 cells/second. Better clamp injection per grid cell to some range [-maxV, maxV] (e.g., 10-20). Since simulation dt max 0.05, displacement per substep max=1 cell. OK.

Dye injection:
- Along segment from prev to current grid point with step size maybe 0.5 cell: for t from 0 to 1 step 0.5/len; inject radially around point with Gaussian kernel radius 1.5 cells. Add to r/g/b arrays.
- Dye amount per point = speed * dyeScale, min.
- Use Gaussian: `factor = exp(-((i-x)^2+(j-y)^2)/(2*sigma^2))`. This is a few cells; manually loop over bounding box.

Potential performance: For each frame, injecting over small kernel (e.g., 5x5=25 cells per point) and maybe 20 points = 500 cells. Fine.

Velocity injection: Add along segment similarly, but simply add velocity at each point; maybe use same kernel. 

Let's implement `splat(x, y, vx, vy, r, g, b, amount)`. It loops over local area and applies Gaussian. For velocity, amount scaled. For dye, color scaled.

Also, to make swirling colors, inject different hue at different times; maybe also periodically auto-inject a small pulse at center if no mouse for a demo? Could be optional. The user wants responsive to dragging. Auto injection could be nice. But maybe keep simple: initial dye via instructions? Let's add a startup hint that disappears. Also maybe on first drag, hide instructions.

Now code structure:
- The single file should have all code inside.
- Use `DOMContentLoaded` or script at end of body.
- Canvas id `fluid`.
- On resize, reinitialize fluid? We can calculate grid size from window dimensions and `targetCells`. We also update `canvas.width`, `canvas.height`. On resize, we lose current simulation; acceptable. We might even set them fixed on load and not update to avoid resets. But responsive full-screen means we'd want resize. We can handle resize by creating new Fluid and resetting mouse positions. Simpler: set once at load and not update on resize, but then full-screen if window size changes maybe black borders. The simplest robust: init on load and on resize recreate. It resets when user resizes; acceptable for a demo.

However, on mobile orientation changes, resets. That's fine.

Now, write the solver functions. Let's design class:

```
class Fluid {
  constructor(nx, ny) {
    this.nx = nx;
    this.ny = ny;
    this.size = (nx + 2) * (ny + 2);
    ...
    this.u = new Float32Array(this.size);
    ...
    this.p = new Float32Array(this.size); // actually project uses p, div; use uPrev/vPrev...
  }
  IX(i,j) { return i + j * (this.nx+2); }
  set_bnd(b, x) { ... }
  lin_solve(b, x, x0, a, c, iter) { ... }
  diffuse(b, x, x0, diff, dt, iter) { ... }
  project(u, v, p, div, iter) { ... }
  advect(b, d, d0, u, v, dt) { ... }
  step(dt, visc, diff, iter) { ... }
  addSplat(...) { ... }
  render(imageData) { ... }
}
```

Wait after velocity advect we project with temp p,div arrays. We can use `uPrev`/`vPrev` as p/div. The function `project(u, v, p, div)` takes arrays u,v to make divergence-free, using p, div as temporaries.

Then step:
```
const sdt = dt; // may be sub
// Velocity step
this.diffuse(1, this.uPrev, this.u, viscosity, sdt, iter);
this.diffuse(2, this.vPrev, this.v, viscosity, sdt, iter);
this.project(this.uPrev, this.vPrev, this.u, this.v, iter); // uPrev,vPrev divergence-free; u,v contain p,div
this.advect(1, this.u, this.uPrev, this.uPrev, this.vPrev, sdt);
this.advect(2, this.v, this.vPrev, this.uPrev, this.vPrev, sdt);
this.project(this.u, this.v, this.uPrev, this.vPrev, iter); // u,v divergence-free; uPrev/vPrev tmp
// Dye step
this.diffuse(0, this.rPrev, this.r, diff, sdt, iter);
this.advect(0, this.r, this.rPrev, this.u, this.v, sdt);
this.diffuse(0, this.gPrev, this.g, diff, sdt, iter);
this.advect(0, this.g, this.gPrev, this.u, this.v, sdt);
this.diffuse(0, this.bPrev, this.b, diff, sdt, iter);
this.advect(0, this.b, this.bPrev, this.u, this.v, sdt);
```
But dye diffuse target array is `this.rPrev`, then advect reads `this.rPrev` and writes `this.r`. That means final dye is in `this.r`. Good. But if we want to decay, multiply `this.r` by decay after advect. Also maybe need to set_bnd for dye (b=0) inside diffuse/advect. Done.

Velocity arrays after final project: `this.u`, `this.v` are divergence-free. Good.

Now angle: `diffuse` parameter b is passed to lin_solve and set_bnd. For velocity x component b=1, y component b=2. For dye b=0. For project, we call set_bnd(0, div) and set_bnd(0,p), then after adjusting set_bnd(1,u) and set_bnd(2,v).

Boundary conditions in `set_bnd`:
For b==1, vertical edges negative. For b==2 horizontal edges negative. For corners average adjacent boundaries. This is standard.

`lin_solve`: loops interior, uses cRecip. To improve performance, we can unroll row pointers. We'll include.

`advect`: 
```
const nx = this.nx, ny = this.ny;
const dt0x = dt * nx; // if u is domain units? But we decided grid units. Use dt.
```
Let's use `dt` direct. But we may want trace step to account for grid resolution? In grid unit, yes dt direct. We'll perhaps multiply dt by a global scaling factor `timeScale` to control speed. Use `const dt0 = dt * 1.0;`. But if we have `dt` small (0.016), movement very little. To make visually fluid move at natural speed, we can set `dt` injected to solver as e.g. `0.5` per frame? Hmm. Let's think: Simulation is discrete grid units, but we don't express velocities physically. If mouse injects velocity `u` equal to grid delta per frame divided by dt (so ~dx*60). Suppose mouse moves 3 cells per frame -> u≈180. With dt=0.016 -> trace back 2.9 cells per step. Over substeps, fine. So direct dt works.

But `dt` from frame time ~0.016. If we use dt directly in diffusion a = dt * diff * nx*ny/50000 ~ 0.016*diff*1. That's okay. We set diff maybe 0.5 -> a=0.008. 20 iterations -> minimal.

Actually, if we set `dt` too small, simulation slow; if we scale dt up artificially (e.g. set timeScale=5), movement faster but may need more substeps. Simpler: in step, compute `const simDt = Math.max(0.001, dt) * timeScale;` and set `timeScale` maybe 4. Run substeps of max `maxSubDt` maybe 0.04. Then each substep dt' = simDt/subSteps. So simulation runs in a time-accelerated manner. With dt=0.016, simDt=0.064, subSteps=2 (since max 0.04), sdt=0.032. Good. Mouse velocity injection uses `dx / realDt` or `dx / simDt`? If we want cursor push independent of `timeScale`, use realDt. But simulation velocity is in grid cells per simulation second. The cursor moved dx over real dt; injection should be dx/simDt? Because the simulation uses simDt as its time unit. If timeScale=4, simDt is larger, so injected velocity smaller (per simulation time) so that it traverses dx in simDt. Actually displacement over one frame: velocity * simDt = dx (if no lag). So injected velocity = dx / simDt. With simDt larger than real dt, velocity is smaller, but since simulation advances faster, the net effect matches. This decouples visual speed from simulation dt. But not necessary; simpler: inject using `dx / realDt * strength`; simulation's `dt` used for advection is real dt (scaled by timeScale). Hmm.

Let's avoid confusion. Choose a fixed simulation approach:
- `dt` per frame = elapsed seconds (real).
- `subSteps` = ceil(dt / maxSimDt). `simDt = dt / subSteps`.
- Pass `simDt` to solver as the `dt`. No timeScale. Mouse injection: `vInj = dx / simDt * strength`? If we inject per substep? We inject once per frame before substeps, not per substep. Then injected velocity should be in cells per simulation second (i.e., over simDt). But u is in cells per simulation time; over one substep (simDt), displacement = u*simDt. To make fluid move by dx in one substep (matching cursor), u should be dx/simDt. Over `subSteps` substeps, the same velocity pushes fluid by dx per substep if kept unchanged? Actually velocity may be updated each substep; but initial impulse pushes. However, after injection, each substep moves dye by u*simDt; if there are subSteps, total displacement could be subSteps*dx. To avoid overpushing when multiple substeps, inject velocity as `dx / dt * strength` (per frame), and solver's dt = real dt maybe split; if velocity persists over substeps, displacement sum = u * dt = dx*subSteps*strength? Wait u = dx/dt * strength. Across subSteps of simDt=dt/subSteps, total displacement if velocity constant = u * dt = dx/dt*strength*dt = strength*dx. So independent of subSteps. Good. So injection based on real dt is correct.

But if we inject each frame and solver runs substeps, injected velocity in simDt terms is `vSim = dx/dt*strength`. Then each substep uses simDt; advect step uses dt=simDt. Over all substeps, displacement = vSim * dt_total = dx*strength. Good. So injection based on real dt.

Thus:
- `dt = now - lastTime`. 
- `subSteps = Math.max(1, Math.ceil(dt / MAX_STEP))`.
- `stepDt = dt / subSteps`.
- call `fluid.step(stepDt, ...)` for each substep. But injection happens only once before substeps. Good.

But wait `fluid.step` includes velocity diffusion etc. If we call it with `stepDt`, the diffusion coefficient uses `stepDt`. That is correct. But the injected velocity from mouse is in per real-second units; within a substep, velocities are per simulation time step? Actually the velocity unit is consistent if u has units cells per simulation time. With stepDt as dt, the advection moves u*stepDt cells. Sum over substeps gives u*dt. If u injected as dx/dt*strength, total displacement = strength*dx. Good.

Now what about the projection time step? Diffusion uses dt. Fine.

Now to improve visual swirl: We can add vorticity confinement if needed. But standard projection creates turbulence.

I think we can code CPU solver with 320x180-ish grid and get good visuals.

Let's write efficient code.

Data types: Use Float32Array for velocity/dye. Use Uint8ClampedArray imageData. Use `Math.min`, `Math.max`.

Detailed code for each function:

`set_bnd`:
```
const nx=this.nx, ny=this.ny, stride=nx+2;
for (let i=1; i<=nx; i++) {
  const bottom = i; // j=0
  const top = i + (ny+1)*stride; // j=ny+1
  if (b==2) {
    this.v? wait x[bottom] = -x[bottom+stride];
  } else {
    x[bottom] = x[bottom+stride];
  }
  if (b==2) {
    x[top] = -x[top-stride];
  } else {
    x[top] = x[top-stride];
  }
}
for (let j=1; j<=ny; j++) {
  const row = j*stride;
  x[row] = b==1 ? -x[row+1] : x[row+1];
  x[row+nx+1] = b==1 ? -x[row+nx] : x[row+nx];
}
// corners
const s = stride;
x[0] = 0.5*(x[1]+x[s]);
x[(ny+1)*s] = 0.5*(x[(ny+1)*s+1] + x[ny*s]);
x[nx+1] = 0.5*(x[nx]+x[nx+1+s]);
x[(ny+1)*s+nx+1] = 0.5*(x[(ny+1)*s+nx] + x[ny*s+nx+1]);
```
Need ensure `b` for v should reverse on horizontal edges (bottom/top), which is b=2. We pass b accordingly.

`lin_solve(b,x,x0,a,c,iter)`:
```
const nx=this.nx, ny=this.ny, stride=nx+2, cRecip=1/c;
for (let k=0;k<iter;k++) {
  for (let j=1;j<=ny;j++) {
    let cur=j*stride;
    let up=cur-stride;
    let down=cur+stride;
    for (let i=1;i<=nx;i++) {
      x[cur+i] = (x0[cur+i] + a*(x[cur+i-1]+x[cur+i+1]+x[up+i]+x[down+i]))*cRecip;
    }
  }
  this.set_bnd(b,x);
}
```
For diffuse, c = 1 + 4*a. Good.

`diffuse(b,x,x0,diff,dt,iter)`:
```
const a = dt * diff * (this.nx*this.ny) / 50000; // domain scaling approx
this.lin_solve(b,x,x0,a,1+4*a,iter);
```
But if a is 0, c=1. Fine.

`project(u,v,p,div,iter)`:
```
const nx=this.nx, ny=this.ny, stride=nx+2;
// compute divergence
for (let j=1;j<=ny;j++) {
  let cur=j*stride;
  for (let i=1;i<=nx;i++) {
    div[cur+i] = -0.5 * (u[cur+i+1] - u[cur+i-1] + v[cur+i+stride] - v[cur+i-stride]);
    p[cur+i] = 0;
  }
}
this.set_bnd(0, div);
this.set_bnd(0, p);
this.lin_solve(0, p, div, 1, 4, iter);
for (let j=1;j<=ny;j++) {
  let cur=j*stride;
  for (let i=1;i<=nx;i++) {
    u[cur+i] -= 0.5 * (p[cur+i+1] - p[cur+i-1]);
    v[cur+i] -= 0.5 * (p[cur+i+stride] - p[cur+i-stride]);
  }
}
this.set_bnd(1, u);
this.set_bnd(2, v);
```

`advect(b,d,d0,u,v,dt)`:
```
const nx=this.nx, ny=this.ny, stride=nx+2;
let dt0x = dt, dt0y = dt; // if grid units
for (let j=1;j<=ny;j++) {
  let cur=j*stride;
  for (let i=1;i<=nx;i++) {
    let x = i - dt0x*u[cur+i];
    let y = j - dt0y*v[cur+i];
    if (x < 0.5) x = 0.5;
    if (x > nx+0.5) x = nx+0.5;
    if (y < 0.5) y = 0.5;
    if (y > ny+0.5) y = ny+0.5;
    let i0 = Math.floor(x);
    let i1 = i0+1;
    let j0 = Math.floor(y);
    let j1 = j0+1;
    let s1 = x - i0;
    let s0 = 1 - s1;
    let t1 = y - j0;
    let t0 = 1 - t1;
    let row0 = j0*stride;
    let row1 = j1*stride;
    d[cur+i] = s0*(t0*d0[row0+i0] + t1*d0[row1+i0]) + s1*(t0*d0[row0+i1] + t1*d0[row1+i1]);
  }
}
this.set_bnd(b,d);
```
Note with `dt` small (0.02), trace back small. Good.

Potential issue: velocity boundary conditions cause velocities at edge to reflect, making fluid bounce. Good.

Performance: advect uses manual indexing. Good.

Now `render(imageData)`:
We need map floats to RGBA. Use background? Could make a subtle black background. Set alpha 255. The pixel data length = (nx+2)*(ny+2)*4? ImageData size should be nx x ny (interior). We create `imageData = ctx.createImageData(nx, ny)`. Render interior cells j=1..ny, i=1..nx.
```
const data = imageData.data;
let idx = 0;
const stride=nx+2;
for (let j=1; j<=ny; j++) {
  const row = j*stride;
  for (let i=1; i<=nx; i++) {
    const id = row+i;
    const r = this.r[id] * 255;
    const g = this.g[id] * 255;
    const b = this.b[id] * 255;
    data[idx] = r>255?255:(r<0?0:r);
    data[idx+1] = same g;
    data[idx+2] = same b;
    data[idx+3] = 255;
    idx += 4;
  }
}
ctx.putImageData(imageData,0,0);
```
Could multiply by exposure or gamma; maybe use `Math.min(255, r*255)`. Since decay, dye likely <1 unless many injects. Use `const exposure = 1.0`; maybe ramp with gamma for darker look. Could also use a tone curve like `value = 255 * (1 - exp(-density*2))` but linear clamp is simpler. Perhaps use `v = 255 * Math.min(1, density)` which saturates quickly. We can use a softcap: `v = 255 * density / (1+density)` to avoid hard clipping. This looks better (HDR-ish). Let's implement: `const f = density / (1 + density);` for each channel. Initial density 0 -> black; density 1 -> 128; density 3 -> 192; high never exceeds 255. This is nice. If we inject bright colors, they'll saturate softly.

Also maybe add subtle background color (deep blue/purple) for empty areas? Since values near 0 produce black. Softcap with density ~0.01 yields near 0. Good.

However, repeated injection may make density ~5-10; softcap gives near 255. After decay, fine.

Now decay: in step, after dye advect, multiply all dye cells by `decay = Math.pow(0.99, dt*60)`? For per simulation dt, if dt=0.016, decay per substep ~0.9984. We can do after each step: `this.scaleArray(this.r, 1 - dt * decayRate)`; where decayRate ~0.5 per second. For dt=0.016 -> factor 0.992. But if subSteps sum dt, apply in step once with dt_total? Simpler: after entire `fluid.step(dt)` (which itself loops substeps), multiply arrays by `Math.pow(0.99, dt*60)`? But if step called multiple times? We'll call once per frame with total dt; inside step, after dye advect, apply decay factor `Math.pow(0.995, dt*60)` maybe. We can include in step after dye advect: `const decay = Math.exp(-dt * dyeFade)`; multiply r,g,b. Similarly velocity damping `velocityDecay = Math.exp(-dt * velocityFade)` to kill momentum.

But the solver internally uses substeps; if we multiply by `dt` inside step once, that's per frame. Good.

Alternatively, each substep decay uses subdt. Let's apply after all substeps with dt to be consistent. In step, loop substeps. After substeps, multiply arrays.

Mouse injection:
In animation loop, if dragging, compute motion and injection. Or use pointer events to update injection per frame. We need to map mouse direction. To avoid injecting at exact grid cell, we can store `lastGridPos` and `gridPos`; in each frame if pointer down, we add a line splat between positions.

Splat function details:
```
addSplat(x, y, vx, vy, r, g, b, size=1.5, amount=1.0) {
  const invSigma2 = 1 / (2*size*size);
  const x0 = Math.max(1, Math.floor(x - size*3));
  const x1 = Math.min(this.nx, Math.ceil(x + size*3));
  const y0 = ...
  const stride=this.nx+2;
  for (let j=y0; j<=y1; j++) {
    const dy = j-y;
    const row = j*stride;
    for (let i=x0; i<=x1; i++) {
      const dx=i-x;
      const w = Math.exp(-(dx*dx+dy*dy)*invSigma2);
      const id=row+i;
      this.u[id] += vx * w;
      this.v[id] += vy * w;
      this.r[id] += r * amount * w;
      this.g[id] += g * amount * w;
      this.b[id] += b * amount * w;
    }
  }
}
```
But this splat adds velocity to current velocity arrays `u`/`v`. If we call before step, the solver will process them. Good.

Color injection line:
```
// mx,my current grid coords; lmx,lmy previous
const dx = mx-lmx, dy = my-lmy;
const dist = Math.hypot(dx,dy);
if (dist > 0) {
  const steps = Math.ceil(dist * 2); // 0.5 cell spacing
  for (let s=0; s<=steps; s++) {
    const t = s/steps;
    const sx = lmx + dx*t;
    const sy = lmy + dy*t;
    // velocity injection proportional to segment
    // We want to spread total force across points or per point? If n points, per-point force maybe stronger. Use total force per frame.
    const fx = (dx/dt) * forceScale; // cells per second maybe
    const fy = (dy/dt) * forceScale;
    const speed = Math.hypot(fx, fy);
    const amount = Math.min(2, speed * dyeScale); // dye amount
    const color = hslToRgb(hue/360, 1.0, 0.55);
    this.addSplat(sx, sy, fx, fy, color.r, color.g, color.b, radius, amount);
  }
}
```
But dividing by `dt` when dt could be 0? Use max dt. Use `const force = (dx / dt) * forceScale`. If dt=0.016 and dx=5 -> force=~5/0.016=312*0.2=62 huge. Wait forceScale small, e.g., 0.05 -> ~15. okay. Then advect in one dt moves 15*0.016=0.24 cells; over 1 frame, small compared to cursor dx=5. Hmm if we want fluid to track cursor better, we might want displacement ~0.2*dx maybe. That's okay: the dye is pushed less than cursor, producing trails.

But if user expects "dragging pushes swirling dye", maybe trail behind is fine.

Clamp injection velocity per cell to maybe 15. Dye amount per cell maybe 0.2.

Actually we can define injection in simulation velocity units differently: Use `forceScale = 0.5`, not divide by dt. So `vx = dx * forceScale` (cells per frame). Then advect moves by `vx*dt` per substep. Since dt small, slow. If no dt division, a drag of 5 cells/frame gives velocity 2.5 cells per frame, but sim time per frame includes maybe multiple substeps; over full dt (1 frame), displacement ~2.5 cells (since u updated each substep? Actually if velocity persists, displacement = u * dt_total = 2.5*1 = 2.5 cells). That still less than cursor dx=5. Trails. Fine. This avoids huge velocities. We can use `vx = dx * forceScale`, no division by dt. With `forceScale` maybe 1.0 to 2.0. But mouse delta can be large (fast movement across screen e.g. 100 cells in one frame) -> u=100-200 causing instability. Need clamp. We can clamp per injection to maxVel maybe 15. Dye injection proportional to clamped speed.

Maybe simpler and stable: 
```
let vx = dx * velocityScale;
let vy = dy * velocityScale;
const speed = Math.hypot(vx,vy);
if (speed > maxVel) { const s=maxVel/speed; vx*=s; vy*=s; }
```
velocityScale maybe 2.0. maxVel maybe 8.0. This means fluid move roughly up to 8 cells per frame. Good.

Dye amount = Math.min(amountMax, speed * dyeScale) maybe 1.0. Use color from hsl.

Then we don't divide by dt. This decouples from frame rate? dx is grid cells per frame; if frame rate lower, dx larger, but clamp prevents insane; sim dt splitting with substeps still moves at clamped speed. It might appear faster at lower fps because same dx compressed. Acceptable.

Add velocity to `this.u` and `this.v` directly (these are current arrays). Add dye to `this.r` etc. The solver's step uses current as initial for diffuse etc. Since injection is before step, good.

Radii: radius ~1.5 cells. Fast move steps ~ 2 cells between splats. Good.

Hue rotation: Use global `hue += 60*dt` maybe; cycle. Saturation 1.0; lightness 0.55.

UI/instructions: Add simple hint overlay: "Click and drag to push colored dye" at top center with white text, fade out after first drag. Since no external resources, use a `<div>` with CSS. Hidden once interacted.

Touch support: Use `pointerdown`, `pointermove`, `pointerup` on canvas. `pointer-events` good. Map via `e.clientX` for all.

Display/hint: Style with `position: absolute; top: 20px; width: 100%; text-align: center; color: rgba(255,255,255,.6); font-family: sans-serif; pointer-events:none; transition: opacity .5s`. 

Now initial demo: Maybe an initial splash of color in center? Could do slow swirling? But if user wants drag interaction, initial blank or with hint. Maybe add a pre-init dye pulse at center with swirl so something visible at start? But the request focuses on dragging. A start hint is enough.

Let's consider using WebGL for better effect. But the CPU version is simpler for a single file with no external shaders. Yet the user asked "real-time fluid simulation" and "swirling colored dye". The CPU solver with 320x180 should produce beautiful swirls. Good.

Potential issues with CPU performance:
- 57k cells per step, with lin_solve loops of 20 iterations for 4 diffuse + 40 for projection. Each iteration loops over all cells: operations ~ (4+4)*57k*20? Let's estimate: velocity diffuse 2 arrays * 20 = 2.3M, project ~2 lin_solve of 40 iterations = 4.6M, dye diffuse 3*20=1.7M, advect ~4 arrays 4*57k = 228k. Total ~9M per step. At 60fps, 540M ops/s. JS can handle maybe ~1-2B ops/s. Should be okay on desktop. Mobile might be slower but still okay at 30fps. We can reduce iterations to 10 for mobile. Maybe adapt iter based on resolution. Use `iter = Math.min(40, Math.max(8, Math.floor(5000000 / (nx*ny))))`? Or simpler `iter = 20`. Actually we can use 16. Let's keep iter = 16 maybe. That reduces ~7M per step. Good.

We can also implement Gauss-Seidel with successive over-relaxation to converge faster; but simple Jacobi enough.

Wait, `lin_solve` with c = 1+4a and a small, Jacobi converges slowly for diffusion but diffusion is small. For projection (Poisson) we need accurate pressure to avoid divergence; 16 Jacobi iterations okay if velocities not huge. Good.

One thing: `dt` inside `fluid.step` after our time scaling may be small per substep (max 0.02). Diffusion `a = dt*diff*nx*ny/50000` with diff maybe 0.001 and dt 0.02 -> a ~0.0002*area/50000 = 0.0002*1 = 0.0002. c=1.0008. Iterations 16 enough. Projection with a=1,c=4 works independent of dt.

Now simulation grid resolution. If we set targetCells=50000, on a 4K screen width 3840, height 2160, aspect 16/9, base sqrt(50000/(16/9))=167.7; nx=298, ny=168. Fine. On small mobile 375x667 aspect 0.56, base sqrt(50000/0.562)=298; nx=168, ny=298. Good. Canvas internal size ~300x300, css full screen; performance okay.

Could allow toggle of resolution? Not needed.

Let's code carefully. We'll create `getGridSize()`:
```
function getGridSize() {
  const target = 45000;
  const aspect = Math.max(window.innerWidth, 1) / Math.max(window.innerHeight, 1);
  const base = Math.sqrt(target / aspect);
  let nx = Math.round(base * aspect);
  let ny = Math.round(base);
  nx = Math.max(32, nx);
  ny = Math.max(32, ny);
  return {nx, ny};
}
```
Use target maybe 50000-60000. For high DPI maybe less. This will re-init on resize.

Now, let's code the complete file.

Implementation steps:
1. HTML skeleton with style and canvas, hint div.
2. Inline script after body (or at end). 
3. Helper functions: `hslToRgb(h,s,l)`, `clamp`.
4. Fluid class.
5. Setup variables: canvas, ctx, imageData, fluid, mouse state, hue.
6. Resize handler: update canvas internal size, recreate fluid. Since fluid reset, okay.
7. Pointer event handlers: down/up/move. Keep `gridPos` arrays. `pointerMove`: if dragging, `currentGridPos` computed; inject line; update `lastPos = currentPos`. 
8. Loop: `requestAnimationFrame`. Compute dt; max dt; substeps; for step `fluid.step(stepDt)`. Render; putImageData. Maybe update hue in loop using real dt.
9. Initialize: resize, requestAnimationFrame.

Potential subtlety: When not dragging, but previously injected velocities still present, fluid continues to evolve, causing swirls. Good.

Dye accumulation: repeated dragging in same spot saturates to white because r,g,b all high. Since color changes, mixing produces rainbow/white. Good.

Velocity damping: apply each step. If not damped, fluid will keep swirling. Gentle damping looks nice.

Boundary conditions: sealed box. Swirls reflect. Good.

Now, ensure performance of giant loops; let's avoid function calls in inner loops, rely on flat arrays. `Fluid` methods use `this.` lookups; okay. For inner loops, we can assign local vars. Good.

Potential bug: `lin_solve` for project uses `a=1, c=4` for Poisson; `cRecip = 0.25`. Good.

Need to think about `diffuse` after dye: We use `this.rPrev` as target. But after velocity step, `this.rPrev` etc. are used as previous for dye. At end of step, dye current arrays are `r,g,b`. Good. The constructor initializes all arrays zero. `this.rPrev` and `this.r` both zero. Step sets diffused and advected correctly.

Let's write `step(dt, viscosity, diffusion, iterations, fade)` function. It loops substeps? No, maybe substep loop in main. But if we want decay based on total dt, step does substeps internally. Let's have main compute `subSteps` and call `fluid.step(totalDt, subSteps)`. So step signature includes `subSteps`. Or we call `fluid.step(stepDt)` per substep. Then decay in each step uses stepDt. If we call per substep, velocity/dye decay per substep applies multiple times across frame with factor `pow(decayPerSubstep, subSteps)` effectively same as `pow(decay, totaldt?)` but not linear if per step factor constant; but small decay approx linear. Simpler: in main loop, just call `fluid.step(stepDt)` for each substep. Each step handles one substep including damping. That yields damping slightly dependent on substep count (higher at low fps). Acceptable.

To make damping time-rate independent, we can call step with total dt and internal substeps. I prefer step internal substeps for simplicity of render. Let's implement `step(dt, viscosity, diff, iter, velocityFade, dyeFade)` and inside compute `const subSteps = Math.max(1, Math.ceil(dt / MAX_SIM_DT)); const sdt = dt / subSteps;` then loop substeps. After each substep, the `u/v` current arrays are divergence-free. We should not decay per substep? We can apply decay after each substep? Approx okay. For consistency, apply after all substeps to total dt: `const vDecay = Math.exp(-dt * velocityFade); const dDecay = Math.exp(-dt * dyeFade);` multiply arrays after loops. This is time-consistent. Good.

Inside step, after substeps loops but before decay? The splat injection occurs before step. Dye is in current arrays. Step diffuses/advects current. Then decay. So final arrays used next frame and render after step. Good.

Now order of step:
```
step(dt) {
  const substeps = ...; const sdt = dt/substeps;
  for (let s=0; s<substeps; s++) {
    // velocity
    this.diffuse(1, this.uPrev, this.u, viscosity, sdt, iter);
    this.diffuse(2, this.vPrev, this.v, viscosity, sdt, iter);
    this.project(this.uPrev, this.vPrev, this.u, this.v, iter);
    this.advect(1, this.u, this.uPrev, this.uPrev, this.vPrev, sdt);
    this.advect(2, this.v, this.vPrev, this.uPrev, this.vPrev, sdt);
    this.project(this.u, this.v, this.uPrev, this.vPrev, iter);
    // dye
    this.diffuse(0, this.rPrev, this.r, diff, sdt, iter);
    this.advect(0, this.r, this.rPrev, this.u, this.v, sdt);
    ... g,b
  }
  // decay
  const vDecay = Math.exp(-dt*velocityFade);
  const dDecay = Math.exp(-dt*dyeFade);
  for each array multiply by decay
}
```
Note: After diffuse velocity into uPrev/vPrev, project makes uPrev/vPrev divergence-free, using u/v as p/div. Then advect uses uPrev/vPrev as velocity. Advect writes to u and v. Then project u/v divergence-free using uPrev/vPrev as p/div. Good.

But there is subtle bug: at beginning of step, `this.u` and `this.v` are current velocities. diffuse reads `this.u` and writes `this.uPrev`. After first diffuse, `this.v` still original current. Good. Then project uses `uPrev,vPrev`. Fine.

Dye: `this.r` is current, diffuse writes `this.rPrev`, advect reads `this.rPrev` and writes `this.r`. Good. But note `this.u`/`this.v` after projection are final velocities; we use them for dye advect. Good.

Decay after substeps: multiply `this.u`, `this.v`, `this.r,g,b`. But `this.uPrev/vPrev` used as temp, no need. However, in render we use `this.r,g,b` only. Good.

Now, we should ensure `this.rPrev/gPrev/bPrev` arrays not aliased; constructor separate. Good.

`set_bnd` modifies boundary values. For dye b=0, edges set equal to interior. That means dye at boundary visible; no flux. Good.

Potential issue: `lin_solve` for diffuse with a small, c=1+4a close to 1, `cRecip` ~0.999. Works.

Now, render: Maybe we want a black background. If dye density low, black. Good. Use softcap `v = 255 * d/(1+d)`. That means if density > 0.1, visible. A small splat of amount 1 at a cell yields ~128 for that cell, blurred. Good.

However, when dye has been diffused, density per cell may be <0.01 and nearly black. To keep visible trails, inject enough. Good.

Now, drawing via `putImageData` and CSS scaling. By default, canvas CSS dimension full screen; internal small. Browser smoothing will bilinearly scale image. That is intended. We can set `image-rendering: auto`. Good.

Now pointer coordinates and grid coords. `const rect = canvas.getBoundingClientRect(); const mx = (e.clientX - rect.left) / rect.width * fluid.nx;` But fluid may be undefined at first; after init fine. For resize, after recreate, old positions irrelevant. On pointerdown, set `lastGridX = mx`, `lastGridY = my`. On pointermove, compute current and line.

Need to handle touch? Pointer events cover touch/mouse.

Maybe hide cursor over canvas? On drag, maybe use `e.pointerId` and `canvas.setPointerCapture(e.pointerId)`.

Now, js code style. Use `const fluidSim = ...`? We'll name `fluid`.

Let's write hslToRgb:
```
function hslToRgb(h, s, l) {
  // h in [0,1]
  let r, g, b;
  if (s === 0) { r=g=b=l; } else {
    const hue2rgb = (p,q,t) => {
      if (t<0) t+=1;
      if (t>1) t-=1;
      if (t<1/6) return p+(q-p)*6*t;
      if (t<1/2) return q;
      if (t<2/3) return p+(q-p)*(2/3-t)*6;
      return p;
    };
    const q = l < 0.5 ? l*(1+s) : l+s-l*s;
    const p = 2*l-q;
    r = hue2rgb(p,q,h+1/3);
    g = hue2rgb(p,q,h);
    b = hue2rgb(p,q,h-1/3);
  }
  return {r,g,b};
}
```

Now params tuning:
- `viscosity = 0.00001` maybe very low.
- `diffusion = 0.00002` maybe low.
- `velocityFade = 0.05` (per second damping of velocities; means e-fold 20 sec). Maybe lower to keep swirl. Let's use `0.02`.
- `dyeFade = 0.15` (dye fades over ~6 sec). Good.
- `iter = 16`.
- `maxSimDt = 0.03`.
- `forceScale = 2.5` or `3.0`.
- `maxVel = 8.0`.
- `dyeScale = 0.15`.
- `splatRadius = 1.2`.
- `hueSpeed` = 40 deg/sec maybe.

But we need to see actual behavior. Since cannot run code in this environment; we must reason and tune.

Potential issue: `advect` with `dt` direct small. If velocity u is 5 and dt=0.03, trace back 0.15 cells, movement minimal. But the simulation updates velocities each frame; the fluid movement may feel slow. However, if there are 60fps and dt=0.016, trace back 0.08 cells per step; over a second, u=5 => move 80 cells. So active swirls traverse domain quickly. Wait, movement per substep is `dt` (seconds) times velocity (cells per second). If u=5, dt=0.016 -> 0.08 cells/frame -> 5 cells/sec. That's modest. The visual speed = velocity in cells/sec. If maxVel=8, dye moves up to 8 cells/sec. Across 300 cells width, takes 37 seconds. That might feel too slow. If we want more responsive, maxVel should be larger, e.g., 40 cells/sec. But then dt*N not: velocity=40, dt=0.016 -> trace 0.64 cells/frame, fine. Across width in 7.5 sec. Good.

But if user drags fast, mouse delta across screen maybe 200 cells in 0.2 sec => 1000 cells/sec. If we clamp velocity to 40, fluid push not as strong but okay. We can set maxVel = 30-50.

So injection velocity should be in cells/sec, not cells/frame. We need to convert dx (cells per frame) to cells/sec: `vx = dx / dt * forceScale`. If dx per frame ~5, dt=0.016, vx =312*forceScale. To get maxVel 40, forceScale~0.128. Fine. We can use `forceScale ≈ 0.15`. If user drags slowly (dx=1), vx=9. Good.

Then dye amount = speed * dyeScale; with speed up to 40, dyeScale maybe 0.05 gives amount up to 2. Good.

So use `dx/dt` but clamp to maxVel. That makes injection physically meaningful and avoid fps dependence (mostly).

Wait, if frame rate drops, dt large, `dx/dt` might be similar because dx per frame larger but dt larger. If frame rate is constant. If one frame skipped, dx accumulates; injection velocity still similar. Good. But we might clamp.

However, in `addSplat` we inject at segment positions over one frame. The total velocity injected in the field is per frame force. The solver uses `dt` per substep; over the total dt, the velocity will advect dye by v*dt = (dx/dt*forceScale)*dt = dx*forceScale. Good. So if forceScale=1, dye moves with cursor; with lower, lags. That's intuitive. So use `vx = dx/dt * forceScale`.

But if dt is tiny due to fast refresh (0.008), dx may be 2 cells; vx=2/0.008=250*forceScale. Clamp to maxVel. Good.

ForceScale maybe 1.0? That would make v = dx/dt, exactly cursor speed in cells/sec. With maxVel maybe 80, most drags not clamped. Let's use `forceScale = 1.0`, `maxVel = 60`. Dye moves at cursor speed, though fluid lags due to inertia/projection.

But if user drags super fast dx=100 in dt=0.016, v=6250, clamp 60. Good.

Dye amount = clamp(speed * dyeScale, 0, maxDye) maybe dyeScale=0.04. With speed=60 -> amount=2.4. Good. Use maxDye=3.

Now `splatSpeed`: line segment from last to current; number of splats = ceil(dist*2). At each splat add the same velocity and amount proportional to speed. If we divide along points, the total added velocity/dye per cell along line could be n times. If n is based on dist, total injected amount grows with speed*dist (area under line). That is appropriate for a brush stroke. It's okay.

But for fast movement, n large; each splat adds maxVel; overlap can amplify. However, splats spaced 0.5 cells apart with radius 1.2 overlap. The local cell might receive ~3-4 splats each at maxVel -> velocity > maxVel. To prevent runaway, after adding splats, the solver will diffuse/ project. Could still be high. We can instead add velocity only once per injection point, but line interpolation at spacing 0.5. Overlap is inherent. It's okay; projection helps.

Maybe add a per-step total velocity cap? Not necessary.

Let's decide to clamp the final velocity arrays? In step after project? We could clamp u/v to something like [-50,50] to prevent overflow. Could add optional clamping in step loop. But not essential.

Alternative advection stability: `advect` clamps x to [0.5, nx+0.5]. If u*dt large, trace x might be far outside, clamped to boundary, causing streaking but stable. With dt small and v clamped 60, dt=0.03 trace 1.8 cells, within nx+0.5. Good.

Now maximum dt in `step`: `maxSimDt = 0.03`. At 60fps dt=0.016 -> 1 substep. If velocities 60, trace 0.96 cells per substep. Good. At low fps dt=0.1 -> 4 substeps of 0.025; trace 1.5 cells per substep. Good.

`advect` uses bilinear interpolation with boundary extended; set_bnd after advect ensures no weird edge.

Need to consider that `dt` passed to step is total frame dt. Step internal uses substeps. Good.

Now initial hint: Let's have a `showHint=true` until pointerdown. Use CSS opacity.

Now code.

Testing mentally:
- Canvas internal resolution set by grid size. CSS full window. 
- First frame dt maybe 0, skip step. Good.
- Mouse down hint hide.
- On move, handle if pointer down. Use `if (!pointer.down) return;` etc.

Coordinates:
```
let rect = canvas.getBoundingClientRect();
const gx = (e.clientX - rect.left) / rect.width * fluid.nx;
const gy = (e.clientY - rect.top) / rect.height * fluid.ny;
```
Use `clientX`. For touch, pointer event clientX works.

Store `lastGx, lastGy`. On pointerdown, set both.

`pointermove`:
```
if (!down) return;
const gx = ..., gy = ...;
const dx = gx - lastGx, dy = gy - lastGy;
const dist = Math.hypot(dx,dy);
if (dist < 0.01) return;
const dt = Math.max(frameDt, 0.001);
// velocity injection
const fx = (dx/dt)*forceScale;
const fy = (dy/dt)*forceScale;
const speed = Math.hypot(fx,fy);
if (speed>maxVel) { const k=maxVel/speed; fx*=k; fy*=k; }
const amount = Math.min(maxDye, speed*dyeScale);
const col = hslToRgb((hue%360)/360,1,0.55);
const steps = Math.max(1, Math.ceil(dist*2));
for (let s=0; s<=steps; s++) {
  const t = s/steps;
  const sx = lastGx + dx*t;
  const sy = lastGy + dy*t;
  fluid.addSplat(sx, sy, fx, fy, col.r, col.g, col.b, splatRadius, amount);
}
lastGx=gx; lastGy=gy;
```
But `amount` per point not divided by steps; maybe if many steps, total dye injected per line length integrated. That's fine.

However, if dist is zero and pointerdown but not moving, we still want dye? Maybe not. If user holds down, no injection. Good.

`frameDt`: We can use global `lastFrameDt` from loop or time between current frame and previous pointer event? Since injection happens once per frame in pointer event, we can use a global `simDt` from last animation frame, but pointer event may fire multiple times per frame. To avoid over-injection, we could store injection per frame and apply in animation loop, not directly. But simpler: process each pointer event; multiple pointermove events per frame may inject more dye than expected but usually one per frame. The dt used can be `Math.max(currentFrameDt, 0.001)` stored globally, updated by loop. That's okay.

Maybe compute `fx = dx * velocityScale` without dt to avoid dt issues and use event timestamps: `const dt = (e.timeStamp - lastTimeStamp)/1000;`. That is frame-independent. Use movement delta per event. Good:
```
let lastTs = 0;
on pointermove:
  const t = e.timeStamp*0.001;
  const dt = t - lastTs;
  lastTs = t;
  if (dt<=0 || dt>0.2) return;
  const fx = (dx/dt)*forceScale; ...
```
Initialize lastTs on pointerdown. This is precise per event. Good. Then ignore frame dt for injection. Loop dt used for simulation. This is robust against multiple events per frame because each event uses event's own dt.

But if multiple moves in same frame, each has dt small; dx small per event; same total. Good.

For simulation, loop computes dt separately.

Use `forceScale` maybe 1.0. Dye amount speed*dyeScale.

Now, we still need to prevent injection when `dt` extremely small (0.001). Use `dt = Math.max(0.001, ...)` and clamp if dt>0.2 then reset.

If `e.movementX` supported, we could use it, but computing via positions is fine.

Potential issue: On pointerdown, `lastGx` set; first move dx is from down to move. Good.

Now about result: It should be visually decent. Could add a simple effect: background slightly tinted based on dye? Black fine.

Let's consider using slightly larger grid for better smoothness. target 50000. Maybe 60000 as `TARGET_CELLS`. Use target = 55000. Good.

Now code the class with careful memory. Let's write it out.

Potential array indexing errors: `stride = nx+2`. For `advect`, `d0[row1+i1]` where row1 = j1*stride. Good. Need ensure j1 <= ny+1; since clamp y to ny+0.5, floor y <= ny, j1 <= ny+1. Good. i1 <= nx+1.

`set_bnd` corners reference `x[(ny+1)*s]` and `x[ny*s]`, etc. Ensure indexes inside. Good.

Now performance of `lin_solve`: It calls `this.set_bnd` per iteration. Inside set_bnd loops edges. Fine.

Now main loop:
```
let lastTime = 0;
let hue = 0;
function frame(time) {
  const now = time * 0.001;
  let dt = now - lastTime;
  lastTime = now;
  if (dt > 0.1) dt = 0.1;
  if (dt <= 0) dt = 0.001;
  hue = (hue + dt * 50) % 360;
  fluid.step(dt, ...);
  fluid.render(ctx.getImageData? Actually have persistent imageData). 
  requestAnimationFrame(frame);
}
```
But note if window resizes and fluid recreated, imageData needs recreation. We'll maintain `imageData` variable in setup and update on resize.

`fluid.render(imageData)` writes into imageData. Then `ctx.putImageData(imageData, 0, 0)`. Good.

Now parameter values. Let's decide final values:
- viscosity = 0.00001
- diffusion = 0.000005 (very low dye diffusion)
- velocityFade = 0.05 per second? Maybe 0.02 for long-lived swirls.
- dyeFade = 0.2 (trail fades over ~5 sec)
- iter = 16
- maxSimDt = 0.03
- forceScale = 1.0
- maxVel = 60
- dyeScale = 0.05
- maxDye = 2.5
- splatRadius = 1.3

But diffusion coefficient `diffusion` used with domain scaling. Let's choose `diff = 0.00005`. With nx*ny/50000 ~1, dt=0.03 => a=1.5e-6. Negligible. Maybe increase diff to `0.0003`. Then a=9e-6; still tiny but over many frames some spread. Actually dye diffusion doesn't need to be obvious; advection does the swirling. Use `diff = 0.0001`.

viscosity similarly `0.00001` or 0.000005. Fine.

Now potential issue: `lin_solve` coefficient c = 1+4a with a tiny; iteration 16 maybe enough for projection independent of a. Good.

Now, should we include a real `u`/`v` boundary reflective but `advect` trace clamp leads to dye and velocity accumulating at edges? Dye boundary Neumann means no dye escapes at edges; dye near edge stays. Velocity boundary bounce: no fluid leaves. Fine.

Let's also consider a "vorticity confinement" to make swirling more interesting. Not necessary but nice. Could be added without too much complexity: compute curl at each interior cell, then add confinement force perpendicular to gradient of |curl|. However, that requires extra arrays and loops. Maybe skip.

Now, one concern is color mixing: If dye colors are pure rgb (e.g., red, green), mixing yields yellow etc. Good.

Now, maybe set CSS `canvas { touch-action: none; }` to prevent scrolling on mobile.

Let's also handle window resize: On resize, recompute grid and fluid; reset mouse positions? This will reset simulation. We may want to debounce to avoid resetting continuously. Use a timeout: when resize occurs, wait 200ms then reinit. Simpler: reinit immediately; but heavy. We can attach `resize` with `setTimeout` debounce. We'll set up initial at load and optional reinit after resize quiet. But for a simple demo, immediate reinit acceptable. However, continuous resize events may lag. Let's debounce:
```
let resizeTimeout;
window.addEventListener('resize', () => {
  clearTimeout(resizeTimeout);
  resizeTimeout = setTimeout(init, 200);
});
```
`init()` sets canvas size and creates fluid,imageData.

Now, during resize with no init, canvas CSS stretches old internal image; after init, new simulation. Fine.

Now, ensure we don't call `ctx.getImageData` each frame, because it's expensive? Actually createImageData once. `putImageData` each frame. Fine.

Let's also set canvas `willReadFrequently: false`? default. Create context: `const ctx = canvas.getContext('2d', { alpha: false });`.

Now final code. Let's craft the entire index.html.

Let's code in our head/notes then write.

We need to be careful with class methods using arrow functions or regular. We'll use regular methods and `this`.

Will the inline script be large? Yes.

Let's produce file via write tool.

Before writing, let's do a final mental simulation for potential bugs:

- `Fluid` arrays length: `(nx+2)*(ny+2)`. Example nx=300, ny=170 => stride=302, rows=172 => 51944. Good.
- `set_bnd` for velocity: bottom edge j=0 row index i. x[i] = b==2 ? -x[i+stride] : x[i+stride]; Wait row j=0 at `i`, row j=1 at `i+stride`. bottom edge (j=0) should mirror j=1. If b==2 (v), reverse vertical component at horizontal walls: v at bottom should be -v[inside]. For top edge row = i + (ny+1)*stride, inside row above is i + ny*stride. Good. For left/right edges: x[row] = b==1 ? -x[row+1] : x[row+1]; row corresponds j*stride + i=0, inside at row+1. For b==1 (u), reverse horizontal component at vertical walls. Good.
- corners: average of two adjacent boundary cells. Example `x[0]` bottom-left: adjacent `x[1]` (bottom edge i=1) and `x[stride]` (left edge j=1). Human formula: `x[0] = 0.5*(x[1] + x[stride])`. In code: `x[0] = 0.5*(x[1] + x[stride]);`. Good. `x[(ny+1)*stride]` bottom? Actually j=ny+1 at left edge: `x[(ny+1)*stride]`. Adjacent cells: `x[(ny+1)*stride + 1]` (top edge inside) and `x[ny*stride]` (left edge inside). Code uses `(ny+1)*s + 1` and `ny*s`. Good. `x[nx+1]` bottom-right: adjacent `x[nx]` and `x[nx+1+stride]`. Code: `x[nx]`, `x[nx+1+s]`. Good. `x[(ny+1)*s+nx+1]` top-right: adjacent `x[(ny+1)*s+nx]` and `x[ny*s+nx+1]`. Code correct.

- `lin_solve` indexing: `let cur = j*stride; let up = cur - stride; let down = cur + stride;`. j loops 1..ny, so up valid >= stride, down <= (ny+1)*stride. i loops 1..nx, cur+i-1 and cur+i+1 valid. Good.

- `project` divergence: uses `u[cur+i+1] - u[cur+i-1] + v[cur+i+stride] - v[cur+i-stride]`. But we need to divide by h^2? For h=1, divergence approx = 0.5*(du/dx + dv/dy) with sign? In Stam code: `div[IX(i,j)] = -0.5*h*(u[IX(i+1,j)]-u[IX(i-1,j)] + v[IX(i,j+1)]-v[IX(i,j-1)]);` and then p solve with `lin_solve(0,p,div,1,4)` (a=1,c=4). Then `u[IX(i,j)] -= 0.5*(p[IX(i+1,j)]-p[IX(i-1,j)])*h;`. If h=1, that matches our code. Good.

- `advect`: We clamp x between 0.5 and nx + 0.5. But if i ranges 1..nx, d0 at row0 and row1 within 0..ny+1. i0 in 0..nx, i1 in 1..nx+1. Good. For y similarly.

Now, since velocities at boundary are reflective, advect might sample from boundary cells. Fine.

Potential issue: `advect` reads `u` and `v` to trace back; these arrays boundary values matter. Good.

The render uses `this.r` etc. If dye density > 0, visible. We multiply by 255*(d/(1+d)). Bound to [0,255] since d>=0. Good.

After decay, multiply dye arrays by dDecay <1, ensuring values decrease. Good.

Now, `addSplat` directly modifies `this.u`, `this.v`, `this.r`, etc. It might write to boundary cells if x near edge. We restrict x0/x1 to 1..nx and y0/y1 to 1..ny. Good. It doesn't touch boundaries. Fine.

Now, the injection line uses `steps` = ceil(dist*2). If dist=0, we skip because dist<0.01. Good.

Now, `hslToRgb` returns values 0..1. Use as dye components. Good.

Now, main loop: `requestAnimationFrame(frame);` At each frame, calls `fluid.step(dt)`. Render. Good.

Now formatting: Write as a single clean file.

Possible improvements: add a small amount of random noise to dye injection or velocity vectors to create more turbulence. For example, add a slight tangential component. We can compute `tx = -dy/dist, ty = dx/dist` and add a rotational component: `fx += tx * swirlStrength; fy += ty * swirlStrength;`. This makes every drag inject vorticity. Use `swirlStrength` maybe speed*0.2. That will produce swirling. Let's do it:
```
const nx = -dy/dist;
const ny = dx/dist;
fx += nx * speed * 0.15;
fy += ny * speed * 0.15;
```
Then clamp to maxVel again. That adds rotational swirl to straight drags. Good.

Now finalize parameters:
```
const CONFIG = {
  targetCells: 55000,
  viscosity: 0.000005,
  diffusion: 0.0001,
  iterations: 16,
  maxSimDt: 0.03,
  velocityFade: 0.03,
  dyeFade: 0.18,
  forceScale: 1.0,
  swirlStrength: 0.2,
  maxVel: 60,
  dyeScale: 0.04,
  maxDye: 2.5,
  splatRadius: 1.3,
  hueSpeed: 45
};
```
Potential problem: diffusion scaling uses `dt * diff * (nx*ny)/50000`. For targetCells=55000, factor ~1.1. Fine.

We should clamp max velocity after adding swirl and before splat. In `addSplat`, velocity added to each cell with Gaussian. If multiple points overlap, local velocity may exceed maxVel. Acceptable.

Now, maybe we should also apply diffusion and projection and render with `imageData` of `canvas.width x canvas.height` (grid size). Since canvas CSS full screen, the putImageData writes pixel data to low-res canvas, then browser scales. Good.

One more: For high-DPI displays, CSS width/height are CSS pixels, while internal grid size maybe smaller, browser scales. Fine.

Let's write the file.

One potential bug: In `fluid.step`, after velocity advect and project, we use `this.u` and `this.v` final for dye advect. However, after the first project, `this.u`/`this.v` contain p/div (temporary). Then advect writes over them. OK. After advect, project writes divergence-free `this.u`,`this.v`; `this.uPrev`,`this.vPrev` now p/div. Good. Then dye advect uses `this.u`,`this.v`. Fine.

But wait: At the start of the next substep, `this.uPrev` and `this.vPrev` are p/div, not old velocities. `diffuse` reads `this.u` (current) and writes `this.uPrev`, overwriting p/div. Fine. So no need `uPrev` to hold old velocity. Good.

Now, accumulative dt inside step for decay: we use `dt` total. Good.

Now, what if `dt` is very small e.g. 0? We clamp to 0.001. `substeps = ceil(0.001/0.03)=1`, sdt=0.001. Good.

Now pointer injection uses event dt which can be very small on high-refresh monitors. Use `dt = Math.max(0.0001, eTime - lastTime)`. Good.

Potential issue: `e.timeStamp` may be relative to performance.now or epoch; difference in ms. We convert to seconds. Good.

Now, let's write.

Let's craft the file content as a single string and write to /workspace/index.html.

Wait, current directory is /workspace. Use relative path `index.html`.

Let's write.

Before final, consider whether to include any UI for resetting? Not necessary.

Let's also consider `meta viewport` for mobile.

Style:
```
html, body { margin:0; padding:0; width:100%; height:100%; overflow:hidden; background:#000; }
#fluid { display:block; width:100%; height:100%; touch-action:none; image-rendering:auto; }
#hint { ... }
```
Use `#hint` with pointer-events none.

Let's produce final file.

  I've created the index.html file. Now I should verify it by perhaps reading key parts or doing a quick syntax check. I can use node to check syntax? The file is HTML with JS; node can't parse html but can extract script? Alternatively just check that file exists. I can also open with headless? No browser. But maybe run `node -c` on script? Could extract. Let's just maybe use bash to check file presence and size, and maybe use Node to parse the JS portion.

I can write a small command to extract script content and run `node --check`. But not necessary. However, to avoid obvious syntax errors, we can run a check. Let's do it.

Use `sed -n '/<script>/,/<\/script>/p' index.html | sed 's/<script>//' | node --check -` maybe. But Node 20 may parse class, const etc. Good.

Alternatively, use `python` or `node` to parse. Let's just run the check.

  I attempted to run the script in vm, but it executes at top level and `document` is null. I provided `document:null`. To syntax-check, I should compile via `new vm.Script` without running, or use `vm.runInNewContext` but script will run. We can use `new Script(m[1])` which compiles but doesn't execute; syntax errors would be caught. Let's do that.

  Great, syntax is valid. Now perhaps also verify `index.html` exists. It does. The file is ready.

Should I make any final modifications? Let's think about possible improvements and bugs.

Potential bug: `addSplat` writes to `this.u`/`this.v` current arrays. After injection, `step` diffuses them. Good.

Potential bug: The `pointermove` event listener is on `window`. If user drags outside window, it still tracks. On `pointerup` also on window. Good.

Potential bug: `canvas.setPointerCapture` may throw if pointerdown not on canvas? It is on canvas. `onPointerUp` checks `canvas.hasPointerCapture`, but `hasPointerCapture` introduced? Modern enough. Add check `if (canvas.hasPointerCapture)`. Good.

Potential bug: On `pointercancel`, release pointer capture. Good.

Potential bug: `hslToRgb` expects h in [0,1]. We pass `(hue % 360)/360`. Good.

Potential bug: tone curve `data[idx] = Math.min(255, 255 * this.r[id] / (1 + this.r[id]))`. Since `this.r` Float32, fine. If negative? Decay and advection should keep non-negative if initial non-negative. Boundary copy may preserve. Good.

Potential bug: large velocities could cause `advect` trace x extremely negative (if dt large and velocity > maxVel due to overlap). We clamp dt per substep to 0.03 and maxVel 60, trace up to 1.8 cells; fine. But if overlapping splats produce local velocity >60 temporarily, trace could be larger for that cell. However pressure/project spreads. But after multiple overlapping, u maybe 200. `advect` trace 200*0.03=6 cells, still within nx+0.5 (300) clamped. Bilinear still stable but may smear. Acceptable.

Potential performance: `linSolve` loops with `this` lookups each iteration maybe okay. Could optimize by storing arrays locally but fine.

Potential issue: On high-DPI / small screen, CSS `canvas` width/height are CSS pixels; internal grid small. Browser smoothing looks nice. `touch-action: none` prevents scroll.

Potential issue: `image-rendering: auto` may cause bilinear smoothing which can look blurry but also smooth. Good.

Potential memory: All arrays Float32. Fine.

Potential issue: We don't handle `pointermove` when down false; okay.

Potential issue: `resize` reinit resets simulation; acceptable.

Potential issue: On some devices, `window.innerWidth` might be 0 at load? We call `init()` immediately; if size is 0, grid min 32. Usually okay. We could wait for DOMContentLoaded, but script is at end of body so `document.getElementById` works. `window.innerWidth` maybe available. Fine.

Potential issue: `requestAnimationFrame(frame)` initial `lastFrameTime=0`, first frame dt will be time since page load (~small) but might be large; we clamp dt>0.2? We set `if(dt>0.2) dt=0.2`. Good.

Potential bug: `hint` uses `#hint` style; on pointer down we set `opacity` to 0. Good.

One possible improvement: Add center swirl init? Not necessary.

Now, let's answer to user with concise note.