Establish GACODE integration with Sauter/Redl and NEO (#550 phases 1-4)

[HongSik-Yun-Fusion]

Implements phases 1-4 of #550: the analytic neoclassical foundation, the shared
GACODE runtime, input.gacode interoperability, and the NEO adapter. Phases 5-8
(IMAS projection, the validation study, the 48224 notebook, database promotion)
stay open as follow-ups under milestone 25.

VAFT had no kinetic/transport external-code boundary — a repo-wide grep for
gacode|GACODE|NEO|TGLF|CGYRO returned nothing — and no neoclassical physics
beyond bootstrap_current_fraction, a 0.3·√β_p heuristic with unsourced
coefficients (#361). This adds both halves, and keeps them as distinct
computational identities: there is no model= dispatcher, because a fitted
formula and a drift-kinetic solve answer related but different questions.

What landed

vaft.formula.neoclassical — twelve pure array/scalar kernels: trapped
fraction, the Sauter Coulomb logarithms and collisionalities, the Spitzer
reference conductivity, and both the Sauter 1999 (with the 2002 erratum) and
Redl 2021 bootstrap and conductivity formulations. Both are provided because at
VEST's trapped fraction the 1999 fit is an extrapolation and the 2021 refit is
not — see the measurement below.

vaft.code.gacode — suite-level runtime shaped like vaft.code.gpec and
disciplined like vaft.code.nubeam. $GACODEHOME with $GACODE_ROOT as a
compatibility fallback, explicit platform resolution, a typed GACODEProfile,
and a pure-Python input.gacode reader/writer. GACODE's own env contract is
set for the subprocess, never redefined, so a tree configured through
shared/bin/gacode_setup keeps working.

vaft.code.gacode.neo — input generation, execution through the launcher
(which is what stamps out.neo.version), parsers for every out.neo.* product,
and a complete NeoOutputs container with a versioned JSON round trip. Nothing
is written to an IDS; that audit is phase 5 and is deliberately still open.

Verification, against a real GACODE build

Built from gafusion/gacode 6357db30 via the new external/gacode/macos.sh.

Check	Result
NEO's shipped reg18 regression case	out.neo.prec = 0.12268957E+02, exact
A VAFT-written input.gacode fed back to NEO	bit-identical out.neo.* to the original run
input.gacode read/write round trip	byte-identical except three unit labels expro itself renamed
VAFT Sauter / Redl vs NEO's own, reg18	4e-9 / 2e-9
VAFT Sauter / Redl vs NEO's own, VEST 48224	9e-9 / 3e-8
Packaged 48224 kinetic state, ODS → NEO	runs end to end, f_trap = 0.73

Both cross-checks run offline against committed fixtures from real runs, so
CI verifies them without GACODE installed. The one test that needs an
installation is gated on $GACODEHOME, in the style test_nubeam_adapter.py
established.

The physics result

NEO carries its own Sauter and Redl implementations and writes both to
out.neo.theory on every run, so the same file gives the analytic models and
the drift-kinetic answer. At VEST's 48224 state, r/a = 0.5, f_trap = 0.73:

• Sauter 1999 overshoots the drift-kinetic <j∥B> by 6.9%
• Redl 2021 overshoots it by 1.8%

At reg18's conventional aspect ratio the two fits agree with each other to
under 1%. That gap growing with the trapped fraction is the justification for
implementing both, and it is asserted as a test rather than left as prose.

Four decisions made explicitly, not guessed

• A rho_tor_norm that is really the sqrt(psi_N) proxy is re-derived, or the
  conversion stops (rho_tor_norm is actually rho_pol: sqrt(psi_norm) written into the toroidal-flux field #276, [Process contract D] Document profile and atomic processing #420). Writing sqrt(psi_N) under GACODE's rho is
  exactly the defect those issues exist to prevent.
• Equilibrium and core_profiles slices are not paired outside a derived
  tolerance; requested, equilibrium and profile times are all recorded.
• A non-positive density or temperature stops the conversion naming the grid
  point, since GACODE takes logarithmic gradients. The packaged 48224 profiles
  reach exactly zero at the boundary, so this fires — rho_max= makes the
  truncation the caller's explicit, recorded decision rather than a silent clip.
• Signs are converted, COCOS 11 → 2. input.gacode runs its toroidal angle
  clockwise from above, and expro reads the field directions from the signs of
  torfluxa and q. GACODE is registered in vaft.data.cocos as COCOS 2 —
  marked unconfirmed, because the index is inferred from expro and from reg18's
  own signs rather than documented upstream.

Worth knowing for anyone reading NEO output elsewhere

pygacode's own read_theory is stale against neo_theory.f90: it reads the
per-species block three-wide, where the writer emits two values per species and
then two trailing scalars — 24 columns against the 23 reg18 actually has. The
layout here is the writer's, pinned by tests at two and three species, which is
what separates the two readings.

Tests

146 new tests across three files, none in the core selection. Local run:
-m core 5067 passed; the new suites 146 passed, 1 skipped (the
$GACODEHOME-gated integration test, which passes when set).

Cold review

A fresh review pass found six issues. All six are fixed in e206409, each with a
regression test; five were reproduced against the real NEO build.

Finding	Evidence
ok=True when NEO rejects its input: NEO logs to out.neo.run and exits 0, and the launcher creates output files before NEO starts	n_theta=16 → rejected by NEO, yet run_neo(check=True) did not raise
A rerun in the same directory returned the previous run's out.neo.* as its own	a good run, then a rejected one → the first run's <j.B> reported
kappa was not required; expro reads an absent tag as zero and NEO returns NaN with nothing logged	reg18 minus kappa → f_trap = nan, ok=True
torfluxa applied the psi storage factor to phi, which is always in Wb → 2π too large for a per-radian ODS	the repo's own phi convention (vaft/data/eqdsk.py)
COCOS 11 signs written verbatim, so NEO saw VEST's field and current both clockwise	expro_locsim.f90:202-203
bcentr read at time index 0 rather than at the converted slice	code path; it matters because VEST's b0 drifts (#325)

The sign fix flips NEO's normalised <j.B> and its b_unit together. NEO
normalises by a signed B_unit, so the physical current (their product) is
unchanged and stays parallel to B, as it must when Ip and Bt are parallel. That
is also why the scalar cross-checks could not catch the mirroring: helicity was
preserved. The VEST fixture is regenerated under the corrected convention, and
VAFT's Sauter and Redl still reproduce NEO's to 1e-8 on it. Phase 5 must apply
b_unit with its sign when dimensionalising.

🤖 Generated with Claude Code