[orbit-6d-03] Fix CPP full-metric Newton Jacobian convergence

[krystophny]

Stage: orbit tracing / 6D CPP
Source language: Fortran
Manuals to read first: DOC/derivations/cp_cpp_derivation.wl, src/orbit_cpp_canonical.f90, src/simple.f90, src/orbit_full.f90, Xiao and Qin 2021 CPC 265 107981 (10.1016/j.cpc.2021.107981)

Goal

Make orbit_model = 5 use the CPP slow-manifold integrator at GC-sized timesteps for Boozer/VMEC fields. After the PR lands, CPP loss results should not change from prompt losses to full confinement merely by increasing npoiper2, provided the requested step remains inside the documented CPP slow-manifold regime.

Files to edit

• src/orbit_cpp_canonical.f90: fix the full-metric CPP symmetric-midpoint Newton Jacobian and convergence test.
• test/tests/test_cpp_canonical.f90: add a finite-difference Jacobian check for the full-metric CPP residual path.
• test/tests/test_cpp6d_vs_gc.f90: add a behavioral guard that CPP does not accept a step with a large post-update residual in the W7-X Boozer path.

Syntax or behavior to implement

Production CPP is selected by:

orbit_model = 5
orbit_coord = 1

This route goes through orbit_cpp_canonical with MODEL_CPP_SYM. CPP is a canonical Pauli particle on the slow manifold. It should use the implicit midpoint solve for slow guiding-center motion, not the CP gyro-resolved step-size rule.

The full-metric Jacobian currently omits the midpoint-coordinate derivative terms for the CPP force gradient:

d/dqmid_k [ 0.5 m dg_ab,j v_a v_b + qc dA_a,j v_a - mu dB_j ]

Those terms include second derivatives of g, A, and B. Implement them analytically if the field API exposes them consistently, or compute the full-metric Jacobian by finite differencing the same residual block used by Newton. Keep the source of truth in one residual routine.

Newton convergence must check the residual at the accepted state. A small correction alone is not enough when the Jacobian is approximate or numerically degraded.

Scaffold

call sym_residual_blk(st, z_trial, fvec)
call sym_jacobian_blk(st, z_trial, fjac)
call solve_linear(fjac, fvec, dz, ierr)
z_new = z_trial - dz
call sym_residual_blk(st, z_new, fnew)
if (maxval(abs(fnew)) < newton_tol) then
   st%z = z_new
   exit
end if
z_trial = z_new

For the full-metric residual, validate sym_jacobian_blk with central differences:

call sym_residual_blk(st, z + eps*ej, fp)
call sym_residual_blk(st, z - eps*ej, fm)
jfd(:, j) = (fp - fm)/(2*eps)

Positive fixtures to add

• test/tests/test_cpp_canonical.f90: full-metric CPP Jacobian agrees with finite differences for a representative Boozer/VMEC state to a tolerance set by the field interpolation accuracy.
• test/tests/test_cpp6d_vs_gc.f90: one CPP Boozer step at a GC-sized timestep exits with ierr == 0 and a small post-update residual.

Negative fixtures to add

• test/tests/test_cpp6d_vs_gc.f90: force or expose a degraded Jacobian path and assert that the step reports non-convergence instead of accepting a large residual.

Observed failure

A W7-X alpha-loss benchmark at SIMPLE v1.5.0-116-g4cbca07 showed strong CPP sensitivity to npoiper2:

10 us, GC npoiper2=512:       confined 1.000
10 us, CPP npoiper2=512:      confined 0.946
10 us, CPP npoiper2=2048:     confined 0.946
10 us, CPP npoiper2=8192:     confined 1.000
10 us, CPP npoiper2=16384:    confined 1.000
10 us, CP  npoiper2=16384:    confined 1.000
10 ms, GC  npoiper2=1024:     confined 0.829
10 ms, CPP npoiper2=1024:     confined 0.318, first loss at 6.8e-9 s

The trend points to the CPP nonlinear solve. Smaller steps hide the defect. CPP itself is meant to run on the slow manifold; the large-step path should not behave like full CP gyro-resolution.

Success criteria

math -script DOC/derivations/cp_cpp_derivation.wl
fo test test_cpp_canonical
fo test test_cpp6d_vs_gc
fo

A focused W7-X benchmark should also show no prompt CPP loss caused solely by the production timestep choice.

Non-goals

• Do not change CP6D gyro-resolved stepping.
• Do not change CPP initialization or normalization unless a failing test proves a mismatch with DOC/derivations/cp_cpp_derivation.wl.
• Do not replace Newton with fixed-point Picard iteration for MODEL_CPP_SYM.
• Do not weaken current conservation or energy checks.

Verification

math -script DOC/derivations/cp_cpp_derivation.wl && fo test test_cpp_canonical && fo test test_cpp6d_vs_gc && fo

[krystophny]

Related: the CPP slow-manifold formalization (#418, repo gitlab.tugraz.at/plasma/proj/cpp-derivation) has a code-facing certificate lane (lean-cpp-14) intended to export checks that guard the CPP Jacobian and the implicit-midpoint residual in src/orbit_cpp_canonical.f90. Once the Hamilton-equation and midpoint statements are proved, those certificates can back the analytic full-metric Jacobian this issue is about.

[krystophny]

Root cause found: the Jacobian is correct; the gap is the midpoint solve at trapped bounces

Investigated the CPP 10 ms loss gap on the reactor-scaled W7-X high-mirror case (B0 = 9 T, a = 1.30 m, rho* ~ 0.023). Summary: the analytic full-metric Jacobian this issue asks to fix is already complete and correct. It is not the cause of the remaining loss discrepancy.

The Jacobian is complete and correct

The requested second-derivative terms are implemented analytically in jacobian_vmec_analytic:

dgrad_dz(k,j) += 0.5*( 0.5 m d2g_ab,kj v^a v^b + qc d2A_a,kj v^a - mu d2B_kj )

with d2g, d2A, d2Bmod supplied analytically by boozer_field_metric_eval (Boozer pullback through splint_vmec_data_d3 and delthe_delphi_BV_d3).

I checked the analytic Jacobian against a central finite difference of the same residual block at the exact iterates where Newton stalls (not only the single point test_cpp_jacobian_fd covers). Across 30 sampled stalls the maximum relative difference is 9.65e-10 (e.g. entry J(4,2) = 3.92e3, analytic vs FD agree to ~1e-9). Replacing the analytic Jacobian with a full finite-difference Jacobian stalls in the same places. So no Jacobian correction, and no quasi-Newton (Broyden/BFGS) substitute, changes the outcome: at those configs there is no root to find.

Where the loss actually comes from

The implicit-midpoint equation has no reachable root at trapped-orbit bounce points at the CPP step (dtOmega_c >> 1). This is a step-size resonance of the residual O(rho) fast mode (Hairer, Lubich, Wang, Numer. Math. 144 (2020) 787; Ascher and Reich, SIAM J. Sci. Comput. 21 (1999) 1045: implicit midpoint does not pin a slow invariant when fast and slow modes are tightly coupled). The 4D guiding center has no perpendicular mode and does not see it; the 6D classical Pauli particle carries it.

Evidence:

• Large dt is fine for passing orbits. At 100 us, npoiper2 = 512, CPP confines 0.970 against GC 0.972, passing and trapped alike, straight through the start of the loss shoulder. So the macrostep size is not the problem per se.
• The gap is trapped-orbit specific and time-cumulative. Over 10 ms the trapped fraction collapses to ~0 while passing survives; trapped orbits bounce many times and each bounce is a chance to stall.
• Finer dt is worse, not better. npoiper2 = 8192 (16x finer) confines 0.174 against npoiper2 = 1024's 0.498: the spurious losses scale with step count (more bounces sampled), opposite to a gyration-resolution deficit.
• Stalls are genuine non-convergence, not the floor: residual ~1e-3 after maxit, not the ~1e-11..1e-13 round-off floor the converging steps reach.

What does not fix it (tried)

• Globalized Newton (Armijo / monotone line search): the monotone constraint cuts off the non-monotone convergence plain Newton needs, making it worse.
• Relative + absolute residual criterion ||F|| <= max(atol, rtol||F0||): cleaner than the plateau acceptance and reaches >= 1e-12 with the COORD_TOK python oracle preserved, but does not change the trapped 10 ms confinement.
• Local dt bisection (sub-stepping) on stall: reduces stalls (~191 -> ~47 per 300 particles at depth 4) but does not eliminate them and is expensive; the bounce can need dt smaller than dt/16.
• Best-effort commit of the stalled iterate: corrupts the orbit (residual ~1e-3 state) and makes confinement far worse (~0.08).

Path forward

This needs an integrator-level treatment, not a Jacobian or solver change:

1. Resolve the bounce: adaptive sub-stepping keyed on |v_par| approaching a turning point, deep enough that the midpoint root exists, with the substep cost confined to the bounce.
2. Or a structure-preserving scheme designed for the unresolved gyration: Xiao and Qin's VSIP2 / BAP2 on the regular Pauli Lagrangian, or a filtered large-step Boris (Hairer, Lubich, Wang 2020; Lubich and Shi, BIT 63 (2023) 44), which hold the modified magnetic moment near-conserved under under-resolution.
3. Independent of the above: a non-converged step must not auto-count as a lost particle. Only an s >= 1 crossing is a physical loss. The current macrostep maps any integrator error to a loss, so numerical stalls inflate the loss fraction; this should be separated (and is now instrumented via diag_counters events for the CPP/CP paths in local diagnostic builds).

The clean Jacobian + relative-residual convergence work landed on feature/simple-6d-port (commit eeee802); this comment records that the Jacobian is not the open problem and reframes the remaining work around the midpoint-at-bounce resonance.

[krystophny]

The trapped-bounce root loss is a fast-slow midpoint resonance, and the guiding center is provably free of it

Building on the root-cause comment above (the implicit-midpoint equation has no reachable root at trapped bounces at the CPP step), here is the mechanism pinned on a minimal model, a clean statement of why the guiding center does not have it, and where this sits relative to the literature.

Minimal model that reproduces it

A two-degree-of-freedom mirror is enough: H = u^2/2 + (1/(2 eps^2))(B(x) q^2 + p^2), B(x) = 1 + x^2, with u the slow parallel velocity, (q,p) the fast perpendicular gyration of frequency Omega(x) = sqrt(B(x))/eps. The implicit-midpoint Newton step needs 1 - L invertible, L = (dt/2) J^{-1} Hess H.

• Away from the bounce, the fast 2x2 block has determinant 1 + (dt*Omega/2)^2 > 0, so det(1 - L) never vanishes and passing orbits always solve, even at dt*Omega >> 1.
• At the trapped bounce the parallel kinetic energy vanishes, the slow block softens, the O(1/eps^2) cross term B'(x) q / eps^2 feeds the fast block in, and det(1 - L) sweeps through zero on a dt*Omega lattice. There 1 - L is singular and Newton stalls. This is exactly the Ascher-Reich tight fast-slow coupling, localized to the bounce: the midpoint averages the fast oscillation out of the residual value while its frequency stays in the Jacobian.

This reproduces the three observations in the root-cause comment: passing orbits are fine at large dt; the loss is trapped-specific; and finer dt is worse, because more macrosteps land across each bounce, so the spurious stall rate scales like p_window / dt and rises (the opposite of a gyro-resolution deficit). Each bounce draws a fresh gyrophase, so a fixed positive fraction stalls.

Why the guiding center is free of it (and this part is machine-checked)

The reduction averages the gyration into mu (a parameter, conserved) and drops (q,p): there is no Omega in the reduced residual. The reduced Hamiltonian H_GC = u^2/2 + mu B(x) has Hessian diag(mu B'', 1), which is eps-independent and O(1). So 1 - (dt/2) J^{-1} Hess_GC is invertible for dt up to an eps-free bound, and there is no lattice to hit: a unique root exists at and around the bounce for every dt up to that O(1) bound. Resonance is a property of a mode the reduction does not carry.

The eps-free solvability of the reduced midpoint is a formal theorem in the companion Lean development (reduced_solvable / large_step_governed_by_slow_flow in gitlab.tugraz.at/plasma/proj/cpp-derivation, machine-checked, the gyrofrequency absent from every hypothesis and conclusion). The minimal model confirms the GC Hessian is genuinely eps-free O(1), so the hypothesis holds for the magnetized case.

Relative to the literature

Covered: the general fast-slow midpoint resonance (Ascher and Reich 1999); the strong-field resonance condition and the (-1)^n parasitic mode as h*omega -> pi, with filtered fixes that hold a modified magnetic moment near-conserved under under-resolution (Hairer, Lubich, Wang 2020; Lubich and Shi 2023); and the large-step structure-preserving Pauli schemes (Xiao and Qin).

Beyond: those analyses suppress parallel |B| variation (no turning points), so the trapped-bounce locus and its det(1 - L) = 0 characterization are not in them; and the rigorous guiding-center resonance-free statement, with an eps-free admissible step tied to a slow-manifold existence proof, is new.

Ways out, with trade-offs

1. Reduce to / project onto the guiding center: provably resonance-free, at the cost of the full-orbit and FLR information the 6D scheme was added for. This is the only one that is currently a theorem.
2. A structure-preserving large-step scheme for the unresolved gyration (Xiao-Qin VSIP2/BAP2 on the regular Pauli Lagrangian, the filtered Boris / large-step integrators): keeps the modified mu near-conserved so the root survives at dt*Omega >> 1.
3. Adaptive sub-stepping keyed on |v_par| at the bounce, deep enough that the midpoint root exists, with the cost confined to the bounce.
4. Independently of the integrator: a non-converged step must not auto-count as a physical loss (only an s >= 1 crossing is one), which removes the inflated trapped loss fraction directly.

[krystophny]

Update: this is an integrator flaw, not a flaw in the slow-manifold theory, and the way out is a mu-preserving step

Following up on the root-cause and the resonance characterization above. Two objects have to be kept apart.

• The continuous slow-manifold theory. The classical Pauli particle has a slow invariant manifold carrying guiding-center dynamics, it is normally stable, and the magnetic moment mu* is an adiabatic invariant conserved to O(eps^{N+1}) through the slow motion INCLUDING trapped bounces. This is sound. In the companion formalization (gitlab.tugraz.at/plasma/proj/cpp-derivation) it is a machine-checked, axiom-clean theorem (cpp_normal_stability, the bounce-inclusive conservation).
• The discrete integrator. The implicit midpoint at dt*Omega >> 1 does NOT preserve that slow invariant at the bounce: it "does not pin the slow invariant" (Ascher-Reich). The det(1 - L) = 0 root loss is the symptom. The scheme tries to ride the manifold and falls off it exactly where mu* conservation is tightest.

So the failure is localized to the midpoint discretization of the 6D system at a large step through a trapped bounce. The theory, the passing orbits, the small-dt full orbit, and the guiding center are all fine.

The criterion

A large-step CPP integrator is resonance-free (the implicit solve has a root and the orbit stays confined to the slow manifold) precisely when it preserves mu* to O(eps^{N+1}) across the step, bounce included. The plain midpoint violates this at the bounce; the guiding center satisfies it trivially (mu* is a parameter); a filtered, projected, or variational scheme that holds mu* near-conserved inherits the safety. This is the bridge between the continuous bounce-inclusive conservation (proved) and the discrete solvability.

The way out

1. Guiding center: provably resonance-free. The reduced Newton solve has a root for an eps-free dt, machine-checked (reduced_midpoint_resonance_free). The gyrophase one would "lose" is unresolved at dt*Omega >> 1 by definition.
2. GC + finite-Larmor-radius reconstruction: integrate the 4D guiding center at large dt carrying the conserved mu*, and dress each guiding center with its Larmor circle for the perpendicular/FLR information. Large steps AND FLR detail, without ever forming the singular 6D midpoint solve. This is issues FLR loss tubes: dress guiding centers with the Larmor circle at wall checks #415 and Adaptive guiding-center to full-orbit switching near the wall for exact footprints #416; the criterion above says this is the provably-safe realization of what the 6D scheme was for.
3. A mu-preserving large-step 6D scheme (filtered Boris; Xiao-Qin VSIP2/BAP2 on the regular Pauli Lagrangian): no slow invariant to break, no resonance. Sound by the conservation criterion.
4. Independently of the integrator: a non-converged bounce step is structural root loss, not a physical loss, and must not be counted as one. Only an s >= 1 crossing is a loss.

Open edge, and why this is new

The existing mu-preservation proofs for filtered and large-step schemes suppress parallel |B| variation, so they do not cover turning points. Certifying a specific practical scheme preserves mu* THROUGH a trapped bounce is therefore open. The continuous foundation it needs, bounce-inclusive mu* conservation, is the new ingredient and is now formalized; the turning-point resonance locus and the machine-checked GC-safe / midpoint-fails dichotomy are not in the literature. The next step is to formalize the criterion itself, turning "use a mu-preserving scheme" from a hope into a theorem.

[krystophny]

Resolution: retire the implicit midpoint, integrate CPP in Cartesian (see #420)

The premise here was to fix the implicit-midpoint Jacobian and convergence. That route is closed. The Jacobian is correct: analytic matches central differences to about 1e-9 at the stall iterates. The stalls are genuine. At a trapped-orbit bounce (v_par -> 0) and a guiding-center-sized step (dt*Omega_c >> 1), the midpoint map has no reachable root, a step-size resonance of the residual O(rho*) perpendicular mode. No Jacobian, line search, or quasi-Newton variant finds a root that does not exist.

Decision: the CPP integrates in Cartesian with the explicit modified Boris (regular Pauli Lagrangian, frozen mu, mirror force, optional rotation-angle filter). No metric, no Christoffel, no axis singularity, no dynamics solve. The implicit-midpoint CPP is deleted, not fixed. The large-step question becomes the explicit-scheme one (the filtered Boris through the turning point), tracked with the Cartesian path in #420.

The caller-side bug is separate and stands: numerical non-convergence was counted as a confinement loss in macrostep. The Cartesian push removes the bounce non-convergence; the rule is that only s(x) >= 1 is a loss.

[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.