# Fluid physics in PolyCFD

Reviewed 2026-09-16 against the current solver and validation code; projection
refinement documentation updated 2026-09-18; static incremental pressure correction
and its validation updated 2026-09-20; Jacobi residual compatibility updated
2026-09-21; Chebyshev diffusion updated 2026-09-26; planned transport selection and algorithm review
linked 2026-09-26 (no runtime change). This guide
uses current method names rather than line numbers. Start with
[architecture.html](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/docs/architecture.html) for the system map and
[API.md](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/docs/API.md) for setup examples.

The [GPU execution requirement](/docs/plans/gpu-resident-processing), agreed 2026-09-26,
requires all future simulation processing on the GPU, including complete reductions,
numerical decisions, geometry and diagnostics, without computational D2H copies.
Current scalar/reduction and control paths still need migration; the equations and
current behavior below do not claim compliance. GPU-01 changes execution placement,
not the governing physics, and requires separate numerical/performance evidence.

## 1. Variables, grid, and governing equations

The constant-density incompressible equations, using kinematic pressure `p = P/rho`, are:

```text
du/dt + (u · grad)u = -grad(p) + nu Laplacian(u) + f
div(u) = 0
```

`P` is physical pressure, `rho` density, `nu` kinematic viscosity, and `f` force
per unit mass. The solver's `PInterior` stores **kinematic pressure**. Converting
it to physical pressure requires a density; force diagnostics account for this
normalization explicitly.

The staggered MAC grid stores pressure at cell centres and each velocity component
at its corresponding face. Solid geometry is represented by an SDF and cut-cell
metrics, not a body-fitted mesh. Apertures and volume fractions weight the discrete
operators. For effective fluid volume `V` and geometric face area `A_face`:

```text
D_alpha(u) = sum_faces(outward_sign × alpha_face × A_face × u_face) / V
G(p) = staggered pressure gradient
A = -D_alpha G
```

The implementation masks solid rows and limits small effective volumes. Fractions
are half precision; velocity and pressure fields are float32. Geometry fractions
are dimensionless: physical face area and cell volume must not be applied twice.
See [MAC_GRID_CONVENTIONS.md](/docs/MAC_GRID_CONVENTIONS) and
[halo-layout.md](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/docs/halo-layout.md).

## 2. One time step

`SimulationSession.Advance` drives the integrator and returns a `StepReport` with actual time/dt, pressure and per-axis diffusion results. It does not download full fields. Source terms receive borrowed `SourceStageContext` views; the integrator owns stage boundary application. Scientific output requests explicit fields independently of live preview. See [runtime contracts](/docs/RUNTIME_CONTRACTS).

[TimeIntegratorGpu.Step](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/src/PolyCfd.Gpu/Integration/TimeIntegratorGpu.cs)
uses this order, with boundary operations between stages:

1. Advect velocity by semi-Lagrangian sampling, optionally with MacCormack correction.
2. Add body-force increments `dt f`.
3. For static geometry, add the previous physical pressure force `-dt G(p^n)`.
4. Solve implicit viscous diffusion, including those increments on its RHS.
5. Apply the immersed-face velocity constraints.
6. Solve a pressure increment and correct velocity, then accumulate physical pressure,
   restore its gauge and refresh physical pressure boundary halos.

The current working tree uses **incremental pressure correction for static geometry**.
[Bounded validation] (unavailable local reference: `../../openfoam/compare/INCREMENTAL_PROJECTION.md`) passes analytical
controls, nine mandatory regressions, the full 128³ cavity and the existing timing gate.
The full frozen enclosure qualification is pending; the original channel protocol
retains its failed steady temporal-ratio verdict. Previous saved comparisons retain
their original builds. For static geometry:

```text
u_adv  = Advect(u^n, dt)
u_rhs  = u_adv + dt f - dt G(p^n)
(I - dt nu L) u_star = u_rhs
A phi = -D_alpha(u_star) / dt        # phi starts at zero
u^(n+1) = u_star - dt G(phi)
p^(n+1) = p^n + phi                 # restore the pressure gauge afterward
```

`ProjectionGpu.ApplyPressurePredictor` and `ProjectIncremental` form that pair.
The physical pressure stays in `PInterior` across timesteps, including fresh-integrator
experiments. One retained cell-sized float buffer saves it during the increment solve;
existing gradient buffers/kernels and pressure solvers are reused. There is one primary
Poisson solve. Constant prescribed pressure uses **zero prescribed increment** at the
boundary; physical nonzero pressure is applied to the predictor and final physical halo.
Wall/inlet pressure-Neumann and periodic types remain paired with the existing velocity
masks. Correction BCs borrow the original device patches and allocate no device indices.
The low-level `Project` API still projects with a full pressure field, as required by
standalone operator tests and the existing moving-geometry formulation.

The negative RHS sign matches PCG's negative Laplacian. Search directions and multigrid
error equations remain homogeneous. Accumulated physical pressure is gauged with the
solver's volume weights in every connected region without a fixed-pressure boundary:
closed or periodic static domains and sealed pockets beside an outlet
([section 4](#pcg-and-multigrid)). Moving geometry continues
using the existing full-pressure projection and per-component gauges described in section 5;
incremental pressure history across changing geometry is not claimed here.

The previous non-incremental method applied the entire pressure force after diffusion.
Its ideal steady planar-channel profile contained `dt * G`, where `G = -dp/dx`, with
relative flow excess `2(h/H)^2 + 12 nu dt/H^2`. The static incremental method balances
pressure and diffusion at steady state; the half-cell spatial/quadrature error remains.
The [narrow-gap investigation] (unavailable local reference: `../../openfoam/compare/ENCLOSURE_GAP_INVESTIGATION.md`)
retains the old-build evidence. Converged pressure and small divergence alone do not
establish momentum accuracy. Backward-Euler diffusion remains first order in time.

The uniform validation `cavity` additionally enables `RefinePressureProjection`.
After a successful finite primary increment solve, it recomputes velocity divergence,
solves a residual increment from zero and corrects velocity again before accumulating
both increments into physical pressure. Both solves must converge; the final native
`maxDiv < 1e-5` criterion is unchanged. The option remains restricted to uniform static
all-fluid domains with homogeneous Neumann/periodic pressure boundaries; see
[residual projection](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/docs/PRESSURE_PROJECTION_REFINEMENT.md). The prior validated cavity
refinement/timing matrix used the non-incremental predictor and is preserved. The
new 128³ trajectory passes the same native and cross-code limits; its full grid/time
refinement matrix has not been rerun.

The static projection drives discrete divergence toward zero to the configured
solver tolerances. It does not promise machine-zero divergence or exact
momentum conservation. Moving immersed geometry has a different discrete
continuity target, described in section 5.

Scientific validation runners reject nonfinite or unconverged diffusion in any
velocity component before publishing a completed step. Pressure convergence
is required for every scenario and scientific verification run.
This acceptance guard consumes `StepReport`; it does not change iteration
tolerances or the numerical update. A finite state alone is insufficient.
Its message includes the pressure solver's bounded failure summary (true
residual, worst cells, solution drift) and any regional imbalance; see
[pressure failure summary](/src/PolyCfd.Gpu/Projection/README#pressure-failure-summary).

## 3. Advection, sources, and diffusion

### Advection

[AdvectionGpu](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/src/PolyCfd.Gpu/Advection/AdvectionGpu.cs) backtraces with RK2 and
interpolates staggered velocities. Limited `MacCormack` is the default for the
operator, time integrator, execution options and packaged CLI cases: forward and
reverse transport, a correction, and limiting to the forward interpolation stencil's
range. Scenario JSON requires an explicit `numerics.advection`; exported cases
write `"macCormack"`. Select `"semiLagrangian"` (or `AdvectionScheme.SemiLagrangian`
in the runtime API) for the simpler single-pass method.

MacCormack usually preserves transported features better, as the Gaussian transport
and backward-step comparisons demonstrate. It adds a reverse advection and correction
pass, and retains four face-buffer sets instead of one. Lower dissipation does not
make it best for every problem, timestep or memory budget. Explicit saved scenarios
and scientific regression configurations retain their recorded scheme.

MacCormack's advection correction and residual pressure refinement do not by
themselves establish second-order time accuracy for the complete split integrator.

The semi-Lagrangian remap introduces numerical diffusion and is not a conservative
finite-volume momentum-flux update. Stability at a large timestep does not imply
accurate transport. Adaptive CFL control is an accuracy/cost control as well as
protection against excessive motion between geometry rebuilds.

**Planned, not implemented:** [TRN-01](/docs/amrex-alignment/SINGLE_LEVEL_TRANSPORT)
will evaluate conservative MAC momentum transport on the current single grid before
S2/S3 transport acceptance. It separates momentum budgets from divergence and energy
checks, and requires demonstrated accuracy benefit on preselected workloads plus
measured cost to reach target accuracy. Conservation alone does not establish lower
error or fix splitting, backward-Euler diffusion or cut-boundary approximations.
AMR-02 data alignment preserves the current advection behavior; this numerical change
has its own design and validation gates. The plan retains useful current methods and
adds one qualified conservative option, with one immutable choice per run and explicit
geometry/topology, timestep and predictor/projection requirements. Initial S3 support
is limited to its qualified conservative method; other choices retain their validated
single-level scope. No new method, default or numerical guarantee is delivered here.
The [algorithm review guide](/docs/amrex-alignment/ALGORITHMS) and
[numerical explainer](/algorithms.html) separate current
equations from proposed ones and connect mathematical claims to required evidence.

Geometry-aware sampling moves departure points out of the solid and excludes
blocked faces from the interpolation stencil. It is always used for moving
geometry and is enabled by default for the resolved static cylinder,
channel-obstacle, and cavity-sphere cases.

Opt-in [stage observation](/docs/RUNTIME_CONTRACTS#optional-solver-stage-observation)
can capture native velocity before/after diffusion and projection to diagnose their separate contributions.

### Sources

[ISourceTermGpu](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/src/PolyCfd.Gpu/SourceTerms/ISourceTermGpu.cs) applies momentum
increments before diffusion. `UniformBodyForceGpu.FromPressureGradient` supplies
the pressure-gradient forcing for Poiseuille flow. A constant forcing in a
periodic channel is represented by a body force, not a discontinuous pressure
boundary value. The rotating fan case has no body-force source: wall motion
drives its flow. The stationary fan case is a pressure-driven STL obstacle case.

### Diffusion and walls

The relative diffusion residual is normalized by the **effective** right-hand side, including
prescribed normal velocities, tangential Dirichlet ghost forcing and immersed-wall terms.
A quiescent field driven by a moving lid therefore has a nonzero physical RHS norm. Normalizing
only by the stored velocity RHS falsely rejected the first cavity-sphere step (relative residual
about 1.54e4); independent double tridiagonal reference tests verify the corrected normalization.
The target remains 1e-6. The current bounded fallback policy below supersedes the historical 200/400 caps. No scientific reference
or acceptance threshold was changed to make the scenario migration pass.

`DiffusionGpu` solves a backward-Euler Helmholtz problem for each velocity component
to a relative L2 residual of `1e-6` within 800 total iterations, for static and moving
steps alike (`DiffusionSolveSettings.Default`, recorded in the effective settings): a
primary phase of at most 400 steps of Chebyshev semi-iteration on undamped Jacobi, then
colored Gauss-Seidel on the same row, auditing the stored velocity against the effective
right-hand side in deterministic double reductions. An input already within the tolerance
returns with 0 iterations. The
operator, exits, reporting convention, budget and limits are in the
[Diffusion README](/src/PolyCfd.Gpu/Diffusion/README);
[diffusion convergence](/docs/ADAPTIVE_DIFFUSION_RELATIVE_TOLERANCE) keeps the rationale and the
original enclosure failures. The budget is not a guarantee for every grid/timestep.
Domain-wall ghost relations remain inline, without GPU velocity halo buffers.

For a tangential component half a cell from a Dirichlet wall:

```text
u_ghost = 2 u_wall - u_interior
```

A normal boundary component lies on the wall and is a prescribed row of the diffusion
system (residual p − u), including zero normal flow for slip/symmetry; the audited result
holds it at its value.
Only faces selected by the velocity patch mask are constrained. The final boundary
application does not repair an unconstrained diffusion solve. Zero-gradient
and periodic boundaries have their own stencil/diagonal treatment. See
[BOUNDARY_CONDITIONS_HALOS.md](/docs/BOUNDARY_CONDITIONS_HALOS).

For resolved static immersed walls, `ImmersedWallCoefficients` reconstruct wall
intercepts along velocity stencil links and place the viscous boundary at the
surface. Coefficients are built once for static geometry and reused by the
iterations. The moving path does not yet rebuild these coefficients for each
geometry change; imposed blocked-face wall velocities are a separate mechanism.
Translated-wall, tilted-wall and circular-Couette tests cover the static coefficient
path. Diffusion multigrid is not implemented. The primary phase is Chebyshev
semi-iteration, then colored Gauss-Seidel; its cost still grows with α/h², and at large
α/h² its stall rule can hand a solve to Gauss-Seidel early
([diffusion limits](/src/PolyCfd.Gpu/Diffusion/README#limits)).

## 4. Pressure solve and cut-cell geometry

### PCG and multigrid

[PcgSolverGpu](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/src/PolyCfd.Gpu/PcgSolverGpu.cs) uses a Jacobi preconditioner;
[PcgMgSolverGpu](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/src/PolyCfd.Gpu/Multigrid/PcgMgSolverGpu.cs) uses a geometric
multigrid V-cycle selected with `numerics.pressure.method="pcgMultigrid"`; its level plan, cycle, smoother,
transfers and footprints are specified in the [Multigrid README](/src/PolyCfd.Gpu/Multigrid/README). Both run
one shared PCG loop ([PcgIteration](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/src/PolyCfd.Gpu/PcgIteration.cs)) with their own preconditioner; the solver
interface, the loop, the regions, the results and the resource budget are specified in the
[Projection README](/src/PolyCfd.Gpu/Projection/README). Both apply
the pressure operator without assembling a sparse matrix and use volume-weighted inner products:

```text
<a, b>_V = sum_i(V_i a_i b_i)
```

The weighted operator `A = V^-1 L` is self-adjoint in `<.,.>_V`. It has one
constant null mode per connected fluid region (cells joined through positive face
apertures, and through periodic boundaries when both faces of an axis are periodic)
that has no fluid face on a fixed-pressure boundary. A sealed pocket is such a
region even when the rest of the domain has an outlet; only exactly disconnected
regions count, so a pocket joined through a tiny positive aperture is one region
with a near-null mode. Each solver's PCG loop owns a
[PressureComponentsGpu](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/src/PolyCfd.Gpu/Projection/PressureComponentsGpu.cs),
built from its operator's own apertures, classification and volume weights and
rebuilt when the boundary topology or the geometry revision changes.

The PCG loop subtracts each unanchored region's volume-weighted mean from the
initial guess, the initial residual, the recursive residual **after every update**
and the final solution; anchored regions are untouched. A start-of-solve projection
alone is insufficient: it cannot remove a mean smaller than half an ulp of the float
residual entries, and PCG then lets that constant grow through the preconditioner
until float cancellation corrupts `A p` (the multigrid solver stalled in 1 of 96
128³ all-periodic fan solves at 1e-6, a two-region fine level under one-region coarse
levels kept as a regression test, and on uniform periodic grids at 1e-7 and 1e-8). With
several regions the projection must be per region: projecting only the global mean in
the loop still stalls sealed-pocket solves. When one
unanchored region holds all fluid volume this is the historical global projection
(`GridBlasGpu.ProjectOutMean`), which also shifts zero-volume solid cells; the
per-region path leaves solid and anchored cells unchanged. Solid-cell pressure has
no stencil coupling and carries no meaning; moving projection resets it to zero.
Up to 16 regions are reduced in one fused pass, deterministic by construction (fixed
grid, shared-memory trees, fixed-order totals). A topology with more regions is
classified on the host at every rebuild (every moving step) and its smaller regions
are reduced from compact cell lists, one thread group per region; that adds one
reduction launch per projection and per-rebuild host transfers (see
[resource estimates](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/docs/RESOURCE_PREFLIGHT.md)). The single-region path inherits
`ProjectOutMean`'s ILGPU reduction, whose atomic cross-group combine can change the
last bits of the mean between runs.

The component removed from the initial residual is the right-hand side's regional
incompatibility. A domain without fixed-pressure faces has always tolerated a net
source: it is still projected and reported as `RhsMean` (the documented policy).
With several unanchored regions that tolerated source is attributed to the largest
one (where open inflow/outflow patches are), so a sealed pocket is charged only its
own mean. The weighted norm of what is removed beyond it, relative to the larger of
`||rhs||_V` and the projected initial residual norm, is reported as
`SolveResult.RegionalImbalance` (also on converged results, through `StepReport`).
The solve is refused as not converged with `InconsistentRhs` when that value exceeds
the topology's refusal limit (and the removed norm exceeds the absolute tolerance),
because removing a real net source, such as inflow into a region sealed from the
outlet, would report a clean convergence with a velocity that is not
divergence-free. The limit is `max(1e-3, 5 sqrt(k V_cell / V_fluid))` for `k`
unanchored regions: when uncorrelated per-cell roundoff dominates the right-hand
side (an incremental projection near steady state), each region's removed mean has
a weighted norm of about `h^1.5 sigma` whatever its size, about
`1/sqrt(fluid cells)` of `||rhs||_V`, which exceeds a fixed 1e-3 below about a
million cells. The limit is 0.010 for one pocket at 64³ and 0.0034 at 128³; the
plate test's fed region measures 0.27 against its limit of 0.059. Roundoff that is
coherent along an axis (an exactly extruded flow) can reach `sqrt(N_axis)` times
the white-noise level and is not covered. Detection scales with the region: a region
holding volume fraction `f` whose net source is comparable to the rest of the
right-hand side measures only about `sqrt(f)`, so small pockets with real sources
are projected silently. Moving projection balances its right-hand side per region
before the solve (section 5), so its solves see only roundoff.

Convergence always refers to the recursive residual. Below a relative tolerance of
about 1e-6 the float true residual `rhs - A p` no longer follows it: measured true
residuals are 4e-7 to 2.3e-6 relative after multigrid solves and up to 1.1e-5 after
480-iteration Jacobi solves at a requested 1e-8. A 1e-8 tolerance therefore does
not deliver 1e-8 accuracy (see also the residual gap handled by
[residual projection](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/docs/PRESSURE_PROJECTION_REFINEMENT.md)). The configured tolerances,
iteration caps and breakdown thresholds are unchanged.

`GpuMgHierarchy` rediscretizes cut-cell geometry from the SDF at coarse levels.
Its default `CoarseThickeningMode.ConstantPhysical` keeps the physical solid
offset fixed rather than keeping the same fraction of each increasingly large
cell. Connectivity, weighted V-cycle symmetry, and frozen fan solves have tests;
they do not replace testing evolving geometry. See
[CUT_CELL_MULTIGRID_DESIGN.md](/docs/CUT_CELL_MULTIGRID_DESIGN).

### Boundary patches and resolved bodies

`PatchBcSet` configures domain faces; `DevicePatchBcSet` uploads that configuration.
Common pairings are no-slip velocity/zero-gradient pressure at walls, specified
inflow velocity/zero-gradient pressure, and zero-gradient outflow velocity/fixed
pressure. Periodic pressure and velocity need consistent opposite-face pairs.
GPU pressure halos support operator/gradient boundary stencils and multigrid;
they are not exclusive to multigrid.

`CutCellBuilder` and `CutCellBuilderGpu` use corner samples for face apertures and
a fitted plane for cell volume fractions. Thin features missed by the corners
retain a centre-based fallback and solid thickening. These measures do not make
sub-cell geometry spatially resolved.

`BlockedOnly` projection preserves partial-face fluid flux and constrains blocked
faces. It is the default for the cylinder, channel-obstacle, and cavity-sphere,
along with immersed wall shear and geometry-aware advection. `NoSlip` projection
closes all cut faces, effectively thickening the body; it remains an available
closure for thin features; a more general sub-cell blade closure remains open.
Geometry-aware advection avoids interpolating blocked-face zeros into the fluid;
immersed wall shear supplies the viscous wall interaction explicitly. These options
work together: a resolved no-slip body needs wall shear as well as geometry-aware
transport. Cylinder force accuracy and grid-resolution limits are documented in
[CYLINDER_FORCES.md](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/docs/CYLINDER_FORCES.md).

Cut-cell validation distinguishes aperture/volume-weighted divergence, velocity
on closed faces, and reconstructed surface-normal velocity at the immersed surface.
A zero closed-face velocity checks the imposed constraint; it does not establish
impermeability at a curved surface. The reconstructed normal velocity has interpolation
error and is not known to converge to zero under refinement. The cylinder currently
bounds it at 0.35 times inflow velocity, independently of its force-reference check.
See `SimulationHelpers.ComputeMaxSurfaceNormalVelocity` and the cylinder parameters.

## 5. Moving geometry and pressure components

`StepWithMovingGeometry` rebuilds cut cells at **t + dt** from the rotated cached
SDF, updates the pressure/multigrid geometry, initializes changed cells, and runs
geometry-aware advection, sources, and diffusion. It imposes wall velocity on
blocked faces before projection and reapplies it afterward.

For rigid rotation:

```text
u_wall(x) = omega × (x - centre)
S_raw = D_alpha(u_wall)
b_raw = -(D_alpha(u_star) - S_raw) / dt
```

The moving projection uses the pressure solver's own regions
(`IPcgSolverGpu.PreparePressureRegions`, section 4), rebuilt for the current
geometry revision before the right-hand side is balanced. For each
region without an open fixed-pressure boundary, it removes the volume-weighted mean
of the RHS and resets solid-cell pressure to zero before and after the solve; the
solvers themselves keep every such region's residual compatible and gauge it. One
global mean is insufficient when a sealed fluid pocket is present, as in the fan
hub. Supplying the regions is a required member of every pressure solver.

The saved target is reconstructed from the compatible RHS:

```text
S_compatible = D_alpha(u_star) + dt b_compatible
maxDivError = max |D_alpha(u^(n+1)) - S_compatible|
```

Consequently, raw `maxDiv` can be large beside a moving wall even when the
continuity residual is small. Nonrotating `IMovingGeometry` implementations
retain the limited volume-change source based on old/current cell volumes;
that fallback is not the rotating fan's current source.

The moving runner stops on a failed/nonfinite pressure solve before publishing
that timestep. Pressure gauge handling and the compatible source restore
numerical stability, but do not establish exact momentum conservation across
cell transitions. See [MOVING_GEOMETRY_DESIGN.md](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/docs/MOVING_GEOMETRY_DESIGN.md) and
[MOVING_FAN_VALIDATION.md](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/docs/MOVING_FAN_VALIDATION.md).

## 6. What validation checks

| Case or layer | Evidence and scope |
| --- | --- |
| Taylor–Green | `initial.profile.planar=true`: exact face velocity and energy decay with `u=sin(x)cos(y) exp(-2 nu t)`, `v=-cos(x)sin(y) exp(-2 nu t)`, `w=0`; fitted numerical viscosity. Both advection schemes have mandatory `verify-all` [analytical checks](/docs/TAYLOR_GREEN_VALIDATION). Every pressure solve must converge. The default 3D case needs a numerical reference and refinement for nonlinear accuracy, not an exponential-decay fit. |
| Exact periodic controls | Fully 3D Beltrami and traveling shear compare native-face velocity and energy with exact formulas; shear fluctuation metrics prevent a uniform carrier from hiding errors. Four short MacCormack cases are mandatory. The diffusion correction removed the original large Beltrami velocity plateau; higher-wave shear failures and remaining convergence limits are retained in [the validation guide](/docs/PERIODIC_EXACT_VALIDATION). |
| Poiseuille | Periodic X/Z, no-slip Y, uniform pressure-gradient forcing; comparison with the transient analytical channel profile. |
| Cavity / cavity-sphere | Lid-driven flow; sphere adds immersed-wall and cut-cell diagnostics. |
| Channel-obstacle | Pressure-driven flow around a sphere or box; distinct from the cylinder inflow/outflow case. |
| Static cylinder | Wake statistics and analytic-surface pressure/viscous Cd/Cl. `--check-force-reference` additionally gates canonical Re=100 force accuracy after settled cycles. |
| Backward step | Uniform-inlet laminar 2:1 expansion at Re_h=100. Mandatory MacCormack check in `verify-all`: conservation, settling, reattachment, velocity profiles, wall shear, and pressure recovery against refined OpenFOAM data. See [BACKWARD_STEP_VALIDATION.md](/docs/BACKWARD_STEP_VALIDATION). |
| Moving fan | Actual 128³ MG/adaptive geometry rebuilds, finite pressure, bounded velocity and compatible continuity error. Mandatory short regression in both `dotnet test` and `verify-all`. |
| Performance | Separate local elapsed-time and post-warm-up solver-time references; configurable slowdown failures in `verify-all`. |

[CylinderSurfaceForces](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/src/PolyCfd.Core/Diagnostics/CylinderSurfaceForces.cs)
integrates `-p n + nu (grad u + grad u^T)n` on the analytic cylinder surface,
using fluid-side field reconstruction. Its output is force divided by density;
Cd/Cl use `0.5 U_inf^2 D Lz`. Cartesian apertures are not physical wall quadrature
areas. `PeriodicForceStatistics` compares complete settled cycles to avoid
point-by-point shedding phase sensitivity. This is a static-cylinder diagnostic,
not a generic moving-body force API. See [CYLINDER_FORCES.md](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/docs/CYLINDER_FORCES.md).

Saved-baseline reproducibility, physical reference agreement, numerical stability,
and runtime are separate checks. A local baseline can reproduce a biased solution.
The coarse cylinder can pass geometry checks while failing force targets; the
refined measured case passed the declared force tolerances. Further spatial
convergence and moving-wall viscous validation remain work. Performance workflow
and measurement limits are in [PERFORMANCE_VALIDATION.md](/docs/PERFORMANCE_VALIDATION).

## 7. Source map

| Topic | Main source |
| --- | --- |
| Time integration and adaptive/moving entry points | [TimeIntegratorGpu.cs](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/src/PolyCfd.Gpu/Integration/TimeIntegratorGpu.cs) |
| Projection and compatible moving target | [ProjectionGpu.cs](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/src/PolyCfd.Gpu/Projection/ProjectionGpu.cs) |
| Divergence/gradient | [DivGradGpu.cs](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/src/PolyCfd.Gpu/DivGradGpu.cs) |
| Static wall coefficients | [ImmersedWallCoefficients.cs](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/src/PolyCfd.Core/Geometry/ImmersedWallCoefficients.cs) |
| Shared execution and scientific observation | [ScenarioRunner.cs](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/src/PolyCfd.Execution/ScenarioRunner.cs), [ScenarioCaseRunner.cs](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/validation/PolyCfd.Validation/Framework/ScenarioCaseRunner.cs) |
| CLI dispatch and performance policy | [Program.cs](https://github.com/hankbeasley/polycfd/blob/18102808d7066ed97d8908d66755a68df8d0516e/src/PolyCfd.Cli/ValidationCommands.cs) |

The separate CPU time-stepping implementation was removed. CPU reference
operators and host diagnostics remain in `PolyCfd.Core`; ILGPU can also run the
current solver kernels on a CPU accelerator when CUDA is unavailable.