# Algorithms and numerical review

**Review companion, 2026-09-26; GPU execution and architecture requirements updated
2026-09-26; current solver descriptions reconciled with the merged source.** Current formulas below were checked against the
linked source. Proposed formulas describe TRN-01/S2/S3 design obligations, not
implemented or scientifically validated capabilities. The scalar example is an
illustration, not the selected production momentum algorithm. Start with the
[algorithm explainer](/algorithms.html) or [transport plan](/docs/amrex-alignment/SINGLE_LEVEL_TRANSPORT).
Current status belongs on the [work board](/docs/WORKBOARD).

## Reading map

Read the physical model, then the current update, then the proposed flux balance.
The distinction matters: choosing a conservative transport operator changes one
part of a coupled incompressible solver, not its whole accuracy or conservation story.

| Question | Detail |
| --- | --- |
| What physical problem is solved? | [Governing equations](#governing-equations) and [MAC locations](#mac-grid-and-control-volumes) |
| What does today's code calculate? | [Sampling and correction](#current-semi-lagrangian-and-maccormack), then [full step](#current-time-step-and-projection) |
| Why does a flux update conserve? | [Signs, units and cancellation](#proposed-shared-flux-update) |
| What could a concrete algorithm look like? | [Illustrative scalar method](#illustrative-regular-grid-scalar-algorithm) |
| What remains to derive? | [Staggered momentum and geometry](#mac-momentum-and-cut-cell-derivation), [S2/S3](#s2-and-s3-coupling) |
| What would establish correctness? | [Review and validation gates](#review-and-validation-gates), with [source map](#source-map) |

Method selection and support rules are defined in [TRN-01](/docs/amrex-alignment/SINGLE_LEVEL_TRANSPORT).
Use one resolved method per run; preserve useful current methods while evaluating
one conservative candidate. A method must state its timestep, stencil, geometry,
boundary and integration requirements. A shared interface must expose required
stages instead of assuming every algorithm is an interchangeable `Advect()` call.
Neither data-model alignment nor a method name establishes numerical compatibility.

## Required execution architecture

**Approved target; not fully implemented:** [GPU-01](/docs/plans/gpu-resident-processing)
requires all runtime numerical and field-derived processing to execute on one GPU.
This includes geometry builds/rebuilds, transport/ghost kernels, pressure topology,
complete reductions, PCG coefficients/convergence, fixed/adaptive timestep decisions,
diagnostics, coarse/fine transfers and regridding. A scalar, flag, count or partial
sum downloaded for CPU calculations or numerical decisions is a computational
readback and is forbidden. There is no silent CPU fallback. GPU-01 acceptance is
required for every supported AMR-02/TRN-01/S2/S3 path.

Host work is configuration parsing, static capability/method binding, allocation
and kernel scheduling from metadata, and publication of completed device-produced
outputs. Numerical dependencies use device flags/masks and device-produced values.
Requested capture and final status are bounded terminal outputs, not routes to host
residual, CFL, diagnostic or regridding computations. Current source is evidence
of behavior, not proof that this stronger execution requirement already passes:
[CflGpu](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/src/PolyCfd.Gpu/Diagnostics/CflGpu.cs) currently computes final limits
from host-returned reductions, the shared
[PcgIteration](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/src/PolyCfd.Gpu/PcgIteration.cs) uses host scalar
coefficients/decisions for both pressure solvers, and
[PressureComponentsGpu](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/src/PolyCfd.Gpu/Projection/PressureComponentsGpu.cs)
still downloads rebuild metadata and compatibility measurements, with host
classification for more than 16 regions. Its per-region mean subtraction already
completes on-device. Diffusion's Chebyshev recurrence and audit finalization also
use host calculations ([implemented contract](/src/PolyCfd.Gpu/Diffusion/README)).
These remaining paths are GPU-01 gaps, not scalar exemptions.

Use [AMR-02's semantic split](/docs/amrex-alignment/DATA_MODEL_ALIGNMENT): `MeshPatch` identifies a mesh
partition; `BoundaryRegion` entries in `DomainBoundaryConditions` describe physical
boundaries, with `DeviceBoundaryRegion`/`DeviceDomainBoundaryConditions` device bindings.
`IEmbeddedBoundarySource`/`SdfEmbeddedBoundarySource` supply geometry inputs;
`EmbeddedBoundaryData`/`DeviceEmbeddedBoundaryData` hold derived metrics. Kernels
borrow `EmbeddedBoundaryMetricsView` plus a compact derived `PatchKernelDescriptor`.
`FieldPlacement`, `IndexBox`, `DomainGeometry`, `LevelLayout` and `LevelFields` retain
their separate roles. Those target names do not rename current source by documentation.
The level coordinator schedules method stages; patch kernels execute the numerical
operations. Field/geometry views identify layout and geometry revisions, time and
stage; ghost validity is refreshed after its producer changes. Retained storage
or an old completed-frame identity does not certify stage freshness.

## Governing equations

**Current physical model:** constant-density incompressible Newtonian flow.
Velocity is `u = (u,v,w)`, density is constant `rho`, kinematic viscosity is `nu`,
body acceleration is `f`, and kinematic pressure is `p = P/rho`:

```text
partial_t u + (u dot grad)u = -grad p + nu Laplacian(u) + f
div u = 0
```

The stored pressure has units `length^2/time^2`, not pascals. For component `a`,
the continuum identity is

```text
div(u_a u) = (u dot grad)u_a + u_a div u.
```

Therefore advective and conservative momentum forms agree when `div u = 0`.
Their discrete versions need not agree: interpolation, discrete divergence,
limiting, boundaries and stage timing affect the identity. Projection does not
automatically turn a nonconservative remap into a shared-flux momentum update.
The [physics guide](/docs/PHYSICS_GUIDE) defines the supported physical scope.

For a scalar concentration `q`, `partial_t q + div(q u) = s` describes a conserved
amount with source density `s`. For momentum take `q = rho u_a`; pressure and
viscosity then enter as forces. Using `q = u_a` instead tracks momentum divided
by the constant density. State which convention a diagnostic uses.

## MAC grid and control volumes

**Current storage:** pressure at cell centers; each velocity component at its
normal face center. With origin zero and spacings `(dx,dy,dz)`:

```text
p(i,j,k): ((i+1/2)dx, (j+1/2)dy, (k+1/2)dz)
u(i,j,k): (i dx,     (j+1/2)dy, (k+1/2)dz)
v(i,j,k): ((i+1/2)dx, j dy,     (k+1/2)dz)
w(i,j,k): ((i+1/2)dx,(j+1/2)dy, k dz)
```

See [MAC conventions](/docs/MAC_GRID_CONVENTIONS). A stored face sample is not
automatically an average over a derived momentum control volume. On a regular
interior grid, the natural dual volume around an `u` sample is shifted half a
cell in X relative to a pressure cell. Its transport surfaces differ from the
pressure cell's surfaces. V/W have corresponding shifts. Physical boundary
dual volumes, periodic duplicate endpoints and cut dual volumes require explicit
weights; summing every stored face with full cell volume double-counts periodic
endpoints. Existing [exact periodic diagnostics](/docs/PERIODIC_EXACT_VALIDATION)
use half endpoint weights, which do not by themselves define cut-cell momentum.

Device velocities currently have dense interior buffers and inline boundary
stencils; pressure also has packed face halos. There are no general GPU velocity
edge/corner ghost arrays. [Halo layout](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/docs/halo-layout.md) describes actual storage.
TRN-01 must specify its stencil before S2 chooses exchange/storage requirements.

## Current semi-Lagrangian and MacCormack

**Implemented:** [AdvectionGpu](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/src/PolyCfd.Gpu/Advection/AdvectionGpu.cs),
methods `Advect`, `AdvectKernelImpl`, `Sample` and `CorrectKernelImpl`.
Let `I_c` be trilinear interpolation of component `c` on its native staggered
lattice, and `I_u` the resulting vector interpolation. The original velocity
`u^n` carries every component. At output location `x`:

```text
x_mid = x - (dt/2) I_u[u^n](x)
x_dep = x - dt I_u[u^n](x_mid)
S_dt[u^n](q)(x) = I_c[q](x_dep)
```

The geometry-aware path additionally projects both `x_mid` and `x_dep` out of
the solid using the SDF. It excludes blocked interpolation samples and renormalizes
the remaining weights. A face classified Closed or with aperture below `0.01`
is treated as blocked and keeps its input value; an entirely rejected interpolation
stencil returns zero. Nonperiodic samples clamp to the domain; periodic samples
wrap. These existing rules are geometry handling, not a finite-volume balance.

`SemiLagrangian` uses `q_forward = S_dt[u^n](q^n)`. `MacCormack` uses:

```text
q_forward  = S_dt[u^n](q^n)
q_backward = S_-dt[u^n](q_forward)       # SAME original carrying velocity
q_trial    = q_forward + (q^n - q_backward)/2
q_new      = clamp(q_trial, q_min, q_max)
```

`q_min/q_max` are the eligible original samples in the forward interpolation
stencil. This clamp bounds the output against that stencil; it does not impose
an equal/opposite transfer between neighboring control volumes. The method
corrects a round-trip interpolation error, with lineage in
[Selle et al., An Unconditionally Stable MacCormack Method](https://physbam.stanford.edu/papers/stanford2006-09.pdf).
That paper's formal results do not prove the order of our geometry-limited,
split pressure/diffusion solver. The word MacCormack also names conservative
schemes elsewhere; the nonconservation statement here concerns this implementation.

The [advection tests](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/test/PolyCfd.Gpu.Tests/Advection/AdvectionSchemeTests.cs)
check uniform fields and Gaussian transport/overshoot. They do not prove general
momentum conservation, cut-cell accuracy or full-solver convergence order.

## Current time step and projection

**Implemented static geometry:** the integrator applies advection, body increments,
the old pressure force, backward-Euler diffusion, immersed-face constraints and
incremental projection. Boundary operations occur between these stages. Denote
the constrained preprojection velocity by `u_star`:

```text
u_adv = Advect(u^n, dt)
u_rhs = u_adv + dt f - dt G(p^n)
(I - dt nu L) u_diff = u_rhs
u_star = ApplyImmersedConstraints(u_diff)
A phi = b,       A = -D_alpha G,       b = -D_alpha(u_star)/dt
u^(n+1) = u_star - dt G(phi)
p^(n+1) = p^n + phi                    # restore gauge and physical pressure BCs
```

`D_alpha` is outward open-area flux divided by effective fluid volume, `G` the
effective staggered correction gradient including solid-face masks and correction
boundary treatment, `L` the viscous discrete Laplacian, and `phi` a pressure
increment. The applied correction skips solid-adjacent faces; the cancellation
requires the pressure operator to match that correction, not merely an unmasked
gradient kernel. The sign check is `D_alpha(u_new) = D_alpha(u_star) + dt A phi`.
For identical operators/BC treatment and residual `r = b - A phi`, this becomes
`D_alpha(u_new) = -dt r`. Reapplied constraints, compatibility projections and
finite precision require checking the actual post-step divergence as well.
Under GPU-01 the residual, norms, compatibility/gauge corrections, stopping decisions
and diagnostic finalization all remain on device; a final status record may report
the completed decision without feeding a host numerical control loop.

The pressure operator limits effective small volumes and masks solid rows;
its weights are not a derivation of conservative momentum volumes. Pressure
inner products use double-accumulated volume weights with float fields and Half
geometry. Constant prescribed physical pressure gives zero prescribed pressure
increment. Closed/periodic pressure components have constant null modes; fixed
pressure anchors only components connected to that boundary. Both pressure solvers
now share a PCG loop that projects every unanchored region's constant mode and
gauges each region, including static disconnected domains. Moving projection also
balances its geometric-source right-hand side per region. Preserve this existing
component compatibility, gauge and imbalance-refusal behavior as topology changes;
do not replace several null modes by one global mean. The
[projection guide](/src/PolyCfd.Gpu/Projection/README) owns these implemented contracts.

**Implemented moving geometry differs:** it rebuilds geometry at `t+dt`, initializes
changed cells, advects, adds sources, diffuses and imposes wall velocity, then uses
full-pressure projection with a geometric source. It does not use the static
old-pressure predictor/increment accumulation. For rigid rotation:

```text
u_wall(x) = omega cross (x - center)
S_raw = D_alpha(u_wall)
b_raw = -(D_alpha(u_star) - S_raw)/dt
b_compatible = component-compatible RHS
S_compatible = D_alpha(u_star) + dt b_compatible
continuity error = D_alpha(u_new) - S_compatible
```

Other motion uses the existing limited old/current-volume source. A small moving
continuity error does not prove exact geometric conservation or momentum balance.
See [moving design](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/docs/MOVING_GEOMETRY_DESIGN.md) and
[moving fan qualification](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/docs/MOVING_FAN_VALIDATION.md). Backward-Euler diffusion
and splitting remain accuracy limitations even if new advection is second order.
AMReX-Hydro's [projection overview](https://amrex-fluids.github.io/amrex-hydro/docs_html/Projections.html)
is useful context; its other velocity/pressure arrangements are not our current solver.

For the implemented linear-solver algorithms behind these equations, continue to
the [shared pressure iteration](/src/PolyCfd.Gpu/Projection/README#one-pcg-loop),
[multigrid preconditioner](/src/PolyCfd.Gpu/Multigrid/README),
[diffusion method and residual definition](/src/PolyCfd.Gpu/Diffusion/README)
and [boundary stencils](/docs/BOUNDARY_CONDITIONS_HALOS). Diffusion currently uses
Chebyshev semi-iteration on undamped Jacobi followed, when needed, by a bounded
colored Gauss-Seidel phase. The [pressure multigrid design](/docs/GPU_MULTIGRID_DESIGN),
[cut-cell multigrid design](/docs/CUT_CELL_MULTIGRID_DESIGN) and
[diffusion history](/docs/ADAPTIVE_DIFFUSION_RELATIVE_TOLERANCE) retain rationale and
labeled historical algorithms. Solver coarsening is not physical AMR, and convergence
of an algebraic solve is separate from the accuracy of its discretization.

## Proposed shared-flux update

**Proposed contract:** choose a fixed orientation for each transport interface
`f`, from volume `L` to `R`, with normal `n_f`. `Q_i = V_i q_i` is the integrated
quantity, and `s_if` is `+1` for the L volume and `-1` for R. Define the signed
amount transported over a step, already including physical open area and time:

```text
I_f = integral_[t_n,t_(n+1)] integral_[open face f] q (u dot n_f) dA dt
Q_i^(n+1) = Q_i^n - sum_f(s_if I_f) + B_i
```

`B_i` is the integrated source/boundary-force contribution appropriate to `q`.
For a scalar concentration, `I_f` has units of scalar amount. For physical momentum
`q=rho u_a`, it has units `mass*length/time`; for `q=u_a`, units `length^4/time`.
If storing flux density instead, specify its units and apply `alpha*A*dt` exactly
once. This fixed-volume equation needs the moving-volume extension below for motion.

Each internal interface contributes `-I_f` to L and `+I_f` to R, so summing all
volumes cancels internal transfers algebraically. Example: initially `Q_L=4`,
`Q_R=1`, with `I_f=+0.3` and no other transfers gives `3.7,1.3`; the total remains
5. A negative `I_f` reverses the transfer. Independently using outgoing `0.3` and
incoming `0.25` would lose `0.05`. Finite precision still requires roundoff-aware
budget tolerances; identical flux data do not guarantee bitwise invariant sums.
Compute the shared transfers, their update and full double-accumulated budget
reduction on GPU. Returning partial sums for CPU finalization violates GPU-01 even
if the main update is device resident.

Constant preservation is an additional condition: for constant `q`, its update
vanishes only when the transporting volume fluxes satisfy the required discrete
continuity and sources/BCs agree. Shared-flux cancellation alone does not prove
this, boundedness, physical energy evolution or accuracy.

## Illustrative regular-grid scalar algorithm

**Illustrative only:** a fully specified small building block for discussion.
It solves `q_t + a q_x = 0` for constant speed `a`, periodic uniform 1D cells of
width `h`, unit cross section, fixed volumes, and cell averages `q_i`. It is not
a 3D momentum, cut-cell, thermal or selected production algorithm.

Use a limited piecewise-linear reconstruction (slopes here have units of `q`):

```text
sigma_i = minmod(q_i - q_(i-1), q_(i+1) - q_i)
minmod(x,y) = sign(x) min(|x|,|y|) if x*y > 0; otherwise 0
q_i^L = q_i - sigma_i/2;    q_i^R = q_i + sigma_i/2
F_(i+1/2)(q) = a q_i^R       if a >= 0
             = a q_(i+1)^L   if a < 0
L_i(q) = -(F_(i+1/2)(q) - F_(i-1/2)(q))/h
```

Advance with the two-stage, second-order strong-stability-preserving Runge-Kutta
formula, reconstructing and filling periodic neighbors separately at each stage:

```text
q_stage = q^n + dt L(q^n)
q^(n+1) = (q^n + q_stage + dt L(q_stage))/2
I_(i+1/2) = (dt/2) [F_(i+1/2)(q^n) + F_(i+1/2)(q_stage)]
```

The last line is the actual integrated flux of this update; storing only its
second-stage flux would give the wrong diagnostic/register balance. For `a>0`,
a forward-Euler substep can be written using `C=a dt/h` as

```text
q_i + dt L_i(q) = (1/2)q_i^L + (1/2-C)q_i^R + C q_(i-1)^R.
```

For `0 <= C <= 1/2`, the coefficients are nonnegative and sum to one; minmod
endpoints stay within their neighboring cell range. Thus each such substep is
bounded by the input's global range, and the two-stage convex combination retains
that bound. Negative `a` gives the symmetric argument. This is a sufficient bound
for this example, not a universal CFL rule. Conservation follows separately from
flux telescoping. Smooth-region spatial reconstruction and time integration are
nominally second order; limiter activation can reduce local order near extrema.

The two-stage form has the SSP Runge-Kutta basis described by
[Gottlieb, Shu and Tadmor (2001)](https://www.math.umd.edu/~tadmor/pub/linear-stability/Gottlieb-Shu-Tadmor.SIREV-01.pdf).
Our explicit boundedness calculation above is for these particular scalar assumptions.
Multidimensional varying velocity, staggered momentum, boundary reconstruction and
small cut cells require new derivations and timestep limits. AMReX-Hydro's
[reconstruction helpers](https://amrex-fluids.github.io/amrex-hydro/docs_html/Utilities.html)
illustrate further choices; this example does not reproduce its Godunov algorithms.

## MAC momentum and cut-cell derivation

**Unresolved production choices:** TRN-01 must derive the following before claiming
conservative MAC transport. A scalar proof cannot supply these missing contracts.

- Derive each component's dual fluid volume and transport surfaces, quadrature,
  face states and advecting velocity. Establish a discrete continuity identity on
  those volumes. Decide whether additional intermediate velocity projection or
  changed stage timing is needed; document its coupling to pressure and viscosity.
- Define pressure and viscous tractions on physical/immersed boundaries and the
  work they do. On an impermeable stationary wall the advective normal transfer
  is zero, but pressure and viscous forces can change total fluid momentum.
- Derive cut dual geometry and physical aperture areas consistently, including
  centroid/normal information needed by reconstruction. Existing pressure-cell
  fractions or arbitrary cell-to-face averages are not a sufficient derivation.
- Choose a small-cell treatment and prove that its redistribution/merging preserves
  the declared balance without crossing a solid. Small volumes can impose severe
  explicit restrictions. [AMReX's embedded-boundary discussion](https://amrex-codes.github.io/amrex/docs_html/EB.html#small-cell-problem-and-redistribution)
  explains the issue; adoption requires MAC-specific validation here.
- Qualify motion separately, including newly exposed/covered fluid and each stage's
  geometry. Do not equate an SDF rebuild, pressure compatibility correction or
  copied neighbor value with a conservative transfer through moving geometry.

For reference, the moving-control-volume momentum balance uses boundary velocity
`w_b`, outward normal `n`, and physical stress `T = -P I + rho nu (grad u + grad u^T)`:

```text
d/dt integral_C(t) rho u dV
  = -integral_boundary rho u [(u-w_b) dot n] dA
    + integral_boundary T n dA + integral_C(t) rho f dV
dV/dt = integral_boundary w_b dot n dA.
```

The second identity is the geometric conservation obligation. A discrete method
must reconcile volume changes with swept geometry and relative fluxes, and preserve
constant states where the physical BCs permit them. At an impermeable moving wall,
`(u-w_b) dot n=0`, yet traction and wall work remain. Pressure-solver coarsening,
physical AMR coarsening, and previous/current motion geometry are distinct objects.

## S2 and S3 coupling

**Planned S2:** for a conservative method, neighboring same-resolution patches must
use the same interface transfer with opposite signs, including all temporal stages, limiter inputs and
periodic ownership. Patch edges are not physical walls. Native field/flux ownership,
ghost updates and boundary stencils must preserve the accepted single-grid method.
Retained interpolation methods instead require equivalent departure-point sampling,
including forward/reverse/limiting stages; they do not acquire a conservative-flux
guarantee. Test partition independence separately for each supported method and geometry.
`DomainBoundaryConditions` applies at the physical domain, while `MeshPatch`
interfaces use device exchange coordinated at the level. Refresh stage-valid
ghosts before each dependent kernel; component connectivity, duplicate ownership
and final reduction decisions stay on GPU.

**Planned S3:** covered coarse values are replaced by appropriately weighted fine
values (average-down); the composite total counts uncovered coarse plus fine data
once. For a scalar, `V_C q_C = sum_children V_f q_f` requires consistent volumes.
MAC velocity restriction requires its own face/dual-volume and divergence contracts.

Average-down does not correct coarse/fine interface transfers. Orient the interface
outward from an uncovered coarse cell toward fine coverage. If that cell used
`I_C` while the neighboring fine cells received net `I_F` over the same area/time,
the correction to its integrated quantity is

```text
Q_C <- Q_C - (I_F - I_C),
I_F = sum_fine_faces sum_fine_substeps I_f.
```

This replaces its coarse loss by the actual fine transfer. When `I_C=0.3` and
`I_F=0.35`, it must lose another `0.05`. If storing an average, divide the correction
by the appropriate coarse volume. The opposite orientation reverses the signs.
Initially equal level timesteps remove substeps, not spatial flux mismatch.
[AMReX's flux-register example](https://amrex-codes.github.io/amrex/docs_html/AmrCore.html#fluxregisters)
provides the scalar model for this synchronization. Staggered momentum and cut
interfaces still need their derived geometry and transfer accounting.

Refluxing does not enforce composite incompressibility, pressure continuity,
viscous coupling or accuracy. S3 must coordinate it with composite pressure and
diffusion so later corrections respect the final budgets. Mesh adaptation also
requires conservative remapping and rebuilt connectivity/coverage. See the
[S2/S3 plan](/docs/amrex-alignment/S2_S3_PLAN); AMReX alignment alone implements none of these operations.
Tagging, compaction, clustering, transfer plans and capacity decisions remain on
GPU, using explicitly validated bounded/preallocated resources. Downloading tags
or partial plans for CPU clustering is prohibited; terminal layout output for
reproducibility is separate from the numerical regridding path.

## Review and validation gates

The review has two parts: **check the equations and assumptions**, then **measure
the implementation against independent references**. Documentation and algebra are
necessary design evidence; they are not a numerical validation run.

| Claim | Mathematical/design check | Required computational evidence |
| --- | --- | --- |
| Linear solve converges | Operator, BCs, compatible RHS, null modes and residual norm | Convergence plus true operator residual; appropriate component gauges |
| Continuity holds | Consistent discrete D/G and actual target, including motion | Native final divergence or moving-target error, after constraints |
| Transport conserves | Shared signed integrated transfers; complete control-volume/source definition | Stagewise double-reduced budgets, physical boundary exchange, roundoff-aware tolerance |
| Constants/bounds survive | Compatible transporting fluxes; stated limiter/stability assumptions | Constant and sharp-profile tests in all supported geometries/BCs |
| Energy is physically correct | Viscous dissipation, boundary work and force work; consistent quadrature | Correct decay/input and numerical dissipation, separate from momentum totals |
| Solution is accurate | Consistency and full integration order, including boundary treatment | Native-face field/phase/amplitude/force errors, independent grid and dt refinement |
| S2/S3 coupling is correct | Unique coverage, flux orientation, restriction, composite operators | Partition equivalence, interface budgets and multilevel/reference convergence |
| GPU processing is complete | No host numerical work or field-derived decision after any readback, including scalar/partial reductions and flags | Zero computational device-to-host calls/bytes across direct copies, indirect APIs and pinned staging; bounded terminal outputs counted separately; no hidden fallback |
| Device reductions preserve numerical checks | Double accumulation, complete device finalization, declared ordering/error behavior, proper weights and nullspaces | Independent reference comparisons for norms, coefficients, CFL and diagnostics; unchanged convergence, conservation and physical acceptance gates |
| Stage data is current | Layout/geometry/time-stage validity, producer dependencies and ghost invalidation | Failure/rebind/regrid cases cannot consume stale fields, metrics or ghosts; owned/borrowed lifetimes remain valid |
| Cost is acceptable | Memory/stage/timestep requirements explicitly stated | Idle-GPU cost to declared error target, retained memory and required output cost |

Freeze cases, error measures, tolerances, observation windows and reference/build
identities before comparing methods. Regular periodic scalar evidence precedes
MAC momentum; then physical boundaries, stationary cuts and separately motion.
Use traveling shear, planar Taylor-Green, Beltrami, qualified nonlinear references,
backward step and cylinder according to [TRN-01 acceptance](/docs/amrex-alignment/SINGLE_LEVEL_TRANSPORT#single-grid-acceptance-matrix).
Do not infer full order from an advection order label or one small divergence value.
GPU reductions may change summation order; analyze the resulting error and compare
against independent references while retaining double accumulation. A readback
removal passes only when the complete device computation preserves the scientific
checks, including nonfinite detection, convergence/failure handling and physical
diagnostics. Do not disable checks, relax tolerances or replace reference results
to make a zero-transfer audit pass. Current host calculations remain a migration
gap until their device replacements satisfy those gates.

Keep qualified failures visible. The [exact periodic report](/docs/PERIODIC_EXACT_VALIDATION)
records that the diffusion correction removed the large original Beltrami velocity
plateau, while high-wave traveling shear still fails its accuracy gate. The
[physics guide](/docs/PHYSICS_GUIDE#2-one-time-step) retains the original channel
steady temporal-ratio failure and pending full enclosure qualification. New
transport must not silently relabel those studies as passing or replace their references.

For runtime work, follow [CUDA verification and suite selection](/docs/COMPLETE_VALIDATION#choose-coverage-for-the-change)
and [test instructions](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/test/README.md). Focused operator tests do not replace
the relevant `verify-all` suite or full coverage for broad acceptance. Report actual
receipts and unrun checks; CPU backend not run (user directive). This documentation
change itself has no new runtime/scientific validation result.

## Source map

| Area | Implemented source and focused checks |
| --- | --- |
| Transport | [AdvectionGpu](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/src/PolyCfd.Gpu/Advection/AdvectionGpu.cs); [scheme tests](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/test/PolyCfd.Gpu.Tests/Advection/AdvectionSchemeTests.cs) |
| Stage order | [TimeIntegratorGpu](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/src/PolyCfd.Gpu/Integration/TimeIntegratorGpu.cs); [stage diagnostics tests](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/test/PolyCfd.Gpu.Tests/Diagnostics/StepStageDiagnosticsTests.cs) |
| Discrete projection | [DivGradGpu](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/src/PolyCfd.Gpu/DivGradGpu.cs), [ProjectionGpu](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/src/PolyCfd.Gpu/Projection/ProjectionGpu.cs); [incremental tests](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/test/PolyCfd.Gpu.Tests/Projection/IncrementalProjectionTests.cs), [operator tests](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/test/PolyCfd.Gpu.Tests/Numerics/Pressure/PressureOperatorGpuTests.cs) |
| Pressure iteration and preconditioning | [Shared PcgIteration](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/src/PolyCfd.Gpu/PcgIteration.cs), [Jacobi solver](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/src/PolyCfd.Gpu/PcgSolverGpu.cs), [multigrid solver](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/src/PolyCfd.Gpu/Multigrid/PcgMgSolverGpu.cs); [projection contract and tests](/src/PolyCfd.Gpu/Projection/README), [multigrid contract and tests](/src/PolyCfd.Gpu/Multigrid/README) |
| Motion and gauges | [PressureComponentsGpu](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/src/PolyCfd.Gpu/Projection/PressureComponentsGpu.cs); [component tests](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/test/PolyCfd.Gpu.Tests/Projection/PressureComponentsTests.cs), [geometric source tests](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/test/PolyCfd.Gpu.Tests/Projection/GeometricSourceProjectionTests.cs) |
| Viscosity and walls | [DiffusionGpu](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/src/PolyCfd.Gpu/Diffusion/DiffusionGpu.cs), [immersed coefficients](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/src/PolyCfd.Gpu/Diffusion/ImmersedWallCoefficientsGpu.cs); [independent residual tests](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/test/PolyCfd.Gpu.Tests/Diffusion/DiffusionResidualAccuracyTests.cs) |
| Geometry | [CutCellBuilderGpu](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/src/PolyCfd.Gpu/Geometry/CutCellBuilderGpu.cs), [DeviceCellGeometry](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/src/PolyCfd.Gpu/Core/DeviceCellGeometry.cs); [geometry contract tests](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/test/PolyCfd.Gpu.Tests/Projection/GeometryContractTests.cs) |

The links identify review entry points, not a claim that these tests were run for
this page. Preserve the device ownership, selective capture, completed-frame and
stage-observer contracts in [runtime contracts](/docs/RUNTIME_CONTRACTS).