# Projection and pressure solve

This folder holds the pressure projection (`ProjectionGpu`) and the pressure regions (`PressureComponentsGpu`). This
README is the implemented contract of the projection stages, the pressure solver interface `IPcgSolverGpu`, the one
PCG loop both solvers share (`PcgIteration`) and its preconditioner seam, the pressure regions, the solve result and
failure summary, the host-sync and allocation budget, and the tests
([contract homes](https://github.com/hankbeasley/polycfd/blob/b7f5b9ba6f52355d9b196b586e4a6c8719d8feea/docs/w4-design/README.md#where-implemented-contracts-live)). The multigrid preconditioner
(levels, V-cycle, smoother, transfers) is specified in the [Multigrid README](/src/PolyCfd.Gpu/Multigrid/README). The equations
and the physical meaning of the regions are in the [physics guide](/docs/PHYSICS_GUIDE#pcg-and-multigrid).

**Status (2026-09-26, SOL-04e):** one PCG loop (`PcgIteration`) serves the Jacobi and the multigrid solver; each
solver class is a shell that owns its operator, geometry and preconditioner. Results are unchanged from the two
loops it replaced, by the frozen-base comparison of the
[runtime contracts](/docs/RUNTIME_CONTRACTS#constraints-on-future-architecture-changes) (evidence in the
[P09.07 row](https://github.com/hankbeasley/polycfd/blob/b7f5b9ba6f52355d9b196b586e4a6c8719d8feea/docs/ARCHITECTURE_IMPLEMENTATION_EVIDENCE.md#ownership-map-and-conditional-decisions)). SOL-02b
changed only the multigrid preconditioner behind the seam ([Multigrid README](/src/PolyCfd.Gpu/Multigrid/README)).

| File | Role |
| --- | --- |
| `ProjectionGpu.cs` | Projection stages: divergence, right-hand side, solve, gradient correction; owns or borrows the solver |
| `PressureComponentsGpu.cs` | Connected pressure regions of the stencil and their constant-mode projections |
| `../IPcgSolverGpu.cs` | The solver interface and `PressureGeometryUpdate` |
| `../PcgIteration.cs` | `PcgIteration` (the loop), `IPressurePreconditioner`, `PressureLevelView` |
| `../PcgSolverGpu.cs` | Jacobi solver shell and `JacobiPreconditioner` |
| `../Multigrid/PcgMgSolverGpu.cs` | Multigrid solver shell; its preconditioner is `MultigridCycle` |
| `../PressureStencil.cs` | Halo fill plus operator: `ApplyCorrection` (homogeneous) and `ApplyPhysical` |
| `../Diagnostics/PressureFailureProbe.cs` | Bounded failure summary at a non-converged exit |
| `PolyCfd.Core/Numerics/IPressureSolve.cs`, `PressureSolveFailure.cs` | `SolverOptions`, `SolveResult`, `PressureSolveFailure` |

## Projection stages

`ProjectionGpu` builds the solver the scenario selects (`PcgSolverGpu`, or `PcgMgSolverGpu` for
`numerics.pressure.method = "pcgMultigrid"`) or borrows an injected `IPcgSolverGpu`; it disposes the solver it created
and, unless `ownsCustomSolver: false`, an injected one.

- **Static non-incremental** (`Project`): velocity boundaries, divergence, right-hand side `div u / dt`, one `Solve`
  for the full pressure, pressure boundaries, gradient and velocity correction. With `refine` (uniform static domains
  with homogeneous Neumann or periodic pressure only) a second solve from zero corrects the remaining divergence, and
  the result reports the failed solve, when there is one, with `with` (`Refinement` holds the correction).
- **Static incremental** (`ApplyPressurePredictor`, then `ProjectIncremental`): the previous pressure enters the
  viscous predictor; one solve for an increment from zero; the increment is accumulated into physical pressure, which
  `GaugePressure` then gauges.
- **Moving geometry** (`ProjectWithGeometricSource`): the continuity target includes the geometric source
  `-(V^{n+1} - V^n)/(dt V^{n+1})` (or the rigid-rotation wall flux). Once the projection has bound a moving geometry
  (`UpdateGeometry`), the right-hand side and the pressure are balanced per region before the solve with the solver's
  `PreparePressureRegions(...).ProjectOutMeans`, so the solve sees only roundoff-sized regional sources, and the
  pressure is balanced again after it.

Geometry updates reach the solver only through `PrepareGeometry` ([physics guide, section 5](/docs/PHYSICS_GUIDE#5-moving-geometry-and-pressure-components)).

## Solver interface

A pressure solver implements `IPcgSolverGpu`; every member is required and none has a default.

- `Solve(rhs, pressure, options, boundaries)` solves `A p = rhs` on interior-sized device buffers: `rhs` is never
  written, `pressure` is the initial guess on entry and the gauged solution on exit.
- `boundaries` is the whole-face `PressureBoundaryConditions` (callers pass `DevicePatchBcSet.Pressure`), so a solver
  cannot read velocity boundaries. Correction equations use its `Homogeneous` form; only physical residuals use the
  prescribed values (`PressureStencil`).
- `PrepareGeometry(PressureGeometryUpdate)` is the only geometry entry. Its resources are borrowed until the next
  update. A solver that cannot use an update throws (multigrid needs the device SDF and builder to rebuild its coarse
  levels); no update is silently ignored, and custom solvers are not discovered by concrete type.
- `CurrentVolumeWeights` is the solver's borrowed actual cell volumes of the current geometry revision (V^{n+1} for
  moving projection), valid until the next `PrepareGeometry`.
- `PreparePressureRegions(boundaries)` returns the solver's borrowed `PressureComponentsGpu`, rebuilt first when the
  boundary topology or the geometry revision changed.
- `GaugePressure(pressure, boundaries)` gives each region without a fixed-pressure boundary zero volume-weighted mean
  and leaves anchored pressure unchanged. Static incremental projection calls it after accumulating pressure.

A solver must report honestly (below), attach the bounded failure summary when it does not converge, stay within the
resource budget and write no files. It must not know the time step, velocity boundaries, the file system or run
directory, host fields or scenario settings, and logs at most warnings.

## One PCG loop

`PcgIteration` is the preconditioned conjugate-gradient loop of both built-in solvers. It owns:

- the five finest-level iteration vectors x, r, z, p, A·p (`PcgIteration.DeviceBytes` = 20 bytes per cell);
- the region cache: one `PressureComponentsGpu`, keyed by (periodic axes, fixed-pressure faces, geometry revision)
  and rebuilt from the finest operator's own apertures, classes and volumes (`GeometryChanged` bumps the revision);
- the recurrence, in the operation order both solvers had before they were merged: r₀ = rhs − A x₀ with the physical
  boundaries; the initial projection and imbalance measure; z = M r, p = z; per iteration A·p (homogeneous halo),
  α = r·z / p·Ap, x += αp, r −= αAp, the per-region projection of r, ‖r‖_V, the stop test, z = M r, β, p = z + βp;
- the stopping rules: absolute `‖r‖_V < AbsTolerance` or relative `‖r‖_V/‖r₀‖_V < RelTolerance`, and the breakdown
  thresholds |p·Ap| < 1e-30 and |r·z| < 1e-30;
- the gauge of x at the start and at the end of every solve;
- the non-converged exit with `PressureFailureProbe`, before the final gauge;
- every `SolveResult` a solve returns (the only construction sites in runtime code).

The solver classes are shells. `PcgSolverGpu` owns its operator, volumes, halo workspace and the
`JacobiPreconditioner` (the inverse diagonal from `MultigridKernels.BuildJacobiDiagonal`, rebuilt on every geometry
update). `PcgMgSolverGpu` owns the hierarchy, its geometry rebuilds and the `MultigridCycle`. Each lends its finest
level to the loop per call as a `PressureLevelView` (stencil, BLAS, volumes, byte classes) and passes the cell volume
it always used for the refusal limit. `PcgIteration` must not know which preconditioner runs, the levels, the boundary
types (descriptors pass through), how topology is found, geometry building, files or scenario settings.

**Preconditioner seam.** `IPressurePreconditioner.Apply(r, z, boundaries)` computes z = M r. M must be linear and the
same fixed map on every call, self-adjoint and positive (semi)definite in `<a, b>_V = Σ V_i a_i b_i` on the
compatible subspace, use homogeneous boundary values only, never write r, allocate nothing and make no host
synchronization. The implementations are `JacobiPreconditioner` (pointwise D⁻¹) and `MultigridCycle` (one V-cycle from
a zero start; its symmetry and positivity limits are in the [Multigrid README](/src/PolyCfd.Gpu/Multigrid/README#v-cycle)).

## Pressure regions

`PressureComponentsGpu` finds the connected fluid regions of the stencil: cells joined through positive face
apertures, and through periodic links when both faces of an axis are periodic; solid cells belong to none. A region
with no fluid face on a fixed-pressure boundary (`Topology` reads the six whole-face descriptors in face order) has a
constant nullspace, also beside an outlet. Only exactly disconnected regions count.

| `Mode` | When | Projection |
| --- | --- | --- |
| `None` | every region is anchored | nothing |
| `Global` | one unanchored region holds all fluid volume | `GridBlasGpu.ProjectOutMean` over the whole grid (also shifts zero-volume solid cells) |
| `Regions` | several regions, or an unanchored region beside an anchored one | each unanchored region's own mean; solid and anchored cells unchanged |

- `Project` is device-only with zero host synchronization: up to `SlotsPerPass` = 16 regions in one fused pass (per
  group partials, then fixed-order totals), more through a per-region device tail. A topology with more than 16
  regions is classified on the host at every rebuild (about 12 bytes per cell to the host;
  [resource estimates](https://github.com/hankbeasley/polycfd/blob/b7f5b9ba6f52355d9b196b586e4a6c8719d8feea/docs/RESOURCE_PREFLIGHT.md)).
- `ProjectAndMeasure` reads one scalar per region (in one copy) and runs once per solve on the initial residual, and once
  more in a failure summary, never in the loop.
  It returns the weighted norm of the removed constants beyond the global source a domain without fixed-pressure faces
  tolerates; that source is attributed to the largest region.
- `ProjectOutMeans` (moving projection) needs the labels and resets solid cells to zero.
- `SolveResult.RegionalImbalance` is that norm relative to max(‖rhs‖_V, ‖r₀‖_V), reported on the `Regions` path and on
  the `Global` path when the boundaries have fixed-pressure faces that no fluid touches (otherwise null). A solve is refused as not converged with
  `InconsistentRhs` when it exceeds `RefusalLimit(cellVolume, fluidVolume)` = max(`IncompatibilityLimit` = 1e-3,
  5·sqrt(k V_cell / V_fluid)) for k unanchored regions and the removed norm exceeds the absolute tolerance; the
  derivation and the detection limits are in the [physics guide](/docs/PHYSICS_GUIDE#pcg-and-multigrid).
- `DeviceWorkspaceBytes(cells)` is the estimator's `pressure-regions` allowance.
- It must not know Krylov state, the preconditioner, iteration counts or the time step.

**Determinism.** The weighted dot products and norms of the loop (`DotWeighted`, `Nrm2Weighted`) are per-block double
partials summed on the host in index order, and the per-region path uses fixed trees and fixed-order totals: anchored
and multi-region solves are reproducible. The `Global` path and `GridBlasGpu.Dot` (the right-hand-side mean, the
failure summary's means) use an ILGPU reduction whose cross-group combine is atomic, so their last bits can vary
between runs, and singular single-region solves inherit that variation.

## Results

`SolveResult` is a `sealed record` built only by `PcgIteration`; callers that combine results use `with` (the refined
projection), and the static integrator step returns the projection's result itself, so every field reaches
`StepReport`. `Converged` means the stopping rule held. `Iterations` counts completed updates, except after an r·z
breakdown (`RZZero`): that iteration's update of x and r has been applied but is not counted. At the cap it equals
`MaxIterations`. `FinalResidual` and `InitialResidual` are the recursive residual norms `‖r‖_V` (after projection);
below a relative tolerance of about 1e-6 the float true residual no longer follows the recursive one
([physics guide](/docs/PHYSICS_GUIDE#pcg-and-multigrid)). `ConvergenceReason` names the rule that held
(`None` when not converged), and `BreakdownReason` is `MaxIterations`, `PAPZero`, `RZZero` or `InconsistentRhs` on
failure. `RhsMean` is set when a domain without fixed-pressure faces has a right-hand-side mean above 1e-10.
`SolverOptions.LogConvergence` is read by no solver; it stays only because `scenario-run.json` serializes it.

### Pressure failure summary

A non-converged solve returns `SolveResult.FailureSummary` (`PressureSolveFailure`); converged results leave it null
and no work is done for it. `PressureFailureProbe` runs at the non-converged exit (maximum iterations, `pAp`/`rz`
breakdown, or a refused regional imbalance), on the finest level and before the final gauge, so constant drift of the
iterate stays visible. It reports:

- the volume-weighted norm of the true residual `rhs − A x` (physical boundaries) after the solver's own per-region
  compatibility projection, comparable with `FinalResidual`, and separately the norm of the constant components that
  projection removed;
- the worst true-residual, right-hand-side and solution cells with (i,j,k), signed value, cell class and volume
  fraction (unweighted maxima; lowest index on ties; NaN is counted, not located; an infinity is located and counted);
- the weighted right-hand-side mean, the weighted solution mean and the weighted solution RMS after each unanchored
  region's mean is removed (a large mean with a small deviation is constant drift);
- non-finite counts of rhs, residual and solution.

It reuses the loop's dead scratch vectors and a `DeviceStatisticsGpu` allocated at the first failure (memory category
`pressure-failure-statistics`). No field crosses to the host: at most `PressureFailureProbe.MaxHostReads` = 17 reads
of a few scalars, the double block partials of two weighted norms (one per 256-cell block, 64 KiB each at 128³), one
sum per region on the per-region path, three 32-byte statistics and nine single elements. `SolveResult.ToString()`
appends the summary as one invariant-culture line of a few hundred characters, and `ScenarioStepAcceptance` includes it
(and any `RegionalImbalance`) in its step-failure exception, whose message `ScenarioRunner` records as the
`status.json` message and the CLI prints; the workbench keeps a bounded prefix of it as the worker's failure reason.
Statistics are global over the finest level; per-region offsets appear only through the gauged RMS. **Runtime code
writes no files and copies no fields**; a full-field snapshot of a failed solve would be a borrowed-buffer observer
attached by a run-directory owner.

## Resource budget

Counted as `DeviceMemoryAudit` counts: explicit application allocations, uploads, downloads and barriers.
`GridBlasGpu` reads its reduction partials and scalars through the backend's audited `CopyToHostRange` into reused host
arrays.

| Scope | Host reads | Allocations and uploads |
| --- | --- | --- |
| PCG iteration | exactly 3: p·Ap, ‖r‖_V, r·z (each one array of double block partials) | none |
| PCG solve, outside the loop | a constant independent of the grid: ‖r₀‖_V, r·z, the right-hand-side mean and norm, the regional measure (at most 8) | none after warm-up |
| V-cycle, `JacobiPreconditioner`, `PressureComponentsGpu.Project` | none (no host synchronization) | none |
| Failure exit | at most 17 (`PressureFailureProbe.MaxHostReads`): scalars, two arrays of double block partials (8 bytes per 256-cell block), one sum per region and single elements; never a field | the statistics workspace, once |
| Any solve | no read larger than one double per reduction block (never a field); at most one explicit barrier per read; no files; no console output | — |

ILGPU-internal temporaries of `CreateReduction` (`GridBlasGpu.Dot` and `ProjectOutMean`: the unanchored-region paths and
the failure summary's means) are outside the audit. The SOL-04 fingerprint probe measured no free-memory drop over 200
warm solves of each pressure case on the SOL-02a build (`output/validation/sol02a-free-memory.json`, 2026-09-26). The
SOL-02b V-cycle calls no ILGPU reduction and allocates only at setup; it is covered by the application-level audit
(`PressureSolverResourceTests`), not by that probe, which was deleted after SOL-02. Region rebuilds (at the first solve of a topology and after every geometry update) are setup work: they
allocate the region workspace once and read back a few scalars.

## Footprints

`PcgIteration.DeviceBytes(cells)` (20c, estimator category `pressure-pcg` for multigrid, part of `pressure-jacobi` for
Jacobi), `JacobiPreconditioner.DeviceBytes(cells)` (4c) and `PcgSolverGpu.DeviceBytes(size)` report the solver bytes
beside their allocations; the multigrid levels and the lazy region and failure workspaces are in the
[Multigrid README footprint table](/src/PolyCfd.Gpu/Multigrid/README#footprints). `ScenarioMemoryEstimator` composes them.

## Testing

Pressure contracts are `[Theory]` rows over the solver kind (`PressureSolverKind.PcgJacobi`, `PcgMultigrid`, through
the production constructors in `PressureSolvers.Create`) or over `bool multigrid`, and judge correctness with an
independent host operator (`PressureSolvers.HostApply`, double precision, apertures and volume fractions of the block),
never with the solver's own residual. The cut-cell rows also judge with Core's CPU `PressureOperator` and check that the
two host operators agree to float roundoff on the same block. They assert no iteration counts; each drift or residual
bound is a small multiple of the value measured for its row. All of them join the protected multigrid selection by name.

| Contract | Tests |
| --- | --- |
| Zero right-hand side: no iteration, exact zero; tight singular (periodic, zero-gradient; 4 seeds) and cut-cell (sphere, STL) solves meet the independent residual; honest cap (`Iterations = MaxIterations`, `FinalResidual` agrees with `FailureSummary.TrueResidualNorm`); a geometry revision changes the operator; the preconditioner is a fixed, linear, symmetric, positive map that never writes r | `Numerics/Pcg/PressureSolverContractTests` |
| Transfer budget after warm-up (0 allocations and uploads; reads ≤ 3 per iteration + 8, plus 17 when capped; at most one explicit barrier per read; largest read ≤ one double per block) | `Numerics/Pcg/PressureSolverResourceTests` |
| Dirichlet manufactured profile | `NonzeroPressureBoundaryTests` |
| Regions: sealed pockets with and without an outlet, balanced pockets, roundoff-sized sums projected, closed-domain net source, a region fed by net inflow refused, a moving topology change | `Multigrid/PressureRegionSolveTests` |
| Region labels and projections against a host union-find | `Projection/PressureComponentsTests` |
| Failure report: no files and no fields, worst cell, cut-cell class and volume, survives the static step and the refined projection | `Multigrid/PressureSolveFailureReportTests` |
| Solver footprints | `Multigrid/PressureSolverFootprintTests` |

Focused command (the multigrid suite's pressure part):

```bash
dotnet test test/PolyCfd.Gpu.Tests -c Release --settings test/gpu.runsettings --filter "(FullyQualifiedName~Pcg|FullyQualifiedName~Pressure|FullyQualifiedName~Projection|FullyQualifiedName~Multigrid)&Backend!=CPU"
```

## Adding a pressure method or preconditioner

1. A new preconditioner implements `IPressurePreconditioner` with the contract above and lives with its owner (a
   shell like `PcgSolverGpu`, or a class like `MultigridCycle`). The loop, regions, stopping rules, failure exit and
   results stay in `PcgIteration`; a change there applies to both solvers at once.
2. Add a solver shell that lends its finest level as a `PressureLevelView`, and a `DeviceBytes` beside its
   allocations for `ScenarioMemoryEstimator` (bump the ModelVersion with every copy:
   [changing the model version](https://github.com/hankbeasley/polycfd/blob/b7f5b9ba6f52355d9b196b586e4a6c8719d8feea/docs/RESOURCE_PREFLIGHT.md#changing-the-model-version)).
3. Add a `PressureSolverKind` row to the contract theories and the resource test; they must pass unchanged. Mark
   method-specific behaviour `Kind=MethodSpecific` and iteration comparisons `Kind=Efficiency`
   ([test kinds](https://github.com/hankbeasley/polycfd/blob/b7f5b9ba6f52355d9b196b586e4a6c8719d8feea/test/README.md#test-kinds)).
4. A third production method also replaces the `bool Multigrid` selection with an enum and one factory; recording the
   solver in `ScenarioEffectiveSettings` is a separate decision (it changes every `scientificId`).
5. Update this README and its pointers; run the multigrid and quick `verify-all` suites
   ([choose coverage](/docs/COMPLETE_VALIDATION#choose-coverage-for-the-change)); never refresh a baseline
   to pass.