CP full orbit performance: path to ASCOT5 speed (phase 'make it fast')

[krystophny]

CP full orbit performance: path to ASCOT5 speed (phase "make it fast")

Tracking issue for the deferred performance work on the Cartesian-chartmap full-orbit
CP/CPP (orbit_cpp_boris). We are in phase "make it run": the full orbit is correct
(reactor W7-X high mirror, 1 ms, 1000 alphas: CP 0.896 vs ASCOT5 full orbit 0.898,
GC 0.905 vs 0.904 -- itpplasma/benchmark-orbit-proxima#2). Do NOT optimise yet; this
issue records how to make it fast later, with the analysis and the trick list.

Measured gap (same hardware, acluster, 32 threads, 1000 orbits)

run	SIMPLE	ASCOT5
GC 1 ms	3.0 s	9.7 s	SIMPLE 3.5x FASTER
GC 100 ms	188 s	655 s	SIMPLE 3.5x FASTER
full orbit 1 ms	409 s	78 s	SIMPLE 5.2x SLOWER

The symplectic GC is faster; the full-orbit CP is ~5x slower. The reason is
architectural, not tuning.

Why the full orbit is ~5x slower

1. SIMPLE inverts, ASCOT reads. ASCOT's field is on a real-space (R,phi,z) grid:
   B(x) is one direct tricubic spline read (B_field_eval_B), no inversion. SIMPLE's
   field is in logical Boozer/chartmap coordinates while the particle is pushed in
   Cartesian, so every B(x) requires solving the Cartesian->logical inverse (a damped
   Newton on the chartmap forward map, several spline evals) and then B = curl A.

2. Three inversions per step vs zero. orbit_timestep_cp_boris does, every step:
   cpp_boris_step -> 1 cart_field (push), and cpp_boris_to_gc -> 2 locate
   calls (orbit_cpp_boris.f90:566,577): one inverts the particle, one inverts the
   Larmor-corrected guiding centre. The 2 extra inversions exist because SIMPLE decides
   the loss on the reconstructed GC and writes the GC (s,theta,phi) every substep.
   ASCOT gets its loss criterion (marker rho>=1) for free from the field read.

3. Finer step. SIMPLE ~38 steps/gyroperiod (npoiper2=16384) vs ASCOT 20
   (gyrotime/20): ~1.9x more steps. Both are 2nd-order Boris.

Net: ~3 inversions + 3 field assemblies/step (SIMPLE) vs ~2 direct reads (ASCOT),
x ~1.9 steps, partly offset because ASCOT's 3D tricubic read is itself heavier per
eval than SIMPLE's chartmap spline -> ~5x.

What ASCOT trades for speed (risks of its coarser interpolation)

ASCOT is fast partly by approximations that SIMPLE's curl-A chartmap avoids. The risks
to keep in mind if we copy them:

• Not field-aligned. A (R,phi,z) grid is not aligned to flux surfaces, so it needs
  many points to resolve the |B| mirror structure and steep edge gradients; tricubic
  error is O(h^4) for B and O(h^2) for grad-B. SIMPLE's chartmap splines in flux
  coordinates, where the field varies slowly along the surface, so fewer DOF reach the
  same accuracy. A coarse (R,phi,z) grid mis-resolves trapping near the edge/X-point.
• div B != 0. ASCOT splines B_R,B_phi,B_z independently, so the interpolated field
  is not exactly divergence-free; SIMPLE's B = curl A is divergence-free by construction.
  Non-zero div B can drive spurious parallel acceleration / energy drift over long
  (1 s) traces. ASCOT mitigates with fine grids.
• Coarser step. 20 vs 38 steps/gyroperiod: Boris phase error ~ (dt*Omega)^2 grows
  faster; over a slowing-down time this is a spurious-radial-transport risk. ASCOT's
  volume-preserving VPA bounds the growth; a plain Boris at 20/gyroperiod is riskier.

At ASCOT's validated grid + step the fields agree with our chartmap to 0.01% in |B|
and 0.01 deg in direction (benchmark orbit_cmp), so these are margins, not errors --
but they are the price of the speed, and they shrink at coarser settings.

What VENUS-LEVIS does

Pfefferle et al., CPC 185 (2014) 3127. Full-orbit (and GC) with a Fourier
reconstruction in the poloidal/toroidal angles + cubic spline in the radial flux
coordinate, transformed to Cartesian (PSLINE), RK4. Fourier-in-angle is exact and
periodic (no field-period seam -- cf. our seam bug), flux-aligned (accurate per DOF),
and the interpolation is built to preserve the RK4 integrator order for fast-ion
accuracy. Lesson: a flux-aligned Fourier-spline field is both accurate and seam-free;
the speed comes from an interpolation tuned to the integrator, not from a coarse
real-space grid.

Plan to reach ASCOT speed (ordered by payoff; for "make it fast" phase)

1. GC reconstruction off the substep loop. Do cpp_boris_to_gc only at the
   output/loss-check cadence (or only when the particle is within ~a few Larmor radii
   of the edge), not every substep. Removes 2 of the 3 inversions/step -- the single
   biggest win (~2/3 of per-step cost). Needs a cheap near-edge proxy (the push already
   locates the particle, so its rho is in hand for the gate).
2. Share the push and to_gc inversions. The push's cart_field and the first
   locate in to_gc invert nearly-adjacent points; warm-seed one from the other or
   reuse the result.
3. Precompute a direct Cartesian field (ASCOT-style) from the chartmap. Build, once
   at setup, a spline of B_cart(x) (and grad|B|, and the loss level set s(x)) on a
   Cartesian or (R,phi,z) grid by evaluating the chartmap curl-A field at build time.
   The hot loop then reads it directly -- no per-step inversion -- while the chartmap
   stays the exact divergence-free source. This is the big architectural lever and
   keeps our accuracy advantage.
4. Carry the logical coordinate as integration state (ASCOT/GORILLA pattern):
   advance u alongside x and source B by the forward map + phi-reduction, eliminating
   the inverse in the common case (fall back to inversion only on resync).
5. Fewer Newton iterations. Better warm start (predict u from the previous step's
   velocity), tighter accept, drop the redundant covariant_basis where evaluate_cart
   already returns the Jacobian.
6. Vectorise over markers (SIMD). ASCOT pushes NSIMD markers per step; our loop is
   one marker per thread. A marker-batched inner loop helps cache + vectorisation.
7. Coarsen to ~20 steps/gyroperiod to match ASCOT once accuracy is locked (1 ms
   already matches at 38; the finer step buys little).
8. Relativistic / volume-preserving push for long-time error control at the coarser
   step (matches ASCOT's VPA), removing a reason to keep 38/gyroperiod.

Items 1-2 are local to orbit_cpp_boris/simple.f90 and should recover a large chunk
of the 5x with no physics change; item 3 is the architectural step to actually reach
ASCOT speed; 6-8 close the remainder.

Acceptance (later)

Full-orbit CP within ~1.5-2x of ASCOT5 wall-time at the same loss accuracy
(CP confined within statistics of ASCOT), on the reactor W7-X benchmark.

Related: #419, #420, #421; benchmark itpplasma/benchmark-orbit-proxima#2.

[krystophny]

Lever 3 (direct Cartesian field, no per-step inversion) is now a detailed design with code-grounded feasibility in #423 -- reuses the libneo batch splines (cubic+quintic, periodic-phi, GPU acc routine seq + resident kernel), samples the exact curl-A field at setup, and recommends splining the vector potential A and taking curl of the spline to keep div B = 0 (better than ASCOT). It also enables a GPU full-orbit push (the chartmap inversion is GPU-hostile). libneo spline asks in itpplasma/libneo#326.

[krystophny]

Benchmark performance plot cleanup landed in itpplasma/benchmark-orbit-proxima@d49d8bf.

run_simple/w7x_high_mirror/cartesian_cp_test/reactor_performance.png is now the canonical performance comparison. It shows available SIMPLE GC timings and ASCOT GC/FO timings, and explicitly marks fresh SIMPLE FO timing as pending until the current 100 ms / 1 s FO reruns finish.

[krystophny]

The W7-X reactor SIMPLE FO 1 s split completed and is recorded in itpplasma/benchmark-orbit-proxima@0671d2b.

Collected result: 1000 markers, confined fraction 0.773. The production comparison now uses SIMPLE GC/FO and ASCOT GC/FO at 1 s.

This does not close the FO performance question: the split chunk logs lack the final timing line, so the performance plot omits the SIMPLE FO 1 s wall-time bar rather than estimating it from partial progress output.

[krystophny]

Automated triage pass (Claude), grounded in the repository at d8b0e01. Verify before acting.

Still relevant, but the issue text points at a file that no longer exists on main. The
architecture it describes does — verbatim — under a new name. Refining, with the naming fixed
and lever 1 turned into a buildable spec.

Naming correction first. src/orbit_cpp_boris.f90 has never existed on main
(git log main -- src/orbit_cpp_boris.f90 is empty); it lives only on
origin/feature/simple-6d-port (and origin/z_archive/feature/simple-6d-port) and was never
merged. 123d26a ("Extract shared GC Gauss implicit-RK core (orbit_rk_core); drop CPP and
symplectic FO", #429, ancestor of main) dropped the CPP model codes from the shared integrator
core, and 36d97ef ("Add full-orbit (FO) Boris pusher (orbit_model=7)", #426, ancestor of main)
landed the production Cartesian full-orbit Boris pusher as src/orbit_fo_boris.f90
(orbit_model = ORBIT_FULL_ORBIT = 7, src/params.f90:74). Every structural claim in the issue
survives the rename: the three-inversions-per-step count is exactly what main does today (line
references below). Please retitle to orbit_fo_boris so the issue is greppable.

Purpose

The full orbit is SIMPLE's own cross-check on the guiding-centre loss fraction, and the org's
benchmark against ASCOT5 (itpplasma/benchmark-orbit-proxima). At ~5x slower than ASCOT5 it is not
usable as a routine confirmation run: the thread above records that the W7-X reactor 1 s FO run had
to be split into chunks, and the split chunk logs lost the final timing line, so the canonical
performance plot still has no SIMPLE FO 1 s bar. Making the FO push cheap turns "we ran it once,
carefully" into "we run it alongside every GC campaign". Lever 1 alone is a pure book-keeping change
with no physics content — the field, the integrator and the loss criterion are untouched — and the
issue's own analysis puts it at ~2/3 of the per-step cost.

Current state

The push and the guiding-centre reconstruction each pay a full Cartesian→logical damped Newton
inversion, and the reconstruction pays two.

src/simple_main.f90:1644 calls orbit_timestep_fo once per microstep (do ktau = 1, ntau_local,
src/simple_main.f90:1628). This is the only call site in the tree. src/simple.f90:216 does, per
microstep:

    call fo_step(fo, status)                          ! src/simple.f90:228
    ...
    call fo_to_gc(fo, s, th, ph, vpar, status)        ! src/simple.f90:240

• fo_step (src/orbit_fo_boris.f90:444) → cart_field (:455) → invert_cart_warm (:330)
  — inversion 1.
• fo_to_gc (:497) → locate(st%x, st%u, …) (:509) → invert_cart_warm (:348)
  — inversion 2;
  then locate(x_gc, u_p, …) (:520) — inversion 3.

Each invert_cart_warm (:73) runs newton_from (:122) with maxit = 30, maxls = 30
backtracking, and one ref_coords%evaluate_cart + one ref_coords%covariant_basis per iteration,
plus a field_at_logical (:266) afterwards that calls covariant_basis again at the converged
point (:279) and fo_eval_field (src/orbit_fo_field.f90:23 → field_can_mod::evaluate).

Two facts make the gate safe, and both are verifiable in the current code:

1. The particle's own rho is already free. fo_step writes it: st%u = u at
   src/orbit_fo_boris.f90:457, where u comes back from the push inversion. st%u(1) is the
   gate variable; no extra work. (Caveat for the reviewer: cart_field is called at the
   half-drift point x + 0.5*dt*v (:454-455), so st%u is the midpoint chart position, not the
   end-of-step one. The offset is O(dt·v) and is far inside the 0.05 gate margin below, but it is
   why the gate threshold must not be pushed up against 1 - GC_PARTICLE_GAP.)
2. z is not consumed every microstep on the FO path. In macrostep
   (src/simple_main.f90:1598, FO branch :1644-1676) the FO branch uses z only for
   locate_linear_lcfs when ierr_orbit == ORBIT_FO_LOSS (:1649), and for the
   z = z_step_start rollback in warning mode (:1660). Wall checking and collisions are
   error stop-ed for ORBIT_FULL_ORBIT at src/simple_main.f90:1168-1173, and classification is
   rejected at src/simple_main.f90:274-275. trace_orbit stores orbit_traj(:, it) = z only at
   macrostep granularity (src/simple_main.f90:1248).

grad|B| is dead code in this module. field_at_logical computes the pseudo-Cartesian
d|B|/dw chain rule and a full inv3(Jw, Jwinv) at :308-314 on every assembly, and no caller
ever reads the result: fo_step passes gradB into cart_field (:455) and the Boris rotation at
:459-464 is Bvec-only; locate (:341) receives it as the local gw (:345,353) and never
uses it — bhat is built from Bw alone (:355-357). The same is true of the Jc
returned to cart_field and locate. Be honest about the size of this: it is one 3x3 inverse and
about a dozen flops per assembly, against a Newton solve with several spline evaluations. It is
worth removing because it is free and unambiguous, not because it moves the 5x.

Reproduction / evidence

The FO fixture is the reactor-scale Landreman-Paul QA VMEC file, downloaded by CMake at configure
time to test/test_data/wout.nc (test/tests/CMakeLists.txt:17-28; it is not in a fresh checkout)
and symlinked into the test binary dir (:37). test/tests/test_fo_boris.f90:37 initialises
exactly from it.

Baseline timing run. simple.x takes an optional config path (src/simple_main.f90:135-137), but
netcdffile and all output files resolve against the working directory, so run from inside the
case directory:

make CONFIG=Fast          # configure step downloads test/test_data/wout.nc
mkdir -p run_fo
ln -sf ../test/test_data/wout.nc run_fo/wout.nc
cat > run_fo/simple.in <<'EOF'
&config
  netcdffile = 'wout.nc'
  isw_field_type = 2
  orbit_model = 7
  orbit_coord = 1
  integmode = 1
  npoiper2 = 16384
  ntestpart = 32
  ntimstep = 100
  trace_time = 1d-4
  sbeg = 0.6d0
  startmode = 1
  deterministic = .True.
  relerr = 1d-13
/
EOF
cd run_fo && time ../build/simple.x

(isw_field_type = 2 is BOOZER, src/magfie.f90:43; it is a deprecated alias for integ_coords
and is normalised at src/params.f90:612-617. netcdffile back-fills both field_input and
coord_input, src/params.f90:594-604, so init_ref_coords gets the chartmap.)

Observed today: three invert_cart_warm calls per accepted microstep (static count from the call
sites above; make it dynamic with the counter in step 1 below). Expected after lever 1: one
inversion per microstep for interior markers, three only within the gate band, and a wall-time
reduction of ~2.5-3x on this case with confined_fraction.dat and times_lost.dat unchanged to the
tolerances in the acceptance list.

Per-step attribution, to confirm the inversion is the cost and not field_can::evaluate:

cd run_fo && perf record -g ../build/simple.x && perf report --stdio | head -40

Expect newton_from, evaluate_cart and covariant_basis to dominate.

Behaviour (BDD)

Scenario: interior marker pays one inversion per microstep
  Given a marker started at sbeg = 0.6 with orbit_model = 7 and npoiper2 = 16384
  And the gate threshold fo_gc_gate_rho = 0.90
  When one microstep of orbit_timestep_fo runs with fo%u(1) < 0.90
  Then the EVT_FO_INVERT counter increases by exactly 1
  And z(1:5) is left at its previous value (not refreshed this microstep)

Scenario: near-edge marker keeps per-microstep loss detection
  Given a marker whose push inversion returns fo%u(1) >= 0.90
  When orbit_timestep_fo runs
  Then fo_to_gc is called that microstep and EVT_FO_INVERT increases by 3
  And a guiding-centre crossing u_gc(1) >= 1 still returns ierr = ORBIT_FO_LOSS on the
      same microstep it occurs, so locate_linear_lcfs receives a z_step_start one
      microstep old

Scenario: macrostep boundary always reports an exact guiding centre
  Given ntau_local microsteps in macrostep it
  When ktau == ntau_local completes
  Then fo_to_gc has run on that microstep and orbit_traj(:, it) holds the
      Larmor-corrected guiding-centre (s, theta_B, phi_B), identical to the ungated
      build to within 1e-12 relative

Scenario: no change in the physics answer
  Given test/test_data/wout.nc, ntestpart = 32, trace_time = 1d-4, deterministic = .True.
  When the gated and ungated builds are both run
  Then confined_fraction.dat agrees exactly (same integer counts)
  And times_lost.dat agrees to 1e-10 relative per marker

Scenario: the field assembly no longer computes an unused grad|B|
  Given fo_step and both locate call sites pass want_grad = .false.
  When field_at_logical runs
  Then it skips lines 308-314 (the d|B|/dw chain rule and inv3(Jw, Jwinv)) and returns
      gradB = 0
  And Bvec, Bmod and bhat are bit-identical to the want_grad = .true. path
  And the energy drift gate in test_fo_boris (|dE/E0| < 1e-3) still passes

Scenario: a misconfigured gate is rejected at startup
  Given fo_gc_gate_rho = 0.97 and GC_PARTICLE_GAP = 0.05
  When validate_orbit_model_config runs with orbit_model = 7
  Then SIMPLE error stops with a message naming both fo_gc_gate_rho and GC_PARTICLE_GAP

Implementation plan

This is the lever-1/2 PR only. Lever 3 is already specified in #423; lever 5 and levers 6-8 need
their own issues (see Scope and risk).

1. src/diag_counters.f90 — add EVT_FO_INVERT = 15 and EVT_FO_GC_EVAL = 16 following the
   existing pattern at :25-42: bump N_EVENT from 14 to 16 (:42), export the two new parameters
   in the public list at :16-21, add their case labels in event_name (:91-118). Note the
   hard constraint documented at :44-46: STRIDE = 16 and STRIDE >= N_EVENT. 16 is exactly at
   the limit — a third new counter requires bumping STRIDE to 32 as well, or count_event writes
   into the next thread's column.
2. src/orbit_fo_boris.f90 — call count_event(EVT_FO_INVERT) once at the top of
   invert_cart_warm (:73). Keep it out of newton_from (:122) so the multi-seed from_cart
   fallback (:82-88) is not double-counted.
3. src/orbit_fo_boris.f90 — add an optional logical, intent(in) :: want_grad to
   field_at_logical (:266), cart_field (:323) and locate (:341), defaulting .true..
   When .false., field_at_logical skips :308-314 and returns gradB = 0; cart_field then
   skips the gradB = rotz(gw, ca, sa) at :336. Pass want_grad = .false. from fo_step
   (:455) and from all three locate call sites (:401, :429/:436, :509, :520) —
   none of them reads gradB, and locate's bhat is built from Bw only (:355-357),
   so this is safe at every site. Do not remove the gradB dummy argument; keep the API for a
   future mirror-force or drift diagnostic.
4. src/params.f90 — add real(dp) :: fo_gc_gate_rho = 0.90_dp and
   integer :: fo_gc_gate_stride = 0 next to the other FO parameters, and add both to the &config
   namelist (:164-183). fo_gc_gate_stride = 0 means "guiding centre only at macrostep
   boundaries and inside the gate band"; = 1 restores today's every-microstep behaviour and is the
   reference for the regression comparison.
5. src/simple.f90 — give orbit_timestep_fo (:216) an added logical, intent(in) :: force_gc
   argument. There is exactly one caller (src/simple_main.f90:1645), so no other call site needs
   touching. After the successful fo_step at :228, evaluate the gate:
   need_gc = force_gc .or. fo%u(1) >= fo_gc_gate_rho .or. fo_gc_gate_stride == 1. Call fo_to_gc
   (:240) and write back z(1:5) (:242-246, :257-259) only when need_gc; otherwise set
   ierr = 0 and leave z untouched. Count EVT_FO_GC_EVAL when it fires. The ORBIT_FO_LOSS and
   ORBIT_FO_NUMERICAL handling at :241-256 is unchanged, as is the fo = fo_start; z = z_start
   rollback at :234-235.
6. src/orbit_fo_boris.f90 + src/simple_main.f90 — GC_PARTICLE_GAP is currently a private
   parameter (src/orbit_fo_boris.f90:44; the module's public :: list is :48 and :50), so the
   startup assertion below will not compile until it is added to the public list. Then, in
   macrostep (:1644-1645) pass force_gc = (ktau == ntau_local) so the macrostep boundary always
   yields an exact z, and keep the z_step_start = z semantics at :1630: when the gate did not
   fire, z_step_start is the last evaluated guiding centre, which is what locate_linear_lcfs
   (:1649) must interpolate from. Because the gate band is entered before any loss can occur (the
   loss test at src/orbit_fo_boris.f90:536 needs u_gc(1) >= 1 and u_p(1) >= 1 - GC_PARTICLE_GAP), z_step_start at a loss is always one microstep old, exactly as today. Assert
   this: fo_gc_gate_rho must be <= 1 - GC_PARTICLE_GAP; error stop in
   validate_orbit_model_config (src/simple_main.f90:267-280, alongside the existing
   orbit_coord and classification checks) otherwise.
7. test/tests/test_fo_gc_gate.f90 (new) — model it on test/tests/test_fo_boris.f90:1-64
   (same wout.nc init via init_field at :37 and get_boozer_coordinates at :39,
   ans_s = ans_tp = amultharm = 5 at :28, npoiper2 = 16384 via init_params at :40). Trace
   one interior marker (s = 0.5, lambda = 0.9, as at :51) and one near-edge marker
   (s = 0.97, lambda = 0.2) for a fixed step count, once with fo_gc_gate_stride = 1 and once
   with 0, and assert: identical (s, theta, phi) at macrostep boundaries to 1e-12;
   diag_counters_total(EVT_FO_INVERT) reduced by at least 2x for the interior marker; identical
   loss step for the near-edge marker. Call diag_counters_init / diag_counters_reset between
   the two runs.
8. test/tests/CMakeLists.txt — register it beside the existing FO block at :371-380:
   add_executable(test_fo_gc_gate.x test_fo_gc_gate.f90),
   target_link_libraries(test_fo_gc_gate.x simple),
   add_test(NAME test_fo_gc_gate COMMAND test_fo_gc_gate.x WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}),
   set_tests_properties(test_fo_gc_gate PROPERTIES LABELS "integration" TIMEOUT 120). It must sit
   inside the same wout.nc-guarded if() block, i.e. before the endif() at :381.
9. DOC/config.md — add a bullet (the file is a 43-line bulleted parameter guide) documenting
   fo_gc_gate_rho and fo_gc_gate_stride, including the fo_gc_gate_stride = 1 reference mode
   and the <= 1 - GC_PARTICLE_GAP constraint.

No libneo change is required for this PR. Lever 3 (#423) does need libneo first — see
itpplasma/libneo#326 — so the ordering is: this PR (independent) → libneo#326 → #423.

Acceptance criteria

• make CONFIG=Fast builds clean.
• make test TEST=test_fo_boris passes unchanged (energy gate |dE/E0| < 1e-3, near-axis
  crossing, no spurious loss).
• make test TEST=test_fo_gc_gate passes.
• make test passes with no new failures.
• make test-regression passes: no golden-record case moves. All 13 cases under
  test/golden_record/*/simple.in leave orbit_model unset, i.e. at the ORBIT_GC = 0 default
  (src/params.f90:73-75), so this must be a no-op — a diff there means the change leaked
  outside the FO branch.
• On the run_fo case above, confined_fraction.dat is byte-identical between
  fo_gc_gate_stride = 0 and = 1, and times_lost.dat agrees to 1e-10 relative.
• Reported EVT_FO_INVERT per accepted microstep drops from ~3 to <1.3 averaged over the
  32-marker run.
• Wall time of the run_fo case drops by at least 2x on the same machine and thread count.
• fo_gc_gate_rho > 1 - GC_PARTICLE_GAP is rejected at startup with a message naming both
  parameters.
• With want_grad = .false. wired in, Bvec/Bmod outputs are bit-identical to the
  want_grad = .true. path (assert this directly in test_fo_gc_gate rather than inferring it
  from the energy gate).

References

• Boris, J. P. (1970), Relativistic plasma simulation — optimization of a hybrid code, Proc. 4th
  Conf. Num. Sim. Plasmas, NRL — the pusher in orbit_fo_boris.
• Qin, H., Zhang, S., Xiao, J., Liu, J., Sun, Y., Tang, W. M. (2013), Why is Boris algorithm so
  good?, Phys. Plasmas 20, 084503 — the phase-error/long-time argument behind the
  "38 vs 20 steps per gyroperiod" trade in the issue.
• He, Y., Sun, Y., Liu, J., Qin, H. (2015), Volume-preserving algorithms for charged particle
  dynamics, J. Comput. Phys. 281, 135 — the volume-preserving family relevant to lever 8, not
  to this PR. (The issue body attributes a VPA integrator to ASCOT5; that attribution is the issue
  author's and is not independently checked here.)
• Varje, J. et al. (2019), High-performance orbit-following code ASCOT5 for Monte Carlo simulations
  in fusion plasmas, arXiv:1908.02482 — the reference for the direct real-space field read and for
  the NSIMD marker batching described in the issue body (levers 3 and 6).
• Pfefferle, D., Cooper, W. A., Graves, J. P., Misev, C. (2014), VENUS-LEVIS and its spline-Fourier
  interpolation of 3D toroidal magnetic field representation for guiding-centre and full-orbit
  simulations of charged energetic particles, Comput. Phys. Commun. 185, 3127 — the
  flux-aligned Fourier-spline alternative discussed in the issue body. The same author's
  arXiv:1412.5464 is already cited in src/orbit_fo_boris.f90:543 for the pseudo-Cartesian
  axis chart.
• Albert, C. G., Kasilov, S. V., Kernbichler, W. (2020), Symplectic integration with non-canonical
  quadrature for guiding-center orbits in magnetic confinement fusion, J. Comput. Phys. 403,
  109065 — SIMPLE's GC path, the one the issue's table puts at 3.5x faster than ASCOT5.

The gate itself, the counters and the want_grad flag have no specific literature reference; they
are implementation details internal to this code.

Scope and risk

Steps 1-9 are one PR, roughly size/M, tier T2: no physics changes, but it touches the
loss-detection cadence, so the locate_linear_lcfs argument in step 6 is the thing a reviewer must
check hardest. The main risk is a marker that reaches u_gc(1) >= 1 without ever tripping
fo_gc_gate_rho on the particle radius — excluded for fo_gc_gate_rho <= 1 - GC_PARTICLE_GAP by
the loss test at src/orbit_fo_boris.f90:536, which is why step 6 makes that an enforced startup
constraint rather than a convention. The second-order risk is the half-step offset noted above:
fo%u is the chart position at x + 0.5*dt*v, not at the end of the step, so the effective gate is
a fraction of a microstep early — which is the conservative direction.

Lever 5 is not in this PR, and the issue's framing of it does not survive contact with the code.
The issue asks to "drop the redundant covariant_basis where evaluate_cart already returns the
Jacobian". ref_coords%evaluate_cart(u, x) returns only the position, so there is nothing to drop
locally. And the obvious in-repo version — hand newton_from's Jc out and reuse it in
field_at_logical — is wrong as stated: newton_from computes covariant_basis(u, Jc) at
:144 before the line search updates u = ut at :182, so on return Jc belongs to the
iterate before the accepted one, and on the rn < tol early return at :140-142 (which can fire on
entry) Jc is never assigned at all. Restructuring so Jc is evaluated at the accepted iterate
costs one extra covariant_basis to save one, i.e. nothing.

The real lever-5 win is one level down, in libneo: chartmap_covariant_basis
(src/coordinates/libneo_coordinates_chartmap_impl.inc) calls
chartmap_eval_cart_der(self, u, vals, dvals), which already computes the position vals
alongside the derivative dvals and then discards vals. Exposing a fused
evaluate_cart_with_basis(u, x, e_cov) on coordinate_system_t would make each Newton iteration
one spline evaluation instead of two (newton_from evaluates evaluate_cart at ut in the line
search at :170 and then covariant_basis at the same point at :144 on the next iteration), and
would remove field_at_logical's repeat call at :279. That is a libneo ask in the shape of
itpplasma/libneo#326 and deserves its own issue, not a line in this PR.

Do not touch the phi seed in invert_warm_newton (:109-115) as part of any "better warm
start": the deliberate use of atan2(x(2), x(1)) instead of the carried Boozer phi is the
field-period-seam fix documented at :105-108, and reinstating the stale phi will re-open the
spurious-loss bug.

Split recommendation. Keep #422 as the tracking issue, but move the remaining levers out so each
is independently implementable: lever 3 already lives in #423 (blocked on itpplasma/libneo#326);
lever 5 (fused evaluate_cart + covariant_basis in libneo, per the analysis above), lever 6
(marker-batched SIMD inner loop), lever 7 (coarsen to ~20 steps/gyroperiod, needs an accuracy study
against test_fo_boris's energy gate at 1 ms and 100 ms), and lever 8 (volume-preserving push)
each want their own issue. Levers 1 and 2 are the content of this spec.
Related: #420 (the Cartesian-integration decision, whose "keep orbit_cpp_boris"
list is also stale), #423, itpplasma/libneo#326.