Audit and cleanup: 6D CP/CPP pushers (correctness, the #417 loss-vs-nonconvergence bug, mu diagnostic, dead code)

[krystophny]

Summary

Audit of the 6D classical-particle (CP) and classical-Pauli-particle (CPP) orbit pushers and their drivers. Most of the math is correct. There is one real correctness bug (structural non-convergence counted as a confinement loss, the caller-side completion of #417), one misleading magnetic-moment diagnostic that likely explains the reported bad mu conservation for the full particle, model-selection guidance to document, and a cleanup list. This issue collects everything so the cleanup can be done from it directly.

What is correct (verified by reading the code)

• orbit_cpp_boris (modified Boris, the explicit large-step CPP pusher): symmetric macro-step (half drift, half mirror kick, exact magnetic rotation, half mirror kick, half drift); mirror force v -= 0.5*dt*(mu/mass)*gradB with correct sign and scaling, applied only when pauli=.true.; exact Boris/Rodrigues rotation; optional rotation-angle filter for large dt*Omega_c; v_perp = 0 slow-manifold seed; energy 0.5*m|v|^2 + mu|B|. Field push is consistent: B_cart = Jc B^ctr (contravariant), grad|B|_cart = Jc^{-T} d|B|/du (covariant). Frozen mu in the pusher vs a recomputed diagnostic mu is a correct, deliberate split.
• orbit_cpp_canonical (implicit midpoint; CP, CPP_SYM, CPP_VAR): Hamiltonian H = |p - qc A|^2/(2m) + mu|B|, qc = charge/(c*ro0); symplectic-midpoint and variational (discrete Euler-Lagrange) residuals; analytic Jacobian (first-derivative form on spline charts to stay smooth through v_par -> 0); mu a fixed parameter. The integrator correctly flags a non-converged step as ierr=3 with an explicit comment that this is a structural root loss and must not be treated as a physical loss.
• Plain Boris (orbit_full boris_step_cart/boris_step_cyl): correct unmodified Boris (drift-kick-drift, exact rotation, Christoffel half-kicks in the curvilinear variant), no mirror force, as intended for the full charged particle.
• Field provider: dBmod is the covariant gradient d|B|/du; the Cartesian push uses Jc^{-T} (correct for a gradient); the canonical models use dBmod directly in flux space (correct, no Jacobian needed there).

Real bug: CPP-midpoint non-convergence is counted as a confinement loss

• orbit_cpp_canonical flags a non-converged trapped-bounce step as ierr=3, with the explicit intent that this structural root loss is "never silently treated as a loss" (see [orbit-6d-03] Fix CPP full-metric Newton Jacobian convergence #417).
• The driver defeats that intent. In trace_orbit/macrostep, every nonzero ierr_orbit is collapsed to exit; the trajectory terminates early and the particle stops being counted as confined from that step, the same outcome as a genuine wall or boundary loss. There is no branch that distinguishes ierr=3 (structural, recoverable) from ierr=2 (radial boundary) or a real loss.
• Consequence: a trapped particle that merely fails the midpoint solve at a bounce (the [orbit-6d-03] Fix CPP full-metric Newton Jacobian convergence #417 resonance) is spuriously counted as an alpha loss, biasing the confined fraction. This is exactly the failure mode the integrator comment warns about, occurring one level up in the caller.
• Fix options: (a) for production large-step trapped tracing, prefer the explicit modified-Boris full-orbit pusher (ORBIT_CP6D_BORIS), which has no convergence floor and is regular through the axis; and/or (b) make the driver distinguish ierr=3 and recover (reduce dt for that macrostep, retry, or re-seed on the slow manifold) instead of terminating the trajectory.

Magnetic-moment diagnostic for the full particle (the reported bad-mu conservation)

• For the full charged particle, mu = m v_perp^2 / (2|B|) is only an adiabatic invariant. At dt*Omega_c >> 1, plain Boris degrades it; that is expected, not a bug, and the mu-preserving remedy is the modified Boris (freeze mu as a parameter and add the mirror force).
• set_mu_from_state is called only at initialization; state%mu is never refreshed after a Boris step. The carried state%mu is therefore stale (trivially constant), while a live diagnostic recomputes mu from the current state. These two are easy to confuse in post-processing.
• The live recomputation evaluates |B| and the field direction at the start-of-step position, while the leapfrog evaluates the field at the half-drift point. That point/time inconsistency can show up as spurious mu drift even when the orbit itself is accurate.
• Recommended: define a single, consistent mu diagnostic (same |B| and field direction as the push, ideally the gyro-averaged value), and either refresh state%mu after each step or document it as init-only. Then determine which of the three explains the observed degradation: genuine large-dt adiabatic breakdown, sampling inconsistency, or the stale state%mu. The first is cured by the modified Boris; the latter two are code fixes.

Model selection and misunderstandings to document

• The 4D Pauli model (ORBIT_PAULI) is literally the guiding center (its own comment marks it "== GC, tautological"). The genuine 6D Pauli particle is ORBIT_CPP6D.
• ORBIT_CPP6D (flux-canonical implicit midpoint) is singular at the magnetic axis and carries the [orbit-6d-03] Fix CPP full-metric Newton Jacobian convergence #417 trapped-bounce non-convergence. ORBIT_CP6D_BORIS (explicit modified Boris) is regular through the axis and has no convergence floor. For large-step, trapped, or near-axis full-orbit tracing, prefer the Boris pusher.
• The large-step CPP scheme that is theoretically well-founded (modified Boris with the mirror force plus the rotation-angle filter; rigorously analyzed for finite time for the ordinary charged particle by Lubich and Shi, arXiv:2205.12895) is already implemented (orbit_cpp_boris) but marked "Experimental". It should be validated against the midpoint on trapped and near-axis cases and promoted, with the recommended production path documented.

Dead code, small mistakes, and cleanup

• ORBIT_PAULI6D (Cartesian research model) error-stops in both the production validation and the main dispatcher; its constant and code path are effectively unreachable in production. Decide: gate clearly as test-only or remove.
• ORBIT_BORIS is the module-local default of FullOrbitState, yet it is rejected by the production validation ("not available in production alpha-loss tracing"). Harmless, but an inconsistent default; set the default to a production model or document why not.
• set_mu_from_state init-only / stale state%mu (above): fix or document.
• cpp_boris_mu computes |B| two ways (the field-returned Bmod for the denominator and Bmod_of(Bvec) for the field-direction normalization). They agree to round-off by construction; collapse to one.
• Pre-existing TODOs worth folding into the cleanup: orbit_symplectic.f90 ("this doesn't work well yet"; three "move coeff_rk_gauss to preprocessing"), diag_newton.f90 ("fix criterion for convergence"), simple.f90 ("Remove side effects"), simple_main.f90 ("Must be called in this order. TODO: Fix"), params.f90 (deprecated-alias handling).
• Test gap: no test exercises orbit_full's boris_step_cart/boris_step_cyl directly (the CP-Boris tests target orbit_cpp_boris). Add coverage if orbit_full is kept, or consolidate.

Suggested acceptance criteria for the cleanup

• Structural non-convergence (ierr=3) never maps to a confinement loss anywhere in the tracing path; it is recovered or retried.
• A single, documented mu diagnostic with consistent field sampling; state%mu either refreshed each step or clearly marked init-only.
• Orbit models are clearly partitioned into production vs test-only; unreachable models are gated or removed.
• The modified-Boris full-orbit/CPP pusher is validated against the implicit midpoint on trapped and near-axis cases, and the recommended production path is documented.

Related: #417 (the integrator-side root-loss flag this completes on the caller side) and #418.

[krystophny]

Update: the Cartesian decision (#420) supersedes the promote-and-recover options here

The CP/CPP production push moves to Cartesian (orbit_cpp_boris, chartmap field source); the guiding center stays curvilinear. This changes two items above:

• The [orbit-6d-03] Fix CPP full-metric Newton Jacobian convergence #417 ierr=3 recovery options (retry with smaller dt, re-seed) are moot for the dynamics. The Cartesian push has no implicit solve, so there is no bounce non-convergence to recover from. What remains is the caller rule: no solver or field-locate failure counts as a loss; only s(x) >= 1 is a loss.
• "Validate the modified Boris against the midpoint, then promote" becomes "the Cartesian Boris is the production CP/CPP path." The midpoint is deleted, so the comparison gate is Cartesian-vs-GC, not Cartesian-vs-midpoint.

The mu diagnostic item is now concrete in #421 (Larmor-corrected guiding-center reconstruction).

[krystophny]

Field-period-seam over-loss fixed; loss-vs-nonconvergence and dead code addressed

Commit b3e4b59 (branch feature/simple-6d-port) resolves the main remaining
correctness item for ORBIT_CP6D_BORIS and removes the dead code this issue flagged.

The over-loss was a seam-localized loss-detection bug, not transport

On the reactor W7-X high-mirror benchmark (1000 alphas, 1e-5 s) CP lost 147 markers
that GC (0 lost) and ASCOT5 full orbit confine. Every spurious loss fires at a
field-period seam (phi_B = k*2pi/nfp, nfp=5): the reconstructed guiding centre
stays within s in [0.49,0.51] for the whole orbit, then is flagged lost in one step
at mid-radius. Two causes in the wedge inversion (libneo's chartmap is seam-clean):

1. The warm guess u_guess(3) is Boozer phi on the global multi-period sheet; across
   a seam to_wedge rotates the point into the next wedge while the guess stays a
   full period away, stalling the Newton. Fix: seed the toroidal coordinate from the
   in-wedge geometric angle atan2(y,x).
2. accept_or_fail accepted any rho >= rho_edge - 1e-3 as the edge regardless of
   residual, so a Newton stalled at the clamped rho=1 faked a loss from mid-radius.
   Fix: accept a clamped-edge inversion as the edge only when its residual is below
   EDGE_FRAC of a radial cell |dx/drho|, else it is a CPB_LOCATE_FAIL (confined).
   Same criterion now applied to the from_cart fallback, so the warm path and the
   fallback share one loss test.

Result: spurious losses 103 -> 5, CP confined 0.853 -> 0.990 (GC 1.000, ASCOT
1.000), genuine ASCOT-matched losses preserved.

Loss vs non-convergence (the #417 caller bug)

Now consistent: a non-converged inversion is CPB_LOCATE_FAIL -> ierr=3 (counted
cpp_lu_fail, confined), never a loss. The ONLY confinement loss is the
Larmor-corrected guiding centre crossing s>=1 in cpp_boris_to_gc
(CPB_LOSS -> ierr=2, cpp_sbound).

Dead code removed

• The out_of_bounds branch in invert_cart_warm (from_cart never returns
  chartmap_from_cyl_err_out_of_bounds; the constant is declared-but-unused in
  libneo).
• The CPB_LOSS branch after cpp_boris_step in orbit_timestep_cp_boris: the step
  only locates the field (CPB_OK/CPB_LOCATE_FAIL); loss is decided solely in
  cpp_boris_to_gc.
• Stale CPB_LOSS doc comments on invert_cart_warm and cart_field.

Still open

• 5 residual spurious losses at seams (0.5%): a marginal first locate gives a
  slightly-off field, hence a bad Larmor vector and a genuine-looking rho>=1 centre.
  Field-assembly robustness at the seam.
• The cpp_boris_mu double-|B| cleanup and the ORBIT_PAULI6D/default-model items
  from this issue are not yet touched.
• 1 ms / 10 ms validation (acluster) that no genuine losses are now under-counted.

Benchmark, figures, data: itpplasma/benchmark-orbit-proxima#2.

[krystophny]

Full-orbit benchmark extended to 1 s (W7-X high mirror, 1000 alphas, s=0.5, collisionless, reactor scale).

Confined fraction at 1 s:

run	confined
ASCOT5 GC	0.775
SIMPLE GC	0.764
ASCOT5 full orbit	0.714
SIMPLE full orbit (Boris, orbit_model=7)	0.709

SIMPLE full orbit matches ASCOT5 full orbit within the 1000-marker statistics (~0.014). Both codes place the full orbit ~0.06 below their own guiding centre, so the finite-orbit deficit is a real prompt loss, consistent across codes rather than a SIMPLE artifact.

Remaining audit item: 43 markers leave in SIMPLE full orbit before 1 ms that SIMPLE GC, ASCOT GC, and ASCOT full orbit all keep past 1 ms. Most are lost by ASCOT full orbit too, only later, so this is a prompt-loss timing difference, not spurious confinement. Checking by rerunning just those markers at the coarse npoiper2=1000 to separate a finite-Larmor prompt loss from a Boris / chartmap-inversion numerical one.

[krystophny]

Found and fixed a loss-counting bug in the full-orbit prompt-loss subset (PR #427).

The 43 markers that FO loses before 1 ms but GC/ASCOT keep are mostly chartmap field-inversion faults miscounted as losses, not physics. trace_orbit exited on any ierr/=0 and recorded times_lost < trace_time, so an FO_FAULT (ierr=3, inversion non-convergence) counted as a particle loss, against the contract that only the guiding-centre s>=1 crossing (ierr=2) is a loss.

Splitting the terminal event over the 43 markers (W7-X high mirror, 1 ms):

npoiper2	confined before	confined after	fo_loss	fo_fault
100	0.000	0.930	3	40
1000	0.233	0.907	4	29
16384	0.302	0.791	9	21

At production npoiper2=16384, 21 of 43 are faults, only 9 genuine s>=1 losses. diag_traj.x (extended for orbit_model=7) traces them: they fault at interior radius (axis, banana tip s=0.5-0.85), never at s=1, with energy conserved to machine precision. After the fix FO matches GC on these markers (0.79-0.91 vs GC 0.95).

Residual: the orbits still terminate at the fault (inversion robustness), tracked as a follow-up.

[krystophny]

Closing as superseded: this targets the abandoned feature/simple-6d-port CP/CPP pusher (orbit_cpp_*.f90, orbit_model 5/6/8), which was explicitly dropped by #429. None of the referenced files/models exist on main.