# Algorithms and numerical review

**Review companion, 2026-09-26.** 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.

## 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/main/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/main/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](https://github.com/hankbeasley/polycfd/blob/main/docs/amrex-alignment/x)
x_dep = x - dt I_u[u^n](https://github.com/hankbeasley/polycfd/blob/main/docs/amrex-alignment/x_mid)
S_dt[u^n](https://github.com/hankbeasley/polycfd/blob/main/docs/amrex-alignment/q)(x) = I_c[q](https://github.com/hankbeasley/polycfd/blob/main/docs/amrex-alignment/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](https://github.com/hankbeasley/polycfd/blob/main/docs/amrex-alignment/q^n)`. `MacCormack` uses:

```text
q_forward  = S_dt[u^n](https://github.com/hankbeasley/polycfd/blob/main/docs/amrex-alignment/q^n)
q_backward = S_-dt[u^n](https://github.com/hankbeasley/polycfd/blob/main/docs/amrex-alignment/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/main/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.

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. Connected static
solves use their weighted gauge; the moving path explicitly handles disconnected
components. New topology must preserve the appropriate component compatibility
and gauge, not replace several null modes by one global mean.

**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/main/docs/MOVING_GEOMETRY_DESIGN.md) and
[moving fan qualification](https://github.com/hankbeasley/polycfd/blob/main/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 linear-solver algorithms behind these equations, continue to the
[pressure multigrid design](/docs/GPU_MULTIGRID_DESIGN),
[cut-cell multigrid operators](/docs/CUT_CELL_MULTIGRID_DESIGN),
[current diffusion iterations and residual definition](/docs/ADAPTIVE_DIFFUSION_RELATIVE_TOLERANCE)
and [boundary stencils](/docs/BOUNDARY_CONDITIONS_HALOS). The multigrid documents
retain labeled historical sketches; their current notes and linked source take
precedence. 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.

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.

**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.

## 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 |
| 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.

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/main/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/main/src/PolyCfd.Gpu/Advection/AdvectionGpu.cs); [scheme tests](https://github.com/hankbeasley/polycfd/blob/main/test/PolyCfd.Gpu.Tests/Advection/AdvectionSchemeTests.cs) |
| Stage order | [TimeIntegratorGpu](https://github.com/hankbeasley/polycfd/blob/main/src/PolyCfd.Gpu/Integration/TimeIntegratorGpu.cs); [stage diagnostics tests](https://github.com/hankbeasley/polycfd/blob/main/test/PolyCfd.Gpu.Tests/Diagnostics/StepStageDiagnosticsTests.cs) |
| Discrete projection | [DivGradGpu](https://github.com/hankbeasley/polycfd/blob/main/src/PolyCfd.Gpu/DivGradGpu.cs), [ProjectionGpu](https://github.com/hankbeasley/polycfd/blob/main/src/PolyCfd.Gpu/Projection/ProjectionGpu.cs); [incremental tests](https://github.com/hankbeasley/polycfd/blob/main/test/PolyCfd.Gpu.Tests/Projection/IncrementalProjectionTests.cs), [operator tests](https://github.com/hankbeasley/polycfd/blob/main/test/PolyCfd.Gpu.Tests/Numerics/Pressure/PressureOperatorGpuTests.cs) |
| Motion and gauges | [PressureComponentsGpu](https://github.com/hankbeasley/polycfd/blob/main/src/PolyCfd.Gpu/Projection/PressureComponentsGpu.cs); [component tests](https://github.com/hankbeasley/polycfd/blob/main/test/PolyCfd.Gpu.Tests/Projection/PressureComponentsTests.cs), [geometric source tests](https://github.com/hankbeasley/polycfd/blob/main/test/PolyCfd.Gpu.Tests/Projection/GeometricSourceProjectionTests.cs) |
| Viscosity and walls | [DiffusionGpu](https://github.com/hankbeasley/polycfd/blob/main/src/PolyCfd.Gpu/Diffusion/DiffusionGpu.cs), [immersed coefficients](https://github.com/hankbeasley/polycfd/blob/main/src/PolyCfd.Gpu/Diffusion/ImmersedWallCoefficientsGpu.cs); [independent residual tests](https://github.com/hankbeasley/polycfd/blob/main/test/PolyCfd.Gpu.Tests/Diffusion/DiffusionResidualAccuracyTests.cs) |
| Geometry | [CutCellBuilderGpu](https://github.com/hankbeasley/polycfd/blob/main/src/PolyCfd.Gpu/Geometry/CutCellBuilderGpu.cs), [DeviceCellGeometry](https://github.com/hankbeasley/polycfd/blob/main/src/PolyCfd.Gpu/Core/DeviceCellGeometry.cs); [geometry contract tests](https://github.com/hankbeasley/polycfd/blob/main/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).