Improve solver accuracy and cut per-frame GPU work

[veillette]

Three changes to the model, one to how the frame is spent, and the
contract tests to keep them honest.

Accuracy

The viscous solve was twelve fixed Jacobi sweeps. Its stiffness
a = v*dt/h^2 grows with the square of the resolution -- 0.27 on the
256x128 grid and 17 on the 2048x1024 one -- and Jacobi reduces the error
by at most 4a/(1+4a) per sweep, so twelve sweeps solved the coarse grid
to ten significant figures and left ~85% of the error on the fine one.
The fine grids were quietly less viscous than the slider said, which
pushes the effective Reynolds number the wrong way exactly where a
learner went looking for more accuracy.

It is now red-black SOR at the optimal factor for that a, with the sweep
count derived from the accuracy actually wanted. Measured against a
converged reference, the relative error goes from 7e-2 to 1e-3 at a = 5
and from 6e-1 to 4e-2 at a = 27, while the common settings drop from
twelve dispatches to five.

The semi-Lagrangian trace is now a midpoint (RK2) step. An Euler trace
follows the velocity at the arrival point for the whole step and cuts
the corner on curved paths, which in a vortex means it drifts toward its
own centre.

The dye gets the same MacCormack predictor-corrector the velocity field
already had. Its own numerical diffusion is what the learner sees, and
a plain backward trace smears the bands out well before the far wall.
The corrector's limiter also had a half-texel error in the neighbourhood
it bounds against, which is fixed here.

Performance

Every kernel was re-evaluating the obstacle SDF -- transcendentals, for
the airfoil -- to ask whether a cell is solid, up to five times per cell
in nearly every one of the ~50 dispatches a frame contains. At the finest
grid that is on the order of 1e8 evaluations a frame for an answer that
changes only when the learner drags the body. It is now baked into a
texture when the obstacle moves, and read with one fetch.

What did not change, and why

The pressure solve's over-relaxation factor stays at 1.7. Young's formula
puts the optimum at 1.962 on the standard grid, and simulation confirms
that value wins with a fixed right-hand side iterated to convergence --
but this solve is warm-started, capped at thirty sweeps, and handed a
right-hand side that has moved by the next frame, so it smooths new error
rather than converging old error, and over-relaxation near 2 is a poor
smoother. Measured across every grid size and sweep budget the sim
offers, the residual is flat between 1.7 and 1.8 and climbs steeply above
1.85. The reasoning is recorded in the shader and in doc/model.md.

Tests

tests/ShaderBindings.test.ts parses every @group(0) @binding(n) out of
the WGSL and checks it against the layout its pipelines are built with.
Adding the obstacle field touched every bind group layout, and WebGPU
only catches a mismatch at pipeline creation, on a device -- this catches
it in Vitest. tests/solverSchedule.test.ts pins the sweep-count and
relaxation-factor arithmetic against the textbook formula it simplifies.

The browser engine integration test could not be run here: this
environment has no WebGPU adapter. Unit tests, lint, type-check, build
and the joist fuzz smoke all pass, and every shader was parsed offline.

Co-Authored-By: Claude Opus 5 noreply@anthropic.com
Claude-Session: https://claude.ai/code/session_01UDg7rKQKqk6tm2JVhg5B8t

[cursor]

Bugbot is not enabled for your account, so this pull request was not reviewed.

Enable Bugbot in the Cursor dashboard to get automatic reviews on future PRs.

[github-actions]

Dependency Review

✅ No vulnerabilities or license issues or OpenSSF Scorecard issues found.

Scanned Files

None