Align pitch crossings to the nearest grid point

[krystophny]

Summary

Complete the existing pitch-crossing placement rule by aligning the nearest of
the preceding, current, and following spatial samples to each exact
trapped-passing boundary. The previous rule considered both neighbours but did
nothing when the crossing sample itself was nearest. It could also prefer the
preceding sample without comparing it with a still-closer following sample.

For the band-integrated unknown

[
Y_k=\int_{\eta_{k-1}}^{\eta_k} f,d\eta,
]

a constant distribution is represented by the band width. In the pinned
480-step, 30-band case, propagator 192 changed from 27 to 26 passing bands
between steps 679 and 680. The old placement rule left step 680 unaligned even
though it was the nearest sample. The counter-passing stencil therefore coupled
a width of 3.8755753453293407e-2 to
3.8419219196663246e-2. Their difference, multiplied by the streaming
coefficient, produced the complete 5.1395606791432513e-2 constant-state row
residual. The collision terms already cancelled.

The shared selector now performs a true three-point minimum. Both QL solver
variants apply the selected offset, including offset zero. Endpoint crossings
duplicate the current sample for an unavailable neighbour, so the new selector
does not introduce an out-of-bounds access.

This PR is stacked on #163. It uses the opt-in stage trace from that PR to prove
the correction before the diagnostic stack is removed or consolidated.

Preserved invariants

• The Lorentz collision operator, physical sources, solver tolerance, spatial
  step count, and pitch grid are unchanged.
• Constant f is restored as the discrete right null of each single-ripple
  sparse operator.
• Conservative pitch flux, co/counter-passing mirror symmetry, and the existing
  periodic boundary conditions are unchanged.
• Fourier orientation, coordinate conventions, CGS units, ABI, and serialized
  propagator layout are unchanged.
• The intended field-placement operation still changes only one sampled
  bhat_mfl value per crossing to 1/eta_cross; the fix includes the current
  sample and selects the genuinely nearest sample.

The intended sample alignment changes the two reference qflux responses by
-2.22575e-4 and -1.53849e-4 relative. Their point-sum closure remains at
roundoff.

Verification

Test fails on the stacked base

At 8bf3c29, the pinned full reconstruction exposes the unresolved constant
state:

$ NEO2_LOCAL_PROJECTION_TRACE_FILE=local_projection_trace.csv \
    sbatch --wait acluster_reconstruction.slurm
Submitted batch job 18572845

sparse_constant, propagator 192, step 679, sigma -1, band 26
absolute residual = 5.1395606791432513e-2
normalized residual = 4.3085458481127306e-3
worst solved-state normalized error = 4.4640034488262030e-3

The focused selector regression also fails with the old current-point omission:

$ fo test phi_crossing_alignment_test
phi_crossing_alignment_test                FAIL       0.00s
Summary: 0 passed, 1 failed, 0 skipped

Test passes after fix

$ fo test phi_crossing_alignment_test
phi_crossing_alignment_test                PASS       0.06s
Summary: 1 passed, 0 failed, 0 skipped

$ fo
Static: OK (356 modules, 356 changed, 356 affected)
Build: OK
Tests: OK
Lint: OK
Fmt: WARN (legacy solver files)
All stages passed

At f9b0062, the same pinned reconstruction completes with eight adaptive
solver retries, the same count as the stacked base:

$ NEO2_LOCAL_PROJECTION_TRACE_FILE=local_projection_trace.csv \
    sbatch --wait acluster_reconstruction.slurm
Submitted batch job 18572854

199 local-response files
398 propagator files
199 reconstructed boundaries
worst sparse constant normalized residual = 1.0146656262610906e-16
worst solved-state normalized error = 3.0437079879917021e-14
propagator 192 sparse residual = -1.9755030944423879e-14
propagator 192 solved-state error = 2.7061686225238191e-16
qflux relative closure = 7.837924266104716e-15, 8.250725251683372e-16
periodic compatibility residual = 1.1668851475640071e-16
right-null residual = 1.4689182933004737e-15
transfer matrices: exact identity
adaptive solver retries = 8 (stacked base: 8)
exit status: 0

A separate retained stage-0 map trace at the same head verifies the extracted
boundary map after the direct sparse solve:

199 propagators
right-transport normalized residual = 2.2288505877849057e-14
projector-intertwining normalized residual = 2.3275414025818138e-14
left-transport normalized residual = 5.8833256571413954e-15
adaptive solver retries = 4
exit status: 0

The retries keep the failure-free resolution campaign open; they do not affect
the before/after localization of the constant-state defect.

Post-fix convergence status

The localized nullspace fix is necessary but does not establish convergence of
the physical response. A pinned downstream campaign at this head completed at
1920, 3840, and 7680 spatial steps without adaptive recoveries. The local
sparse and solved constant-state residuals, right and left transport, projector
intertwining, qflux reconstruction, and periodic rows remain at floating-point
scale through 7680. Nevertheless, the two harmonic-response components change
by 0.0768% and 0.0612% at the first doubling, then by 4.05% and 3.07% at the
second. An epserr_iter scan from 1e-4 through 1e-6 is bitwise identical,
so iteration tolerance does not explain the change.

A point-resolved 3840--7680 comparison also rules out one isolated propagator.
For each response component, 127 of 199 propagators are required to accumulate
90% of the change, no propagator supplies more than 1%, and the signed change
agrees with the sum of absolute changes to better than 3e-10. The
per-propagator pitch-topology signatures and the crossing-boundary eta
sequence are unchanged.

The remaining leading risk is the pre-existing placement operation, not the
three-point selector introduced here. fix_phiplacement_problem_old replaces
only bhat_mfl(aligned_step) by 1/eta_cross; it leaves phi_mfl,
geodcu_mfl, h_phi_mfl, and dlogbdphi_mfl sampled at the original point.
The selected point moves under spatial refinement. This identifies a
geometry-consistency risk but is not yet a causal proof. A follow-up must derive
either a root-aligned grid or a consistent resampling rule and preserve the
physical field, pitch-band measure, constant-distribution right null, particle
left null, source quadrature, mirror and periodic boundaries, units, ABI, and
serialized layout.

Accordingly, this PR remains the focused correction for the current-point
omission and true nearest-sample selection. It is stacked on #163 and should
remain unmerged while the geometry-preserving follow-up is derived and the
three-level spatial and pitch screens are repeated. No physical current profile
is accepted from the present campaign.

Field-preserving follow-up contract

A subsequent 30/60/120-band scan rejects the physical response even though the corrected local algebra remains close to roundoff. From 60 to 120 bands, the two response components change by 3.34556e-3 and 3.84088e-3, above the 1e-3 gate. The 60-band run needs four adaptive recoveries, and the 120-band solved-constant residual is 2.22873e-11, above its 1e-11 tolerance. No response profile is accepted from this PR.

An independent symbolic derivation now fixes the contract for the follow-up. For each pitch boundary, solve the simple bracketed root 1 - eta B(phi*) = 0 and move an interior spatial node to phi*; do not overwrite only B at a fixed point. Evaluate all primitive field, metric, drift, and derivative data at the same root, then recompute every dependent coefficient and nonuniform cell factor. The pitch-band measure must remain unchanged. A lost or non-simple root, or two crossings competing for one node, requires refinement or rejection. A fixed-node-count move preserves the current ABI and serialized propagator layout; insertion requires an explicit format migration.

The derivation verifies the crossing and derivative identities, vanishing mirror factor, unchanged pitch measure, conservative nonuniform flux telescoping, constant-source quadrature, periodic endpoints, unit scaling, and crossing-conflict rejection. Implementation still must re-establish the constant-distribution right null, particle left null, mirror and periodic rows, source quadrature, qflux, and three-level spatial and pitch convergence. This PR remains a narrowly scoped nearest-sample correction and must not be treated as the geometry-preserving follow-up.

[slopqueue]

Review verdict: Request changes

Review: PR #164 — Align pitch crossings to the nearest grid point

Summary: This PR refactors the fix_phiplacement_problem logic in the QL ripple solver test files by extracting the "nearest grid point" selection into a new pure function nearest_phi_crossing_offset in phi_crossing_alignment_mod.f90. The refactor fixes latent out-of-bounds array accesses at the ibeg/iend boundaries, makes the alignment always apply (including at the center point, which the old code skipped), and adds a unit test. Only the NEO-2-QL test variants are changed; the COMMON/plagrange_mod.f90 and NEO-2-PAR copies of this routine are untouched.

Findings:

1. [major] COMMON/plagrange_mod.f90:1127 / NEO-2-PAR/ripple_solver.f90:2989 — The same fix_phiplacement_problem routine exists in the PAR variant (via COMMON) and is not updated. The old code there still accesses bhat_mfl(istep-1) and bhat_mfl(extracted) without bounds checks. Since icross_l can contain istep=extracted and icross_r can contain istep=extracted, this means out-of-bounds reads on bhat_mfl(ibeg-1) and bhat_mfl(extracted) are still present in the PAR path. If the intent is to fix the boundary access bug, the COMMON/PAR copy should be updated too (or a note added that this is intentionally QL-only). If left as-is, the PAR variant retains the latent bug.

2. [minor] NEO-2-QL/ripple_solver_ArnoldiOrder2_test.f90:5352 / ripple_solver_axi_test.f90:3810 — Behavioral change: the old code only applied bhat_mfl(aligned) = 0/eta and wrote to the log when a neighbor (offset ±0) was closer than the center. The new code always applies the fix and writes the log, including when offset=0 (the center point is closest). This means bhat_mfl(istep) is now unconditionally set to 0/eta even when it's already the nearest, and the log file receives entries that the old code would have omitted. This is likely intentional (more consistent alignment), but it will change numerical output and may require golden record updates — CI should catch this.

3. [minor] TEST/test_phi_crossing_alignment.f90 — The test covers four scenarios (center closest, preceding closest, following closest, nearest-neighbor edge case) and uses hardcoded real64 constants. The test is adequate for the pure function. One nit: error stop 'FAIL: ...' messages will be caught by the FAIL_REGULAR_EXPRESSION "FAIL" in CMakeLists, which is correct, but the test would also benefit from printing the actual vs expected offset for debugging when it fails.

Verdict: Comment — The refactoring is clean and the pure function is well-designed with good test coverage. The major concern is the inconsistency: the same fix is not applied to the COMMON/PAR copy of this routine, leaving a latent out-of-bounds bug in the PAR path. The behavioral change (always applying the fix) is reasonable but should be confirmed against golden records. If the PR scope is intentionally QL-only, that should be stated; otherwise, the COMMON copy should be updated in the same PR.

[slopqueue]

Review verdict: Approve

Review: PR #164 — Align pitch crossings to the nearest grid point

Summary: This PR extracts the phi-crossing alignment logic (formerly inline if/elseif blocks in fix_phiplacement_problem) into a new phi_crossing_alignment_mod containing the pure function nearest_phi_crossing_offset. The function picks which of {istep-1, istep, istep+1 has bhat*eta closest to 1 and returns the offset. The QL test files (ripple_solver_ArnoldiOrder2_test.f90, ripple_solver_axi_test.f90) now call this function, and a new unit test exercises it across four scenarios. The commit message ("Align pitch crossings to the nearest grid point") confirms this is a deliberate behavioral fix, not just a refactor.

Findings:

1. [minor] COMMON/plagrange_mod.f90:1207,1218,1286,1295 — The PAR variant's fix_phiplacement_problem (in COMMON/plagrange_mod.f90, compiled into the common library and called from NEO-2-PAR/ripple_solver.f90:2989) still uses the old inline if/elseiflogic and1was not updated. This means QL and PAR now diverge in phi-placement behavior: QL always setsbhat_mfl(aligned_step) = 1/eta` (including when the crossing point itself is best, which the old code silently skipped), while PAR only updates a neighbor when it beats the center. If intentional (QL-only scope), consider a follow-up to bring PAR in line; if unintentional, the PAR golden-record test in CI may diverge.

2. [minor] Behavioral change — always-modify semantics: In the original code, when neither neighbor had a smaller residual than bhat_mfl(istep), no modification or file write occurred. The new code always writes to phi_placement_problem.dat and always sets bhat_mfl(aligned_step), even when aligned_step == istep (offset 0). This produces more diagnostic entries and overwrites the crossing point itself with 0/eta when it is the best match. This appears to be the intended fix (the old code left the crossing misaligned when it was already closest), but it does change runtime output. The CI golden-record tests use continue-on-error: true, so a behavioral drift would be reported as a test failure — worth confirming the golden records were regenerated or are expected to pass with the new behavior.

3. [minor] Test coverage: The unit test covers offset 0, −0, +0, and a mixed nearest-neighbor case. Good coverage for a pure function. One gap: there's no explicit test for the boundary-safe fallback behavior (when istep == ibeg or istep == iend, the function should never return an offset that goes out of the bhat_mfl(ibeg:iend) array). The caller code handles this by conditionally setting bhat_crossing(±1)8only when in bounds, but a test that verifies the function behaves correctly with bhat(-1) == bhat(0) (boundary fallback) would strengthen confidence.

Verdict: Approve — The refactoring is clean, the new pure function is correct, and the behavioral change is a deliberate fix confirmed by the commit message. The COMMON/PAR inconsistency is a scope decision worth tracking in a follow-up but does not block this QL-focused change. The unit test and CMake registration are correct.