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:

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.