OpenPhysics · veillette · Aug 8, 2026 · Aug 8, 2026
@@ -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 |

@@ -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 |

@@ -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

@@ -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,