CP/CPP full orbit: source field from the existing libneo neo_spline_field (splined A, B=curl A, direct read, GPU-ready)

[krystophny]

CP/CPP full orbit: optional direct real-space field spline (no per-step inversion, GPU-ready)

A concrete realization of the #422 lever 3 ("precompute a direct Cartesian field").
Add an ALTERNATIVE field representation for the Cartesian full-orbit pusher
(orbit_cpp_boris) that reads B from a real-space (R, phi, z) spline instead of
inverting Cartesian->logical every step. Keep the chartmap as the exact build-time
source. The current B = curl A-from-chartmap path stays the default and the accuracy
reference; this is opt-in and benchmarked against it.

Why

• Speed (CP full orbit performance: path to ASCOT5 speed (phase 'make it fast') #422): the full orbit is ~5x slower than ASCOT5 almost entirely because
  every field eval inverts Cartesian->logical (a damped Newton, several spline evals),
  plus the GC reconstruction does 2 more inversions/step. ASCOT reads three direct
  tricubic splines on an (R,phi,z) grid -- no inversion.
• GPU: confirmed asymmetry in the code. libneo's batch-spline EVALUATE is already
  GPU-callable (!$acc routine seq at ~11 sites in batch_interpolate_3d.f90, a
  resident device kernel, managed-memory bug fixed). The chartmap inversion is
  GPU-HOSTILE: data-dependent Newton trip counts + line search (orbit_cpp_boris.f90
  newton_from maxit=30), select type on the hot path, procedure pointers and
  class-dispatched metric/christoffel -- orbit_cpp_chartmap_metric.f90:16 says
  outright "NOT GPU-portable". So a direct spline read is the thing that makes a GPU
  full-orbit push possible at all.

What to build (reuses most of the substrate)

1. Setup-time sampler. Loop a regular (R, phi, z) grid over one field period; at each
   node invert once to logical and call the SAME field_at_logical the hot loop uses
   today, storing the Cartesian field. The inversion is paid once per node at setup,
   never per step -- that is the whole win. The sampler is essentially today's
   cart_field run over a grid.
2. Real-space field type. Clone src/field/field_splined.f90 (which today splines in
   FLUX coords, so it still needs the inversion) onto an (R, phi, z) grid, built with
   construct_batch_splines_3d (libneo), periodic = [.false., .true., .false.] (phi
   periodic over 2pi/nfp, R/z natural -- inherits the seam-clean periodic-phi support
   directly; use cylindrical, NOT Cartesian, so phi is an axis).
3. Pusher branch. cart_field and the loss test in orbit_cpp_boris get a fast path
   that reads the real-space spline instead of to_wedge -> invert_cart_warm -> field_at_logical. Gate as a new orbit_model or a field-representation flag; the
   chartmap stays the source.

What to store, two options (benchmark both):

• (A) spline B_cart directly (ASCOT-style): batch of (B_x, B_y, B_z, s) and read
  grad|B| via evaluate_batch_splines_3d_der. Simplest; but an independently splined
  B is NOT divergence-free (ASCOT lives with this), an energy-drift risk over long traces.
• (B) spline the vector potential A_cart, take B = curl of the spline (via the
  analytic spline derivatives _der). This keeps div B = 0 exactly -- better than
  ASCOT -- while staying a direct, inversion-free, GPU-friendly read. grad|B| then needs
  A's 2nd derivatives (_der2, available). Recommended primary; A is smooth so it
  splines cleanly.

Spline the radial label s(x) too, so the loss test (u_gc>=1 today) becomes a spline
read (GPU-friendly), exactly as ASCOT splines psi.

Orders and the LCFS kink

libneo batch splines already support cubic and quintic (order in [3,5],
batch_interpolate_3d.f90:96); a quintic real-space field is a one-line change. The open
problem is the separatrix kink: a rectangular (R,phi,z) box has corners with rho>1
where the chartmap is undefined and field_at_logical returns the edge-clamped field, so
a grid line crossing the LCFS sees a C1 kink. libneo has ONLY natural/periodic BCs --
no monotone/Steffen/limited/one-sided variant -- and a quintic rings worse at the kink
than a cubic. Mitigations, by effort:

• gate on the splined s(x) so no used cell touches rho>1 (ASCOT approach, cleanest);
• keep the box just inside the LCFS;
• a kink-robust spline variant in libneo (one-sided/limited near the boundary) -- a
  separate libneo task (see the linked libneo issue). The user is fine allowing
  cubic + quintic now and adding a kink-robust variant later.

Risks (the depth)

• Memory is the first constraint. libneo's batch layout stores a dense (order+1)^3
  coefficient block per cell: a 4-quantity QUINTIC at 200x64x200 is ~190 GB
  (prohibitive); cubic ~36 GB; ASCOT's compact form (one set per node, reconstruct on
  the fly) is far lighter. So: cubic + a modest grid + store only (B_x,B_y,B_z,s) (grad
  via _der) + per-field-period storage; and/or add a compact-storage option to libneo
  (linked issue). Size the grid in the design before coding.
• div B != 0 for option (A); option (B) (curl of splined A) avoids it -- prefer (B).
• LCFS-edge accuracy (the kink) lands exactly where losses are decided; benchmark it.
• phi-periodicity: use (R,phi,z), not Cartesian, to keep periodic-phi.

Benchmark plan ("how bad does it get")

On the reactor W7-X high-mirror case (benchmark-orbit-proxima#2), compare the real-space
spline CP against the exact chartmap-inversion CP and ASCOT5, at 1 ms / 100 ms:

• speed (target: within ~1.5-2x of ASCOT, vs ~5x today);
• confined fraction / loss times (target: within statistics of the chartmap-inversion CP);
• energy/mu conservation over the trace (the div-B / kink test);
• cubic vs quintic, option (A) vs (B), grid resolution sweep -> the accuracy-vs-cost curve.

Reuse vs new

Reuse: the whole batch-spline engine (cubic+quintic, periodic-phi, value/der/der2, GPU
acc routine seq + resident kernel, managed-memory fix); field_splined.f90 as the
template; field_at_logical as the per-node sampler; the loss-as-level-set pattern.
New: the setup sampler, the real-space field type, the pusher branch, and (optional,
libneo) a kink-robust BC + compact storage.

Effort: medium. No physics change -- the exact curl-A field is sampled. Related:
#422 (perf umbrella), #419/#420/#421; libneo spline-infra issue (linked);
benchmark itpplasma/benchmark-orbit-proxima#2.

[krystophny]

libneo-side spline-infrastructure asks (kink-robust BC variant + compact coefficient storage for the real-space field grid) tracked in itpplasma/libneo#326.

[krystophny]

Correction (do not invent a new representation -- libneo already has it). neo_field_mesh + neo_spline_field (libneo/src/field/field_mesh.f90, spline_field.f90) already implement exactly this:

• field_mesh_init_with_field(limits, field, n_points, is_periodic) samples ANY field_t's A and B onto a regular 3D grid (frame-agnostic: (R,phi,z) periodic-phi OR Cartesian via is_periodic+limits; a degenerate axis gives 2D/EQDSK);
• spline_field_init splines A1,A2,A3 AND B1,B2,B3;
• runtime compute_afield/compute_bfield are direct reads (no inversion), and compute_afield_derivatives gives dA/dx -> B = curl A (divergence-free).
  Coils are already a source (biotsavart.f90/neo_il, coil_tools.f90 writes vector_potential_il_fourier/Anvac_fourier); field_divB0 and field_mesh3d/field_c_mod are older lab-frame splined fields.

So this issue is NOT a new representation -- it is: (1) wrap the chartmap/Boozer field as a Cartesian/cylindrical field_t (compute_abfield = invert+field_at_logical, paid once per node at build); (2) field_mesh_init_with_field over an (R,phi,z) box + spline_field_init; (3) in orbit_cpp_boris read B from the spline_field_t instead of the per-step inversion. The cylindrical/Cartesian/2D 'variants' are just mesh frames, not new types. Remaining genuinely-new libneo work is only #326 (kink-robust BC + compact storage), attached to neo_spline_field's spline backend. No 'fieldmap' naming needed.

[krystophny]

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

Refine — the premise is alive on main, but the file the issue names is not, and half the libneo argument needs correcting. Three corrections before the spec:

1. orbit_cpp_boris.f90 does not exist on main. It only ever lived on origin/feature/simple-6d-port (git log --all -- src/orbit_cpp_boris.f90 → 703404f, 86b2b2f, b3e4b59, …; git branch -a --contains 703404f → feature/simple-6d-port and z_archive/feature/simple-6d-port only). 123d26a (Extract shared GC Gauss implicit-RK core (orbit_rk_core); drop CPP and symplectic FO #429, "Extract shared GC Gauss implicit-RK core; drop CPP and symplectic FO") removed the MODEL_CPP / MODEL_FO_SYMPL codes and their dispatch shells; it did not delete a CPP source file because none was on main. The production Cartesian full-orbit pusher on main is src/orbit_fo_boris.f90, orbit_model = ORBIT_FULL_ORBIT = 7 (src/params.f90:74), landed by 36d97ef (Add full-orbit (FO) Boris pusher (orbit_model=7) #426). orbit_cpp_chartmap_metric.f90 likewise does not exist on main, so the "class-dispatched metric/christoffel, NOT GPU-portable" half of the GPU-hostility argument is void — there is no metric/Christoffel machinery left. Everything below is written against orbit_fo_boris.
2. The core premise survives intact. The per-step Cartesian→logical inversion is still there verbatim at the new path (see Current state). The work is unimplemented: rg 'neo_spline_field|spline_field_t|field_mesh_t' src/ test/ in SIMPLE returns nothing.
3. The "Correction" comment on the issue is right that libneo already has the representation, and it is right that option (B) exists — but it misses that the backend is not GPU-callable. spline_field_init_acurl (option B, B = curl A) is already implemented including the cylindrical form. However neo_spline_field is built on the non-batch SplineData3D, which has zero OpenACC annotations. Details below.

Purpose

The Cartesian full-orbit Boris pusher is the gyro-resolved reference for SIMPLE's guiding-centre losses and the direct comparison point against ASCOT5. Per #422 it is roughly 5x slower than ASCOT5 on the same case, and the cause is architectural: the field lives in logical Boozer/chartmap coordinates while the particle is pushed in Cartesian, so every field evaluation solves a Cartesian→logical inverse. Sampling the same exact curl A field once per grid node at setup and reading it back from a real-space (R, φ, Z) spline removes the inversion from the hot loop entirely, with no physics change. It is also the precondition for a GPU full-orbit push: a data-dependent damped Newton with a line search cannot be offloaded; a spline read can.

Current state

The inversion, per step. fo_step (src/orbit_fo_boris.f90:444) calls cart_field once (:455), which is:

  subroutine cart_field(x, u_guess, Bvec, Bmod, gradB, u_out, status)
    ...
    call to_wedge(x, xw, ca, sa)
    call invert_cart_warm(xw, u_guess, u, status)
    if (status /= FO_OK) return
    u_out = u
    call field_at_logical(u, Bw, Bmod, gw, Jc, status)

(src/orbit_fo_boris.f90:323-337)

invert_cart_warm (:73-98) calls invert_warm_newton (:80), whose newton_from (:122-187) has integer, parameter :: maxit = 30, maxls = 30 (:126); on FO_LOCATE_FAIL it falls back to the chartmap's multi-seed cs%from_cart inside a select type (:82-88). Both are GPU-hostile: data-dependent trip counts, a line search, and class dispatch on the hot path.

Two more inversions per output step. fo_to_gc (:497-537) calls locate twice — once for the particle (:509), once for the Larmor-corrected guiding centre (:521) — and each locate is another to_wedge → invert_cart_warm → field_at_logical (:347-348, :353). The loss test is a level set on the result:
if (u_gc(1) >= 1.0_dp .and. u_p(1) >= 1.0_dp - GC_PARTICLE_GAP) status = FO_LOSS (:536, with GC_PARTICLE_GAP = 0.05_dp at :44).

The push only needs Bvec. In fo_step, gradB is declared (:447), filled by cart_field (:455) and then never read — the Boris rotation uses Bvec alone (:458-464). This kills the issue's _der2 concern: option (B) (spline A, take B = curl A from first derivatives) is sufficient for the push. grad|B| would only be needed by a future GC-drift consumer.

libneo — what already exists (paths relative to the libneo checkout, fe3f76e):

subroutine spline_field_init_acurl(self, field_mesh, coordinate_system, order)
    ! Spline only A; reconstruct B = curl A from the A-spline. The stored B
    ! arrays are ignored, so div B = 0 holds to spline precision regardless of
    ! grid resolution.

(src/field/spline_field.f90:43-46). compute_bfield dispatches on the b_from_curl flag (src/field/spline_field.f90:114-117) into curl_a (:123-144), whose cylindrical branch (:133-138) is
B(1) = dA_dx(3,2)/R - dA_dx(2,3), B(2) = dA_dx(1,3) - dA_dx(3,1), B(3) = (A(2) + R*dA_dx(2,1) - dA_dx(1,2))/R
— i.e. it expects physical cylindrical components (A_R, A_φ, A_Z) with x = (R, φ, Z). field_mesh_init_with_field (src/field/field_mesh.f90:64-113) samples any field_t onto the grid via field%compute_abfield. So option (B) is done in libneo and needs no new representation.

libneo — the gap. neo_spline_field uses use interpolate, only: SplineData3D, construct_splines_3d / evaluate_splines_3d, evaluate_splines_3d_der (src/field/spline_field.f90:4-5). That is the non-batch backend (SplineData3D type at src/interpolate/interpolate.f90:98-105, constructor construct_splines_3d at :455). grep -c 'acc ' src/interpolate/interpolate.f90 returns 0; the 11 !$acc routine seq sites the issue cites are all in src/interpolate/batch_interpolate_3d.f90 (grep -c 'acc routine seq' → 11). So reusing neo_spline_field buys the CPU speedup and div B = 0 immediately, but not GPU-callability; that needs a batch-backed variant of spline_field_t (a libneo follow-up, not a blocker for the CPU win).

Not to be confused with field_splined. src/field/field_splined.f90 splines in scaled logical coordinates (r = sqrt(s), see the header at :5-7) and therefore still needs the inversion. It does already use BatchSplineData3D (:28, :44, :195), so SIMPLE is familiar with the batch API.

Reproduction / evidence

The motivating measurement is in #422 (SIMPLE full orbit ~5x slower than ASCOT5 on the same case). The in-repo gate and profiling handle:

cd /path/to/SIMPLE
make CONFIG=Fast                      # Fast is a valid build type, CMakeLists.txt:65
make test TEST=test_fo_boris          # existing gate, must stay green

test/tests/test_fo_boris.f90 runs the reactor-scale QA VMEC field on wout.nc — downloaded at configure time by file(DOWNLOAD ...) in test/tests/CMakeLists.txt:23-28 from the simsopt wout_LandremanPaul2021_QA_reactorScale_lowres_reference.nc, and symlinked into the test binary dir at test/tests/CMakeLists.txt:37. (test/test_data/get_test_data.sh is a separate interactive legacy script; it is not what the build uses.) The test calls init_params(norb, 2, 4, 3.5e6_dp, 16384, 1, 1.0d-13) (test_fo_boris.f90:40) and traces three markers (:51-53)

• passing [0.5, 0.5, 0.2, 1.0, 0.9]
• trapped [0.5, 0.5, 0.2, 1.0, 0.2]
• near-axis [0.04, 0.5, 0.2, 1.0, 0.7]

asserting Emax < 1.0e-3_dp (:149).

To confirm the inversion dominates before committing to the work, build make CONFIG=Profile (also a valid build type, CMakeLists.txt:65) and perf record a test_fo_boris run; the expectation is that time sits in newton_from / evaluate_cart / covariant_basis, not in the Boris rotation. Do this first — it is the cheap falsifier for the whole issue.

Verification recipe for the feature (this is a performance/representation change, not a bug fix): the new path must reproduce the chartmap-inversion path to within spline error on the same markers.

• Same three markers, same dtaumin, same number of steps; final (s, θ_B, φ_B) from fo_to_gc agreeing with the reference path to < 1e-3 in s after 1e4 steps at the recommended grid, and |dE/E0| < 1e-3 unchanged. The Boris rotation is exact for any B, so an energy regression means a broken field assembly, not interpolation error.

Behaviour (BDD)

Scenario: Direct real-space field reproduces the chartmap-inversion field
  Given the reactor-scale wout.nc fixture used by test/tests/test_fo_boris.f90
  And a cylindrical (R, phi, Z) field mesh over one field period with
      n_points = [96, 48, 96] and is_periodic = [.false., .true., .false.]
  When B is evaluated at 1000 random Cartesian points with u(1) < 0.95 by both
      cart_field (inversion path) and the new spline read
  Then the relative difference in |B| is below 1e-4
  And the relative difference in each Cartesian component of B is below 1e-3

Scenario: The splined-A field is divergence-free
  Given the same mesh built through spline_field_init_acurl with COORD_CYLINDRICAL
  When div B is evaluated by centred finite differences with h = 1e-4 * (R_max - R_min)
      at 200 interior points
  Then |div B| * L / |B| is below 1e-8, where L = R_max - R_min

Scenario: Full-orbit trace matches the inversion path
  Given orbit_model = 7 and fo_field_repr = 1
  When the passing marker [0.5, 0.5, 0.2, 1.0, 0.9] is traced for 10000 steps
  Then |dE/E0| stays below 1e-3
  And the final guiding-centre s from fo_to_gc agrees with the fo_field_repr = 0 path
      to within 1e-3

Scenario: No node of a used cell lies outside the plasma
  Given the splined radial label s(x) on the same mesh
  When the mesh box is chosen by the build-time LCFS bound
  Then every grid node used by an accepted interpolation cell has s <= 1.0
  And a marker whose splined s crosses 1.0 is flagged FO_LOSS by the same
      GC_PARTICLE_GAP criterion as src/orbit_fo_boris.f90:536

Implementation plan

Steps 1-6 are SIMPLE-only and land in one PR — no libneo change is required, because spline_field_init_acurl already exists. Steps 7-8 are libneo follow-ups and must not block.

1. New file src/field/field_chartmap_cart.f90 — module field_chartmap_cart, type chartmap_cart_field_t, extends(field_t) (libneo module neo_field_base; field_t is abstract with exactly three deferred bindings compute_afield, compute_bfield, compute_abfield, libneo/src/field/field_base.f90:5-10). Implement all three with x = (R, phi, Z):

   • convert (R, φ, Z) → Cartesian, then call to_wedge + invert_cart_warm (currently private in orbit_fo_boris, see the note below);
   • for B: field_at_logical (src/orbit_fo_boris.f90:266-316) already returns the wedge-frame Cartesian Bvec; rotate back with rotz and then into physical cylindrical (B_R, B_φ, B_Z);
   • for A: field_at_logical does not return A. Use fo_eval_field (src/orbit_fo_field.f90:23-53), which returns Acov — covariant components in logical coordinates. Convert with the chartmap Jacobian: Jc = ∂x/∂u from ref_coords%covariant_basis(u, Jc) (used the same way at src/orbit_fo_boris.f90:94-95), then A_cart_k = Acov_i * (Jc^{-1})^i_k, then rotate to physical cylindrical A_R = A_x cosφ + A_y sinφ, A_φ = -A_x sinφ + A_y cosφ, A_Z = A_z. Physical components are what libneo's curl_a expects — the 1/R and A(2) terms at libneo/src/field/spline_field.f90:136,138 are the physical-component cylindrical curl. Getting this wrong is the single most likely source of a wrong B; the BDD scenario 1 tolerance is the check.
   • compute_abfield has no status argument and field_mesh_init_with_field has no error path (libneo/src/field/field_mesh.f90:97-112). A node outside the chart must therefore return a finite clamped/extrapolated value, never error stop. Record which nodes were clamped for step 4.

   This is build-time only — the inversion cost is paid once per node here. To reach to_wedge / invert_cart_warm / field_at_logical, note the module is private (src/orbit_fo_boris.f90:27) with a narrow public list (:48, :50). Either widen that list, or (preferred, smaller blast radius) move the locate machinery into a new src/orbit_fo_locate.f90 that both modules use.

2. Same file: the s-label sampler. field_t has no radial-label member, so carry a separate type(SplineData3D) :: s_spline built directly with construct_splines_3d (libneo/src/interpolate/interpolate.f90:455), filled from u(1)**2 at each node using the same inversion (this is the rho -> s convention used at src/orbit_fo_boris.f90:522). For off-chart nodes store a monotone continuation s = 1 + d/L from the boundary distance so the level set stays continuous.

3. New file src/field/field_realspace.f90 — module field_realspace, holding type(spline_field_t) :: fld (libneo module neo_spline_field), the s_spline, the box limits, nfp, and a build_realspace_field(limits, n_points, order). Both libneo entry points are type-bound, so the calls are:

   type(field_mesh_t) :: mesh
   type(chartmap_cart_field_t) :: src
   call mesh%field_mesh_init_with_field(limits, src, n_points, &
                                        is_periodic=[.false., .true., .false.])
   call self%fld%spline_field_init_acurl(mesh, COORD_CYLINDRICAL, order)

   limits is real(dp), dimension(3,2) with limits(:,1) the minima and limits(:,2) the maxima (libneo/src/field/field_mesh.f90:68, :80-82). COORD_CYLINDRICAL = 1 is a public parameter of neo_spline_field (libneo/src/field/spline_field.f90:9). Use spline_field_init_acurl, not spline_field_init — option (B) in the issue, and it splines only A1,A2,A3 so it costs half the storage.

   Add realspace_bfield_cart(x_cart, Bvec, s, status): wrap phi into [0, 2*pi/nfp) (the field period is already computed at src/orbit_fo_boris.f90:227-236 from cs%num_field_periods), call self%fld%compute_bfield, then rotate the physical cylindrical (B_R, B_φ, B_Z) into the global Cartesian frame.

4. Box sizing (do this before coding the pusher branch). Choose limits from the LCFS bounding box, shrunk so that every used cell has s <= 1 at all eight nodes. SplineData3D allocates coeff(0:order(1), 0:order(2), 0:order(3), n1, n2, n3) (libneo/src/interpolate/interpolate.f90:474) — a dense (order+1)^3 block per grid node, so memory is n_quantities * (order+1)^3 * n1*n2*n3 * 8 bytes. With 4 quantities (A1, A2, A3 + s):

   • [96, 48, 96] cubic → ~0.91 GB
   • [128, 64, 128] cubic → ~2.1 GB
   • [96, 48, 96] quintic → ~3.1 GB
   • [128, 64, 128] quintic → ~7.2 GB

   Start at [96, 48, 96] cubic and sweep upward. Quintic is a memory decision, not a capability one — construct_splines_3d takes order(3) and allocates accordingly; the explicit 3 <= order <= 5 guard is in the batch path (libneo/src/interpolate/batch_interpolate_3d.f90:95-98), so validate the order in SIMPLE rather than relying on libneo to reject a bad value on the non-batch path.

5. Edit src/orbit_fo_boris.f90 — add a module-level logical :: use_realspace_field = .false., set from the new namelist key. In cart_field (:323), branch at the top: if the flag is set, call realspace_bfield_cart and return, skipping to_wedge / invert_cart_warm / field_at_logical entirely. In fo_to_gc (:497), replace the two locate calls (:509, :521) with a splined-s read plus bhat = B/|B| from the same read, preserving the GC_PARTICLE_GAP rejection at :536 unchanged. Keep the u_guess / st%u plumbing so the inversion path stays selectable and byte-identical when the flag is off.

6. Edit src/params.f90 — add fo_field_repr to the &config namelist (declared across src/params.f90:164-...; orbit_model appears at :168), values 0 (chartmap inversion, default) and 1 (real-space spline), plus fo_field_nR / fo_field_nphi / fo_field_nZ and fo_field_order. Document in DOC/config.md. Edit src/CMakeLists.txt to compile field/field_chartmap_cart.f90 and field/field_realspace.f90 before orbit_fo_boris.f90. New test test/tests/test_fo_realspace_field.f90 implementing all four BDD scenarios, registered in test/tests/CMakeLists.txt next to the test_fo_boris block (test/tests/CMakeLists.txt:371-380), with LABELS "integration" and a TIMEOUT well above the existing 120 (the setup sampler is the slow part).

7. libneo follow-up (separate PR, non-blocking): batch-backed spline field. In libneo/src/field/spline_field.f90, add a variant of spline_field_t backed by BatchSplineData3D / construct_batch_splines_3d / evaluate_batch_splines_3d_der so the 11 !$acc routine seq evaluators in libneo/src/interpolate/batch_interpolate_3d.f90 become reachable. Only then is a GPU full-orbit push possible; the CPU win in steps 1-6 does not depend on it.

8. libneo follow-up: batch_interpolate_3d: kink-robust BC variant + compact coefficient storage libneo#326 (kink-robust BC + compact coefficient storage, as noted in the first issue comment). Step 4's "keep every used cell inside the LCFS" makes the separatrix kink a non-issue for the first landing, so Finish GVEC straight-fieldline support end-to-end and add PEST #326 is an accuracy/memory improvement, not a prerequisite.

Acceptance criteria

• make CONFIG=Fast builds with the two new source files wired into src/CMakeLists.txt.
• make test TEST=test_fo_realspace_field passes, covering all four BDD scenarios.
• make test TEST=test_fo_boris still passes with fo_field_repr = 0 (default), with identical trajectories.
• make test-regression shows no change in any golden-record case — the default path is untouched, so any diff is a bug.
• With fo_field_repr = 1 on the wout.nc fixture at [96, 48, 96] cubic: max |B_spline - B_inversion| / |B| < 1e-3 over 1000 interior samples with u(1) < 0.95.
• |div B| * L / |B| < 1e-8 by centred finite differences. This must hold by construction with spline_field_init_acurl; a failure means the covariant-to-physical A conversion (step 1) or the component rotation is wrong.
• |dE/E0| < 1e-3 over 10000 steps for the passing, trapped and near-axis markers of test_fo_boris.
• Measured field-eval time per step drops by at least 3x versus fo_field_repr = 0 on the same fixture (measure with CONFIG=Profile + perf, same recipe as the pre-work falsifier).
• make test-smoke passes.

References

• Boris, J.P. (1970), Relativistic plasma simulation — optimization of a hybrid code, Proc. 4th Conf. Num. Sim. Plasmas, NRL — the pusher already 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, doi:10.1063/1.4818428 — volume preservation, which is why energy stays bounded for any B, splined or not.
• Hirvijoki, E. et al. (2014), ASCOT: Solving the kinetic equation of minority particle species in tokamak plasmas, Comput. Phys. Commun. 185, 1310, doi:10.1016/j.cpc.2014.01.014 — the ASCOT-family (R, φ, Z) direct-read field design this issue mirrors.
• Pfefferlé, D., Cooper, W.A., Graves, J.P., Misev, C., arXiv:1412.5464 — already cited in src/orbit_fo_boris.f90:543-544 for the pseudo-Cartesian near-axis chart; unchanged by this issue.
• Lekien, F. and Marsden, J. (2005), Tricubic interpolation in three dimensions, Int. J. Numer. Meth. Engng 63, 455 — the compact one-set-per-node storage that libneo#326 asks for, in contrast to libneo's dense per-node (order+1)^3 block.
• Steffen, M. (1990), A simple method for monotonic interpolation in one dimension, Astron. Astrophys. 239, 443 — the limiter referenced by libneo#326 for the LCFS kink.

Note on ASCOT5 specifics: the claims in the issue body that ASCOT5 splines B components independently, accepts div B ≠ 0, and reads its rho >= 1 loss criterion off the same spline are the issue author's, and are not verified here — the triage pass could not confirm a citable ASCOT5 code paper backing them. Treat them as motivation, not as a specification input. Nothing in steps 1-8 depends on them: spline_field_init_acurl gives div B = 0 regardless of what ASCOT5 does.

Scope and risk

One SIMPLE PR for steps 1-6 (new opt-in code path, default behaviour untouched), plus two independent libneo follow-ups. Main risks:

• The covariant-to-physical A conversion in step 1. This is new arithmetic with no existing reference in SIMPLE, and a sign or 1/R error produces a plausible-looking but wrong B. BDD scenario 1 against cart_field is the gate; do that before wiring the pusher.
• Memory (step 4): the dense per-node coefficient layout is the binding constraint. Size the grid first.
• The LCFS edge, where losses are decided — mitigated by keeping every used cell inside s = 1.

Because the default fo_field_repr = 0 path is unmodified, golden records cannot move; if they do, the branch in cart_field was placed wrong.

Related, as named in the issue and its comments: #422 (perf umbrella, the ~5x measurement), #419/#420/#421, itpplasma/libneo#326 (kink-robust BC + compact storage), itpplasma/benchmark-orbit-proxima#2 (the W7-X high-mirror benchmark case).

Suggested split, if this is too large for one PR: (a) steps 1-2 plus a unit test that chartmap_cart_field_t reproduces cart_field node-by-node — this is where the A-conversion risk lives and it is worth isolating; (b) steps 3-6, the pusher branch and the benchmark; (c) libneo step 7, the batch-backed spline_field_t, which is what actually unlocks GPU.

Also please update the issue title and body to say orbit_fo_boris, since orbit_cpp_boris does not exist on main.