[spectre-05] Add symplectic field-line Poincare app for SPECTRE equilibria

[krystophny]

Stage: SIMPLE / field-line diagnostics (milestone M1)
Source language: Fortran 2008
Manuals to read first: itpplasma/SPOIC spoic.f90:106-141 (trace — the working symplectic field-line scheme against SPEC data: Newton on the implicit radial update, A_theta as canonical momentum, zeta as time); SPECTRE fortran_src/poincare_mod.F90 and spectre/field_tracer.py (the reference tracers to match); DOC/coordinates-and-fields.md.
Depends on: itpplasma/libneo#371

Goal

SIMPLE ships a standalone app spectre_poincare.x producing Poincare sections and rotational-transform profiles for a SPECTRE equilibrium by symplectic field-line integration, per volume. Field lines never cross MRxMHD interfaces (B.n = 0 on both sides), so no crossing logic is involved; this milestone validates the whole libneo evaluation stack against SPECTRE's own tracer before any guiding-center work.

Files to edit

• app/spectre_poincare.f90: NEW. Namelist input: SPECTRE file, volume list, seeds per volume, turns, section plane zeta0; output one poincare_volN.dat per volume (columns: turn index, s, theta, R, Z).
• test/tests/test_spectre_poincare.py: NEW. Runs the app on the libneo fixture and compares against committed SPECTRE-tracer reference sections.
• DOC/coordinates-and-fields.md: document the SPECTRE field input and its coordinate conventions (mandatory per repo policy).
• src/CMakeLists.txt / CMakeLists.txt: link LIBNEO::spectre, add app target.

Behavior to implement

Field-line flow as canonical system with zeta as time, conjugate pair (A_theta, theta), radial coordinate solved implicitly (semi-implicit symplectic Euler as in SPOIC, second-order midpoint optional):

! Newton solve for r at fixed theta: A_theta(r_new) - A_theta_old = h * dH/dtheta
fun = f%dAth(1)*(f%Ath - Athold) + h*(f%dAzt(2)*f%dAth(1) - f%dAzt(1)*f%dAth(2))
r = r - fun/dfun
theta = theta + h*f%dAzt(1)/f%dAth(1); zeta = zeta + h

Seeds placed on zeta = zeta0, radially spread inside each volume, off-axis offset in volume 1.

Acceptance scenarios (BDD)

• Given the fixture equilibrium, when tracing 16 seeds per volume for 500 toroidal turns, then no field line leaves its volume (|s| <= 1 throughout) — direct numerical witness of B.n = 0.
• Given the same seeds traced by SPECTRE's Poincare tracer (reference data committed), then section point clouds agree: for regular surfaces, radial spread of matched surfaces within 1e-6; island chains reproduced with matching O-point positions and widths within 5%.
• Given a seed on a regular surface, then the fitted rotational transform matches SPECTRE's write_transform output for that surface to 1e-6.
• Given step-halving, then section points converge at the order of the chosen scheme (no first-order pollution from the Newton tolerance).

Success criteria

make && make test TEST=spectre_poincare

Non-goals

• No guiding-center orbits, no interface crossing, no losses.
• No chaos classification beyond visual/reference comparison (SIMPLE classifiers arrive with guiding-center milestones).

Verification

make test TEST=spectre_poincare

[krystophny]

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

Keep open: the app ships, but two of the four acceptance scenarios are met at 3-4 orders of magnitude looser tolerance than specified, and the island criterion is not implemented at all

Purpose

The deliverables listed under "Files to edit" all exist on main, so the construction half of this milestone is done. What is not done is the validation half, which is the stated point of the milestone ("this milestone validates the whole libneo evaluation stack against SPECTRE's own tracer before any guiding-center work"). The committed test passes, but it passes against tolerances that do not test what the acceptance scenarios ask for. This comment narrows the issue to the remaining validation gap plus one stale-documentation fix.

Current state

What is present and verified:

• app/spectre_poincare.f90 (268 lines), added in 9fd7d6c ("Add symplectic field-line Poincare app for SPECTRE equilibria ([spectre-05] Add symplectic field-line Poincare app for SPECTRE equilibria #437)", 2026-07-11, ancestor of origin/main). Reads the &poincare namelist — spectre_file, zeta0, nturns, nsteps_per_period, nseeds_per_volume, volumes, seed_file, spline_n_{s,theta,zeta}, spline_order (app/spectre_poincare.f90:32-34) — and writes poincare_vol<N>.dat with columns seed, section, s, theta, R, Z plus iota_vol<N>.dat (app/spectre_poincare.f90:11-12, :165-166). Default seeds sit at theta = 0.2 with s spread over [-0.85, 0.85] (:21-23, :131-141); no seed lands on s = 0, so the volume-1 off-axis requirement holds.
• src/field_line_spectre.f90 (301 lines). field_line_step (:102-133) advances (theta, psi) by an implicit midpoint map, solved by solve_midpoint / midpoint_residual / midpoint_jacobian (:135-207) with an endpoint radius recovered in recover_radius (:209-230). This upgrade from the original semi-implicit Euler landed in fc20806 ("Use implicit midpoint for SPECTRE field lines (Use implicit midpoint for SPECTRE field lines #483)", 2026-07-13, also on main).
  Note: the midpoint Jacobian is not analytic. midpoint_jacobian (:184-207) builds both columns by central finite differences with delta = fd_scale * ..., fd_scale = epsilon(1.0_dp)**(1/3) (:30); the module's own comment at :99-101 says "Numerical Jacobians use only values and first derivatives of the vector potential." Only recover_radius (:228) uses an exact derivative (dAth(1)).
• test/tests/test_spectre_poincare.py (202 lines), registered at test/tests/CMakeLists.txt:852-863 against ${libneo_SOURCE_DIR}/test/spectre/data/spectre_test.h5 and test/tests/data/poincare.npz, label integration;python, timeout 120 s.
• Build wiring: CMakeLists.txt:522-535 defines spectre_poincare.x and links simple. make test TEST=spectre_poincare resolves through Makefile:24 (ctest -R "$(TEST)"), so the success-criteria command is real.
• DOC section 6.6 "SPECTRE Field-Line Coordinates" at DOC/coordinates-and-fields.md:589-625.

Where the test falls short of the acceptance scenarios:

Issue scenario	Issue tolerance	What the test actually asserts
Confinement, 16 seeds/volume, 500 turns	|s| <= 1 throughout	nseeds=4, nturns=60 (test_spectre_poincare.py:100)
Section agreement on regular surfaces	radial spread 1e-6	median nearest-neighbour distance < NN_MEDIAN_TOL = 5.0e-3 m (:26, :138)
Island chains: O-point positions and widths within 5%	5%	no such check exists — rg -i 'island|o-point' over the test and DOC returns nothing
Rotational transform vs SPECTRE write_transform	1e-6	iota anywhere inside the reference min/max band widened by 10% of its own width (:146-154)

The iota band is not a tight check: computed from poincare.npz, it spans [0.29011, 0.29376], [0.28709, 0.29730], [0.28299, 0.30092] for volumes 1-3, i.e. a 1.2%-6.2% acceptance window against a specified 1e-6 (~3e-6 relative). Nor is it per-surface: the app seeds at theta = 0.2 over s in [-0.85, 0.85] while the reference seeds_vol0 is theta = 0 over s in [-0.5, 0.85], so in the axis volume the app's innermost seed lies radially inside anything the reference sampled and is being compared to a band it was never meant to be in.

The 5e-3 m section tolerance is documented in the test itself as a floor, not a target: "The reference sections floor at ~3e-3 m relative to our tracer (SPECTRE uses a different integrator and radial resolution)" (:24-26). That is an honest note, but it means the 1e-6 cross-code agreement the issue asked for has never been demonstrated, only bypassed.

Separately, the DOC prose is now wrong on two counts. DOC/coordinates-and-fields.md:610-614 still describes "a semi-implicit symplectic Euler map: a Newton solve on s for the implicit momentum update (analytic Jacobian from the ss and st second derivatives, |F| tolerance 1e-13), then explicit theta and zeta updates". Since fc20806 the map is implicit midpoint, and the Jacobian was never analytic in the current code.

Reproduction / evidence

cd build && ctest --test-dir test -R test_spectre_poincare --output-on-failure -V

The test passes. Its own stdout is the witness for the gap: BDD2 vol<N>: outer nn median = ... prints a residual on the order of 1e-3 m against a 5e-3 m gate, and BDD3 prints the app and reference iota ranges without ever differencing them per surface.

To see that the specified comparison is not being made, inspect the reference seeds:

python3 -c "import numpy as np; z=np.load('test/tests/data/poincare.npz'); print(z['seeds_vol0']); print(z['seeds_per_vol'], z['num_transits'])"
# seeds at theta = 0, s in [-0.5, 0.85]; 4 seeds; 200 transits

The app is never run on those seeds, even though seed_file makes it trivial to do so.

Behaviour (BDD)

• Given the fixture equilibrium and the exact seed set stored in poincare.npz (seeds_vol<i>, written to a seed_file), when spectre_poincare.x traces them for the same num_transits as the reference, then each app section curve matches the reference curve for the same seed, surface by surface, and the test asserts and prints the achieved residual against a tolerance justified in a comment that names the limiting error source.
• Given the same matched run, when the fitted iota is compared to iota_vol<i> seed by seed, then the per-seed difference is asserted against a stated tolerance, not membership in a band spanning all seeds.
• Given the fixture equilibrium, when 16 seeds per volume are traced for 500 toroidal turns, then no seed terminates and max|s| < 1 — this scenario needs no reference data and can be gated on wall time separately from the comparison test.
• Given the fixture equilibrium contains an island chain (or, if it does not, given that fact recorded explicitly), when sections are compared, then either O-point positions and widths agree within 5%, or the issue is amended to drop the island criterion with a one-line reason.

Implementation plan

1. In test/tests/test_spectre_poincare.py, add a matched-seed comparison: dump ref["seeds_vol<i>"] to a seeds.dat in the volume-index convention the app expects (volume_seeds at app/spectre_poincare.f90:225-243 matches on ivol == lvol - 1, i.e. 0-based in the file), run with seed_file= and nturns equal to ref["num_transits"], and compare per seed rather than per cloud. The seed_file path is already exercised by the existing BDD5 block (:182-191), so this is reuse, not new plumbing.
2. Replace the band membership check for iota with a per-seed difference against iota_vol<i>. Pick the tolerance from the matched run's measured residual, and put the measured number in the comment.
3. Determine empirically what the matched-seed residual actually is. If it lands near 1e-6, tighten NN_MEDIAN_TOL accordingly and the scenario is genuinely satisfied. If it plateaus at 1e-3 because the reference tracer's own accuracy dominates, say so in the test and amend the issue's tolerance with that measurement as the justification — do not leave a 5e-3 gate standing against a 1e-6 requirement with no recorded reason.
4. Add the 16-seed/500-turn confinement run as its own registered test (or an opt-in slow label), so the specified confinement witness exists without inflating the 120 s comparison test.
5. Resolve the island criterion: either implement the O-point check, or record in the issue that spectre_test.h5 has no island chain to check and strike the criterion.
6. Fix DOC/coordinates-and-fields.md:610-614 to describe the implicit midpoint map and the finite-difference Jacobian (fd_scale = eps**(1/3), central differences in s and theta_delta), keeping the 1e-13 residual tolerance, which is still correct (src/field_line_spectre.f90:29).

Acceptance criteria

• make test TEST=spectre_poincare passes, and the test compares app-vs-reference sections on the reference's own seeds, per surface.
• Every numeric tolerance in the test is either the one written in this issue or accompanied by a measured justification for the value used.
• A confinement scenario at 16 seeds per volume and 500 turns exists and passes.
• The island criterion is implemented or explicitly retired in this issue with a reason.
• DOC/coordinates-and-fields.md section 6.6 describes the integrator that is actually in src/field_line_spectre.f90.

References

• app/spectre_poincare.f90, src/field_line_spectre.f90, test/tests/test_spectre_poincare.py, test/tests/data/poincare.npz, DOC/coordinates-and-fields.md:589-625, CMakeLists.txt:522-535, test/tests/CMakeLists.txt:852-863.
• Commits 9fd7d6c (initial app, [spectre-05] Add symplectic field-line Poincare app for SPECTRE equilibria #437) and fc20806 (implicit midpoint, Use implicit midpoint for SPECTRE field lines #483), both ancestors of origin/main.
• No specific literature reference; the tolerance question here is an implementation detail internal to the code and its reference fixture.

Scope and risk

size/S, test and documentation only. No change to app/spectre_poincare.f90 or src/field_line_spectre.f90 is anticipated unless step 3 reveals a real disagreement with the reference tracer, which is exactly the outcome this milestone was written to detect. Risk of the alternative — closing now — is that the SPECTRE evaluation stack is treated as validated to 1e-6 by the downstream guiding-center milestones (#438-#442) when it has only been checked to 5e-3 m and a few percent in iota.