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 preconditioned conjugate gradients with Jacobi or geometric multigrid preconditioning. Cut-cell weights define the relevant inner product and nullspace handling. Multigrid's coarse correction grids accelerate the linear solve; they are separate from physical AMR levels.

Diffusion solves the backward-Euler equations with weighted 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 pressure multigrid design, cut-cell multigrid operators, current diffusion algorithm and residual definition, and boundary stencils. The multigrid documents label their historical sketches; use their current notes and linked source when reviewing implemented behavior.

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.