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); + } + } + }); +});