CelerisLab/research/repair-2026-10-10/repair-numerical.md

# Numerical repairs — 2026-10-10

Owner `/root/numerical`. User-authorized repairs, not jet/AMR implementation. Scope: numerical CUDA operator/boundary/step headers, `lbm/stepper.py`, permanent actual-operator regressions. No generated repository config headers, PTX, simulation.py, environment or GPU state modified by this agent. Historical audit scripts/results are retained as evidence, not rewritten to conceal their original failures.

## Reference and design review before edits

Read repository architecture guard. Applicable physics/acceptance: advertised D2Q9/D3Q19 SRT/TRT/MRT local mass/momentum/equilibrium invariants, low-Mach nonuniform closed-channel FP32 reference, same-time storage conversion, moving impermeable wall; no production-Re extrapolation. Code path mapped from `LBMStepper.step` to collision dispatch, force-region postpass, EsoPull sync, curved-wall prepass and in-place step. Historical `b110591` explicitly introduced force_region as an operator-split correction; `6e3756c` introduced initialization-compatible EsoPull semantics. Those intentional choices are retained.

Reference set (local checked-in sources):

1. `ref/lbmpy/src/lbmpy/forcemodels.py`, Guo lines 273–300: source correction is `M^-1 (I-S/2) M Fin`, with half-force equilibrium velocity. A single SRT scalar is not the general MRT/TRT correction.
2. Same file, EDM lines 493–510: equilibrium difference at `u+F/rho` minus equilibrium at `u`. For the repository's quadratic equilibrium, raw Guo source at midpoint `u+F/(2rho)` is exactly this difference. It preserves the existing post-step force-density interface without claiming equivalence to collision-integrated global Guo at all orders.
3. `ref/lbmpy/src/lbmpy/boundaries/boundaryconditions.py`, UBB lines 909–935: fluid-to-wall outgoing direction maps to opposite reflected population with **negative** wall-velocity correction; compressible convention multiplies by density. Also inspected NoSlipLinearBouzidi for interpolation branch context.
4. `ref/FluidX3D/src/kernel.cpp`, load/store around 1333–1343 and streaming section: paired EsoPull layout is a parity-dependent mapping. Same-time conversion must invert the read addresses; advancing store is a different operation. Actual CelerisLab initialization and caller time-count evidence determine which stage is used here.
5. `ref/LBM_boundary_methods_for_stable_high_Re.md` and current source descriptors informed scope: existing west inlet is a ghost-source scheme; only its algebraic shifted-weight restoration is repaired. No general stability claim or boundary-method replacement.

## Dispositions

| Finding | Status | Repair and preserved design |
|---|---|---|
| NUM01 global TRT/MRT force | confirmed / fixed | TRT source split into even/odd rates. MRT adds half source before linear relaxation and half after, exactly implementing the model-specific relaxation matrix. Existing half-force velocity and model rates remain. D2 entire vector checked against independent moment-matrix oracle. |
| NUM02 force-region underforcing | confirmed / fixed | Keep sparse post-step split impulse and public fx/fy force-density semantics. Remove spurious collision scalar. Midpoint source equals exact equilibrium difference, including second moments. Actual ForceRegionKernel regression checks full vector with non-unit rho, nonzero velocity, opposite signed forces and both layouts/parities. Overlapping actual fluid footprints are rejected by state owner's compact-list validation to avoid duplicate writers. |
| NUM03 full D3 entrypoint compile | confirmed / runtime capability rejected; implementation deferred | Architecture/root add `validate_runtime()` rejection for D3 Simulation. Legal schema and local D3 operator tests remain possible. Isolated complete D3 entrypoint still fails 18 dimension-mismatch errors in aux/sync. This round does **not** claim D3 solver support or complete dimension-qualified boundaries. |
| NUM03b shifted D3 MRT | confirmed / local operator fixed | D3 projection uses `g-feq` in the same representation. Removed inconsistent physical-unshift/shift around compile-time shifted equilibrium. Equilibrium/force local regression covers shift on/off. |
| NUM04 same-t EsoPull write | confirmed / fixed | Store uses exactly `esopull_src_cell` and `esopull_slot`, inverse to read at identical count, never the advancing store mapping. Nonuniform exact-address inverse at t=0,1,2,3. |
| NUM05 sync collision/shift | confirmed / fixed | Encode only subtracts weights when needed and encodes at same count. No unsolicited collision of patch/halo. Actual decode/encode preserves arbitrary nonequilibrium populations within FP32 conversion roundoff. |
| NUM06 moving wall | confirmed / fixed | Fluid→solid direction gets negative moving-wall reflection correction. Both paths use local physical density consistent with compressible equilibrium; EsoPull obtains density from pre-streaming donors. q<1/2, q≥1/2 and fallback retained. No omega API sign change. |
| NUM07 EsoPull y-wall stage/race | confirmed / fixed | New read-only `PrepareEsoPullWalls` snapshots 3 donors × 2 rows before any curved/source/in-place mutation. Main step reads compact scratch. t=0 uses initialization-aware helper; later t uses previous spatial postcollision values. Pure double-buffer helpers no longer use timestep as a storage-mode sentinel. |
| Shifted inlet density closure | confirmed / fixed | Both descriptor-weight sums equal 1; replace 5/6 and 19/18. Preserve existing west ghost-source boundary schemes, timing and damping. |
| LES model-family differences | design retained | No unsupported claim that closed-form and iterative LES are interchangeable, no change to LES definition. |

## Phase and interface details

`load_physical_node_esopull(k,t)` reconstructs the next pull state at the completed count. Its spatial meaning is not the same-cell postcollision state reconstructed by `load_physical_dir_pre_streaming`. Decode/encode and snapshot restore retain the existing next-pull interpretation. No clock advancement is hidden in encoding.

`EsoPullStep` now takes `const float* wall_values` immediately before `t`. `LBMStepper` allocates `2*nx*3*sizeof(float)` scratch and launches the prepass on the same stream before `_launch_curved`. This costs one small boundary launch and `24*nx` bytes; no performance-neutrality claim is made. Separating preparation from writes removes the added wall cross-thread read/write hazard; it is not a blanket race proof for all existing inlet/outlet schemes.

Force-region updates remain split and distinct from global collision-integrated Guo. Target momentum per cell per step is F. Macroscopic interpretation with time-dependent forces must respect each scheme's stage; no new global half-force convention or external control semantics is introduced. D2 Guo collision full-vector oracle checks more than momentum rescaling.

## Verification

Permanent tests: `tests/unit/test_numerical_operators.py` and two `tests/unit/cpp/*_regression.cpp` fixtures. All C++ uses actual repository headers/kernels with CPU CUDA-qualifier shims; no copied production arithmetic replaces the invoked operators. Python contains independent D2 matrix/TRT references. Optional nvcc tests copy the kernel tree and edit **only temporary headers**, then emit temporary PTX.

Command:

```text
PYTHONPATH=/home/frank14f/CelerisLab/src /home/frank14f/anaconda3/envs/pycuda_3_10/bin/python -m pytest tests/unit/test_numerical_operators.py -q
```

Coverage:

- 12 local dispatch cases: both dimensions × shift on/off × SRT/TRT/MRT, each at omega .7/1.25/1.8, non-unit density, nonzero velocity and nonequilibrium. Mass/momentum increments, D2 full-vector relaxation, equilibrium invariance and inlet closure.
- 12 semantic cases: D2 models × shift × no-slip/free-slip. Same-time nonuniform encoding, actual no-edit sync, actual wall prepass at initialization and both later parities, moving-wall all directions/branches with rho1.17, actual force-region kernel, and four actual closed-channel steps compared at identical postcollision spatial phase between double-buffer and EsoPull. x periodic, no inlet/outlet donors; this is a deterministic CPU schedule check, not GPU concurrency validation.
- 12 isolated D2 entrypoint compile cases: collision × shift × FP32/FP16S, covering both entrypoints and both wall branches. FP16S is compile coverage only in this agent's work.

Initial 28-test run passed (before expanding semantic cases to all collision models); final expanded result is recorded below. FP32 summation tolerance permits four unit-density epsilons (5e-7) for D3 19-term accumulation; this is separated from per-population vector tolerances. Original bugs exceed tolerances by orders of magnitude.

GPU execution is owned by root, to prevent source/cache/device competition. Add root's observed GPU outcomes separately; no GPU or long-flow claims are inferred from these local tests. High-Re stability, curved-body force/torque convergence, inlet/outlet in-place race audit, LES combinations, complete D3 and performance remain beyond the demonstrated surface.

### Cross-review and GPU handoff

State owner independently accepted same-count inverse addresses, encode-without-collision, and staging before curved/in-place writes on the same stream. I reviewed `field.py` snapshot/restore and `common/checkpoint.py` and reran `tests/unit/test_state_regressions.py`: **21 passed**. Field-only snapshots decode at save count and encode at current count without resetting time. The t=0 initialization special pre-streaming mapping means restoring an initial field at later time promises restored next-pull field, not reproduction of the original initialization step. Full checkpoint restores raw storage plus original count and is the restart mechanism.

Expanded semantic CPU run: **12 passed**, covering all SRT/TRT/MRT, shifted/unshifted, no-slip/free-slip combinations. Together with 12 local-operator cases and 12 isolated compilation cases, all **36** permanent numerical cases have been exercised successfully across the recorded focused runs.

Added `tests/integration/test_numerical_gpu.py` for root's coordinated execution: 18 bounded FP32 cases (4 nonuniform closed-channel layout comparisons over four steps, 4 actual t=0 initialization/staged-wall cases, 6 global-force momentum cases, 4 force-region total impulse cases). Syntax checked; this agent has not run GPU. The layout comparison constructs identical **spatial postcollision** input at completed count 1, rather than incorrectly comparing DB postcollision diagnostics with EsoPull next-pull diagnostics. Actual t=0 initialization is checked separately. No FP16 runtime or quantization tolerance is silently substituted.