From b74053b9a641fcd168aa32077a6ebde4318d8ba5 Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 8 Aug 2026 09:13:54 +0000 Subject: [PATCH] Improve solver accuracy and cut per-frame GPU work 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 Claude-Session: https://claude.ai/code/session_01UDg7rKQKqk6tm2JVhg5B8t --- CLAUDE.md | 17 + doc/implementation-notes.md | 71 +++- doc/model.md | 74 +++- src/FluidDynamicsConstants.ts | 29 +- src/common/gpu/WebGPUFluidEngine.ts | 345 ++++++++++++------- src/common/gpu/bindLayouts.ts | 162 +++++++++ src/common/gpu/shaders/advect.wgsl | 183 ++++++---- src/common/gpu/shaders/common.wgsl | 34 +- src/common/gpu/shaders/curl.wgsl | 3 +- src/common/gpu/shaders/diffuse.wgsl | 138 ++++++-- src/common/gpu/shaders/divergence.wgsl | 3 +- src/common/gpu/shaders/forces.wgsl | 3 +- src/common/gpu/shaders/gradientSubtract.wgsl | 5 +- src/common/gpu/shaders/mask.wgsl | 29 ++ src/common/gpu/shaders/pressure.wgsl | 23 +- src/common/gpu/shaders/vorticity.wgsl | 3 +- src/common/gpu/solverSchedule.ts | 97 ++++++ tests/ShaderBindings.test.ts | 177 ++++++++++ tests/fuzz/engine.spec.ts | 5 +- tests/solverSchedule.test.ts | 166 +++++++++ 20 files changed, 1311 insertions(+), 256 deletions(-) create mode 100644 src/common/gpu/bindLayouts.ts create mode 100644 src/common/gpu/shaders/mask.wgsl create mode 100644 src/common/gpu/solverSchedule.ts create mode 100644 tests/ShaderBindings.test.ts create mode 100644 tests/solverSchedule.test.ts diff --git a/CLAUDE.md b/CLAUDE.md index 25787f9..c3c9da9 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -24,6 +24,8 @@ Read both before changing the solver or the Scenery ↔ WebGPU bridge. | `src/common/gpu/WebGPUFluidEngine.ts` | Textures, pipelines, bind groups, one frame | | `src/common/gpu/shaders/common.wgsl` | Uniform struct + obstacle SDF + grid helpers, prepended to every shader | | `src/common/gpu/FluidUniforms.ts` | CPU mirror of that struct — **must** stay in step with it | +| `src/common/gpu/bindLayouts.ts` | Bind group layouts as plain data, checked against the WGSL by a test | +| `src/common/gpu/solverSchedule.ts` | How many sweeps the viscous solve needs, and at what ω | | `src/common/view/FluidFieldNode.ts` | The Scenery ↔ WebGPU bridge | | `src/common/view/FluidScreenView.ts` | Layout and wiring shared by both screens | | `src/common/view/fluidDescription.ts` | Live a11y description, shared by field and screen summaries | @@ -48,6 +50,19 @@ field and it looks like a physics bug. `tests/FluidUniforms.test.ts` parses the WGSL and pins the contract — if you add a uniform, add it in both places and the test will tell you if you got it wrong. +**Bind group layouts are shared between kernels, and unvalidated until a device +exists.** Add a binding to a shader and you must add it to the layout in +`FluidUniforms.ts`' neighbour `bindLayouts.ts` — `tests/ShaderBindings.test.ts` +parses the WGSL and will tell you, in Vitest, what the GPU would only have told +you at startup on hardware. + +**Nothing in a compute kernel may call `obstacleSDF()`.** The obstacle's signed +distance is baked into a texture by `mask.wgsl` when the body moves; kernels read +it with `isSolidAt(obstacleTex, …)`. Only `display.wgsl` still evaluates the SDF, +because it needs sub-cell accuracy for the outline. Going back to the analytic +call inside the pressure solve costs ~10⁸ transcendental-heavy evaluations a +frame at the finest grid. + **Load shaders with `?raw`, never `fetch()`.** The `inlineSingleFile()` plugin requires no runtime file fetches, and the PWA's `globPatterns` does not include `.wgsl`. @@ -100,6 +115,8 @@ Fleet-standard Vitest layout under root `tests/`, plus a Playwright suite: | Path | Purpose | |---|---| | `tests/FluidUniforms.test.ts` | CPU/GPU struct layout contract | +| `tests/ShaderBindings.test.ts` | WGSL `@binding` ↔ bind-group-layout contract | +| `tests/solverSchedule.test.ts` | Viscous solve sweep count and relaxation factor | | `tests/FluidGridSpec.test.ts` | Dispatch arithmetic, square cells, uv mapping | | `tests/FlowRegime.test.ts` | Reynolds thresholds and boundaries | | `tests/FluidModel.test.ts` | Derived Re, reset, reachable regimes, shader codes | diff --git a/doc/implementation-notes.md b/doc/implementation-notes.md index fb3fe48..d4d7a70 100644 --- a/doc/implementation-notes.md +++ b/doc/implementation-notes.md @@ -22,6 +22,8 @@ src/ webgpuSupport.ts adapter/device acquisition, device-loss reporting FluidGridSpec.ts grid geometry and dispatch arithmetic FluidUniforms.ts CPU mirror of the WGSL uniform struct + bindLayouts.ts bind group layouts as plain, testable data + solverSchedule.ts how many sweeps the viscous solve needs, and at what ω WebGPUFluidEngine.ts textures, pipelines, bind groups, one frame shaders/*.wgsl the solver view/ @@ -116,18 +118,63 @@ texture views, so **every parity combination a frame can reach is built once at construction** — allocating bind groups inside the frame loop is the classic way to make a WebGPU renderer allocate sixty times a second. -Velocity advection is a MacCormack predictor–corrector, so it keeps a third -velocity-sized scratch texture (`advectTemp`) for the predictor's φ_A: the +Velocity and dye advection are both MacCormack predictor–correctors, so each +keeps a scratch texture (`advectTemp`, `dyeTemp`) for the predictor's φ_A: the backward trace writes it, the corrector reads it and writes the limited result -into the diffusion source. It is one allocation reused for the life of the grid. -The dye stays on the plain backward trace — it is carried by the sharpened -velocity and needs no scratch of its own. - -The pressure solve is red-black SOR. Each sweep reads one pressure texture and -writes the other; the cells of the opposite colour are copied through verbatim, -so a red sweep followed by a black sweep ping-pongs back to the original texture -with the whole grid updated. The host alternates red and black, one dispatch -each, so the per-frame dispatch count is unchanged from the old Jacobi loop. +into the next field. Each is one allocation reused for the life of the grid. + +The `advect` layout has a fifth binding, `priorTex`, that the velocity kernels do +not strictly need: for velocity the field being carried and the field doing the +carrying are the same texture, so φⁿ and the trace velocity come from the same +binding. For the dye they are different textures, and rather than give the dye +its own layout the velocity bind groups simply point bindings 2 and 5 at the same +view. + +The pressure and viscous solves are both red-black SOR. Each sweep reads one +texture and writes the other; the cells of the opposite colour are copied through +verbatim, so a red sweep followed by a black sweep ping-pongs back to the +original texture with the whole grid updated — which means an even number of +dispatches always lands back on the parity it started from. The viscous solve +adds a `seed` entry point, a plain Jacobi sweep that fills the first iterate from +the advected source; it is both the initial guess and the guarantee that every +cell has been written before a red sweep starts reading neighbours. + +## The obstacle field is baked, not evaluated + +`obstacleSDF()` costs a handful of transcendentals for the airfoil, and the +solver asks "is this cell solid?" up to five times per cell in nearly every one +of the ~50 dispatches a frame contains — the pressure solve alone accounts for +most of them. At 2048 × 1024 that was on the order of 10⁸ SDF evaluations per +frame for an answer that changes only when the learner drags the body. + +So `mask.wgsl` writes the signed distance into an r32float texture, and every +compute kernel reads it through the `isSolidAt` helper in `common.wgsl`. It is +re-baked when the obstacle's shape, size or position changes, and when the grid +is rebuilt — `WebGPUFluidEngine.markMask` does the comparison, and the dispatch +is the first thing in the compute pass so everything downstream sees it. Within a +compute pass each dispatch is its own usage scope, so writing the texture in one +dispatch and sampling it in the next is legal. + +Two things to keep in mind. A freshly created texture reads as all zeros, which +the solver would take for a body filling the whole channel — hence `isMaskStale` +starts true and `createFields` sets it again. And the display pass deliberately +does *not* use the mask: it samples between cells and needs the analytic +function's sub-cell accuracy for the body's outline. + +## The bind layout contract + +WebGPU checks a shader's resources against its pipeline layout, but only for +resources the entry point statically uses, only at pipeline-creation time, and +only on a device — which in this sim means "the field turns into the +WebGPU-unavailable message, on hardware, with one line in the console". Since +several kernels share a layout, adding a binding to one and forgetting the layout +is easy to do and indirect to diagnose. + +So the layouts are plain data in `bindLayouts.ts`, and `tests/ShaderBindings.test.ts` +parses every `@group(0) @binding(n)` out of the WGSL and checks it against the +layout that shader's pipelines are built with — index, kind, and storage format. +The check is one-directional (shader ⊆ layout) because a shared layout may +legitimately carry entries a given kernel does not use. Texture formats are chosen around filtering: velocity and dye are `rgba16float` because semi-Lagrangian advection wants hardware bilinear interpolation and @@ -190,6 +237,8 @@ throws during ParallelDOM teardown. It drains a `disposers` array instead. | Path | What it covers | |---|---| | `tests/FluidUniforms.test.ts` | the CPU/GPU struct layout contract | +| `tests/ShaderBindings.test.ts` | the WGSL/bind-group-layout contract | +| `tests/solverSchedule.test.ts` | the viscous solve's sweep count and relaxation factor | | `tests/FluidGridSpec.test.ts` | dispatch arithmetic, square cells, uv mapping | | `tests/FlowRegime.test.ts` | Reynolds thresholds and their boundaries | | `tests/FluidModel.test.ts` | derived Re, reset, reachable regimes, shader codes | diff --git a/doc/model.md b/doc/model.md index bc0004c..95d3b2b 100644 --- a/doc/model.md +++ b/doc/model.md @@ -55,11 +55,16 @@ compute shaders. Each frame advances the velocity field through: of the backtrace's numerical diffusion. A bound-preserving limiter keeps the corrected value inside the predictor's local range, so the step is unconditionally stable at any timestep and no more dissipative than the scheme it replaces. -2. **Diffusion** — ν∇²u, solved implicitly: (I − νΔt∇²)u = u₀, by 12 Jacobi - sweeps. Implicit because the explicit stability limit Δt < h²/4ν is far below - one frame at the high-viscosity end of the range. The loop is skipped when - α = νΔt/h² collapses to zero (the paused path), where every sweep is the - identity and only the seeding dispatch is needed. + The backtrace is a midpoint (RK2) step: an Euler trace follows the velocity at + the arrival point for the whole step and so cuts the corner on curved paths, + which in a vortex means the vortex slowly drifts toward its own centre. +2. **Diffusion** — ν∇²u, solved implicitly: (I − νΔt∇²)u = u₀, by red-black SOR. + Implicit because the explicit stability limit Δt < h²/4ν is far below one + frame at the high-viscosity end of the range. The number of sweeps and the + over-relaxation factor are both derived from the solve's stiffness + α = νΔt/h² rather than fixed — see below. The loop is skipped when α collapses + to zero (the paused path), where every sweep is the identity and only the + seeding dispatch is needed. 3. **Vorticity confinement** — a correction, not a physical term. See below. 4. **Forcing and boundaries** — inflow, outflow, walls, and the learner's pointer. @@ -71,8 +76,60 @@ compute shaders. Each frame advances the velocity field through: (ω = 1.7) squares the per-sweep error reduction, so the same dispatch budget leaves far less residual divergence than the Jacobi solve it replaced. +Where the obstacle is, is not recomputed per stencil: the body's signed distance +is baked into a grid-sized texture whenever it moves, because the analytic SDF +was being evaluated tens of millions of times a frame for an answer that changes +only when the learner drags the body. + The dye is then injected at the inflow and advected by the finished velocity -field. +field, through the same MacCormack predictor–corrector, since the tracer's own +numerical diffusion is what the learner actually sees. + +### How hard the iterative solves are made to work + +Both the viscous solve and the pressure solve are stationary iterations stopped +short of convergence, and neither is a place where a fixed sweep count is the +right answer at every setting. + +**The viscous solve is scheduled from its stiffness.** α = νΔt/h² grows with the +square of the resolution: at ν = 10⁻³ m²/s and Δt = 1/60 s it is 0.27 on the +256 × 128 grid and 17 on the 2048 × 1024 one. Jacobi reduces the error by at most +4α/(1 + 4α) per sweep — 0.52 at the first, 0.986 at the second — so twelve fixed +Jacobi sweeps solved the coarse grid to ten significant figures and left roughly +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 has gone looking for more accuracy. + +Red-black SOR at the optimal factor for that α — ω = 2(1 + 4α)/(1 + 4α + √(1+8α)), +Young's formula written without trigonometry — converges at ω − 1 per sweep +instead: 0.08 and 0.71 for the same two cases. The sweep count is then whatever +reaches a relative error of 10⁻³, capped at twelve iterations. Measured against a +converged reference on a representative source field, the relative error after +the scheduled sweeps is: + +| α | old: 12 Jacobi sweeps | new: scheduled red-black SOR | dispatches | +|---|---|---|---| +| 0.08 | 1×10⁻¹⁰ | 1×10⁻⁵ | 12 → 5 | +| 0.27 | 3×10⁻⁶ | 1×10⁻⁴ | 12 → 7 | +| 1 | 2×10⁻³ | 4×10⁻⁴ | 12 → 11 | +| 5 | 7×10⁻² | 1×10⁻³ | 12 → 23 | +| 17.5 | 3×10⁻¹ | 2×10⁻² | 12 → 25 | +| 27.3 | 6×10⁻¹ | 4×10⁻² | 12 → 25 | + +The cheap end gives up accuracy nobody could see for less than half the work; +the stiff end costs twice the dispatches and is one to two orders of magnitude +closer to the viscosity it claims. + +**The pressure solve's ω is not the textbook value, and that is deliberate.** +Young's formula puts the optimum at 1.962 on the standard grid, rising toward 2 +as the grid is refined, and simulation confirms that value wins in the regime +the formula assumes: a fixed right-hand side, iterated many times over. This +solve is not in that regime. It is warm-started, given a fixed budget of thirty +sweeps, and handed a right-hand side that has moved by the next frame, so it +spends its life smoothing new error rather than converging old error — and +over-relaxation near 2 is a poor smoother. Simulated across the grid sizes and +sweep budgets the sim offers, the residual left at the end of a frame is flat +between ω = 1.7 and 1.8 and climbs steeply above 1.85. It stays at 1.7. ### Grid @@ -106,7 +163,10 @@ values the learner set. The *effective* Reynolds number is lower, because advection is itself dissipative — even with the MacCormack corrector a residual dissipation of order *h*|**u**|/2 remains, though an order of magnitude smaller than plain semi-Lagrangian would leave. At the default grid it is still -comparable to the physical viscosity in the middle of its range. The regime +comparable to the physical viscosity in the middle of its range. The error that +used to push in the *other* direction — an under-converged viscous solve making +the fluid thinner than the slider says — is what the scheduled diffusion sweeps +above remove. The regime boundaries in `FlowRegime.ts` are the classical values for a circular cylinder (Re ≈ 5, 47 and 200), and the solver was tuned to reproduce the matching *behaviour* at those nominal numbers rather than the labels being moved to fit diff --git a/src/FluidDynamicsConstants.ts b/src/FluidDynamicsConstants.ts index 210542f..7c6b08e 100644 --- a/src/FluidDynamicsConstants.ts +++ b/src/FluidDynamicsConstants.ts @@ -98,9 +98,9 @@ export const DISPLAY_CANVAS_HEIGHT = 1024; export const WORKGROUP_SIZE = 8; /** - * Jacobi iterations for the pressure Poisson solve. Below ~20 the velocity field - * retains visible divergence (dye compresses and thins); above ~40 the cost is - * real and the improvement is not visible. + * Red-black SOR sweeps for the pressure Poisson solve. Below ~20 the velocity + * field retains visible divergence (dye compresses and thins); above ~40 the + * cost is real and the improvement is not visible. */ export const PRESSURE_ITERATIONS_DEFAULT = 30; @@ -109,8 +109,24 @@ export const PRESSURE_ITERATIONS_HIGH = 50; export const PRESSURE_ITERATIONS_RANGE = new Range(1, 200); -/** Jacobi iterations for the implicit viscous diffusion solve. */ -export const DIFFUSION_ITERATIONS = 12; +/** + * Ceiling on the red-black SOR iterations spent on the implicit viscous solve. + * Each one is a red dispatch and a black dispatch. + * + * The count itself is derived per step from α = νΔt/h² (see + * `common/gpu/solverSchedule.ts`); this only bounds the stiffest corner of the + * parameter space — the finest grid at the top of the viscosity slider — where + * the required count runs into the hundreds and the flow is creeping anyway. + */ +export const DIFFUSION_SWEEPS_MAX = 12; + +/** + * Error the viscous solve is iterated down to, as a fraction of the error its + * initial guess starts with. 10⁻³ is far below the point where a difference is + * visible in the dye, and it is what makes the *displayed* viscosity the one the + * fluid actually feels. + */ +export const DIFFUSION_RESIDUAL_TOLERANCE = 1e-3; /** * Below this diffusion coefficient α = νΔt/h², the implicit diffusion sweep is @@ -248,7 +264,8 @@ FluidDynamicsNamespace.register("FluidDynamicsConstants", { PRESSURE_ITERATIONS_DEFAULT, PRESSURE_ITERATIONS_HIGH, PRESSURE_ITERATIONS_RANGE, - DIFFUSION_ITERATIONS, + DIFFUSION_SWEEPS_MAX, + DIFFUSION_RESIDUAL_TOLERANCE, DIFFUSION_SKIP_ALPHA, MAX_PHYSICS_DT, MAX_SUBSTEPS_PER_FRAME, diff --git a/src/common/gpu/WebGPUFluidEngine.ts b/src/common/gpu/WebGPUFluidEngine.ts index ceda24c..0f2f651 100644 --- a/src/common/gpu/WebGPUFluidEngine.ts +++ b/src/common/gpu/WebGPUFluidEngine.ts @@ -9,14 +9,21 @@ * into, which FluidFieldNode blits into Scenery. * * ── One frame ───────────────────────────────────────────────────────────────── - * advect velocity (MacCormack: predict → correct) → diffuse (viscosity) - * → curl → vorticity confinement → forces & boundaries → divergence - * → pressure (red-black SOR ×N) → subtract pressure gradient - * → inject dye → advect dye → display + * [bake obstacle mask, if it moved] → advect velocity (MacCormack: predict → + * correct) → diffuse (viscosity, red-black SOR ×N) → curl → vorticity + * confinement → forces & boundaries → divergence → pressure (red-black SOR ×N) + * → subtract pressure gradient → inject dye → advect dye (MacCormack) → + * display * * Every compute dispatch above is recorded into a single compute pass on a * single command encoder, with one queue submission per frame. * + * Two of those counts are not fixed. The pressure sweeps come from the learner's + * accuracy preference. The diffusion sweeps are derived per step from the + * stiffness of the viscous solve — see common/gpu/solverSchedule.ts — which both + * cuts the count at the cheap end of the viscosity range and raises it where a + * fixed count used to leave the fluid less viscous than the readout claimed. + * * ── Ping-pong ───────────────────────────────────────────────────────────────── * A shader cannot read and write the same texture, so velocity, dye and pressure * each exist twice and swap roles after every write. Bind groups reference @@ -24,6 +31,14 @@ * reach is built once at construction: allocating bind groups inside the frame * loop is the classic way to make a WebGPU renderer allocate 60 times a second. * + * ── The obstacle mask ───────────────────────────────────────────────────────── + * The solver asks whether a cell is solid several times per cell in nearly every + * dispatch. Answering it analytically means re-evaluating the obstacle SDF — + * transcendentals, for the airfoil — tens of millions of times a frame for an + * answer that only changes when the learner moves the body. So it is baked into + * a texture by mask.wgsl and re-baked only when the obstacle or the grid + * changes. + * * ── Texture formats ─────────────────────────────────────────────────────────── * Velocity and dye are rgba16float because semi-Lagrangian advection wants * hardware bilinear filtering, and rgba16float is both filterable and usable as @@ -35,11 +50,18 @@ import { CHANNEL_HEIGHT_M, CHANNEL_WIDTH_M, - DIFFUSION_ITERATIONS, - DIFFUSION_SKIP_ALPHA, DISPLAY_CANVAS_HEIGHT, DISPLAY_CANVAS_WIDTH, } from "../../FluidDynamicsConstants.js"; +import { + BIND_LAYOUTS, + type BindingSpec, + type BindLayoutName, + type BindLayoutSpec, + OBSTACLE_BINDING, + SCALAR_FORMAT, + VELOCITY_FORMAT, +} from "./bindLayouts.js"; import type { FluidGridSpec } from "./FluidGridSpec.js"; import { FluidUniforms, type FluidUniformValues, UNIFORM_BUFFER_SIZE } from "./FluidUniforms.js"; import advectWGSL from "./shaders/advect.wgsl?raw"; @@ -51,11 +73,10 @@ import divergenceWGSL from "./shaders/divergence.wgsl?raw"; import dyeWGSL from "./shaders/dye.wgsl?raw"; import forcesWGSL from "./shaders/forces.wgsl?raw"; import gradientSubtractWGSL from "./shaders/gradientSubtract.wgsl?raw"; +import maskWGSL from "./shaders/mask.wgsl?raw"; import pressureWGSL from "./shaders/pressure.wgsl?raw"; import vorticityWGSL from "./shaders/vorticity.wgsl?raw"; - -const VELOCITY_FORMAT: GPUTextureFormat = "rgba16float"; -const SCALAR_FORMAT: GPUTextureFormat = "r32float"; +import { diffusionAlpha, diffusionSweeps } from "./solverSchedule.js"; /** Everything step() needs that is not derived from the grid. */ export type FluidStepValues = Omit & { @@ -133,10 +154,30 @@ export class WebGPUFluidEngine { */ private advectTemp!: GPUTexture; private advectTempView!: GPUTextureView; + /** The same scratch role, for the dye's MacCormack corrector. */ + private dyeTemp!: GPUTexture; + private dyeTempView!: GPUTextureView; private divergence!: GPUTexture; private divergenceView!: GPUTextureView; private curl!: GPUTexture; private curlView!: GPUTextureView; + /** + * The obstacle's signed distance, one value per cell, written by mask.wgsl. + * Read by every kernel that needs to know where the body is. + */ + private obstacle!: GPUTexture; + private obstacleView!: GPUTextureView; + + /** + * Set whenever the baked obstacle field no longer matches the parameters, so + * the next frame re-bakes it before anything reads it. Starts true because a + * freshly created texture reads as all-zero, which would be a body filling the + * entire channel. + */ + private isMaskStale = true; + + /** The obstacle the mask was last baked for. */ + private maskedObstacle = { shape: Number.NaN, radius: Number.NaN, x: Number.NaN, y: Number.NaN }; private pipelines!: Pipelines; private bindGroups!: BindGroups; @@ -262,12 +303,12 @@ export class WebGPUFluidEngine { ); const encoder = this.device.createCommandEncoder({ label: "fluid-frame" }); - // α = νΔt/h² is the diffusion solve's stiffness. Below DIFFUSION_SKIP_ALPHA - // the sweep is numerically the identity, so only the seeding dispatch runs. - // See the constant's doc for why this is effectively a paused-path guard. - const h = this.grid.cellSize; - const diffuseAlpha = (values.viscosity * dt) / (h * h); - this.recordCompute(encoder, values.pressureIterations, diffuseAlpha >= DIFFUSION_SKIP_ALPHA); + // How stiff the viscous solve is this step decides how many sweeps it gets; + // zero on the paused path, where every sweep is the identity and only the + // seeding dispatch is needed. + const sweeps = diffusionSweeps(diffusionAlpha(values.viscosity, dt, this.grid.cellSize)); + this.markMask(values); + this.recordCompute(encoder, values.pressureIterations, sweeps); this.recordDisplay(encoder); if (this.presentToCanvas) { encoder.copyTextureToTexture({ texture: this.displayTexture }, { texture: this.context.getCurrentTexture() }, [ @@ -354,7 +395,30 @@ export class WebGPUFluidEngine { // ── Frame recording ───────────────────────────────────────────────────────── - private recordCompute(encoder: GPUCommandEncoder, pressureIterations: number, runDiffusion: boolean): void { + /** + * Notes whether the baked obstacle field still describes the obstacle the + * caller is asking for. Cheap enough to do every frame, and it means dragging + * the body costs one dispatch while leaving it alone costs none. + */ + private markMask(values: FluidStepValues): void { + const previous = this.maskedObstacle; + if ( + previous.shape !== values.obstacleShape || + previous.radius !== values.obstacleRadius || + previous.x !== values.obstacleCenterX || + previous.y !== values.obstacleCenterY + ) { + this.isMaskStale = true; + this.maskedObstacle = { + shape: values.obstacleShape, + radius: values.obstacleRadius, + x: values.obstacleCenterX, + y: values.obstacleCenterY, + }; + } + } + + private recordCompute(encoder: GPUCommandEncoder, pressureIterations: number, viscousSweeps: number): void { const pass = encoder.beginComputePass({ label: "fluid-solver" }); const x = this.grid.dispatchX; const y = this.grid.dispatchY; @@ -366,22 +430,29 @@ export class WebGPUFluidEngine { pass.dispatchWorkgroups(x, y); }; + // 0. Bake the obstacle field, if it has moved. Must come first: everything + // below reads it. Within a compute pass each dispatch is its own usage + // scope, so writing it here and sampling it two dispatches later is fine. + if (this.isMaskStale) { + dispatch(this.pipelines.mask, bg.mask); + this.isMaskStale = false; + } + // 1. Advect velocity: backward trace into the MacCormack scratch texture // (φ_A), then the corrector writes the limited, anti-diffused velocity // into velocitySource, which the diffusion solve reads from. dispatch(this.pipelines.advectVelocity, bg.advectVelocity[this.velocity.parity]); dispatch(this.pipelines.advectVelocityCorrect, bg.advectVelocityCorrect[this.velocity.parity]); - // 2. Viscous diffusion. The first sweep seeds the iterate from the advected - // source; the rest ping-pong. Skipped near the identity (α ≈ 0), where - // the seeding sweep alone already holds the unchanged field. - dispatch(this.pipelines.diffuse, bg.diffuseSeed); + // 2. Viscous diffusion. The seeding dispatch fills the iterate from the + // advected source; the red-black SOR sweeps then alternate colours, one + // dispatch each, so a red/black pair is one full Gauss–Seidel iteration + // and an even count lands back on the parity it started from. + dispatch(this.pipelines.diffuseSeed, bg.diffuseSeed); this.velocity.parity = 0; - if (runDiffusion) { - for (let i = 1; i < DIFFUSION_ITERATIONS; i++) { - dispatch(this.pipelines.diffuse, bg.diffuse[this.velocity.parity]); - this.velocity.swap(); - } + for (let i = 0; i < viscousSweeps * 2; i++) { + dispatch(i % 2 === 0 ? this.pipelines.diffuseRed : this.pipelines.diffuseBlack, bg.diffuse[this.velocity.parity]); + this.velocity.swap(); } // 3-4. Restore the small-scale vorticity advection dissipated. @@ -407,10 +478,14 @@ export class WebGPUFluidEngine { dispatch(this.pipelines.gradientSubtract, bg.gradientSubtract[this.velocity.parity][this.pressure.parity]); this.velocity.swap(); - // 9-10. Dye: inject at the inflow, then carry it with the finished velocity. + // 9-11. Dye: inject at the inflow, then carry it with the finished velocity + // — through the same MacCormack predictor–corrector the velocity gets, + // because the dye is what the learner actually looks at and a plain + // backward trace smears the bands out well before the far wall. dispatch(this.pipelines.injectDye, bg.injectDye[this.dye.parity]); this.dye.swap(); dispatch(this.pipelines.advectDye, bg.advectDye[this.velocity.parity][this.dye.parity]); + dispatch(this.pipelines.advectDyeCorrect, bg.advectDyeCorrect[this.velocity.parity][this.dye.parity]); this.dye.swap(); pass.end(); @@ -462,6 +537,25 @@ export class WebGPUFluidEngine { }); this.advectTempView = this.advectTemp.createView(); + this.dyeTemp = device.createTexture({ + label: "dye-temp", + size: [grid.width, grid.height], + format: VELOCITY_FORMAT, + usage: storageAndSample, + }); + this.dyeTempView = this.dyeTemp.createView(); + + this.obstacle = device.createTexture({ + label: "obstacle-field", + size: [grid.width, grid.height], + format: SCALAR_FORMAT, + usage: storageAndSample, + }); + this.obstacleView = this.obstacle.createView(); + // A new texture reads as zero, which the solver would take for "solid + // everywhere". Nothing may read it before mask.wgsl has filled it in. + this.isMaskStale = true; + this.divergence = device.createTexture({ label: "divergence", size: [grid.width, grid.height], @@ -487,6 +581,8 @@ export class WebGPUFluidEngine { this.pressure.destroy(); this.velocitySource.destroy(); this.advectTemp.destroy(); + this.dyeTemp.destroy(); + this.obstacle.destroy(); this.divergence.destroy(); this.curl.destroy(); } @@ -496,6 +592,8 @@ export class WebGPUFluidEngine { const layouts = this.pipelines.layouts; const uniform: GPUBindGroupEntry = { binding: 0, resource: { buffer: this.uniformBuffer } }; const sampler: GPUBindGroupEntry = { binding: 1, resource: this.linearSampler }; + /** Bound to every compute layout, whether or not the kernel reads it. */ + const obstacle: GPUBindGroupEntry = { binding: OBSTACLE_BINDING, resource: this.obstacleView }; const velocityViews = this.velocity.views; const dyeViews = this.dye.views; @@ -508,6 +606,8 @@ export class WebGPUFluidEngine { const perParity = (build: (parity: 0 | 1) => T): [T, T] => [build(0), build(1)]; return { + mask: group(layouts.mask, [uniform, { binding: 1, resource: this.obstacleView }]), + // Predicts φ_A: backward-advects velocity by itself into the scratch texture. advectVelocity: perParity((p) => group(layouts.advect, [ @@ -516,12 +616,14 @@ export class WebGPUFluidEngine { { binding: 2, resource: velocityViews[p] }, { binding: 3, resource: velocityViews[p] }, { binding: 4, resource: this.advectTempView }, + { binding: 5, resource: velocityViews[p] }, + obstacle, ]), ), - // MacCormack corrector: reads φ^n + the trace velocity from velocity, φ_A - // from the scratch texture, and writes the corrected velocity to the - // diffusion source. velocityTex and sourceTex share the advect layout. + // MacCormack corrector: traces with φⁿ (binding 2), reads φ_A from the + // scratch texture (binding 3) and φⁿ itself again (binding 5), and writes + // the corrected velocity to the diffusion source. advectVelocityCorrect: perParity((p) => group(layouts.advect, [ uniform, @@ -529,15 +631,18 @@ export class WebGPUFluidEngine { { binding: 2, resource: velocityViews[p] }, { binding: 3, resource: this.advectTempView }, { binding: 4, resource: this.velocitySourceView }, + { binding: 5, resource: velocityViews[p] }, + obstacle, ]), ), - // First diffusion sweep: iterate seeded from the source itself. + // Seeding dispatch: iterate seeded from the advected source itself. diffuseSeed: group(layouts.twoInRGBA, [ uniform, { binding: 1, resource: this.velocitySourceView }, { binding: 2, resource: this.velocitySourceView }, { binding: 3, resource: velocityViews[0] }, + obstacle, ]), diffuse: perParity((p) => @@ -546,6 +651,7 @@ export class WebGPUFluidEngine { { binding: 1, resource: this.velocitySourceView }, { binding: 2, resource: velocityViews[p] }, { binding: 3, resource: velocityViews[p === 0 ? 1 : 0] }, + obstacle, ]), ), @@ -554,6 +660,7 @@ export class WebGPUFluidEngine { uniform, { binding: 1, resource: velocityViews[p] }, { binding: 2, resource: this.curlView }, + obstacle, ]), ), @@ -563,6 +670,7 @@ export class WebGPUFluidEngine { { binding: 1, resource: velocityViews[p] }, { binding: 2, resource: this.curlView }, { binding: 3, resource: velocityViews[p === 0 ? 1 : 0] }, + obstacle, ]), ), @@ -571,6 +679,7 @@ export class WebGPUFluidEngine { uniform, { binding: 1, resource: velocityViews[p] }, { binding: 2, resource: velocityViews[p === 0 ? 1 : 0] }, + obstacle, ]), ), @@ -579,6 +688,7 @@ export class WebGPUFluidEngine { uniform, { binding: 1, resource: velocityViews[p] }, { binding: 2, resource: this.divergenceView }, + obstacle, ]), ), @@ -588,6 +698,7 @@ export class WebGPUFluidEngine { { binding: 1, resource: pressureViews[p] }, { binding: 2, resource: this.divergenceView }, { binding: 3, resource: pressureViews[p === 0 ? 1 : 0] }, + obstacle, ]), ), @@ -598,6 +709,7 @@ export class WebGPUFluidEngine { { binding: 1, resource: velocityViews[v] }, { binding: 2, resource: pressureViews[p] }, { binding: 3, resource: velocityViews[v === 0 ? 1 : 0] }, + obstacle, ]), ), ), @@ -607,9 +719,11 @@ export class WebGPUFluidEngine { uniform, { binding: 1, resource: dyeViews[d] }, { binding: 2, resource: dyeViews[d === 0 ? 1 : 0] }, + obstacle, ]), ), + // Dye predictor: traced by the finished velocity, into the dye scratch. advectDye: perParity((v) => perParity((d) => group(layouts.advect, [ @@ -617,7 +731,25 @@ export class WebGPUFluidEngine { sampler, { binding: 2, resource: velocityViews[v] }, { binding: 3, resource: dyeViews[d] }, + { binding: 4, resource: this.dyeTempView }, + { binding: 5, resource: dyeViews[d] }, + obstacle, + ]), + ), + ), + + // Dye corrector: unlike the velocity's, the field being carried and the + // field doing the carrying are different textures — hence binding 5. + advectDyeCorrect: perParity((v) => + perParity((d) => + group(layouts.advect, [ + uniform, + sampler, + { binding: 2, resource: velocityViews[v] }, + { binding: 3, resource: this.dyeTempView }, { binding: 4, resource: dyeViews[d === 0 ? 1 : 0] }, + { binding: 5, resource: dyeViews[d] }, + obstacle, ]), ), ), @@ -657,22 +789,18 @@ export class WebGPUFluidEngine { // ── Pipeline construction ───────────────────────────────────────────────────── -type Layouts = { - readonly advect: GPUBindGroupLayout; - readonly twoInRGBA: GPUBindGroupLayout; - readonly oneInRGBA: GPUBindGroupLayout; - readonly oneInScalar: GPUBindGroupLayout; - readonly mixedRGBA: GPUBindGroupLayout; - readonly twoInScalar: GPUBindGroupLayout; - readonly display: GPUBindGroupLayout; -}; +type Layouts = { readonly [K in BindLayoutName]: GPUBindGroupLayout }; type Pipelines = { readonly layouts: Layouts; + readonly mask: GPUComputePipeline; readonly advectVelocity: GPUComputePipeline; readonly advectVelocityCorrect: GPUComputePipeline; readonly advectDye: GPUComputePipeline; - readonly diffuse: GPUComputePipeline; + readonly advectDyeCorrect: GPUComputePipeline; + readonly diffuseSeed: GPUComputePipeline; + readonly diffuseRed: GPUComputePipeline; + readonly diffuseBlack: GPUComputePipeline; readonly curl: GPUComputePipeline; readonly vorticity: GPUComputePipeline; readonly forces: GPUComputePipeline; @@ -685,6 +813,7 @@ type Pipelines = { }; type BindGroups = { + readonly mask: GPUBindGroup; readonly advectVelocity: [GPUBindGroup, GPUBindGroup]; readonly advectVelocityCorrect: [GPUBindGroup, GPUBindGroup]; readonly diffuseSeed: GPUBindGroup; @@ -697,97 +826,45 @@ type BindGroups = { readonly gradientSubtract: [[GPUBindGroup, GPUBindGroup], [GPUBindGroup, GPUBindGroup]]; readonly injectDye: [GPUBindGroup, GPUBindGroup]; readonly advectDye: [[GPUBindGroup, GPUBindGroup], [GPUBindGroup, GPUBindGroup]]; + readonly advectDyeCorrect: [[GPUBindGroup, GPUBindGroup], [GPUBindGroup, GPUBindGroup]]; readonly display: [ [[GPUBindGroup, GPUBindGroup], [GPUBindGroup, GPUBindGroup]], [[GPUBindGroup, GPUBindGroup], [GPUBindGroup, GPUBindGroup]], ]; }; -const COMPUTE = GPUShaderStage.COMPUTE; -const FRAGMENT = GPUShaderStage.FRAGMENT; - -const uniformEntry = (visibility: number): GPUBindGroupLayoutEntry => ({ - binding: 0, - visibility, - buffer: { type: "uniform" }, -}); - -const filterable = (binding: number, visibility: number): GPUBindGroupLayoutEntry => ({ - binding, - visibility, - texture: { sampleType: "float" }, -}); +/** The half of a layout entry that says what kind of resource is bound. */ +function layoutResource(binding: BindingSpec): Omit { + if (binding.kind === "uniform") { + return { buffer: { type: "uniform" } }; + } + if (binding.kind === "sampler") { + return { sampler: { type: "filtering" } }; + } + if (binding.kind === "texture") { + return { texture: { sampleType: binding.sampleType } }; + } + return { storageTexture: { access: "write-only", format: binding.format } }; +} -/** - * r32float is not guaranteed filterable, so pressure, divergence and curl must - * be declared unfilterable and read with textureLoad. - */ -const unfilterable = (binding: number, visibility: number): GPUBindGroupLayoutEntry => ({ - binding, - visibility, - texture: { sampleType: "unfilterable-float" }, -}); - -const storageOut = (binding: number, format: GPUTextureFormat): GPUBindGroupLayoutEntry => ({ - binding, - visibility: COMPUTE, - storageTexture: { access: "write-only", format }, -}); +/** Turns one of the plain-data layout specs in bindLayouts.ts into a real one. */ +function createLayout(device: GPUDevice, spec: BindLayoutSpec): GPUBindGroupLayout { + const visibility = spec.stage === "compute" ? GPUShaderStage.COMPUTE : GPUShaderStage.FRAGMENT; + + return device.createBindGroupLayout({ + label: spec.label, + entries: Object.entries(spec.bindings).map(([index, binding]) => ({ + binding: Number(index), + visibility, + ...layoutResource(binding), + })), + }); +} function createPipelines(device: GPUDevice, canvasFormat: GPUTextureFormat): Pipelines { - const layouts: Layouts = { - advect: device.createBindGroupLayout({ - label: "advect", - entries: [ - uniformEntry(COMPUTE), - { binding: 1, visibility: COMPUTE, sampler: { type: "filtering" } }, - filterable(2, COMPUTE), - filterable(3, COMPUTE), - storageOut(4, VELOCITY_FORMAT), - ], - }), - twoInRGBA: device.createBindGroupLayout({ - label: "two-in-rgba", - entries: [uniformEntry(COMPUTE), filterable(1, COMPUTE), filterable(2, COMPUTE), storageOut(3, VELOCITY_FORMAT)], - }), - oneInRGBA: device.createBindGroupLayout({ - label: "one-in-rgba", - entries: [uniformEntry(COMPUTE), filterable(1, COMPUTE), storageOut(2, VELOCITY_FORMAT)], - }), - oneInScalar: device.createBindGroupLayout({ - label: "one-in-scalar", - entries: [uniformEntry(COMPUTE), filterable(1, COMPUTE), storageOut(2, SCALAR_FORMAT)], - }), - mixedRGBA: device.createBindGroupLayout({ - label: "mixed-rgba", - entries: [ - uniformEntry(COMPUTE), - filterable(1, COMPUTE), - unfilterable(2, COMPUTE), - storageOut(3, VELOCITY_FORMAT), - ], - }), - twoInScalar: device.createBindGroupLayout({ - label: "two-in-scalar", - entries: [ - uniformEntry(COMPUTE), - unfilterable(1, COMPUTE), - unfilterable(2, COMPUTE), - storageOut(3, SCALAR_FORMAT), - ], - }), - display: device.createBindGroupLayout({ - label: "display", - entries: [ - uniformEntry(FRAGMENT), - { binding: 1, visibility: FRAGMENT, sampler: { type: "filtering" } }, - filterable(2, FRAGMENT), - filterable(3, FRAGMENT), - unfilterable(4, FRAGMENT), - unfilterable(5, FRAGMENT), - ], - }), - }; + const layouts = Object.fromEntries( + Object.entries(BIND_LAYOUTS).map(([name, spec]) => [name, createLayout(device, spec)]), + ) as Layouts; // WGSL has no include mechanism, so the shared struct and helpers are // concatenated ahead of every shader. @@ -805,11 +882,14 @@ function createPipelines(device: GPUDevice, canvasFormat: GPUTextureFormat): Pip const advectLayout = device.createPipelineLayout({ bindGroupLayouts: [layouts.advect] }); const pressureModule = module("pressure", pressureWGSL); const pressureLayout = device.createPipelineLayout({ bindGroupLayouts: [layouts.twoInScalar] }); + const diffuseModule = module("diffuse", diffuseWGSL); + const diffuseLayout = device.createPipelineLayout({ bindGroupLayouts: [layouts.twoInRGBA] }); const displayModule = module("display", displayWGSL); return { layouts, + mask: compute("mask", layouts.mask, maskWGSL, "main"), advectVelocity: device.createComputePipeline({ label: "advect-velocity", layout: advectLayout, @@ -825,7 +905,26 @@ function createPipelines(device: GPUDevice, canvasFormat: GPUTextureFormat): Pip layout: advectLayout, compute: { module: advectModule, entryPoint: "advectDye" }, }), - diffuse: compute("diffuse", layouts.twoInRGBA, diffuseWGSL, "main"), + advectDyeCorrect: device.createComputePipeline({ + label: "advect-dye-correct", + layout: advectLayout, + compute: { module: advectModule, entryPoint: "advectDyeCorrect" }, + }), + diffuseSeed: device.createComputePipeline({ + label: "diffuse-seed", + layout: diffuseLayout, + compute: { module: diffuseModule, entryPoint: "seed" }, + }), + diffuseRed: device.createComputePipeline({ + label: "diffuse-red", + layout: diffuseLayout, + compute: { module: diffuseModule, entryPoint: "solveRed" }, + }), + diffuseBlack: device.createComputePipeline({ + label: "diffuse-black", + layout: diffuseLayout, + compute: { module: diffuseModule, entryPoint: "solveBlack" }, + }), curl: compute("curl", layouts.oneInScalar, curlWGSL, "main"), vorticity: compute("vorticity", layouts.mixedRGBA, vorticityWGSL, "main"), forces: compute("forces", layouts.oneInRGBA, forcesWGSL, "main"), diff --git a/src/common/gpu/bindLayouts.ts b/src/common/gpu/bindLayouts.ts new file mode 100644 index 0000000..f95cf45 --- /dev/null +++ b/src/common/gpu/bindLayouts.ts @@ -0,0 +1,162 @@ +/** + * bindLayouts.ts + * + * The bind group layouts the solver's pipelines are built from, as plain data. + * + * WebGPU checks a bind *group* against its layout, and a shader's declared + * resources against the pipeline layout — but only for resources the entry point + * statically uses, and only at pipeline-creation time, on a device. That is a + * long way from the code being edited: adding a binding to a kernel and + * forgetting to add it to the layout it shares with five other kernels produces + * a validation error at startup, on hardware, and nothing at all in the type + * checker. + * + * So the layouts live here as data rather than as `device.createBindGroupLayout` + * calls buried in the engine, and `tests/ShaderBindings.test.ts` parses the WGSL + * to check that every `@group(0) @binding(n)` a shader declares exists in the + * layout its pipelines use, with a compatible kind and storage format. That test + * runs in Vitest with no GPU. + * + * Several kernels share a layout, and a layout may legitimately carry entries a + * given kernel does not use — the obstacle field is bound to every compute + * layout, but the dye injection kernel has no use for it. The check is therefore + * one-directional: shader ⊆ layout. + */ + +/** Texture formats the fields use, shared with WebGPUFluidEngine. */ +export const VELOCITY_FORMAT: GPUTextureFormat = "rgba16float"; +export const SCALAR_FORMAT: GPUTextureFormat = "r32float"; + +/** + * Binding index of the baked obstacle field, in every compute layout that has + * one. A single high index rather than "the next free one per layout", so the + * number means the same thing in every kernel that reads it. + */ +export const OBSTACLE_BINDING = 6; + +/** What a binding holds, in the terms both WebGPU and WGSL can be checked in. */ +export type BindingSpec = + | { readonly kind: "uniform" } + | { readonly kind: "sampler" } + | { readonly kind: "texture"; readonly sampleType: "float" | "unfilterable-float" } + | { readonly kind: "storageTexture"; readonly format: GPUTextureFormat }; + +export type BindLayoutSpec = { + readonly label: string; + /** Which shader stage the layout's entries are visible to. */ + readonly stage: "compute" | "fragment"; + readonly bindings: Readonly>; +}; + +const uniform: BindingSpec = { kind: "uniform" }; +const sampler: BindingSpec = { kind: "sampler" }; + +/** A texture read through the filtering sampler, so it must be filterable. */ +const filterable: BindingSpec = { kind: "texture", sampleType: "float" }; + +/** + * r32float is not guaranteed filterable, so pressure, divergence, curl and the + * obstacle field are declared unfilterable and read only with textureLoad. + */ +const unfilterable: BindingSpec = { kind: "texture", sampleType: "unfilterable-float" }; + +const velocityOut: BindingSpec = { kind: "storageTexture", format: VELOCITY_FORMAT }; +const scalarOut: BindingSpec = { kind: "storageTexture", format: SCALAR_FORMAT }; + +export const BIND_LAYOUTS = { + /** Bakes the obstacle SDF. The one kernel that does not read the result. */ + mask: { + label: "mask", + stage: "compute", + bindings: { 0: uniform, 1: scalarOut }, + }, + + /** + * Semi-Lagrangian transport. 2 is the field the trace follows, 3 the field + * being carried, 5 the pre-advection field the MacCormack corrector needs. + */ + advect: { + label: "advect", + stage: "compute", + bindings: { + 0: uniform, + 1: sampler, + 2: filterable, + 3: filterable, + 4: velocityOut, + 5: filterable, + [OBSTACLE_BINDING]: unfilterable, + }, + }, + + /** Diffusion: the fixed source, the current iterate, the next one. */ + twoInRGBA: { + label: "two-in-rgba", + stage: "compute", + bindings: { 0: uniform, 1: filterable, 2: filterable, 3: velocityOut, [OBSTACLE_BINDING]: unfilterable }, + }, + + oneInRGBA: { + label: "one-in-rgba", + stage: "compute", + bindings: { 0: uniform, 1: filterable, 2: velocityOut, [OBSTACLE_BINDING]: unfilterable }, + }, + + oneInScalar: { + label: "one-in-scalar", + stage: "compute", + bindings: { 0: uniform, 1: filterable, 2: scalarOut, [OBSTACLE_BINDING]: unfilterable }, + }, + + mixedRGBA: { + label: "mixed-rgba", + stage: "compute", + bindings: { 0: uniform, 1: filterable, 2: unfilterable, 3: velocityOut, [OBSTACLE_BINDING]: unfilterable }, + }, + + twoInScalar: { + label: "two-in-scalar", + stage: "compute", + bindings: { 0: uniform, 1: unfilterable, 2: unfilterable, 3: scalarOut, [OBSTACLE_BINDING]: unfilterable }, + }, + + display: { + label: "display", + stage: "fragment", + bindings: { 0: uniform, 1: sampler, 2: filterable, 3: filterable, 4: unfilterable, 5: unfilterable }, + }, +} as const satisfies Record; + +export type BindLayoutName = keyof typeof BIND_LAYOUTS; + +export const BIND_LAYOUT_NAMES = Object.keys(BIND_LAYOUTS) as BindLayoutName[]; + +/** + * What a layout binds at an index, or undefined if it binds nothing there. + * + * `as const satisfies` above narrows each layout's bindings to a literal record, + * which is what makes a typo in a layout a compile error — but it also means the + * records cannot be indexed by a computed number. This is the one place that + * widening happens. + */ +export function layoutBinding(name: BindLayoutName, binding: number): BindingSpec | undefined { + return (BIND_LAYOUTS[name].bindings as Readonly>)[binding]; +} + +/** + * Which layout each shader's pipelines are created with. Every shader file must + * appear here, so a new one cannot quietly escape the binding check. + */ +export const SHADER_LAYOUTS = { + "mask.wgsl": "mask", + "advect.wgsl": "advect", + "diffuse.wgsl": "twoInRGBA", + "curl.wgsl": "oneInScalar", + "vorticity.wgsl": "mixedRGBA", + "forces.wgsl": "oneInRGBA", + "divergence.wgsl": "oneInScalar", + "pressure.wgsl": "twoInScalar", + "gradientSubtract.wgsl": "mixedRGBA", + "dye.wgsl": "oneInRGBA", + "display.wgsl": "display", +} as const satisfies Record; diff --git a/src/common/gpu/shaders/advect.wgsl b/src/common/gpu/shaders/advect.wgsl index a48885e..abb0248 100644 --- a/src/common/gpu/shaders/advect.wgsl +++ b/src/common/gpu/shaders/advect.wgsl @@ -5,90 +5,121 @@ // the velocity field and take whatever was there. Unconditionally stable at any // timestep, at the cost of being diffusive. // -// The backward trace alone (the predictor, advectVelocity) behaves like an extra -// viscosity of order h·|u|/2 — enough to smear the shear layer off the obstacle -// and weaken the Kármán street. So the velocity field is advected with a -// MacCormack predictor–corrector (Selle, Fedkiw, Lanson, Molemaker & Bridson, -// 2008): the predictor's backward step is corrected by an anti-diffusive term -// built from a forward step, cancelling most of the numerical diffusion. A -// bound-preserving limiter clamps the result to the predictor's local range so -// the correction cannot create new extrema and destabilise the flow. +// Two things are done about that cost, and both apply to the velocity field and +// to the dye: // -// The dye is a passive tracer and stays on the plain backward trace: it is -// carried by the (now sharper) velocity field, so it sharpens for free, and the -// dye texture needs no extra storage. +// 1. The trace itself is a midpoint (RK2) step rather than a single Euler step. +// An Euler trace follows the velocity at the arrival point for the whole +// step, which cuts the corner on every curved path — and a vortex is nothing +// but curved paths. Sampling the velocity once more at the halfway point +// makes the trace second order in Δt for the cost of one extra fetch, and +// the vortices of the Kármán street stop drifting toward their own centres. // -// The bilinear filtering in textureSampleLevel is what makes the backtrace -// second-order in space; it is also why the velocity and dye textures are -// rgba16float (filterable) rather than rg32float (not, without an optional -// feature). +// 2. The plain backward trace behaves like an extra viscosity of order h·|u|/2 +// — enough to smear the shear layer off the obstacle and wash the dye bands +// out well before they reach the far end of the channel. So both fields are +// advected with a MacCormack predictor–corrector (Selle, Fedkiw, Lanson, +// Molemaker & Bridson, 2008): the predictor's backward step is corrected by +// an anti-diffusive term built from a forward step, cancelling most of the +// numerical diffusion. A bound-preserving limiter clamps the result to the +// predictor's local range, so the correction cannot create new extrema — +// which for the dye also means it can neither go negative nor invent a +// colour that was not already in the neighbourhood. +// +// The bilinear filtering in textureSampleLevel is what makes the trace second +// order in space; it is also why the velocity and dye textures are rgba16float +// (filterable) rather than rg32float (not, without an optional feature). +// +// ── Bindings ────────────────────────────────────────────────────────────────── +// velocityTex the field the trace follows — always a velocity field +// sourceTex the field being carried: φⁿ in a predictor, φ_A in a corrector +// priorTex φⁿ again, in a corrector; unused by the predictors +// outTex the result +// +// Velocity advection binds velocityTex and priorTex to the same texture, since +// the field being carried is the one doing the carrying. Dye advection does not. @group(0) @binding(0) var u : SimUniforms; @group(0) @binding(1) var linearSampler : sampler; @group(0) @binding(2) var velocityTex : texture_2d; @group(0) @binding(3) var sourceTex : texture_2d; @group(0) @binding(4) var outTex : texture_storage_2d; +@group(0) @binding(5) var priorTex : texture_2d; +@group(0) @binding(6) var obstacleTex : texture_2d; -// Where the fluid now at `uv` came from one timestep ago. -// Velocity is in m/s and uv is dimensionless, so the displacement is divided by -// the domain size to convert metres to uv. +// Velocity is in m/s and uv is dimensionless, so a displacement in metres has to +// be divided by the domain size to become a displacement in uv. +fn traceOffset() -> vec2 { + return vec2(u.dt, u.dt) / u.domainSize; +} + +fn traceVelocity(uv: vec2) -> vec2 { + return textureSampleLevel(velocityTex, linearSampler, uv, 0.0).xy; +} + +// Where the fluid now at `uv` came from one timestep ago, by the midpoint rule. fn backtrace(uv: vec2) -> vec2 { - let velocity = textureSampleLevel(velocityTex, linearSampler, uv, 0.0).xy; - return uv - velocity * u.dt / u.domainSize; + let offset = traceOffset(); + let midpoint = uv - 0.5 * traceVelocity(uv) * offset; + return uv - traceVelocity(midpoint) * offset; +} + +// The same step taken forwards, which is what makes the MacCormack error +// estimate below an estimate of the *backward* step's error. +fn forwardTrace(uv: vec2) -> vec2 { + let offset = traceOffset(); + let midpoint = uv + 0.5 * traceVelocity(uv) * offset; + return uv + traceVelocity(midpoint) * offset; +} + +fn isSolidHere(cell: vec2) -> bool { + return isSolidAt(obstacleTex, vec2(cell), u); } +// ── Predictors ──────────────────────────────────────────────────────────────── + @compute @workgroup_size(8, 8) fn advectVelocity(@builtin(global_invocation_id) id: vec3) { let cell = id.xy; - if (cell.x >= u32(u.gridSize.x) || cell.y >= u32(u.gridSize.y)) { + if (!inGrid(cell, u)) { return; } - let uv = cellToUV(cell, u); - // No-slip: the body does not move, so neither does the fluid inside it. // Enforced here as well as in the dedicated pass, so the backtrace never // picks up a stale velocity from inside the obstacle. - if (isSolid(uvToMetres(uv, u), u)) { + if (isSolidHere(cell)) { textureStore(outTex, cell, vec4(0.0, 0.0, 0.0, 1.0)); return; } - let source = textureSampleLevel(sourceTex, linearSampler, backtrace(uv), 0.0); + let source = textureSampleLevel(sourceTex, linearSampler, backtrace(cellToUV(cell, u)), 0.0); textureStore(outTex, cell, vec4(source.xy, 0.0, 1.0)); } @compute @workgroup_size(8, 8) fn advectDye(@builtin(global_invocation_id) id: vec3) { let cell = id.xy; - if (cell.x >= u32(u.gridSize.x) || cell.y >= u32(u.gridSize.y)) { + if (!inGrid(cell, u)) { return; } - let uv = cellToUV(cell, u); - // Dye inside the body would be visible through the obstacle overlay's edges. - if (isSolid(uvToMetres(uv, u), u)) { + if (isSolidHere(cell)) { textureStore(outTex, cell, vec4(0.0, 0.0, 0.0, 0.0)); return; } - // dyeDissipation is the fraction remaining after one second, so the per-step - // factor is that raised to dt — framerate-independent fading. - let decay = pow(u.dyeDissipation, u.dt); - let source = textureSampleLevel(sourceTex, linearSampler, backtrace(uv), 0.0); - textureStore(outTex, cell, source * decay); + // No dissipation here: the corrector applies it once, to the finished value. + let source = textureSampleLevel(sourceTex, linearSampler, backtrace(cellToUV(cell, u)), 0.0); + textureStore(outTex, cell, source); } -// ── MacCormack corrector for the velocity field ─────────────────────────────── -// -// Bindings: velocityTex holds the pre-advection field φ^n (and is sampled for the -// trace velocity); sourceTex holds the predictor result φ_A (advectTemp); -// outTex receives the corrected, limited velocity (velocitySource). +// ── MacCormack corrector ────────────────────────────────────────────────────── // -// φ_A = backward advect of φ^n (the predictor, above) +// φ_A = backward advect of φⁿ (the predictor, above) // φ_D = forward sample of φ_A at x + u·dt (the reverse trace) -// φ_raw = φ_A + ½(φ^n(x) − φ_D) (anti-diffusive blend) +// φ_raw = φ_A + ½(φⁿ(x) − φ_D) (anti-diffusive blend) // φ = clamp φ_raw to [min, max] of φ_A around the departure point // // The clamp is what makes the scheme unconditionally stable: without it the @@ -96,46 +127,66 @@ fn advectDye(@builtin(global_invocation_id) id: vec3) { // corrected value never leaves the range the predictor already produced, so the // whole step is no less stable than plain semi-Lagrangian. -fn forwardTrace(uv: vec2) -> vec2 { - let velocity = textureSampleLevel(velocityTex, linearSampler, uv, 0.0).xy; - return uv + velocity * u.dt / u.domainSize; +fn maccormack(cell: vec2) -> vec4 { + let uv = cellToUV(cell, u); + + let phiN = textureLoad(priorTex, vec2(cell), 0); + let phiA = textureSampleLevel(sourceTex, linearSampler, uv, 0.0); + let phiD = textureSampleLevel(sourceTex, linearSampler, forwardTrace(uv), 0.0); + + let raw = phiA + 0.5 * (phiN - phiD); + + // Bound the corrected value to the range of φ_A over the four texels whose + // bilinear interpolation produced the predictor's value at the departure + // point. A value outside this range is an overshoot the limiter removes. + // + // The −0.5 is the half-texel offset between a uv coordinate and the texel + // *centres* the hardware interpolates between: without it the four cells + // sampled here are the wrong ones whenever the departure point falls in the + // far half of its cell, and the limiter bounds the correction against a + // neighbourhood the predictor never looked at. + let base = clampCell(vec2(floor(backtrace(uv) * u.gridSize - vec2(0.5, 0.5))), u); + let p00 = textureLoad(sourceTex, clampCell(base + vec2(0, 0), u), 0); + let p10 = textureLoad(sourceTex, clampCell(base + vec2(1, 0), u), 0); + let p01 = textureLoad(sourceTex, clampCell(base + vec2(0, 1), u), 0); + let p11 = textureLoad(sourceTex, clampCell(base + vec2(1, 1), u), 0); + + return clamp(raw, min(min(p00, p10), min(p01, p11)), max(max(p00, p10), max(p01, p11))); } @compute @workgroup_size(8, 8) fn advectVelocityCorrect(@builtin(global_invocation_id) id: vec3) { let cell = id.xy; - if (cell.x >= u32(u.gridSize.x) || cell.y >= u32(u.gridSize.y)) { + if (!inGrid(cell, u)) { return; } // No-slip is enforced by the predictor too; repeated here so a correction can // never put velocity back inside the body. - if (isSolid(uvToMetres(cellToUV(cell, u), u), u)) { + if (isSolidHere(cell)) { textureStore(outTex, cell, vec4(0.0, 0.0, 0.0, 1.0)); return; } - let c = vec2(cell); - let uv = cellToUV(cell, u); + textureStore(outTex, cell, vec4(maccormack(cell).xy, 0.0, 1.0)); +} - let phiN = textureLoad(velocityTex, c, 0).xy; - let phiA = textureSampleLevel(sourceTex, linearSampler, uv, 0.0).xy; - let phiD = textureSampleLevel(sourceTex, linearSampler, forwardTrace(uv), 0.0).xy; +@compute @workgroup_size(8, 8) +fn advectDyeCorrect(@builtin(global_invocation_id) id: vec3) { + let cell = id.xy; + if (!inGrid(cell, u)) { + return; + } - let raw = phiA + 0.5 * (phiN - phiD); + if (isSolidHere(cell)) { + textureStore(outTex, cell, vec4(0.0, 0.0, 0.0, 0.0)); + return; + } - // Bound the corrected value to the range of φ_A over the four texels whose - // bilinear interpolation produced the predictor's value at the departure - // point. A value outside this range is an overshoot the limiter removes. - let backUV = backtrace(uv); - let base = clampCell(vec2(floor(backUV * u.gridSize)), u); - let p00 = textureLoad(sourceTex, clampCell(base + vec2(0, 0), u), 0).xy; - let p10 = textureLoad(sourceTex, clampCell(base + vec2(1, 0), u), 0).xy; - let p01 = textureLoad(sourceTex, clampCell(base + vec2(0, 1), u), 0).xy; - let p11 = textureLoad(sourceTex, clampCell(base + vec2(1, 1), u), 0).xy; - let lo = min(min(p00, p10), min(p01, p11)); - let hi = max(max(p00, p10), max(p01, p11)); - - let corrected = clamp(raw, lo, hi); - textureStore(outTex, cell, vec4(corrected, 0.0, 1.0)); + // dyeDissipation is the fraction remaining after one second, so the per-step + // factor is that raised to dt — framerate-independent fading. Applied once, + // here, so the predictor's output is a clean φ_A for the corrector to bound + // against rather than a field that has already faded by a different amount. + let decay = pow(u.dyeDissipation, u.dt); + textureStore(outTex, cell, maccormack(cell) * decay); } diff --git a/src/common/gpu/shaders/common.wgsl b/src/common/gpu/shaders/common.wgsl index 40c1b99..1a6100a 100644 --- a/src/common/gpu/shaders/common.wgsl +++ b/src/common/gpu/shaders/common.wgsl @@ -125,10 +125,6 @@ fn obstacleSDF(p: vec2, uniforms: SimUniforms) -> f32 { return 1.0e9; } -fn isSolid(p: vec2, uniforms: SimUniforms) -> bool { - return obstacleSDF(p, uniforms) <= 0.0; -} - // ── Grid helpers ────────────────────────────────────────────────────────────── // True when the invocation is inside the grid. Dispatch sizes are rounded up to @@ -150,12 +146,32 @@ fn uvToCell(uv: vec2, uniforms: SimUniforms) -> vec2 { return clampCell(vec2(floor(uv * uniforms.gridSize)), uniforms); } -fn isSolidCell(cell: vec2, uniforms: SimUniforms) -> bool { - let clamped = clampCell(cell, uniforms); - return isSolid(uvToMetres(cellToUV(vec2(clamped), uniforms), uniforms), uniforms); -} - // Cell edge length in metres, the h in every finite difference below. fn cellSize(uniforms: SimUniforms) -> f32 { return uniforms.domainSize.x / uniforms.gridSize.x; } + +// ── The baked obstacle field ────────────────────────────────────────────────── +// +// obstacleSDF() above is not cheap — the airfoil alone costs a sqrt, a sin, a +// cos and a quartic — and the solver asks "is this cell solid?" up to five times +// per cell in nearly every one of the ~50 dispatches that make up a frame. The +// pressure solve alone accounts for most of them. +// +// So the distance is evaluated once per cell by mask.wgsl, whenever the obstacle +// actually changes, and every compute kernel reads that texture instead. The +// analytic function is still used by the display pass, which samples between +// cells and needs sub-cell accuracy for the body's outline. +// +// The texture is a parameter rather than a global here because the binding is +// declared by each kernel that needs one, and display.wgsl has none. + +fn obstacleDistanceAt(field: texture_2d, cell: vec2, uniforms: SimUniforms) -> f32 { + return textureLoad(field, clampCell(cell, uniforms), 0).x; +} + +// Clamping before the lookup gives a cell outside the grid the solidity of the +// edge cell nearest it, which is what the boundary stencils below rely on. +fn isSolidAt(field: texture_2d, cell: vec2, uniforms: SimUniforms) -> bool { + return obstacleDistanceAt(field, cell, uniforms) <= 0.0; +} diff --git a/src/common/gpu/shaders/curl.wgsl b/src/common/gpu/shaders/curl.wgsl index e4f1235..6fa433c 100644 --- a/src/common/gpu/shaders/curl.wgsl +++ b/src/common/gpu/shaders/curl.wgsl @@ -8,11 +8,12 @@ @group(0) @binding(0) var u : SimUniforms; @group(0) @binding(1) var velocityTex : texture_2d; @group(0) @binding(2) var outTex : texture_storage_2d; +@group(0) @binding(6) var obstacleTex : texture_2d; fn velocityAt(cell: vec2) -> vec2 { // Solid cells read as stationary, so the shear layer forms at the body's // surface rather than at the first fluid cell outside it. - if (isSolidCell(cell, u)) { + if (isSolidAt(obstacleTex, cell, u)) { return vec2(0.0, 0.0); } return textureLoad(velocityTex, clampCell(cell, u), 0).xy; diff --git a/src/common/gpu/shaders/diffuse.wgsl b/src/common/gpu/shaders/diffuse.wgsl index 54d841d..441677b 100644 --- a/src/common/gpu/shaders/diffuse.wgsl +++ b/src/common/gpu/shaders/diffuse.wgsl @@ -1,48 +1,144 @@ -// diffuse.wgsl — implicit viscous diffusion, one Jacobi sweep. +// diffuse.wgsl — implicit viscous diffusion, one red-black SOR sweep. // // The viscous term ν∇²u is solved implicitly: (I − νΔt∇²)u_new = u_old. Solving // it explicitly would need Δt < h²/(4ν), which at high viscosity is far smaller -// than a frame. The implicit form is unconditionally stable and reduces to the -// same Jacobi iteration used for pressure: +// than a frame. The implicit form is unconditionally stable and reduces to a +// linear system with one number in it: // -// u_new = (u_old + α·Σ neighbours) / (1 + 4α), α = νΔt/h² +// (1 + 4α)·u_c − α·Σ u_n = u_source, α = νΔt/h² // -// This is the term the Reynolds-number slider actually moves. At low viscosity α -// is tiny and the sweep is nearly the identity, which is correct — there the -// numerical diffusion of the advection step dominates, and that ceiling on the -// achievable Reynolds number is documented in doc/model.md. +// This is the term the Reynolds-number slider actually moves, so how well it is +// solved decides whether the fluid is as viscous as the readout claims. +// +// ── Why not Jacobi ──────────────────────────────────────────────────────────── +// α grows with the square of the resolution: at ν = 10⁻³ m²/s and Δt = 1/60 s it +// is 0.27 on the 256 × 128 grid and 17 on the 2048 × 1024 one. Jacobi's error +// reduction is 4α/(1 + 4α) per sweep — 0.52 at the first and 0.986 at the second +// — so a fixed twelve sweeps solved the coarse grid comfortably and barely +// touched the fine one. The fine grids were quietly *less* viscous than asked, +// in the direction that inflates the effective Reynolds number. +// +// Red-black SOR fixes both ends at once. Ordering the grid like a checkerboard +// makes Gauss–Seidel run in parallel — the red sweep updates red cells from +// their black neighbours, the black sweep updates black cells from the red ones +// just written — and over-relaxing past the Gauss–Seidel value squares the error +// reduction again. At the optimal factor the rate becomes ω − 1: 0.08 and 0.71 +// for the same two cases. The host then picks the sweep count from α, so the +// easy end of the range costs less than it used to and the stiff end is right. +// +// ── Boundaries ──────────────────────────────────────────────────────────────── +// Solid cells are held at zero every sweep, so a neighbour inside the body reads +// as stationary fluid: no-slip on the obstacle. clampCell gives the channel +// walls a zero-gradient (free-slip) condition for free. @group(0) @binding(0) var u : SimUniforms; @group(0) @binding(1) var sourceTex : texture_2d; @group(0) @binding(2) var previousTex : texture_2d; @group(0) @binding(3) var outTex : texture_storage_2d; +@group(0) @binding(6) var obstacleTex : texture_2d; -@compute @workgroup_size(8, 8) -fn main(@builtin(global_invocation_id) id: vec3) { +/** + * Optimal over-relaxation factor for this α. + * + * The Jacobi spectral radius of the system above is ρ = 4α/(1 + 4α), and + * Young's formula ω = 2/(1 + √(1 − ρ²)) simplifies — since 1 − ρ² is exactly + * (1 + 8α)/(1 + 4α)² — to the closed form below: one square root, no + * trigonometry, and ω stays inside (1, 2) for every α. + * + * Mirrored by `diffusionOmega` in common/gpu/solverSchedule.ts, which the host + * uses to choose how many sweeps this rate needs. The two must agree. + */ +fn sorOmega(alpha: f32) -> f32 { + let diagonal = 1.0 + 4.0 * alpha; + return 2.0 * diagonal / (diagonal + sqrt(1.0 + 8.0 * alpha)); +} + +/** α = νΔt/h², the stiffness of the system being solved. */ +fn stiffness() -> f32 { + let h = cellSize(u); + return u.viscosity * u.dt / (h * h); +} + +/** The four-neighbour sum of the current iterate, with clamped (Neumann) edges. */ +fn neighbourSum(c: vec2) -> vec2 { + return + textureLoad(previousTex, clampCell(c - vec2(1, 0), u), 0).xy + + textureLoad(previousTex, clampCell(c + vec2(1, 0), u), 0).xy + + textureLoad(previousTex, clampCell(c - vec2(0, 1), u), 0).xy + + textureLoad(previousTex, clampCell(c + vec2(0, 1), u), 0).xy; +} + +/** + * One sweep. `redPass` selects which colour is updated; cells of the other + * colour are copied through verbatim, so a red sweep followed by a black one + * ping-pongs back to the original texture with the whole grid updated. + */ +fn sweep(id: vec3, redPass: bool) { let cell = id.xy; if (!inGrid(cell, u)) { return; } - if (isSolidCell(vec2(cell), u)) { + let c = vec2(cell); + + if (isSolidAt(obstacleTex, c, u)) { textureStore(outTex, cell, vec4(0.0, 0.0, 0.0, 1.0)); return; } + let centre = textureLoad(previousTex, c, 0).xy; + + let isRed = ((c.x + c.y) & 1) == 0; + if (isRed != redPass) { + // Not this sweep's colour: carry the previous value through unchanged. + textureStore(outTex, cell, vec4(centre, 0.0, 1.0)); + return; + } + + // `previousTex` is the last iterate; `sourceTex` is the pre-diffusion field, + // which stays fixed across the whole solve. + let alpha = stiffness(); + let original = textureLoad(sourceTex, c, 0).xy; + + // Gauss–Seidel value, then over-relaxation past it. + let gs = (original + alpha * neighbourSum(c)) / (1.0 + 4.0 * alpha); + let relaxed = centre + sorOmega(alpha) * (gs - centre); + + textureStore(outTex, cell, vec4(relaxed, 0.0, 1.0)); +} + +/** + * The seeding dispatch: a plain Jacobi sweep from the advected field, which both + * writes every cell — giving the sweeps above a complete iterate to start from — + * and is already a better initial guess than the source alone. + */ +@compute @workgroup_size(8, 8) +fn seed(@builtin(global_invocation_id) id: vec3) { + let cell = id.xy; + if (!inGrid(cell, u)) { + return; + } + let c = vec2(cell); - let h = cellSize(u); - let alpha = u.viscosity * u.dt / (h * h); - // `previousTex` is the last Jacobi iterate; `sourceTex` is the pre-diffusion - // field, which stays fixed across the whole solve. - let neighbours = - textureLoad(previousTex, clampCell(c - vec2(1, 0), u), 0).xy + - textureLoad(previousTex, clampCell(c + vec2(1, 0), u), 0).xy + - textureLoad(previousTex, clampCell(c - vec2(0, 1), u), 0).xy + - textureLoad(previousTex, clampCell(c + vec2(0, 1), u), 0).xy; + if (isSolidAt(obstacleTex, c, u)) { + textureStore(outTex, cell, vec4(0.0, 0.0, 0.0, 1.0)); + return; + } + let alpha = stiffness(); let original = textureLoad(sourceTex, c, 0).xy; - let result = (original + alpha * neighbours) / (1.0 + 4.0 * alpha); + let result = (original + alpha * neighbourSum(c)) / (1.0 + 4.0 * alpha); textureStore(outTex, cell, vec4(result, 0.0, 1.0)); } + +@compute @workgroup_size(8, 8) +fn solveRed(@builtin(global_invocation_id) id: vec3) { + sweep(id, true); +} + +@compute @workgroup_size(8, 8) +fn solveBlack(@builtin(global_invocation_id) id: vec3) { + sweep(id, false); +} diff --git a/src/common/gpu/shaders/divergence.wgsl b/src/common/gpu/shaders/divergence.wgsl index 6afc132..567c988 100644 --- a/src/common/gpu/shaders/divergence.wgsl +++ b/src/common/gpu/shaders/divergence.wgsl @@ -12,9 +12,10 @@ @group(0) @binding(0) var u : SimUniforms; @group(0) @binding(1) var velocityTex : texture_2d; @group(0) @binding(2) var outTex : texture_storage_2d; +@group(0) @binding(6) var obstacleTex : texture_2d; fn velocityAt(cell: vec2) -> vec2 { - if (isSolidCell(cell, u)) { + if (isSolidAt(obstacleTex, cell, u)) { return vec2(0.0, 0.0); } return textureLoad(velocityTex, clampCell(cell, u), 0).xy; diff --git a/src/common/gpu/shaders/forces.wgsl b/src/common/gpu/shaders/forces.wgsl index 5e9a9a4..1d8bc9c 100644 --- a/src/common/gpu/shaders/forces.wgsl +++ b/src/common/gpu/shaders/forces.wgsl @@ -16,6 +16,7 @@ @group(0) @binding(0) var u : SimUniforms; @group(0) @binding(1) var velocityTex : texture_2d; @group(0) @binding(2) var outTex : texture_storage_2d; +@group(0) @binding(6) var obstacleTex : texture_2d; // Width of the inflow and outflow boundary strips, in cells. const BOUNDARY_CELLS: i32 = 2; @@ -34,7 +35,7 @@ fn main(@builtin(global_invocation_id) id: vec3) { let width = i32(u.gridSize.x); let height = i32(u.gridSize.y); - if (isSolidCell(c, u)) { + if (isSolidAt(obstacleTex, c, u)) { textureStore(outTex, cell, vec4(0.0, 0.0, 0.0, 1.0)); return; } diff --git a/src/common/gpu/shaders/gradientSubtract.wgsl b/src/common/gpu/shaders/gradientSubtract.wgsl index 32cb253..79fc94a 100644 --- a/src/common/gpu/shaders/gradientSubtract.wgsl +++ b/src/common/gpu/shaders/gradientSubtract.wgsl @@ -13,12 +13,13 @@ @group(0) @binding(1) var velocityTex : texture_2d; @group(0) @binding(2) var pressureTex : texture_2d; @group(0) @binding(3) var outTex : texture_storage_2d; +@group(0) @binding(6) var obstacleTex : texture_2d; fn pressureAt(neighbour: vec2, centre: f32) -> f32 { if (neighbour.x >= i32(u.gridSize.x)) { return 0.0; } - if (isSolidCell(neighbour, u)) { + if (isSolidAt(obstacleTex, neighbour, u)) { return centre; } let clamped = clampCell(neighbour, u); @@ -37,7 +38,7 @@ fn main(@builtin(global_invocation_id) id: vec3) { let c = vec2(cell); - if (isSolidCell(c, u)) { + if (isSolidAt(obstacleTex, c, u)) { textureStore(outTex, cell, vec4(0.0, 0.0, 0.0, 1.0)); return; } diff --git a/src/common/gpu/shaders/mask.wgsl b/src/common/gpu/shaders/mask.wgsl new file mode 100644 index 0000000..0fc617b --- /dev/null +++ b/src/common/gpu/shaders/mask.wgsl @@ -0,0 +1,29 @@ +// mask.wgsl — bakes the obstacle's signed distance into a grid-sized texture. +// +// Everything downstream needs to know which cells are solid, and several of the +// stencils need it for four neighbours as well as the centre. Evaluating the +// analytic SDF each time costs transcendentals — for the airfoil, several — and +// it is the same answer every time until the learner moves or resizes the body. +// +// So this kernel runs once whenever the obstacle (or the grid) changes, and the +// rest of the solver reads the result with a single texture fetch. See the +// helpers at the bottom of common.wgsl, and the dirty-tracking in +// WebGPUFluidEngine.recordCompute. +// +// The stored value is the distance in metres, positive outside the body, not a +// 0/1 flag: it is no more expensive to keep, and it leaves the door open to +// weighting the boundary by how much of a cell the body actually covers. + +@group(0) @binding(0) var u : SimUniforms; +@group(0) @binding(1) var outTex : texture_storage_2d; + +@compute @workgroup_size(8, 8) +fn main(@builtin(global_invocation_id) id: vec3) { + let cell = id.xy; + if (!inGrid(cell, u)) { + return; + } + + let distance = obstacleSDF(uvToMetres(cellToUV(cell, u), u), u); + textureStore(outTex, cell, vec4(distance, 0.0, 0.0, 1.0)); +} diff --git a/src/common/gpu/shaders/pressure.wgsl b/src/common/gpu/shaders/pressure.wgsl index 37b6028..bb114a3 100644 --- a/src/common/gpu/shaders/pressure.wgsl +++ b/src/common/gpu/shaders/pressure.wgsl @@ -32,15 +32,28 @@ // instead of piling up against a closed boundary. // Over-relaxation factor ω ∈ (1, 2). Above the Gauss–Seidel value of 1, larger -// ω converges faster but approaches the unstable edge at 2. 1.7 is conservative -// across every grid resolution the Lab screen offers and stays well clear of the -// instability that the collocated-grid checkerboard mode could otherwise amplify. +// ω converges faster but approaches the unstable edge at 2. +// +// It is worth knowing why this is not the textbook value. For the 5-point +// Laplacian on a W × H grid, Young's formula puts the optimal factor at 1.962 +// on the standard grid, rising toward 2 as the grid is refined — and measured +// against a converged reference, that value *does* win, but only in the regime +// the formula assumes: a fixed right-hand side iterated many times over. +// +// This solve is not in that regime. It gets a warm start, a fixed budget of +// thirty sweeps, and a right-hand side that has moved by the next frame, so it +// spends its life smoothing new error rather than converging old error — and +// over-relaxation near 2 is a poor smoother. Simulated over the sweep budgets +// and grid sizes this sim actually offers, the residual left at the end of a +// frame is flat between 1.7 and 1.8 and climbs steeply above 1.85. So the +// conservative value stands, on measurement rather than on caution. const PRESSURE_SOR_OMEGA: f32 = 1.7; @group(0) @binding(0) var u : SimUniforms; @group(0) @binding(1) var pressureTex : texture_2d; @group(0) @binding(2) var divergenceTex : texture_2d; @group(0) @binding(3) var outTex : texture_storage_2d; +@group(0) @binding(6) var obstacleTex : texture_2d; // Pressure at a neighbour, applying the boundary conditions. fn pressureAt(neighbour: vec2, centre: f32) -> f32 { @@ -50,7 +63,7 @@ fn pressureAt(neighbour: vec2, centre: f32) -> f32 { } // Solid neighbour (obstacle, or a cell past the wall after clamping): // reflect the centre value, giving zero normal gradient. - if (isSolidCell(neighbour, u)) { + if (isSolidAt(obstacleTex, neighbour, u)) { return centre; } let clamped = clampCell(neighbour, u); @@ -72,7 +85,7 @@ fn solveColour(id: vec3, redPass: bool) { // Pressure is undefined inside a solid; hold it at zero so the reflected // boundary values above stay well behaved. - if (isSolidCell(c, u)) { + if (isSolidAt(obstacleTex, c, u)) { textureStore(outTex, cell, vec4(0.0, 0.0, 0.0, 1.0)); return; } diff --git a/src/common/gpu/shaders/vorticity.wgsl b/src/common/gpu/shaders/vorticity.wgsl index 95c723b..44975ce 100644 --- a/src/common/gpu/shaders/vorticity.wgsl +++ b/src/common/gpu/shaders/vorticity.wgsl @@ -17,6 +17,7 @@ @group(0) @binding(1) var velocityTex : texture_2d; @group(0) @binding(2) var curlTex : texture_2d; @group(0) @binding(3) var outTex : texture_storage_2d; +@group(0) @binding(6) var obstacleTex : texture_2d; fn curlAt(cell: vec2) -> f32 { return textureLoad(curlTex, clampCell(cell, u), 0).x; @@ -31,7 +32,7 @@ fn main(@builtin(global_invocation_id) id: vec3) { let velocity = textureLoad(velocityTex, vec2(cell), 0).xy; - if (isSolidCell(vec2(cell), u)) { + if (isSolidAt(obstacleTex, vec2(cell), u)) { textureStore(outTex, cell, vec4(0.0, 0.0, 0.0, 1.0)); return; } diff --git a/src/common/gpu/solverSchedule.ts b/src/common/gpu/solverSchedule.ts new file mode 100644 index 0000000..e7e3d57 --- /dev/null +++ b/src/common/gpu/solverSchedule.ts @@ -0,0 +1,97 @@ +/** + * solverSchedule.ts + * + * How hard the iterative solves have to work this step. + * + * Both the viscous diffusion solve and the pressure solve are stationary + * iterations, and both have a well-known optimal over-relaxation factor. Rather + * than pinning a sweep count and a relaxation factor to whatever looked right at + * one grid resolution, these functions derive them from the problem the solver + * has actually been handed. + * + * ── Why this matters for the diffusion solve ────────────────────────────────── + * The implicit viscous solve is (I − νΔt∇²)u = u₀, whose stiffness is the single + * number α = νΔt/h². A plain Jacobi sweep reduces the error by at most + * ρ_J = 4α/(1 + 4α) per sweep, and α grows with the *square* of the resolution: + * at ν = 10⁻³ m²/s and Δt = 1/60 s it is 0.27 on the 256 × 128 grid but 17 on the + * 2048 × 1024 one. Twelve Jacobi sweeps leave 0.8 % of the error at the first and + * 84 % of it at the second — so at the fine resolutions the fluid was quietly + * less viscous than the learner asked for, in the direction that inflates the + * effective Reynolds number. + * + * Red-black SOR with the optimal factor for that α converges at ω − 1 per sweep + * instead, which is 0.08 and 0.71 for the same two cases. The sweep count then + * follows from the accuracy actually wanted, so the cheap end of the range gets + * cheaper and the stiff end gets correct. + * + * No GPU or scenery dependency: this is arithmetic, and it is unit-tested. + */ + +import { + DIFFUSION_RESIDUAL_TOLERANCE, + DIFFUSION_SKIP_ALPHA, + DIFFUSION_SWEEPS_MAX, +} from "../../FluidDynamicsConstants.js"; + +/** + * α = νΔt/h², the stiffness of the implicit viscous solve. + * + * Zero when the sim is paused (Δt = 0), which is the case the skip guard below + * exists for. + */ +export function diffusionAlpha(viscosity: number, dt: number, cellSize: number): number { + return (viscosity * dt) / (cellSize * cellSize); +} + +/** + * Optimal SOR factor for the diffusion system, for a given α. + * + * The Jacobi iteration matrix of (1 + 4α)x_c − αΣx_n = f has spectral radius + * ρ_J = 4α/(1 + 4α), and Young's formula gives ω = 2/(1 + √(1 − ρ_J²)). Both + * radicals simplify: 1 − ρ_J² = (1 + 8α)/(1 + 4α)², so + * + * ω = 2(1 + 4α) / (1 + 4α + √(1 + 8α)) + * + * which needs one square root and no trigonometry. **The same closed form is + * evaluated in shaders/diffuse.wgsl** — the two must agree, because the sweep + * count below is chosen for the rate this ω produces. + * + * ω → 1 as α → 0 (the system becomes the identity and Gauss–Seidel is already + * exact) and → 2 as α → ∞, staying inside the convergent range at every α. + */ +export function diffusionOmega(alpha: number): number { + const diagonal = 1 + 4 * alpha; + return (2 * diagonal) / (diagonal + Math.sqrt(1 + 8 * alpha)); +} + +/** + * Error remaining after one red-black SOR iteration (one red sweep plus one + * black sweep) of the diffusion solve. This is ω − 1, the asymptotic rate of SOR + * at its optimal factor. + */ +export function diffusionConvergenceRate(alpha: number): number { + return diffusionOmega(alpha) - 1; +} + +/** + * Red-black SOR iterations to run on the viscous solve this step — each one is + * a red dispatch and a black dispatch. + * + * Chosen as the fewest that drive the error below DIFFUSION_RESIDUAL_TOLERANCE, + * capped so that the stiffest corner of the parameter space (the finest grid at + * the top of the viscosity slider) cannot turn into an unbounded dispatch count. + * + * Returns 0 below DIFFUSION_SKIP_ALPHA, where every sweep is the identity to + * within float precision and only the seeding dispatch is needed. + */ +export function diffusionSweeps(alpha: number): number { + if (!(alpha >= DIFFUSION_SKIP_ALPHA)) { + return 0; + } + const rate = diffusionConvergenceRate(alpha); + if (rate <= 0) { + return 1; + } + const exact = Math.log(DIFFUSION_RESIDUAL_TOLERANCE) / Math.log(rate); + return Math.min(DIFFUSION_SWEEPS_MAX, Math.max(1, Math.ceil(exact))); +} diff --git a/tests/ShaderBindings.test.ts b/tests/ShaderBindings.test.ts new file mode 100644 index 0000000..c68b479 --- /dev/null +++ b/tests/ShaderBindings.test.ts @@ -0,0 +1,177 @@ +/** + * Tests for the WGSL/bind-group-layout contract. + * + * A compute kernel's resources are checked against its pipeline layout by + * WebGPU, but only at pipeline-creation time and only on a real device — which + * in this sim means "the whole field turns into the WebGPU-unavailable message, + * on hardware, with a validation error in the console". Since several kernels + * share a layout, adding a binding to one of them and forgetting the layout is + * an easy mistake with a very indirect symptom. + * + * So the layouts are declared as plain data in bindLayouts.ts and these tests + * parse the shader sources to check that every binding a kernel declares exists + * in the layout its pipelines are built with, with a compatible kind and, for + * storage textures, the same format. No GPU required. + * + * The check is one-directional — shader ⊆ layout — because a shared layout may + * legitimately carry entries a given kernel has no use for. The obstacle field + * is bound to every compute layout, but the dye injection kernel never asks + * where the body is. + */ + +import { describe, expect, it } from "vitest"; +import { + BIND_LAYOUT_NAMES, + BIND_LAYOUTS, + type BindingSpec, + layoutBinding, + OBSTACLE_BINDING, + SCALAR_FORMAT, + SHADER_LAYOUTS, +} from "../src/common/gpu/bindLayouts.js"; + +/** + * The preamble, which declares no bindings of its own and is concatenated ahead + * of all the others rather than compiled as a shader. + */ +const PREAMBLE = "common.wgsl"; + +/** Every shader in the folder, keyed by file name, so a new one cannot hide. */ +const SHADER_SOURCES = Object.fromEntries( + Object.entries( + import.meta.glob("../src/common/gpu/shaders/*.wgsl", { query: "?raw", import: "default", eager: true }), + ) + .map(([path, source]) => [path.split("/").pop() ?? path, source as string] as const) + .filter(([file]) => file !== PREAMBLE), +); + +type Declaration = { + readonly binding: number; + readonly group: number; + readonly name: string; + readonly spec: BindingSpec; +}; + +/** + * Parses the `@group(g) @binding(b) var name : type;` declarations out of + * a shader, mapping each to the layout entry it requires. + * + * WGSL cannot express the difference between a filterable and an unfilterable + * texture — both are `texture_2d` — so a plain texture is only checked to + * be a texture. Storage textures carry their format in the type and are checked + * exactly, which is the half that actually drifts. + */ +function declarations(source: string): Declaration[] { + const pattern = /@group\((\d+)\)\s*@binding\((\d+)\)\s*var(?:<([^>]*)>)?\s+(\w+)\s*:\s*([^;]+);/g; + const found: Declaration[] = []; + + for (const match of source.matchAll(pattern)) { + const [, group, binding, addressSpace, name, rawType] = match; + const type = (rawType ?? "").trim(); + + let spec: BindingSpec; + if ((addressSpace ?? "").trim() === "uniform") { + spec = { kind: "uniform" }; + } else if (type === "sampler") { + spec = { kind: "sampler" }; + } else { + const storage = type.match(/^texture_storage_2d<\s*(\w+)\s*,/); + if (storage) { + spec = { kind: "storageTexture", format: storage[1] as GPUTextureFormat }; + } else { + expect(type, `unrecognised binding type for ${name}`).toMatch(/^texture_2d { + it("covers every shader file with a layout", () => { + expect(Object.keys(SHADER_SOURCES).sort()).toEqual(Object.keys(SHADER_LAYOUTS).sort()); + }); + + it("finds the declarations at all — every kernel binds at least a uniform and an output", () => { + // Without this, a regex that stopped matching would make every check below + // pass by looping over nothing. + for (const [file, source] of Object.entries(SHADER_SOURCES)) { + const found = declarations(source); + expect(found.length, `${file} parsed to no bindings`).toBeGreaterThanOrEqual(2); + expect( + found.map((declaration) => declaration.spec.kind), + file, + ).toContain("uniform"); + } + }); + + it("declares every resource in group 0, which is the only group any pipeline binds", () => { + for (const [file, source] of Object.entries(SHADER_SOURCES)) { + for (const declaration of declarations(source)) { + expect(declaration.group, `${file}: ${declaration.name}`).toBe(0); + } + } + }); + + it("declares no binding index twice within one shader", () => { + for (const [file, source] of Object.entries(SHADER_SOURCES)) { + const indices = declarations(source).map((declaration) => declaration.binding); + expect(new Set(indices).size, `${file} reuses a binding index`).toBe(indices.length); + } + }); + + it("matches every declared binding to an entry of the same kind in its layout", () => { + for (const [file, layoutName] of Object.entries(SHADER_LAYOUTS)) { + const source = SHADER_SOURCES[file]; + expect(source, `no source for ${file}`).toBeDefined(); + + for (const declaration of declarations(source ?? "")) { + const entry = layoutBinding(layoutName, declaration.binding); + expect( + entry, + `${file}: binding ${declaration.binding} (${declaration.name}) is not in layout ${layoutName}`, + ).toBeDefined(); + if (entry === undefined) { + continue; + } + + expect(entry.kind, `${file}: ${declaration.name} is a ${declaration.spec.kind}`).toBe(declaration.spec.kind); + + if (declaration.spec.kind === "storageTexture" && entry.kind === "storageTexture") { + expect(entry.format, `${file}: ${declaration.name} writes a different format than its layout declares`).toBe( + declaration.spec.format, + ); + } + } + } + }); + + it("puts the obstacle field at the same binding in every kernel that reads it", () => { + for (const [file, source] of Object.entries(SHADER_SOURCES)) { + const obstacle = declarations(source).find((declaration) => declaration.name === "obstacleTex"); + if (obstacle !== undefined) { + expect(obstacle.binding, `${file} binds the obstacle field somewhere else`).toBe(OBSTACLE_BINDING); + expect(obstacle.spec.kind).toBe("texture"); + } + } + }); + + it("offers the obstacle field to every compute layout except the one that writes it", () => { + for (const name of BIND_LAYOUT_NAMES) { + if (name === "mask" || BIND_LAYOUTS[name].stage !== "compute") { + continue; + } + expect(layoutBinding(name, OBSTACLE_BINDING), `${name} cannot see the obstacle`).toEqual({ + kind: "texture", + sampleType: "unfilterable-float", + }); + } + }); + + it("writes the obstacle field as the scalar format the readers are declared with", () => { + expect(BIND_LAYOUTS.mask.bindings[1]).toEqual({ kind: "storageTexture", format: SCALAR_FORMAT }); + }); +}); diff --git a/tests/fuzz/engine.spec.ts b/tests/fuzz/engine.spec.ts index 9ed863e..cd32f39 100644 --- a/tests/fuzz/engine.spec.ts +++ b/tests/fuzz/engine.spec.ts @@ -101,8 +101,9 @@ async function render(page: Page, format: string, steps: number, overrides: Reco test.describe("WebGPU fluid engine", () => { // The solver runs hundreds of steps per case; the default 30 s is not enough - // on a software rasterizer. The MacCormack velocity corrector adds a dispatch - // per frame, so this is budgeted generously rather than tuned to the wire. + // on a software rasterizer. The per-frame dispatch count now varies with the + // viscosity (the viscous solve schedules its own sweeps), so this is budgeted + // for the stiff end of the range rather than tuned to the wire. test.setTimeout(360_000); test("dye is carried downstream and around a cylinder", async ({ page }) => { diff --git a/tests/solverSchedule.test.ts b/tests/solverSchedule.test.ts new file mode 100644 index 0000000..370950e --- /dev/null +++ b/tests/solverSchedule.test.ts @@ -0,0 +1,166 @@ +/** + * Tests for how hard the iterative solves are made to work. + * + * The viscous solve's sweep count and over-relaxation factor are derived from + * α = νΔt/h² rather than pinned, which is what stops the fine grids from being + * quietly less viscous than the Reynolds-number readout claims. Getting that + * arithmetic wrong is invisible — the sim still runs, it is just solving a + * different problem than the one on the label — so it is pinned here. + * + * The `diffusionOmega` formula is also evaluated, in the same closed form, by + * shaders/diffuse.wgsl. These tests check the CPU half against the textbook + * definition it is a simplification of; the WGSL half is a transcription of the + * same three lines. + */ + +import { describe, expect, it } from "vitest"; +import { FluidGridSpec } from "../src/common/gpu/FluidGridSpec.js"; +import { + diffusionAlpha, + diffusionConvergenceRate, + diffusionOmega, + diffusionSweeps, +} from "../src/common/gpu/solverSchedule.js"; +import { + DIFFUSION_RESIDUAL_TOLERANCE, + DIFFUSION_SKIP_ALPHA, + DIFFUSION_SWEEPS_MAX, + MAX_PHYSICS_DT, + VISCOSITY_DEFAULT, + VISCOSITY_RANGE, +} from "../src/FluidDynamicsConstants.js"; + +/** Young's formula, written out, for the Jacobi radius of the viscous system. */ +function textbookOmega(alpha: number): number { + const jacobiRadius = (4 * alpha) / (1 + 4 * alpha); + return 2 / (1 + Math.sqrt(1 - jacobiRadius * jacobiRadius)); +} + +const ALPHAS = [0, 1e-6, 1e-3, 0.01, 0.08, 0.273, 1, 5, 17.5, 27.3, 200, 1747]; + +describe("diffusionAlpha", () => { + it("is νΔt/h²", () => { + expect(diffusionAlpha(1e-3, 1 / 60, 1 / 128)).toBeCloseTo((1e-3 * (1 / 60)) / (1 / 128) ** 2, 9); + }); + + it("grows with the square of the resolution, which is why a fixed sweep count could not work", () => { + const standard = diffusionAlpha( + VISCOSITY_DEFAULT, + MAX_PHYSICS_DT, + FluidGridSpec.forResolution("standard").cellSize, + ); + const ultra = diffusionAlpha(VISCOSITY_DEFAULT, MAX_PHYSICS_DT, FluidGridSpec.forResolution("ultraFine").cellSize); + + // Eight times the linear resolution, sixty-four times the stiffness. + expect(ultra / standard).toBeCloseTo(64, 6); + }); + + it("is zero when the sim is paused", () => { + expect(diffusionAlpha(VISCOSITY_DEFAULT, 0, 1 / 128)).toBe(0); + }); +}); + +describe("diffusionOmega", () => { + it("agrees with Young's formula at every stiffness", () => { + for (const alpha of ALPHAS) { + expect(diffusionOmega(alpha), `α = ${alpha}`).toBeCloseTo(textbookOmega(alpha), 9); + } + }); + + it("stays inside the convergent range (1, 2)", () => { + for (const alpha of ALPHAS) { + expect(diffusionOmega(alpha), `α = ${alpha}`).toBeGreaterThanOrEqual(1); + expect(diffusionOmega(alpha), `α = ${alpha}`).toBeLessThan(2); + } + }); + + it("collapses to Gauss–Seidel when the system is the identity", () => { + expect(diffusionOmega(0)).toBe(1); + }); + + it("rises with stiffness", () => { + const sorted = [...ALPHAS].sort((a, b) => a - b); + for (let i = 1; i < sorted.length; i++) { + expect(diffusionOmega(sorted[i] ?? 0)).toBeGreaterThanOrEqual(diffusionOmega(sorted[i - 1] ?? 0)); + } + }); +}); + +describe("diffusionConvergenceRate", () => { + it("beats the Jacobi rate it replaced, by more the stiffer the system is", () => { + for (const alpha of [0.273, 1, 17.5, 27.3]) { + const jacobi = (4 * alpha) / (1 + 4 * alpha); + expect(diffusionConvergenceRate(alpha), `α = ${alpha}`).toBeLessThan(jacobi); + } + + // The case that motivated the change: the default viscosity on the finest + // grid, where twelve Jacobi sweeps removed almost none of the error. + const alpha = diffusionAlpha(VISCOSITY_DEFAULT, MAX_PHYSICS_DT, FluidGridSpec.forResolution("ultraFine").cellSize); + const jacobiAfter12 = ((4 * alpha) / (1 + 4 * alpha)) ** 12; + const sorAfter12 = diffusionConvergenceRate(alpha) ** 12; + + expect(jacobiAfter12).toBeGreaterThan(0.5); + expect(sorAfter12).toBeLessThan(0.05); + }); +}); + +describe("diffusionSweeps", () => { + it("skips the solve entirely when every sweep would be the identity", () => { + expect(diffusionSweeps(0)).toBe(0); + expect(diffusionSweeps(DIFFUSION_SKIP_ALPHA / 2)).toBe(0); + expect(diffusionSweeps(Number.NaN)).toBe(0); + }); + + it("runs at least one sweep as soon as the solve does anything", () => { + expect(diffusionSweeps(DIFFUSION_SKIP_ALPHA)).toBeGreaterThanOrEqual(1); + }); + + it("never exceeds the cap, however stiff the system gets", () => { + for (const alpha of ALPHAS) { + expect(diffusionSweeps(alpha), `α = ${alpha}`).toBeLessThanOrEqual(DIFFUSION_SWEEPS_MAX); + } + expect(diffusionSweeps(1e9)).toBe(DIFFUSION_SWEEPS_MAX); + }); + + it("asks for enough sweeps to reach the tolerance whenever the cap allows it", () => { + for (const alpha of ALPHAS) { + const sweeps = diffusionSweeps(alpha); + if (sweeps === 0 || sweeps === DIFFUSION_SWEEPS_MAX) { + continue; + } + expect(diffusionConvergenceRate(alpha) ** sweeps, `α = ${alpha}`).toBeLessThanOrEqual( + DIFFUSION_RESIDUAL_TOLERANCE, + ); + // …and not one sweep more than it needs. + expect(diffusionConvergenceRate(alpha) ** (sweeps - 1), `α = ${alpha} is not overspent`).toBeGreaterThan( + DIFFUSION_RESIDUAL_TOLERANCE, + ); + } + }); + + it("asks for more sweeps as the system stiffens", () => { + const sorted = [...ALPHAS].sort((a, b) => a - b); + for (let i = 1; i < sorted.length; i++) { + expect(diffusionSweeps(sorted[i] ?? 0)).toBeGreaterThanOrEqual(diffusionSweeps(sorted[i - 1] ?? 0)); + } + }); + + it("costs less than the twelve fixed sweeps it replaced at the settings both screens open on", () => { + const cellSize = FluidGridSpec.forResolution("standard").cellSize; + const sweeps = diffusionSweeps(diffusionAlpha(VISCOSITY_DEFAULT, MAX_PHYSICS_DT, cellSize)); + + // Each sweep is a red dispatch and a black one, plus the seeding dispatch. + expect(2 * sweeps + 1).toBeLessThan(12); + }); + + it("stays bounded across the whole parameter space the sliders can reach", () => { + for (const resolution of ["standard", "fine", "veryFine", "ultraFine"] as const) { + const cellSize = FluidGridSpec.forResolution(resolution).cellSize; + for (const viscosity of [VISCOSITY_RANGE.min, VISCOSITY_DEFAULT, VISCOSITY_RANGE.max]) { + const sweeps = diffusionSweeps(diffusionAlpha(viscosity, MAX_PHYSICS_DT, cellSize)); + expect(sweeps, `${resolution} at ν = ${viscosity}`).toBeGreaterThanOrEqual(1); + expect(sweeps, `${resolution} at ν = ${viscosity}`).toBeLessThanOrEqual(DIFFUSION_SWEEPS_MAX); + } + } + }); +});