From physical laws
to discrete updates.

∂u/∂t + (u · ∇)u = −∇p + ν∇²u + f
∇ · u = 0
p = physical pressure / ρ

Here u is velocity [m/s], ρ is constant density [kg/m³], p is kinematic pressure [m²/s²], ν is kinematic viscosity [m²/s], and f is acceleration [m/s²]. Every term in the momentum equation has acceleration units.

∇ · (u ⊗ u) = (u · ∇)u + u(∇ · u)

The two transport forms agree for a continuously divergence-free field. Discrete interpolation, divergence and flux operators need their own compatibility argument; the continuum identity does not prove that the implemented updates agree.

A cell-centered scalar uses a cell control volume. U momentum uses a control volume around an X face, shifted by half a cell; V and W have their own volumes. On a uniform periodic grid their regular volumes equal the cell volume, but their locations differ. Boundary and cut volumes need explicit definitions.

u_adv = Advect(uⁿ, Δt)
u_rhs = u_adv + Δt f − Δt G(pⁿ)
(I − Δt ν L)u_diff = u_rhs
u* = ApplyImmersedConstraints(u_diff)
Aφ = −D(u*) / Δt,    A = −DG
uⁿ⁺¹ = u* − Δt Gφ
pⁿ⁺¹ = pⁿ + φ

D is the implemented divergence, G the effective correction gradient including solid-face masks and correction boundary treatment, L the viscous operator, and φ the pressure increment. These equations abbreviate boundary forcing and imposed face constraints. The pressure operator must match the correction actually applied, including skipped solid-adjacent faces.

D(uⁿ⁺¹) = D(u*) − Δt DGφ
          = D(u*) + Δt Aφ ≈ 0

This cancellation requires the same compatible D and G in the solve and correction. Remaining divergence depends on solver error, boundary treatment and geometry. Closed/periodic connected components need a compatible right-hand side and a pressure gauge; a gauge fixes pressure's arbitrary constant, not velocity error.

Backward-Euler diffusion is first order in time. A higher-order transport step alone does not make this complete split method second order. Optional residual projection has its own restricted scope. Moving geometry uses a separate full-pressure formulation and a geometry-dependent continuity target.

The pressure equation uses one shared preconditioned conjugate-gradient loop with Jacobi or geometric multigrid preconditioning. It handles the weighted compatibility and gauge of each unanchored pressure region, including static disconnected domains. Multigrid's coarse correction grids accelerate the linear solve; they are separate from physical AMR levels.

Diffusion solves the backward-Euler equations with Chebyshev semi-iteration on undamped Jacobi and a bounded colored Gauss–Seidel fallback. Their stopping criteria assess the discrete equation residual, including boundary forcing. A converged algebraic solve does not establish that the spatial or temporal discretization is accurate.

Continue into the implemented shared pressure iteration, multigrid preconditioner, diffusion algorithm and residual definition, and boundary stencils. The pressure multigrid design, cut-cell multigrid design and diffusion history keep the rationale and labeled historical algorithms.

x_mid = x − (Δt/2) I[uⁿ](x)
x_dep = x − Δt I[uⁿ](x_mid)
F = S_Δt[uⁿ](qⁿ) = I[qⁿ](x_dep)
B = S_−Δt[uⁿ](F)
q_adv = clamp(F + ½(qⁿ − B), q_min, q_max)

q is one face-velocity component; I interpolates its staggered samples, and S is the backtrace/interpolation operation. Semi-Lagrangian returns F. MacCormack uses the final expression. Both passes use the original carrying velocity uⁿ, including the reverse pass; the reverse transported field is F.

The bounds are the forward interpolation stencil's valid minimum and maximum. Geometry-aware sampling adjusts midpoint and departure locations out of solids and excludes blocked samples. Those modifications and limiting affect accuracy; the regular-grid formula alone does not qualify cut cells.

There is no shared momentum flux subtracted from one volume and added to the next. For an interpolation matrix W, preserving a constant means W1 = 1; preserving a weighted total would also need vᵀW = vᵀ, where v contains control-volume weights. The first condition alone does not imply the second. Limiting makes the overall operation nonlinear as well.

Q_K = V_K q̄_K
T_f = ∫[tⁿ,tⁿ⁺¹] ∫[face f] F · n_f dA dt
Q_Kⁿ⁺¹ = Q_Kⁿ − Σ_f s_Kf T_f + S_K

For this fixed-volume balance, q̄ is a volume average, Q an amount, and S the integrated source. Choose one orientation n_f for each shared face. The incidence sign s is +1 when it points outward from K, and −1 when inward. T includes area and time exactly once.

For scalar advection, F = q u. In three dimensions, T has units [q] m³, matching Vq. For density-normalized momentum q = u_i, Q has units m⁴/s; multiply by constant density for physical momentum [kg m/s]. Pressure and viscous stress are additional transfers, not included in this advective flux.

ΔQ_A = −T_f,    ΔQ_B = +T_f
Δ(Q_A + Q_B) = 0

Floating-point accumulation still needs a declared, scale-aware tolerance. At physical boundaries only the interior contribution exists, so the domain total changes by boundary exchange and sources.

Example only: a fixed 1D regular grid, constant speed a, periodic boundaries, cell-average scalar q, spacing h. This is not the selected 3D momentum algorithm.

σ_i = minmod(q_i − q_(i−1), q_(i+1) − q_i)
q_L = q_i + σ_i/2
q_R = q_(i+1) − σ_(i+1)/2
F_(i+½) = max(a,0) q_L + min(a,0) q_R
L(q)_i = −(F_(i+½) − F_(i−½)) / h
q¹ = qⁿ + Δt L(qⁿ)
qⁿ⁺¹ = ½qⁿ + ½[q¹ + Δt L(q¹)]

Minmod returns the smaller magnitude with the common sign, or zero when signs differ. Reconstruct fluxes again at q¹. The two-stage SSP Runge–Kutta update uses the average of both stage fluxes for its integrated transfer. A cautious sufficient bound for this specific constant-speed, minmod example is |a|Δt/h ≤ 1/2.

Second-order behavior concerns smooth resolved regions; limiting reduces order near extrema and sharp profiles. Neither this bound nor this reconstruction can be copied unchanged into multidimensional cut-cell MAC momentum. Derive and test that method's own conditions.

1. Derive each component's volume, face area and location, including physical boundaries.
2. Define reconstruction of momentum and normal carrying velocity at every transport face and time stage.
3. Show that the carrying fluxes obey the required continuity relation on those volumes; a small primal-cell divergence alone does not prove this.
4. Specify any predicted velocity or additional projection, and how it couples to diffusion, pressure history and sources.
5. Derive momentum and kinetic-energy balances separately. Conservative momentum can still have excessive energy dissipation or energy growth.

Only then choose the limiter, temporal scheme, ghost width and stable timestep. The existing outer workflow can be shared where mathematically valid; extra numerical stages must be explicit.

V_fluid = κ V_cell
A_open = α A_face
Δq̄ = (−Σ outward transfers + sources) / V_fluid

κ and α are dimensionless fractions. Small V_fluid can impose a restrictive explicit timestep. Redistribution or merging is a candidate treatment, but its conservation, boundedness and wall behavior need derivation and tests. Current geometry weights and volume floors cannot automatically become the new physical momentum measure.

d/dt ∫[Ω(t)] q dV
  + ∫[∂Ω(t)] q (u − w) · n dA = ∫[Ω(t)] s dV

This moving-volume scalar identity uses boundary velocity w and source density s. For momentum, pressure and viscous tractions must also be included. Fluid volume change, swept transfers, newly exposed samples and wall work must agree. Treat current and previous geometry distinctly; a static derivation does not qualify motion.

S2: fill the chosen method's stencil at every needed stage, then use one agreed transfer on the shared interface with opposite signs. Interpolation methods instead need valid departure-point sampling across patch boundaries. Demonstrate partition equivalence for each enabled method.

q̄_C = (Σ_fine V_f q̄_f) / V_C
V_C = Σ_fine V_f

This average-down equation applies to compatible nested scalar volumes. Staggered fields need placement-aware restriction. Covered coarse cells are excluded from composite physical totals. Inconsistent cut volumes invalidate the simple identity.

Orientation: uncovered coarse volume → fine region
coarse update used:  −T_C
fine update received: +T_F
coarse correction:  +(T_C − T_F)

T_F sums fine transfers over the same physical face and time interval. After correction the coarse update is −T_F, so the interface contributions cancel. Averaging covered coarse values and refluxing the adjacent uncovered coarse volume do different jobs.

Synchronized timesteps do not remove the spatial mismatch. Refluxing also does not enforce the incompressibility constraint: compatible transfers, composite projection, diffusion and geometry remain separate S3 work. The equations above illustrate the accounting; MAC momentum interfaces still need their specific construction.

GPU fields → GPU partial reduction → GPU final result
           → GPU coefficient / decision → GPU update

Requested output: GPU final product → external I/O

A PCG dot product must remain on-device through its final double-precision sum and coefficient calculation. A CFL reduction must feed a device-computed admissible timestep. Component labeling, pressure gauges, geometry rebuilding, diagnostics and future regridding follow the same rule. No scalar or setup-only readback exemption is implied.

Host code may schedule launches using configuration metadata. GPU-resident flags govern numerical convergence and dependent work. A completed-operation status may be published, but polling numerical partials on the host to decide further iterations is prohibited. Explicit final measurements and captures never feed CPU-computed values back into the solver.

Current gaps: weighted reduction partials, the shared PCG recurrence, diffusion weights and audits, and CFL host arithmetic still involve CPU numerical work. Component mean subtraction now completes on-device, but topology rebuild decisions and imbalance measurement still read back data; larger topologies also use host classification. The formulas elsewhere on this page describe current numerics, not compliance with the new execution gate.