Compare EQDSK->Boozer output to GPEC Boozer coordinates

[krystophny]

Migrated from archived itpplasma/efit_to_boozer#7. Cross-check libneo's Boozer coordinates against GPEC. Original: itpplasma/efit_to_boozer#7

[krystophny]

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

Still relevant. The EQDSK→Boozer path moved into libneo and grew a real test suite, but every existing check is self-consistency (structure, signs, Ampère's law); there is no cross-code validation of the Boozer angle itself. Spec below.

Purpose

efit_to_boozer builds Boozer coordinates from an axisymmetric g-file by field-line integration in symmetry-flux coordinates, then Fourier-transforms with a Jacobian weight (app/efit_to_boozer.f90:135-137). The resulting .bc files and chartmaps feed SIMPLE, NEO-2, NEO-RT and rabe, where a wrong θ_B origin, a wrong handedness, or a wrong Jacobian weight silently corrupts every resonance condition m θ_B + n φ_B. The only way to prove the transform is right is to compare against an independent implementation on the same g-file. GPEC/DCON's EQUIL module is a natural referee: it constructs flux coordinates from a g-file by a different route (Grad-Shafranov flux-surface reconstruction with a prescribed Jacobian selected by jac_type), and the group already runs it (itpplasma/rmp_torque).

Current state

libneo side. app/efit_to_boozer.f90 writes three .bc variants plus flux_functions.dat:

! app/efit_to_boozer.f90:135-137
      calE(it)=exp(cmplx(0.d0,theta+G/q, kind=cdp))
      aJb=R*(Br**2+Bp**2+Bz**2)/(C_norm*Bp)
      calEm_times_Jb(it)=cmplx(aJb,0.d0, kind=cdp)

so θ_B = θ_symflux + G/q, and aJb = R·|B|²/(C_norm·B_p) is the weight carried into the Fourier sums (it enters sqrtg00 at app/efit_to_boozer.f90:149 as calEm_times_Jb/B_oft**2). Confirming that this weight is the correct Boozer measure is part of what this issue is meant to establish — do not assume it. flux_functions.dat (app/efit_to_boozer.f90:239-243) carries R_beg, r, q, psi_pol, psi_tor_vac, psi_tor, i.e. the q profile and the poloidal↔toroidal flux map — exactly the two things GPEC also publishes. B amplitudes are written as Bmn_c(m)*1d-4 (app/efit_to_boozer.f90:216), confirming the internal field is CGS Gauss.

The Python front ends are python/libneo/eqdsk_to_boozer_chartmap.py:116 (convert_eqdsk_to_chartmap) and python/efit_to_boozer/boozer.py:100 (get_B0_of_s_theta_boozer(stor, num_theta)), which returns one periodic CubicSpline of |B|(θ_B) per surface — the ideal comparison object. Its |B| comes straight from _efit_to_boozer.efit_to_boozer.magdata (python/efit_to_boozer/boozer.py:116-135), so it is in Gauss, unconverted.

Tests today. python/tests/test_eqdsk_to_boozer_chartmap.py (143 lines) checks structure (nfp==1 at :56, nzeta==1 at :57), Bmod>0 (:65), plausible axis R (:73), |I/G| < 0.5 (:89), I(edge) ≈ 0.2·Ip in G·cm (Ampère, :102-105), and sign consistency of A_phi with torflux (:138). All of these pass for a transform with a shifted or reflected θ_B. Nothing pins the angle.

GPEC side, already in the org. itpplasma/rmp_torque has GPEC decks and ASCII readers:

• src/rmp_torque/templates/gpec/equil.in:6 — jac_type = 'boozer'. This is the only Boozer deck in the org. The template comment del.B ~ B_p**power_bp * B**power_b / R**power_r (:7, :9) documents that jac_type overrides the individual powers.
• run/AUG/33353_2900_EQH/GPEC/equil.in — a real, exercised EQUIL deck with out_eq_1d/out_eq_2d/out_2d = .true. (:38, :42, :46), mpsi = 128, mtheta = 256, psilow = 0.01, psihigh = 0.995. Note it runs jac_type = 'pest' (:6), not 'boozer' — useful as a working set of output/grid flags, but the coordinate choice must be changed for this issue.
• src/rmp_torque/gpec_ascii.py:158-216 read_equilibrium_grid() — parses idcon_equil.out: a Flux functions: table (6 columns; col 0 = ψ_N, col 1 = ψ_tor,N, col 3 = q) and a 2D functions: table (8 columns; col 0 = ψ_N, col 1 = θ, col 7 = |B|), with rectangularity, strict-monotonicity, θ-closes-at-2π and |B|-periodicity validation. Returns GpecEquilibriumGrid(psi, theta, q, b_magnitude) with the duplicated θ endpoint dropped.
• src/rmp_torque/gpec_ascii.py:290-321 read_flux_map() — returns GpecFluxMap(psi_normalised, psitor_normalised) only; it does not return q. q is available only from the col-3 interpolation inside read_equilibrium_grid.

So the reference reader exists but lives downstream. libneo is upstream of rmp_torque and must not depend on it — port, do not import.

Units are an open question, not a settled fact: gpec_ascii.py carries no unit annotation anywhere, and rmp_torque commits no idcon_equil.out fixture (find rmp_torque -name 'idcon*' is empty). DCON is conventionally SI, so col 7 is expected in Tesla, but the first thing the comparison must do is verify that against the libneo Gauss values rather than assume it — a silent 1e4 is exactly the kind of error this issue exists to catch.

Reproduction / evidence

Produce the libneo side today (nothing external needed). Note that python/tests/field_divB0.inp contains relative paths ('../../test/resources/input_efit_file.dat', '../../test/resources/convexwall.dat'), so it cannot simply be copied to an arbitrary directory — either run from python/tests or rewrite those two lines:

cd /path/to/libneo && make                        # cmake --preset default && cmake --build
mkdir -p /tmp/e2b && cd /tmp/e2b
LIBNEO=/path/to/libneo
cp "$LIBNEO/python/tests/efit_to_boozer.inp" .    # 3600/500/500/10000/2000/12/1e10
sed -e "s|'\.\./\.\./test/resources/|'$LIBNEO/test/resources/|g" \
    "$LIBNEO/python/tests/field_divB0.inp" > field_divB0.inp
"$LIBNEO/build/efit_to_boozer.x"
head -8 flux_functions.dat        # R_beg, r, q, psi_pol, psi_tor_vac, psi_tor
head -8 fromefit.bc

Observed: fromefit.bc, fromefit_neo_lhs.bc, fromefit_neo_rhs.bc, flux_functions.dat, box_size_axis.dat. Expected but missing: any statement of how q(s) and |B|(θ_B) compare to a second code. rg -ril gpec over libneo hits only matlab/GPEC_interface/, whose last commit is 2020-06-25 ("Modified templates for GPEC v1.4") — no reader, no comparison, no fixture.

Verification recipe for the finished feature, on the repo's own fixture test/resources/input_efit_file.dat:

1. Run GPEC EQUIL with eq_type='efit', eq_filename=<input_efit_file.dat>, jac_type='boozer' (the run deck ships 'pest' — this must be changed), mpsi=64, mtheta=128, psilow=0.01, psihigh=0.98, out_eq_2d=.true., out_2d=.true.. Start from rmp_torque/run/AUG/33353_2900_EQH/GPEC/equil.in and edit those keys. Run via dcon (rmp_torque/run/run_gpec.sh: module load gpec && dcon && gpec && pentrc; only the dcon stage is needed here). This yields idcon_equil.out.
2. Confirm the |B| unit first: compare mean(|B|) on a mid-radius surface between GPEC col 7 and libneo (Gauss) and assert the ratio is within 10% of either 1 or 1e4; record which. Everything below assumes GPEC Tesla, libneo Gauss ×1e-4.
3. q(s_tor) from flux_functions.dat col 3 vs the GPEC flux table col 3, both interpolated onto s_tor ∈ [0.1, 0.9]: target agreement 1e-3 relative.
4. |B|(θ_B) from get_B0_of_s_theta_boozer vs GPEC 2D col 7 at s_tor = 0.25, 0.50, 0.75, after fitting one global (Δθ, σ) with σ = ±1 for handedness: target max|ΔB|/⟨B⟩ < 2e-3 at s ≤ 0.75.
5. The fitted Δθ must be the same on all three surfaces to within 5e-3 rad — an s-dependent Δθ means the two codes are not using the same poloidal angle at all, only the same |B| envelope. This is the load-bearing check; it is the one that a shifted-angle bug cannot survive.

The numeric tolerances in steps 3-4 are targets chosen a priori, not measured values — nobody has run this comparison yet. On the first real run, record the achieved numbers in the fixture README and in doc/ExtraDocuments/Conventions.md, then set the test thresholds to those numbers with modest headroom. Do not silently loosen a tolerance to make the test pass; a miss is the finding.

Behaviour (BDD)

Scenario: GPEC idcon_equil.out is parsed into flux functions and a Boozer-angle B grid
  Given the fixture test/resources/gpec/idcon_equil_input_efit_file.out
  When read_gpec_equilibrium_grid() is called on it
  Then the returned psi grid is strictly increasing and len(theta) == mtheta
  And theta spans [0, 2*pi) after the duplicated endpoint is dropped
  And b_magnitude has shape (len(psi), len(theta)) and is strictly positive

Scenario: a non-rectangular 2D table is rejected
  Given an idcon_equil.out whose "2D functions:" block has one row deleted
  When read_gpec_equilibrium_grid() is called on it
  Then ValueError is raised with a message containing "non-rectangular"

Scenario: the |B| unit is identified before any comparison
  Given libneo |B| in Gauss at s_tor = 0.50 and the GPEC |B| row at the nearest surface
  When the ratio mean(B_libneo) / mean(B_gpec) is computed
  Then that ratio is within 10 percent of 1.0e4
  And the test reports the detected unit explicitly

Scenario: q profiles from libneo and GPEC agree on the same g-file
  Given flux_functions.dat produced by efit_to_boozer.x on test/resources/input_efit_file.dat
  And the GPEC flux table from the same g-file with jac_type='boozer'
  When both q profiles are interpolated onto s_tor = linspace(0.1, 0.9, 40)
  Then max relative difference in q is below 1.0e-3

Scenario: |B| on the Boozer angle grid agrees after a single rigid angle alignment
  Given the libneo |B|(theta_B) splines from get_B0_of_s_theta_boozer at s_tor in {0.25, 0.50, 0.75}
  And the GPEC |B| rows at the nearest GPEC surfaces
  When a single (delta_theta, sign) pair is fitted by minimising the L2 residual at s_tor = 0.50
  And that same (delta_theta, sign) is applied at s_tor = 0.25 and 0.75
  Then max|B_libneo - B_gpec| / mean(B_gpec) < 2.0e-3 on every surface
  And the fitted delta_theta varies by less than 5.0e-3 rad across the three surfaces

Scenario: Boozer harmonic magnitudes agree without any angle alignment
  Given |B|(theta_B) from both codes at s_tor = 0.50
  When |b_m| = |(1/2pi) * oint B exp(-i m theta) dtheta| is evaluated for m = 0..6
  Then |b_0| agrees to 1.0e-4 relative and |b_1|..|b_4| agree to 2.0e-3 relative

Scenario: the comparison is skipped, not failed, without external data
  Given LIBNEO_GPEC_FIXTURE is unset and the committed fixture is absent
  When pytest runs python/tests/test_gpec_boozer_comparison.py
  Then the test reports skipped and the suite exit status is 0

Implementation plan

1. Create python/libneo/gpec_equil.py. Port the parsing and validation of read_equilibrium_grid (rmp_torque/src/rmp_torque/gpec_ascii.py:158-216) and read_flux_map (:290-321), keeping every check: strict ψ and θ monotonicity, len(grid) == len(psi)*len(theta), θ closing at 2π to atol=1e-7, |B| periodic across the seam. Expose frozen dataclasses GpecFluxFunctions(psi_poloidal, psi_toroidal, q) — note read_flux_map upstream returns only the two flux columns, so take q from flux col 3 as read_equilibrium_grid does — and GpecBoozerGrid(psi_poloidal, psi_toroidal, theta, q, b_magnitude), the second resolving ψ_tor per 2-D surface by interpolating the flux table so consumers work in s_tor directly. Module docstring states provenance (itpplasma/rmp_torque), the column layout, and that the |B| unit is unannotated upstream and is checked at comparison time. numpy only, no new dependency.
2. Create python/libneo/compare_gpec_boozer.py. Pure functions, no I/O: harmonic_magnitudes(theta, values, mmax); fit_rigid_angle_alignment(theta_ref, b_ref, spline_test) returning (delta_theta, sign, residual) by brute-force scan over 360 offsets × 2 signs refined with scipy.optimize.minimize_scalar; compare_profiles(s, q_a, q_b) returning max relative difference; detect_field_unit_scale(b_libneo, b_gpec) returning the nearest of {1.0, 1.0e4} plus the raw ratio. Keeping the numerics here makes them unit-testable without any GPEC data.
3. Create python/tests/test_gpec_equil_reader.py. Unit-test step 1 against a small synthetic idcon_equil.out written inside the test, plus negative cases: non-rectangular grid, non-monotone ψ, θ not closing at 2π, non-periodic |B|. The "Flux functions:\n" … "2D functions:\n" string literal at rmp_torque/tests/test_gpec.py:96-102 shows the section headers and the 6-column flux row format, but it stops at the 2D functions: header with no grid rows, so it does not exercise the 2-D path — the 8-column rows must be written fresh here (≥3 ψ × ≥3 θ so the rectangularity check is not vacuous).
4. Add fixture test/resources/gpec/. Commit equil.in (the deck from step 1 of the verification recipe), README.md recording GPEC version, host, date, exact command and the measured (Δθ, σ, unit scale), and a trimmed idcon_equil_input_efit_file.out reduced to ~16 surfaces × 128 angles so CI can run without GPEC. Trim by whole-surface row selection only — do not resample, or the reader's rectangularity check becomes vacuous. There is no committed idcon_equil.out anywhere in the org to copy, so this file must come from a real GPEC run.
5. Create python/tests/test_gpec_boozer_comparison.py. Implement the scenarios above. Guard with pytest.importorskip("scipy") and a try: import _efit_to_boozer / except ImportError: pytest.skip(...) block (the pattern at python/tests/test_eqdsk_to_boozer_chartmap.py:27-31, alongside the pytest.importorskip("netCDF4") at :24), plus an os.environ.get("LIBNEO_GPEC_FIXTURE") override pointing at a full GPEC run directory (mirroring GPEC_ITER_FIXTURE at rmp_torque/tests/test_gpec.py:142-144). Run efit_to_boozer via libneo.eqdsk_to_boozer_chartmap._write_inp (python/libneo/eqdsk_to_boozer_chartmap.py:50) and _write_field_divB0_inp (:68) inside a tmp_path so cwd is not polluted and the relative-path problem above does not recur, then call efit_to_boozer.boozer.get_B0_of_s_theta_boozer.
6. Extend doc/ExtraDocuments/Conventions.md. Add a GPEC/DCON row to the kernel table (:54-58) and the sine/cosine table (:71-75), and a short subsection: flux label ψ_N (normalised poloidal flux) vs the ITP s = Ψ_tor/Ψ_tor(sep), the jac_type → Jacobian mapping, and the measured (Δθ, σ) between GPEC and libneo once step 5 reports it. This is the durable output of the issue — the number belongs in the docs, not only in a test.
7. CI wiring. .github/workflows/main.yml:101 already runs python -m pytest python/tests -vv -s; the new tests are picked up automatically. No new package is needed in pyproject.toml — scipy is already in both the dev and docs extras.
8. If a systematic discrepancy is found, fix in app/efit_to_boozer.f90 (the weight/θ_B definition at :135-137) or src/efit_to_boozer/spline_and_interpolate_magdata.f90, and record the resolved convention in efit_to_boozer: tested, clean version with clear COCOS convention #360, which owns the COCOS statement for this path. Do not change .bc output signs without re-running the downstream gate — ci/downstreams gates SIMPLE, NEO-2, MEPHIT, KAMEL and rabe as blocking, with NEO-RT report-only.

Acceptance criteria

• python -m pytest python/tests/test_gpec_equil_reader.py -q passes with no external data, and each negative case raises ValueError.
• python -m pytest python/tests/test_gpec_boozer_comparison.py -q passes against the committed trimmed fixture, and skips cleanly (exit 0) when _efit_to_boozer is not built or the fixture is absent.
• The detected GPEC/libneo |B| unit scale is asserted and printed; it is 1e4 (GPEC Tesla, libneo Gauss) or the discrepancy is documented.
• q agreement ≤ 1e-3 relative on s_tor ∈ [0.1, 0.9] against GPEC on test/resources/input_efit_file.dat.
• |B|(θ_B) agreement ≤ 2e-3 of mean |B| at s_tor = 0.25, 0.50, 0.75 under one rigid (Δθ, σ).
• The fitted Δθ is s-independent to 5e-3 rad; the fitted σ is reported explicitly in the test output.
• |b_0| agrees to 1e-4 relative, |b_1..4| to 2e-3 relative, with no angle alignment applied.
• The achieved numbers for the four criteria above are written into test/resources/gpec/README.md; any threshold that had to be loosened from the target is called out in the PR description with a physical reason.
• doc/ExtraDocuments/Conventions.md contains a GPEC/DCON row in both tables and the measured (Δθ, σ).
• make && make test green; python -m pytest python/tests test/python -q shows no new failures; downstream gate green.

References

• A. H. Boozer, Plasma equilibrium with rational magnetic surfaces, Phys. Fluids 24, 1999-2003 (1981) — the Boozer-coordinate construction; the Boozer Jacobian is ∝ 1/B².
• A. H. Boozer, Establishment of magnetic coordinates for a given magnetic field, Phys. Fluids 25, 520-521 (1982), doi:10.1063/1.863765.
• J.-K. Park, A. H. Boozer, A. H. Glasser, Computation of three-dimensional tokamak and spherical torus equilibria, Phys. Plasmas 14, 052110 (2007), doi:10.1063/1.2732170 — the IPEC/GPEC formulation.
• A. H. Glasser, The direct criterion of Newcomb for the ideal MHD stability of an axisymmetric toroidal plasma, Phys. Plasmas 23, 072505 (2016), doi:10.1063/1.4958328 — DCON, which owns equil.in and idcon_equil.out.
• O. Sauter, S. Yu. Medvedev, Tokamak coordinate conventions: COCOS, Comput. Phys. Commun. 184, 293-302 (2013), doi:10.1016/j.cpc.2012.09.010 — the convention framework efit_to_boozer: tested, clean version with clear COCOS convention #360 must pin down before the sign comparison is meaningful.
• GPEC output documentation: http://princetonuniversity.github.io/GPEC/outputs.html — the authoritative statement of the idcon_equil.out column layout and units; check it before trusting the ported column indices, which were inferred from rmp_torque code rather than from a spec.
• booz_xform (hiddenSymmetries/booz_xform) covers the VMEC→Boozer route and is already wired into libneo via python/libneo/booz_xform_to_boozer_chartmap.py and python/libneo/bc_to_booz_xform.py — a second, independent angle reference if the GPEC run stays blocked.

Scope and risk

Two PRs: (1) reader + alignment helpers + unit tests, which needs no external data and can land immediately; (2) the committed GPEC fixture plus the comparison test and tolerances, which is blocked on someone with a GPEC installation running EQUIL with jac_type='boozer' on test/resources/input_efit_file.dat. That dependency is the main schedule risk — nothing in libneo or rmp_torque contains a usable idcon_equil.out today. Code risk is confined to new Python files and one doc file; no Fortran changes unless step 8 fires, in which case any .bc sign change must go through the downstream gate.

Related: #360 (clean, COCOS-explicit EQDSK→Boozer path — this issue is its external validation and should land after it), #363 (automatic COCOS-3 check on EQDSK input), #309 and itpplasma/rabe#92 (the analogous VMEC→Boozer benchmark gate, where spline-vs-spectral evaluation of |B| extrema was the failure mode — reuse the lesson: evaluate |B| spectrally, not from splines, when tightening tolerances).