# Velocity diffusion This folder solves the implicit viscous step of every time step: for each MAC velocity component, ```text (I − ν·dt·∇²) u_new = u_old ``` This README is the implemented contract for the operator, the solve method, its exits, the settings, the resource budget, the kernel rules and the tests ([contract homes](https://github.com/hankbeasley/polycfd/blob/b7f5b9ba6f52355d9b196b586e4a6c8719d8feea/docs/w4-design/README.md#where-implemented-contracts-live)). The [physics guide](/docs/PHYSICS_GUIDE#diffusion-and-walls) and [boundary guide](/docs/BOUNDARY_CONDITIONS_HALOS#8-bc-aware-diffusion-halo-free-stencil-computation) derive the ghost relations; [diffusion convergence](/docs/ADAPTIVE_DIFFUSION_RELATIVE_TOLERANCE) keeps the history of the fallback, budget and primary-method decisions. **Status (2026-09-26, SOL-02a):** primary method Chebyshev semi-iteration on undamped Jacobi, then colored Gauss-Seidel, with deterministic double audits. It replaced weighted Jacobi (ω = 0.8) behind the same seam ([solver seams](/docs/RUNTIME_CONTRACTS#extension-seams)); the quick, multigrid and full `verify-all` suites pass with it ([complete validation](/docs/COMPLETE_VALIDATION#chebyshev-solvers-full-suite--2026-09-26)). Whether the prescribed-flow enclosure campaign moves to this method is the user's decision ([enclosure plan](https://github.com/hankbeasley/polycfd/blob/b7f5b9ba6f52355d9b196b586e4a6c8719d8feea/docs/plans/prescribed-flow-enclosure.md)). | File | Role | | --- | --- | | `DiffusionGpu.cs` | Row evaluator, kernels, solve phases, result assembly, scratch bytes | | `DiffusionBcDescriptor.cs` | `FaceBc` and `VelocityDiffusionBcs`: per-face boundary type and value, normal-face index lists | | `ImmersedWallCoefficientsGpu.cs` | Device copy of the static `ImmersedWallCoefficients` (mask, diagonal, right-hand side) | | `src/PolyCfd.Core/Numerics/DiffusionSolveSettings.cs` | The recorded budget and tolerance | | `src/PolyCfd.Core/Numerics/DiffusionResult.cs` | The typed result | | `src/PolyCfd.Core/Numerics/Chebyshev.cs`, `ChebyshevStallMonitor.cs` | Host double semi-iteration weights and the stall rule (no device code) | ## Operator One row per stored unknown of a component, with c_e = α/h_e² and α = ν·dt, both float: ```text diag·u − Σ c_e·u_e = b, diag = 1 + 2α(1/hx² + 1/hy² + 1/hz²) ``` `EvaluateRow` is the only encoding of the row and its folds; every kernel calls it: - **Tangential Dirichlet face** (no-slip, fixed value, moving wall): ghost 2·u_wall − u, so the neighbour sum gains 2·u_wall·c and the diagonal gains c. - **Neumann face** (zero gradient, slip, symmetry): ghost u, so the diagonal loses c. - **Periodic face:** the wrapped neighbour. On the component's own (staggered) axis both end faces store the same physical face, so the wrap skips the duplicate: the low neighbour of face 0 is face N − 1 and the high neighbour of face N is face 1. - **Prescribed normal face:** a normal face whose patch holds normal velocity (the Dirichlet types, slip, symmetry) keeps its patch value throughout the solve. A partial velocity mask holds only its selected faces; unselected faces are ordinary rows. - **Immersed walls** (`WallCoefficients`, static geometry only): a held face keeps its wall velocity; every wall link drops its neighbour and adds the coefficient's diagonal and right-hand side, times α. The **effective right-hand side** b_eff is b plus the Dirichlet forcing and the immersed-wall right-hand side. The relative residual that every exit uses is ```text ‖b_eff − A u‖₂ / (‖b_eff‖₂ + 1e-12) ``` per component, where a prescribed row contributes p − u and a held immersed face 0. In the norm of b_eff a prescribed row contributes p and a held face its b. Rows and residuals accumulate in double; fields and coefficients stay float, and the correction is scaled in float. `ComputeResidual` (internal) writes the signed r = b_eff − A u of one component through the production row with no boundary application. It is the operator-test seam and a building block for a future diffusion multigrid. ## Method, exits and result `DiffuseVelocity(state, boundaries, ν, dt, settings)` runs, in order: 1. **Zero viscosity:** nothing to solve; the result has 0 iterations, converged, exit `Skipped`. 2. **Initial pass:** copies b = u_old (also the initial guess) and, in one fused launch per component, reduces the effective-RHS norm and the initial residual ‖b_eff − A b‖ (the audit of b). When all three components are already below the tolerance the solve returns b with 0 iterations and exit `AlreadyConverged`, reporting that audit. 3. **`RunPrimary`:** Chebyshev semi-iteration on undamped Jacobi, within the primary limit: the private constant `PrimaryIterationLimit = 400`, or `MaxIterations` when smaller. With J(x) the undamped Jacobi update of the row, ```text x₁ = J(x₀), x_{k+1} = x_{k−1} + w_{k+1}·(J(x_k) − x_{k−1}), w₁ = 1, w₂ = 2/(2 − ρ²), w_{k+1} = 1/(1 − ρ²w_k/4) ``` with the weights in host double (`Chebyshev.SemiIterationWeight`) passed as a float per step. The first step writes x₁ into the iterate scratch and never reads it; later steps overwrite x_{k−1} in place in the existing ping-pong pair, so there are no new buffers. Prescribed normal faces follow the recurrence with J = p (exactly p at step 1 and, in exact arithmetic, at every odd degree); held wall faces keep their value. Every fourth step (4, 8, ...) also audits its **input** x_k, whose degree is odd. Exits: | Exit | Rule | Then | | --- | --- | --- | | `Converged` | an audit below the tolerance for all three components | returns the **audited** x_k (that step's output is discarded), so the reported residual describes the stored field | | `NonFinite` | an audit that is not finite | restores b for all three components, then the fallback | | `Stalled` | `ChebyshevStallMonitor` (below) | keeps, per component, the better of x_k and b (the initial pass audited b), then the fallback | | `Limit` | the primary limit reached | the fallback continues from x_limit | A primary limit of 0 (the phase-coverage overload) or a bound that is not below 1 in float (S ≳ 2^24) gives exit `Skipped`, and the fallback runs alone. 4. **`RunFallback`** (unconverged primary only): colored Gauss-Seidel passes with the remaining budget, on the same row. Two colors suffice except for an odd periodic cell count, which uses three (each odd periodic axis alternates 0/1 and gives its last cell 2; duplicated staggered endpoints share their physical color, except the two directly linked endpoints of a single periodic cell). Each color is a separate ordered launch, and one full cycle counts as one pass. At every total iteration count that is a multiple of four, and at the cap, the velocity boundary values are applied and the **stored** field is audited against the original b_eff. When no budget remains, the stored field is audited once, so a short or non-multiple-of-four cap reports its final state. 5. **Final `ApplyVelocityBcs`**, then one result assembly. It stays load-bearing: after a `Limit` exit the even-degree iterate does not hold prescribed faces exactly until the fallback or this call sets them. Holding prescribed rows at p every sweep would make it a provable no-op, after which `DiffusionGpu` would no longer need `IBoundaryKernels` (the proposed SOL-05, with the bound-derived exits). **Spectral radius bound** (`SpectralRadiusBound`, from the float coefficients the kernels use): with c_e = α/h_e² and S = 2(c_x + c_y + c_z), ρ = S/(1 + S); isotropic, 6d/(1 + 6d) with d = ν·dt/h² (0.9545 on the 128³ fan, d ≈ 3.5). Proof: every link of an active row has off-diagonal weight o at most its diagonal share g and at most c (ordinary link g = o = c; Dirichlet face g = 2c, o = 0; Neumann face g = o = 0; immersed-wall link g = c·h/d ≥ c/2, o = 0), so the row sum of |D⁻¹N| is Σo/(1 + Σg) ≤ S/(1 + S) < 1, attained by six ordinary links; prescribed rows (J = p) and held rows are zero rows. The spectrum is real because ordinary links are symmetric, so D⁻¹N is similar to D^(−1/2)·N·D^(−1/2) on the active unknowns; links into prescribed faces, duplicated staggered endpoints and links out of solid-centred active unknowns are one-way and add only eigenvalue 0. The one nonsymmetric case, a periodic wrap into a solid-centred active unknown, is guarded by the exits and the fallback (limits). **Stall rule** (`ChebyshevStallMonitor`, the experiment's validated heuristic): with cosh θ = 1/ρ, it starts once the audited degree k satisfies kθ ≥ 2 and after two audits in that regime; a component still at or above the tolerance whose residual exceeds √(T_{k−8}/T_k)·max(r_{k−4}, r_{k−8}) at two consecutive audits has stalled (the float storage floor, or a spectrum the bound does not cover). The proposed SOL-05 replaces it with exits derived from the bound, a change to the monitor class alone. `DiffusionResult` reports `PrimaryIterations` (steps launched, including a discarded final step), `FallbackIterations`, `Iterations` (their sum; residual-only audits are not counted), `PrimaryExit` (`Converged`, `Limit`, `Skipped`, `AlreadyConverged`, `Stalled`, `NonFinite`) and per component `Converged`, `Iterations` and `RelativeResidual`. All components share the iteration count and the convergence decision. The validation case table's `D_Iters` column is `Iterations`. A result above the tolerance is a failure that step acceptance refuses; no convergence is guaranteed for every grid and time step. ## Settings `DiffusionSolveSettings` (Core) holds exactly the recorded fields: `MaxIterations` (800, primary steps plus fallback passes) and `RelativeTolerance` (1e-6). `DiffusionSolveSettings.Default` is used by both integrator call sites, and `ScenarioResolver` records it as `DiffusionMaxIterations`/`DiffusionRelativeTolerance` in `ScenarioEffectiveSettings`. The settings argument is required, so there is one way to call the solver. The split between primary steps and fallback passes and the check interval are the method's own policy, pinned by the build hashes in `scientificId`, not settings. Tests reach an internal overload with an explicit primary limit (0 runs the fallback alone). `DiffusionGpu` must not know scenario settings, stage order, step acceptance, estimator categories, files or the console. ## Resources - **Device bytes:** `DiffusionGpu.ScratchBytes(faces, maxFace) = 8·faces + 4·maxFace`: three right-hand-side buffers and one work buffer holding a largest-component partials region followed by the three iterate buffers, allocated once per grid size and retained. A grid whose work buffer would exceed one allocation (int.MaxValue floats, more than about 812³ cells) keeps the region alone, max(MaxFace, 4G) floats, and gives each iterate its own buffer: the same bytes, since 4G ≤ MaxFace under CUDA groups on any grid that large. `ScenarioMemoryEstimator` reports it as `diffusion-scratch`, and `DiffusionResourceTests` requires the first solve to retain exactly these bytes in either layout. The fallback's Gauss-Seidel and audit kernels and the residual kernel compile on first use and allocate nothing. The Chebyshev recurrence adds no buffer. - **Audit partials.** The step, audit and initial-pass kernels reduce squared residuals per launch group in double (a fixed tree over the actual group), and each audit reads all three components' partials in **one** copy, summed on the host in index order: G = Σ_c groups(c) doubles for the groups `KernelLaunchHelpers.Make3D` launches (2G for the initial pass: norm and residual), read through a reinterpreted slice at offset 0 of the work buffer. The partials region is max(MaxFace, 2G, 4G − faces) floats; the initial pass may extend into the iterate buffers, which are free then. Under CUDA groups (32×4×2) and at least two cells per axis (the validator's minimum) each component needs at most faces_c/6 groups, so 2G ≤ MaxFace and 4G < faces: the region is exactly MaxFace and the CUDA bytes above are exact (128³: G = 27,008, 211 KiB). The ILGPU CPU accelerator's smaller groups (about 16×1×1) can enlarge the region on thin grids, an assumption the estimator states. `EnsureScratch` refuses a group above 256 threads or a partials layout that does not fit, and checks the 8-byte alignment of the double view off CUDA; it never truncates. - **Host reads:** one for the initial pass and one per audit: every fourth iteration, at a cap that is not a multiple of four, and once more on the stored field when the primary phase ends unconverged at the cap (with the cap a multiple of four, its last in-loop audit described the iterate before). A solve of I iterations therefore reads 1 + ⌈I/4⌉ times, plus one in that last case; at most 2 + ⌊I/4⌋. The stall exit reuses the initial pass's audit of b. Exit decisions are deterministic: the partial sums have a fixed order. - **After warm-up** a solve makes no application allocation and no upload, and exactly those device-to-host copies; the largest is the initial pass's 2G doubles, far smaller than a field under CUDA groups (`DiffusionResourceTests`). The counts cover what `DeviceMemoryAudit` sees: explicit application allocations and transfers. - **No field-sized download, no files, no console output.** ## Kernel rules Every new or rewritten diffusion kernel follows these rules: - **No nested integer conditionals.** ILGPU 1.5.3 in a C# Release CUDA build silently narrows a nested integer join with a boolean-derived branch to bool ([bug report](https://github.com/hankbeasley/polycfd/blob/b7f5b9ba6f52355d9b196b586e4a6c8719d8feea/docs/ILGPU_PERIODIC_INDEX_BUG_REPORT.md)). Every integer selection in `EvaluateRow`, the initial-pass, Chebyshev step, audit, Gauss-Seidel and residual kernels and the helpers they inline (`MacComponent.IndexAlong`/`Extent`, `AxisColor`, `EffectiveRhs`, `WriteGroupSum`, the immersed-wall mask load) is a separate single-level statement stored to a local; the periodic wrap keeps the report's workaround verbatim, the row kind is two booleans, not an integer, and the step's first-step and audit switches are bit tests of one flags argument. The report's workaround and controls pass only for the shapes it tested, so the Release CUDA operator tests below remain the evidence. - **Correct for any launch group.** Kernels take their extents from `MacComponent` and bounds-check every thread. The group reduction (`WriteGroupSum`) runs over the actual launch group of `KernelLaunchHelpers.Make3D` (the CPU accelerator uses about 16 threads per group), never the CUDA 32×4×2: shared memory holds the 256-thread maximum, the tree follows the linear thread index and the real group size, and out-of-range threads contribute 0 instead of returning early, so every thread reaches the barriers. - **CPU coverage is reported, not claimed.** Kernels must compile for the CPU accelerator; CPU test passes are not run (user directive), and results say "CPU backend not run (user directive)". - **Independent host checks run through the production step kernels,** at least one per boundary type, not only through `ComputeResidual`: a miscompiled row shared by the step and the audit could otherwise report convergence of the wrong operator. ## Limits - **Large diffusion numbers.** Chebyshev's step count grows like √(1 + 6d) instead of Jacobi's 1 + 6d, but it still grows: the moving-wall rows at α/h² = 8 converge in 80 steps at the test's relative tolerance 2e-7 (`DiffusionBoundaryForcingTests`), and d ≳ 230 exceeds the 400-step reach at 1e-6 and relies on Gauss-Seidel. A diffusion multigrid would replace `RunPrimary`. - **The stall rule can leave Chebyshev early at large d.** Circular Couette at d = 115 (relative tolerance 1e-7) stalls at step 96 of its first solve; Gauss-Seidel then converges most steps within the 800 budget, but one of the first eight ends unconverged at the cap (reported 1.008e-7, equal to the stored field's residual), and the steady state then returns `AlreadyConverged` (`ChebyshevDiffusionMethodTests`). Such a solve is reported honestly and step acceptance refuses it; convergence at large d is claimed only after the proposed bound-based exits (SOL-05). - **Near-unit fields hit the float storage floor:** at α/h² = 6.4 the stored-field residual stalls near 1.12–1.16e-6, above 1e-6; the stall rule hands the solve to Gauss-Seidel, which ends honestly unconverged at the cap. - **A periodic wrap into a solid-centred active unknown** is a one-way link that can make the iteration matrix nonsymmetric, outside the bound's real-spectrum argument. The stall and non-finite exits and the fallback guard it; there is no test of that geometry yet. - **Float coefficients.** The operator is defined by the float coefficients of `EvaluateRow`; on CUDA the immersed-wall folds `diagonal + α·Diag` and `b + α·Rhs` compile to single-rounding fused multiply-adds (a code-generation property; the CPU accelerator need not contract them). At Couette d = 115 the same host row with two roundings is 1.0–3.2e-8 relative higher, 10–32% of that test's 1e-7 tolerance, and would judge every converged step unconverged: convergence at that tolerance and diffusion number lies inside the float rounding of the operator's coefficients. Independent host checks near that level fold the same way as the device (`ChebyshevDiffusionMethodTests`, CUDA only). - Immersed-wall coefficients are built once for static geometry; the moving path does not rebuild them. - **Adaptive steps cap the diffusion number** at `CflGpu.MaxDiffusionCfl` = 3.5 (ν·dt/h² on the smallest spacing), which bounds the implicit solve's cost. Raising it changes every adaptive trajectory, a separate scientific decision; fixed steps are not capped (the enclosure runs d = 6.4). ## Testing ```bash dotnet test test/PolyCfd.Gpu.Tests -c Release --settings test/gpu.runsettings --filter "FullyQualifiedName~Diffusion&Backend!=CPU" dotnet test test/PolyCfd.Tests -c Release --filter "FullyQualifiedName~Chebyshev" # host Chebyshev math, no GPU ``` Tests follow the [test kinds](https://github.com/hankbeasley/polycfd/blob/b7f5b9ba6f52355d9b196b586e4a6c8719d8feea/test/README.md#test-kinds): untagged tests are contracts that any correct method passes; `Kind=MethodSpecific` tests encode the current method and are the only ones rewritten or deleted when it changes. `test/PolyCfd.Gpu.Tests/Diffusion/DiffusionReference.cs` holds the independent host double row (with the row diagonal), a Thomas solver, field fills and the boundary fake that leaves the audited field untouched. | Contract | Tests | | --- | --- | | Operator through `ComputeResidual` against the host row | `DiffusionOperatorTests` (tangential Dirichlet and Neumann ghosts; prescribed normal faces with a partial mask; effective right-hand side with domain and immersed forcing), `DiffusionPeriodicBcTests.PeriodicRowResidualMatchesStencilWithoutBoundaryCopy`, `ImmersedWallDiffusionTests.TiltedWall_LinearShearProfile_IsExactOnInteriorRows` (float roundoff on rows whose wall links are intersected; an O(h) bound on fallback-link rows) | | Independent references through the step kernels | Thomas solutions (`DiffusionBoundaryForcingTests`, `DiffusionNormalEndpointTests`), backward-Euler Fourier modes (`DiffusionPeriodicBcTests`), steady profiles (`DiffusionWallBcTests`, `DiffusionNeumannBcTests`, `ImmersedWallDiffusionTests`) | | Honest reporting | `DiffusionResidualAccuracyTests.HighAlphaSolveMeetsIndependentResidualTarget`: converged exactly when the independent residual is below 1e-6, never beyond the cap, the reported residual is the stored field's (converged or not, so a converged solve returns the audited field), and an unconverged solve ends at the cap | | Initial pass | `DiffusionInitialPassTests`: an input within the tolerance returns unchanged with 0 iterations and exit `AlreadyConverged`, reporting its initial residual; the same input iterates under a tighter tolerance | | Phases | fallback-only rows of `DiffusionBoundaryForcingTests` (primary limit 0: exit `Skipped`, all iterations fallback; 264 passes at α/h² = 8 and the test's relative tolerance 2e-7) | | Resources | `DiffusionResourceTests`: exact retained bytes in both scratch layouts (the separate-iterate layout, forced through an internal limit, gives bit-identical results); after warm-up no allocation or upload and exactly the reads above through a primary, a fallback and a primary-at-cap solve, the largest being the initial pass's partials | | Current method (`Kind=MethodSpecific`) | `HighAlphaCappedBroadbandSolveReportsTrueResidual`, `NearUnitFieldsHandTheFloatFloorToTheFallback` (exit `Stalled` or `Limit`), `ChebyshevDiffusionMethodTests` (ρ = S/(1 + S); an underestimated ρ ends converged or honest and finite through the stall exit; at a stall a component whose audited iterate is worse than b gets b back while another keeps its iterate; a non-finite input takes the `NonFinite` exit and never claims convergence; the Couette d = 115 log, CUDA only, with the reported residual equal to an independent host row with the device's fused folds and to the stored field's, converged exactly when that host residual is below the tolerance) | | Host math (CPU, `test/PolyCfd.Tests`) | `ChebyshevTests`: the Saad schedule and the Golub–Varga weights give the Chebyshev polynomial and the same polynomial, weight and `LogT` edges, and the stall rule on synthetic residual histories | **Known gaps:** no test distinguishes the `NonFinite` exit's restore of b from continuing with the non-finite iterate (a finite b never gives a non-finite iterate on these rows, so the test's non-finite value is in b itself), and no test covers a periodic wrap into a solid-centred active unknown (limits). Tests whose fully qualified name (namespace, class or method) contains a multigrid-suite substring (`Boundary`, `MovingWall`, ...) join the multigrid selection, whose snapshot is protected ([suite selection](https://github.com/hankbeasley/polycfd/blob/b7f5b9ba6f52355d9b196b586e4a6c8719d8feea/test/README.md#runner-settings-and-suite-selection)); choose new names deliberately. ## Changing the diffusion method 1. Replace `RunPrimary` and its step kernel. Build the kernel on `EvaluateRow`; do not copy the row. 2. Keep `DiffusionSolveSettings` and the recorded settings shape. Method constants stay private. 3. Extend `DiffusionPrimaryExit` for new exits; report only through `DiffusionResult`. 4. Run the contract tests unchanged. Rewrite or delete only `Kind=MethodSpecific` tests, and keep one independent host reference per boundary type through the new step kernel. 5. If device buffers change, update `ScratchBytes` and bump the estimator ModelVersion with every copy ([changing the model version](https://github.com/hankbeasley/polycfd/blob/b7f5b9ba6f52355d9b196b586e4a6c8719d8feea/docs/RESOURCE_PREFLIGHT.md#changing-the-model-version)). 6. Update this README and the pointers to it, then run the `verify-all` coverage a numerical solver change needs ([choose coverage](/docs/COMPLETE_VALIDATION#choose-coverage-for-the-change)); never refresh a baseline to pass.