Track .bc retirement: verify chartmap and booz_xform conventions

[krystophny]

Depends on: itpplasma/libneo#343, itpplasma/libneo#344, itpplasma/libneo#345

Goal

SIMPLE already reads Boozer chartmaps (src/field/field_boozer_chartmap.f90). Verify that chartmap and booz_xform inputs produced from the same source data preserve SIMPLE's .bc observables within documented converter/interpolation tolerance, then deprecate .bc input.

Context

Shared converters live in libneo: EQDSK -> Boozer chartmap (itpplasma/libneo#343), .bc -> booz_xform (itpplasma/libneo#344), and a same-source converter/sign-convention gate (itpplasma/libneo#345). Prefer routing through libneo's existing chartmap/booz_xform reader over writing a new parser. Fortran codes can wrap the libneo reader via f90wrap or read the NetCDF directly; Python reads it directly.

Different Boozer transforms need not match to floating-point accuracy. This issue requires consistent signs, mode indices, and field values for converted inputs that originate from the same equilibrium data.

What to do

• Add a reference tokamak case run from a chartmap and from a booz_xform-derived chartmap; assert agreement with the .bc-based result to a documented tolerance appropriate for the same-source conversion path.
• Include convention checks for booz_xform cos(m theta - n zeta) / sin(m theta - n zeta) harmonics and asymmetric coefficients.
• Once green and the Same-source Boozer convention test for .bc retirement libneo#345 gate passes, mark .bc input deprecated.

Non-goals

• Removing .bc before the same-source convention/tolerance gate is green.
• Requiring independently generated Boozer coordinates to match .bc results to floating-point accuracy.

[krystophny]

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

Still open. The libneo converters this issue depends on exist on libneo main; SIMPLE's side does not. There is a first cut on the unmerged branch origin/test/bc-retirement (96dc56f, "test: .bc retirement field-level gate (#430)") — it is not an ancestor of origin/main, so it closes nothing, and as written it cannot observe the conventions this issue exists to gate. Spec below.

Purpose

.bc (Strumberger/Boozer-coefficient) files are the legacy interchange format across ITPplasma codes. SIMPLE itself no longer reads them — field_from_file (src/field.F90:69-140) accepts SPECTRE files, .nc (VMEC wout / Boozer chartmap), and coils*, and error stops otherwise; grep -rn '\.bc' src/ docs/ DOC/ returns nothing. So this is not a code-removal task in SIMPLE but a trust task: before the org drops .bc, show that a chartmap derived from a .bc through the libneo converters reproduces the same geometry and field as the .bc itself, with signs, mode indices and radial staggering correct. A wrong toroidal-angle sign silently mirrors the equilibrium and yields plausible-looking but wrong loss fractions for non-axisymmetric configurations.

Current state

libneo — converters and producer-side sign gates are done.

• python/libneo/bc_to_booz_xform.py — .bc -> boozmn NetCDF (introduced in cdfb006; off-by-one radial-grid fix in 0bc2d60).
• python/libneo/booz_xform_to_boozer_chartmap.py — boozmn -> chartmap.
• python/libneo/eqdsk_to_boozer_chartmap.py (introduced in cdfb006; covariant B_theta fix in fe3f76e, Correct integmode mode list in configs and docs #399).
• Asymmetric phase coefficients preserved in boozmn: f43f491 (Structure-preserving collision integrator to complement the symplectic orbit push #357). Shared Fortran boozmn_reader: 94fd53b (Add coil-clearance geometry helpers #364).

All five commits are ancestors of libneo origin/main. (The dependency issues in the body are quoted as #343/#344/#345; none of these commit subjects carries those numbers, so treat the issue-number mapping as unconfirmed and check it on GitHub before closing anything by reference.)

Critically — and this changes the plan the obvious way — python/tests/test_bc_to_booz_xform.py already gates the mode-index and harmonic-sign conventions on non-axisymmetric fixtures:

• test_symmetric_bc_writes_modes_with_booz_xform_signs (:167) — synthetic nper=4 .bc, asserts ixn == n*nfp and that bmnc/pmns reproduce cos/sin(m*theta - ixn*zeta).
• test_asymmetric_bc_writes_phase_coefficients_with_booz_xform_signs (:228) — synthetic nper=3 lasym .bc, asserts lasym__logical__ == 1, pmns == vmns*2pi/nfp, pmnc == vmnc*2pi/nfp, and the full asymmetric sum.
• test_asymmetric_vmec_bc_to_boozmn_matches_booz_xform_mode_signs (:312) — downloads an asymmetric VMEC wout, runs the real booz_xform library and libneo's converter on it, and asserts the two mode sets (ixm, ixn) agree with no sign-flipped modes.

Do not rewrite these. Two of the four scenarios a naive reading of this issue suggests are already covered.

The one convention that is genuinely unpinned is the sign of p. .bc defines its phase column, verbatim from the fixture header, as (phib-phi)*nper/twopi. libneo carries it straight through:

python/libneo/bc_to_booz_xform.py:58-59

# vmns in .bc = (phib-phi)*nfp/(2*pi); pmns_b = phib-phi in radians
pmns[k] = np.array(bc.vmns[k], dtype=float) * TWOPI / bc.nper

so the Fourier sum of pmns_b equals phi_B - phi_cyl. But the consumer reconstructs the cylindrical angle with the opposite sign, python/libneo/booz_xform_to_boozer_chartmap.py:245:

phi_cyl = zeta_geom[None, None, :] + p

Given the .bc column definition, exactly one of these two lines is wrong: either :58 needs a minus (if booz_xform's pmns_b is defined as phi_cyl - zeta_B), or :245 needs a minus (if it is defined as zeta_B - phi_cyl). The existing libneo cross-check against the real booz_xform library compares mode sets only — it never compares pmns_b values — so nothing in either repo currently decides this. That comparison is the missing gate, and it is a small addition to a test that already has both files in hand.

SIMPLE — nothing exercises .bc. Chartmap coverage on main stops at the VMEC source: test/tests/test_booz_xform_chartmap.py validates wout -> boozmn -> chartmap in three layers (spectral convergence, golden cross-check against export_boozer_chartmap_tool.x, end-to-end confined fractions with tol = 1.0/npart at :341), registered at test/tests/CMakeLists.txt:611-619 under if(SIMPLE_HAVE_BOOZ_XFORM).

SIMPLE vendors its own diverged copy of the converter. tools/booz_xform_to_boozer_chartmap.py (325 lines) is a fork of libneo's python/libneo/booz_xform_to_boozer_chartmap.py (309 lines), and test_booz_xform_chartmap.py:44 imports it via sys.path.insert(0, str(TOOLS_DIR)), not from libneo. Both copies have the same + p reconstruction (libneo :245, SIMPLE tools/…:271), so the hazard is present in both. The issue body says to prefer routing through libneo rather than a second parser, so decide explicitly which copy the new gate binds; if the vendored fork stays, it needs the same fix.

Branch origin/test/bc-retirement (96dc56f) adds test/tests/test_bc_retirement.py (4 pytest functions), fixtures test/test_data/{circ.bc,boozmn_circ.nc,chartmap_circ.nc}, three .gitignore allowlist entries, and an add_test(NAME test_bc_retirement …) block with LABELS "unit;python", TIMEOUT 120. Two structural limits:

1. It compares only Bmod (1e-4 relative) and the m=0 amplitude (1e-2), both against a Python re-summation of the .bc harmonics. No phi_cyl, no geometry, no SIMPLE run — so neither the sign question above nor any orbit observable is gated.
2. The fixture is axisymmetric. circ.bc's header row m0b n0b nsurf nper is 18 0 63 1. Verified across every .bc in the org checkouts — libneo/python/tests/circ.bc, NEO-RT/examples/circ.bc, NEO-2/DOC/example_{bg,pert}_field.bc, BOOZER_MAGFIE/tokamak_circ/tok-synch2-n0.bc, rmp_torque/tests/fixtures/{iter,mastu,demo}/equilibrium.bc — all eight have n0b = 0 and nper = 1. With every n = 0 the toroidal-phase sign is multiplied by zero and is unobservable; the test passes either way.

The fix for (2) is not to hand-write a stellarator .bc. libneo.boozer.BoozerFile can produce one from a wout (convert_vmec_to_boozer then write, python/libneo/boozer.py:803), which is exactly what libneo's own test does at python/tests/test_bc_to_booz_xform.py:337-341. SIMPLE's test_booz_xform_chartmap.py already downloads a non-axisymmetric QA wout (test/test_data/wout.nc, WOUT_URL at :38). Round-tripping that wout to .bc and back gives a non-axisymmetric fixture and a same-source reference run, which is precisely what the issue body asks for.

Separately, libneo/src/field/bc_file.f90:95 hard-errors on non-stellarator-symmetric .bc ("non-stellarator-symmetric .bc files (10 columns) are not supported") while the Python path handles lasym (bc_to_booz_xform.py:50-79). Any asymmetric check must therefore go through Python today.

Reproduction / evidence

Today, on main:

cd /path/to/SIMPLE
grep -rn '\.bc' src/ docs/ DOC/          # no hits: SIMPLE has no .bc input path
grep -rln '\.bc' test/                   # no hits: no test takes a .bc
make CONFIG=Fast build
make CONFIG=Fast test TEST=test_booz_xform_chartmap   # passes; source is wout.nc, never a .bc

Observed: no test in SIMPLE takes a .bc as input. Expected after this issue: make CONFIG=Fast test TEST=test_bc_retirement gates the .bc -> boozmn -> chartmap pipeline including the toroidal-phase sign, on a non-axisymmetric case.

Minimal reproduction of the open sign question (libneo checkout, no SIMPLE build needed):

cd /path/to/libneo/python
python - <<'PY'
import numpy as np, netCDF4
from libneo.boozer import BoozerFile
from libneo.bc_to_booz_xform import convert_bc_to_boozmn
from booz_xform import Booz_xform

wout = "wout.nc"   # any non-axisymmetric wout, e.g. the QA one SIMPLE downloads

bx = Booz_xform(); bx.read_wout(wout, True); bx.mboz = 6; bx.nboz = 8
bx.run(); bx.write_boozmn("boozmn_ref.nc")

bc = BoozerFile(filename=""); bc.convert_vmec_to_boozer(wout, uv_grid_multiplicator=1)
bc.write("round.bc"); convert_bc_to_boozmn("round.bc", "boozmn_libneo.nc")

def pmns(f):
    with netCDF4.Dataset(f) as ds:
        return (np.asarray(ds.variables["ixm_b"][:]),
                np.asarray(ds.variables["ixn_b"][:]),
                np.asarray(ds.variables["pmns_b"][:]))

mr, nr, pr = pmns("boozmn_ref.nc")
ml, nl, pl = pmns("boozmn_libneo.nc")
ref = {(m, n): k for k, (m, n) in enumerate(zip(mr, nr))}
for k, (m, n) in enumerate(zip(ml, nl)):
    if (m, n) in ref and abs(pr[-1][ref[(m, n)]]) > 1e-8:
        print(m, n, "ratio", pl[-1][k] / pr[-1][ref[(m, n)]])
PY

A ratio of -1 means bc_to_booz_xform.py:58 has the wrong sign; +1 means booz_xform_to_boozer_chartmap.py:245 does. Record whichever it is; it decides step 1.

Behaviour (BDD)

Scenario: pmns sign is pinned against the booz_xform library
  Given a non-axisymmetric VMEC wout with nfp > 1
  And boozmn_ref.nc written by booz_xform.Booz_xform.write_boozmn
  And boozmn_libneo.nc written by libneo convert_bc_to_boozmn from the .bc
      produced by BoozerFile.convert_vmec_to_boozer(wout).write()
  When pmns_b is compared mode by mode on the outermost common surface
       for every mode where |pmns_b(ref)| > 1e-8
  Then the ratio pmns_libneo/pmns_ref is +1 within 1e-6 for every such mode
  And no mode has ratio -1

Scenario: cylindrical angle round-trips through the chartmap
  Given the same wout and a chartmap produced from the .bc-derived boozmn
        with nrho=64, ntheta=64, nzeta=64
  When phi_cyl = arctan2(y, x) is read from the chartmap at 200 random grid
       points and compared with phib - v evaluated directly from the .bc
       harmonics at the same (s, theta_B, zeta_B)
  Then the maximum absolute difference is below 1e-6 rad
  And inverting the sign at booz_xform_to_boozer_chartmap.py:245 makes this
      scenario fail with a difference above 1e-3 rad

Scenario: bc-sourced chartmap reproduces Bmod from the .bc harmonics
  Given test/test_data/circ.bc (m0b=18, n0b=0, nsurf=63, nper=1)
  And chartmap_circ.nc built via bc_to_booz_xform then
      booz_xform_to_boozer_chartmap with nrho=64, ntheta=128, nzeta=1
  When Bmod is sampled at 30 random (s, theta_B) with 0.20 < s < 0.90
  Then the relative deviation from direct .bc Fourier summation is below 1e-4
  And min(Bmod) over the whole chartmap is greater than 0

Scenario: orbits agree between VMEC-sourced and bc-sourced chartmaps
  Given test/test_data/wout.nc and the .bc round-tripped from it
  When SIMPLE is run with ntestpart=32, trace_time=1d-2, isw_field_type=2,
       integmode=1, deterministic=.True., sbeg=0.3d0, once with
       netcdffile='wout.nc' and once with
       field_input=coord_input='boozer_chartmap.nc' built from that .bc
  Then max |total_conf(vmec) - total_conf(chartmap)| over the time grid is
       at most 1/32
  And the same holds for pass_conf and trap_conf

Implementation plan

1. libneo — decide the p sign (blocking PR). Run the MRE above. Then extend the existing python/tests/test_bc_to_booz_xform.py:312 test (test_asymmetric_vmec_bc_to_boozmn_matches_booz_xform_mode_signs) — it already generates both boozmn files and _read_boozmn (:60) already loads pmns — with the per-mode pmns_b ratio assertion of scenario 1. Whichever of bc_to_booz_xform.py:58 or booz_xform_to_boozer_chartmap.py:245 is wrong, fix it, and state the resolved convention explicitly in the bc_to_booz_xform.py module docstring (:1-:21, next to the existing .bc vmnc/vmns (phase) = (phib-phi)*nper/(2*pi) line) and in booz_xform_to_boozer_chartmap.py above :245. Do not add a new test file and do not re-test ixn, lasym, or the cos(m theta - n zeta) harmonic signs — :167, :228 and :312 already cover those on nper=3/nper=4 fixtures.
2. libneo — geometric round-trip. Add one test asserting scenario 2 in the same file: convert to a chartmap and compare arctan2(y, x) against direct .bc evaluation. Include the mutation check (flip the sign, assert the test fails) in the PR description, not in the committed test.
3. SIMPLE — land the field-level gate. Cherry-pick 96dc56f from origin/test/bc-retirement: test/tests/test_bc_retirement.py, fixtures test/test_data/{circ.bc,boozmn_circ.nc,chartmap_circ.nc}, and the three .gitignore allowlist entries. The branch's CMake hunk targets an older test/tests/CMakeLists.txt; on current main place the add_test(NAME test_bc_retirement …) block after the endif() at line 619. Keep LABELS "unit;python" and TIMEOUT 120 so it runs in make test-fast (which excludes only slow|regression|performance|scalability). Keep the 1e-4 Bmod and 1e-2 m=0 tolerances. Regenerate the committed fixtures after step 1 if the sign changed.
4. SIMPLE — add the non-axisymmetric orbit gate. Add test/tests/test_bc_retirement_orbits.py, modelled on check_e2e_orbits in test/tests/test_booz_xform_chartmap.py:299. Reuse confined_metrics, load_confined_fraction, run_cmd, download_if_missing from test/tests/boozer_chartmap_artifacts.py (defined at :59, :52, :26). Reuse the QA wout via WOUT/WOUT_URL (test_booz_xform_chartmap.py:37-41); assert nfp > 1 on it so the test fails loudly if the fixture is ever swapped for an axisymmetric one. Build the .bc with libneo.boozer.BoozerFile.convert_vmec_to_boozer + .write(), convert with python -m libneo.bc_to_booz_xform and python -m libneo.booz_xform_to_boozer_chartmap … --nrho --ntheta --nzeta, and use the two namelists at test_booz_xform_chartmap.py:254 (SIMPLE_IN_VMEC) and :271 (SIMPLE_IN_CHARTMAP) unchanged as templates. Tolerance 1.0/npart on the confined fractions, matching :341.
5. SIMPLE — register the slow test. Add add_test(NAME test_bc_retirement_orbits …) with LABELS "slow;python" and TIMEOUT 7200, inside the existing if(SIMPLE_HAVE_BOOZ_XFORM) block at test/tests/CMakeLists.txt:611-619, using ${BOOZER_CHARTMAP_PYTHON} like its neighbours.
6. SIMPLE — resolve the vendored converter. Either delete tools/booz_xform_to_boozer_chartmap.py and repoint test_booz_xform_chartmap.py:44 at libneo.booz_xform_to_boozer_chartmap, or apply the step-1 sign fix to tools/…:271 as well. Do not leave the two copies disagreeing.
7. SIMPLE — document the retirement. Add a "Legacy .bc input" subsection to docs/boozer-chartmap-schema.rst, next to the angle-reversal table at lines 44-59, stating that SIMPLE reads no .bc, that the supported route is .bc -> bc_to_booz_xform -> booz_xform_to_boozer_chartmap -> field_input/coord_input, and spelling out the p/v sign convention resolved in step 1.

Ordering: libneo PR (steps 1-2) merges first; SIMPLE PR (steps 3-7) pins it with make LIBNEO_REF=<sha>.

Acceptance criteria

• libneo pytest python/tests/test_bc_to_booz_xform.py passes, including a new per-mode pmns_b comparison against booz_xform's own output, and that assertion demonstrably fails when the sign at bc_to_booz_xform.py:58 is inverted.
• The resolved p sign is stated in both libneo module docstrings, and bc_to_booz_xform.py:58 and booz_xform_to_boozer_chartmap.py:245 are mutually consistent with the .bc (phib-phi) column definition.
• make CONFIG=Fast test TEST=test_bc_retirement passes and completes under 120 s.
• make CONFIG=Fast test TEST=test_bc_retirement_orbits passes when booz_xform is installed, and its fixture asserts nfp > 1.
• Bmod from the .bc-derived chartmap matches direct .bc Fourier summation to below 1e-4 relative for 0.20 < s < 0.90.
• max |confined_fraction| difference between the .bc-derived chartmap run and the VMEC-direct run of the same equilibrium is at most 1/ntestpart for total_conf, pass_conf and trap_conf.
• tools/booz_xform_to_boozer_chartmap.py is either removed or carries the same sign convention as libneo's copy.
• make test-regression shows no change in the boozer, boozer_ncsx and test_tokamak golden records.
• docs/boozer-chartmap-schema.rst states the .bc phase and v-sign convention explicitly.

References

• Boozer, A. H., Plasma equilibrium with rational magnetic surfaces, Phys. Fluids 24, 1999 (1981) — definition of Boozer coordinates.
• Albert, C. G., Kasilov, S. V., Kernbichler, W., Symplectic integration with non-canonical quadrature for guiding-center orbits in magnetic confinement devices, J. Comput. Phys. 403, 109065 (2020) — SIMPLE's integrators, i.e. what the confined-fraction observable in step 4 exercises.
• booz_xform (https://github.com/hiddenSymmetries/booz_xform) is the reference implementation for the boozmn layout used here. Its variables ixm_b, ixn_b (with n pre-multiplied by nfp), bmnc_b and pmns_b are consumed by libneo/python/libneo/bc_to_booz_xform.py and read back in python/tests/test_bc_to_booz_xform.py:60; the ixn_b = n*nfp convention is already pinned against the library itself at :312. Its pmns_b sign is the quantity step 1 must pin — read it from the library's output, not from prose.
• The rest of this issue is an implementation detail internal to the ITPplasma codes; no specific literature reference applies to the .bc phase-column convention, which is fixed by the (phib-phi)*nper/twopi header in the .bc files themselves.

Scope and risk

Two PRs (libneo sign resolution + round-trip test, then SIMPLE gate + docs), roughly size/M each. The main risk is step 1: if the pmns sign is wrong today, fixing it changes every chartmap produced from a .bc through libneo, and the committed boozmn_circ.nc/chartmap_circ.nc fixtures on origin/test/bc-retirement must be regenerated before the cherry-pick in step 3. Run libneo's downstream-gate workflow (.github/workflows/downstream-gate.yml) and SIMPLE's make test-smoke before merging. Note also that the axisymmetric fixtures cannot detect the change at all, which is exactly why step 4 must use a non-axisymmetric case.

Open question for the issue author, not resolvable from the checkouts: whether .bc stays a supported source format for NEO-2 / NEO-RT / rmp_torque while SIMPLE drops it. All three still ship .bc fixtures.