VMEC toroidal-flux sign is hardcoded for left-handed (signgs=-1) equilibria across the ingestion path; positive-Jacobian VMEC is silently mis-signed

[krystophny]

Problem

The VMEC ingestion path assumes a left-handed equilibrium (signgs = -1) and hardcodes torflux = -phi_vmec[-1]/(2*pi) without ever reading signgs from the wout/boozmn file. For a positive-Jacobian / right-handed equilibrium (signgs = +1), the stored toroidal flux has the opposite sign, so the hardcoded negation yields a toroidal flux of the wrong sign, which mirror-flips iota and A_phi relative to the written geometry. No test covers a signgs = +1 equilibrium.

This is a codebase-wide convention rather than a defect specific to one file: the booz_xform chartmap converter faithfully mirrors libneo's own native VMEC reader, which makes the same assumption.

Where

• libneo:src/vmecinm_m.f90 lines 75-86 (phi = -phi/(2*pi), comment ties this to left-handed VMEC / signgs = -1).
• libneo:src/spline_vmec_data.f90 (compute_boozer_data, sA_phi construction).
• SIMPLE:tools/booz_xform_to_boozer_chartmap.py lines 211-218 (torflux_si = -float(d["phi"][-1])/TWOPI, comment asserts phi_b[-1] > 0 without checking signgs).
• libneo:doc/ExtraDocuments/Conventions.md (toroidal-flux sign convention).

Evidence

• vmecinm_m.f90 applies phi = -phi/(2*pi); the preceding comment explicitly ties the negation to left-handed VMEC with signgs = -1.
• The booz_xform converter applies the same negation (-phi_b[-1]/(2*pi)) with a comment asserting phi_b[-1] > 0, again without inspecting signgs.
• signgs is not read from any wout/boozmn file anywhere in the ingestion path; it appears only in explanatory comments. A signgs = +1 wout fed through the native path would be mis-signed identically to the converter.
• VMEC has supported positive-Jacobian equilibria (wout writes signgs) since v6.20, so signgs = +1 inputs are possible in practice.

Why it matters

When format conversion is consolidated into a single libneo Boozer representation, the in-memory model should be unambiguous about the toroidal-flux sign convention. A positive-Jacobian VMEC equilibrium currently passes through silently with a flipped flux sign, producing a self-consistent but mirror-reflected field. Surfacing this as a documented, checked limitation (read and validate signgs, or reject unsupported sign conventions at the boundary) avoids silent mis-signing in the consolidated path.

Do not fix yet - tracking only

This is a latent, codebase-wide limitation recorded for triage during converter consolidation, not a converter-specific regression. Do not implement a fix in response to this issue; it is a tracking record only.

[krystophny]

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

Still valid at fe3f76e. Nothing has changed since the issue was filed, and the affected surface has grown: the booz_xform converter that the issue attributes to SIMPLE now also lives inside libneo (python/libneo/booz_xform_to_boozer_chartmap.py, src/field/boozmn_reader.f90), so the same unconditional negation exists in three places in this repo. The issue says "tracking only"; the spec below is what to implement when converter consolidation lands.

Purpose

VMEC's signgs is the sign of the Jacobian of (s, theta, zeta); it is -1 for the usual VMEC output and +1 for a positive-Jacobian equilibrium. libneo currently folds the assumption signgs = -1 into a hardcoded sign flip of the toroidal flux and never reads signgs, so a signgs = +1 wout produces a field whose direction is reversed relative to the geometry it was built from. Because the result is internally self-consistent (|B|, flux surfaces and |iota| all still look right), nothing downstream will flag it; only the sign of parallel motion and of any direction-sensitive quantity silently inverts.

Current state

src/vmec_support/vmecinm_m.f90:77-82 is the single place where the VMEC toroidal flux enters libneo:

call nc_get(ncid, 'phi', phi)
! VMEC is left-handed (phi_v = -cylindrical angle, signgs = -1), so the wout
! stores a negative enclosed toroidal flux. SIMPLE works with a positive
! internal psi_tor = -phi_vmec/(2*pi). This is the toroidal flux on the
! poloidal covariant component, A_theta = torflux*s (see spline_vmec_data),
! and the reader recovers iota = -dA_phi/dA_theta = +iota_vmec.
phi = -phi/(2*pi)

The comment (added by 34533c9, "doc: explain the VMEC/SIMPLE toroidal-flux sign convention", comments only, no behaviour change) names signgs = -1 as the justification, but signgs is never read. grep -rn signgs over the whole tree returns exactly five hits — README.md:101, python/libneo/nemec.py:16, python/libneo/nemec.py:322 (the NEMEC→VMEC writer, which hardcodes "signgs": -1), test/python/test_nemec2vmec.py:131, and the comment above. No reader anywhere consumes it.

Downstream of that flip:

• src/vmec_support/spline_vmec_data.f90:72-74 builds sA_phi(k+1,:) = -torflux*splcoe(k,:); interpolate_vector_potential (declared line 477) sets A_theta = torflux*s, dA_theta_ds = torflux (lines 486-487).
• src/vmec_support/spline_vmec_data.f90:1028-1029: Bctrvr_vartheta = -dA_phi_ds/sqg, Bctrvr_varphi = dA_theta_ds/sqg. sqg comes from metric_tensor_symflux (line 1038), which sets sqg = sqgV*cjac with cjac = 1/(1+dl_dt) (lines 1053-1054), and sqgV comes from metric_tensor_vmec (line 1209): sqgV = R*(dR_dt*dZ_ds - dR_ds*dZ_dt). That is the signed VMEC Jacobian — the quantity whose sign signgs reports — and no absolute value is taken anywhere on the path.

So for a signgs = +1 file sqg flips sign while torflux does not, and both contravariant components come out reversed relative to the written rmnc/zmns/iotaf.

One claim in the issue is wrong. iota is not mirror-flipped. aiota = -dA_phi_ds/dA_theta_ds (spline_vmec_data.f90:675), and since dA_phi_ds = -torflux*iota_vmec while dA_theta_ds = torflux, the torflux sign cancels exactly: aiota = +iota_vmec for either sign. The damage is confined to the field components, not to iota.

The same unconditional flip is duplicated in the Boozer ingestion path:

• src/field/boozmn_reader.f90:59 — torflux_si = -b%phi(b%ns)/TWOPI
• python/libneo/booz_xform_to_boozer_chartmap.py:207 — torflux_si = -float(d["phi"][-1]) / TWOPI
• python/libneo/eqdsk_to_boozer_chartmap.py:17 — module docstring torflux = -psitor_max (negative: left-handed SIMPLE convention), matching the code at lines 211-213

and read_boozmn (src/field/boozmn_file.f90:42-70) reads no handedness marker at all.

One more latent assumption in the same family: python/libneo/vmec.py:115 and :140 gate on float(self.phi[-1]) > 0.0 and silently fall back to linear index mapping when phiedge < 0, rather than erroring.

Footgun: byte-identical orphan copies src/vmecinm_m.f90 and src/spline_vmec_data.f90 still exist at the src/ top level (verified identical with diff) and are not in any CMakeLists — only src/vmec_support/CMakeLists.txt:4-5 builds them. Edit the src/vmec_support/ copies. (34533c9 predates the move and patched the then-current top-level files.)

Reproduction / evidence

The test fixture is downloaded, not committed: test/scripts/setup_vmec_wout.py fetches simsopt's wout_LandremanPaul2021_QA_reactorScale_lowres_reference.nc to build/test/wout.nc (ctest target setup_vmec_wout, test/CMakeLists.txt:366-369).

A signgs = +1 wout is produced by reflecting the equilibrium in the midplane, Z -> -Z: rmnc unchanged, zmns -> -zmns, zaxis_cs -> -zaxis_cs, lmns -> -lmns, iotaf/iotas -> -iota, phi unchanged, signgs -> +1. Since R is unchanged and Z is negated, sqgV = R*(dR_dt*dZ_ds - dR_ds*dZ_dt) changes sign, which is exactly what signgs = +1 must mean. Note this is the mirror-image equilibrium, not a relabelling of theta: a pure theta -> -theta relabelling would additionally require xn -> -xn to stay representable in the m*theta - n*zeta phase convention, which is why the reflection is the cleaner fixture. zaxis_cs must be negated too because vmecinm_m.f90:99 reads it for the near-axis limit.

Executables land in the build root (CMakeLists.txt:83 sets CMAKE_RUNTIME_OUTPUT_DIRECTORY to ${CMAKE_CURRENT_BINARY_DIR}), so the binary is build/vmec_to_chartmap.x.

make                                    # cmake --preset default && cmake --build --preset default
cd build/test && ctest -R setup_vmec_wout
python3 - <<'EOF'
import shutil, netCDF4, numpy as np
shutil.copy("wout.nc", "wout_rh.nc")
with netCDF4.Dataset("wout_rh.nc", "a") as ds:
    for v in ("zmns", "lmns", "zaxis_cs", "iotaf", "iotas"):
        if v in ds.variables:
            ds.variables[v][:] = -ds.variables[v][:]
    if "signgs" in ds.variables:
        ds.variables["signgs"][...] = 1
    else:
        raise SystemExit("fixture has no signgs variable; adapt the recipe")
EOF
../vmec_to_chartmap.x wout.nc    lh.chartmap.nc 24 33 34
../vmec_to_chartmap.x wout_rh.nc rh.chartmap.nc 24 33 34
ncdump -h lh.chartmap.nc | grep torflux
ncdump -h rh.chartmap.nc | grep torflux

Observed: both chartmaps carry the same torflux global attribute. The attribute is written unconditionally at src/boozer_chartmap.f90:286 from vector_potentail_mod's torflux (imported at line 36), which vmecin set with the same hardcoded -phi/(2*pi) in both runs. Expected once fixed: torflux must have opposite signs between the two files, so that B^varphi = torflux/sqg — whose sqg did flip — is unchanged. As written, the rh chartmap encodes a reversed field.

The minimal generalisation that reproduces today's behaviour exactly at signgs = -1 is

torflux = signgs * phi_edge / (2*pi)

which equals -phi_edge/(2*pi) for signgs = -1 (bit-for-bit the current path) and flips for signgs = +1.

Behaviour (BDD)

Scenario: left-handed wout keeps today's behaviour bit-for-bit
  Given build/test/wout.nc with signgs = -1
  When vmecin reads it and torflux is taken from vector_potentail_mod
  Then torflux equals -phi(kparb)/(2*pi) * 1e4 * 1e4 (Gauss cm^2, fac_b0=1d4
       and fac_r0=1d2 from vmecinm_m.f90:45) to within 1e-14 relative
  And test_boozer_converter_vs_simple still matches its hard-coded SIMPLE reference

Scenario: midplane-reflected wout is recognised as positive-Jacobian
  Given build/test/wout_rh.nc built by the recipe above (signgs = +1)
  When vmec_field is evaluated at the same (s=0.5, theta, varphi) on both files
       for theta in {0, pi/3, 2*pi/3, pi} and varphi in {0, pi/8}
  Then sqrt(g) of the reflected file equals minus that of the original
       to within 1e-10 relative
  And torflux of the reflected file equals minus that of the original
  And Bctrvr_varphi agrees between the two to within 1e-10 relative
  And Bctrvr_vartheta of the reflected file equals minus that of the original
       to within 1e-10 relative

Scenario: the fixture is self-consistent before it is trusted
  Given build/test/wout_rh.nc
  When sqrt(g) is evaluated at s=0.5, theta=0, varphi=0
  Then its sign equals the signgs value stored in the file
  And the same check on build/test/wout.nc yields signgs = -1

Scenario: wout without signgs is accepted with the documented default
  Given a wout.nc from which the signgs variable has been removed
  When vmecin reads it
  Then torflux is computed with signgs = -1 and one line naming the assumed
       convention is printed

Scenario: an out-of-range signgs is rejected at the boundary
  Given a wout.nc whose signgs variable holds 0
  When vmecin reads it
  Then the run stops with an error naming signgs and the file

Implementation plan

1. src/vmec_support/canonical_coordinates_mod.f90 — add integer :: signgs = -1 to module new_vmec_stuff_mod (declared line 11; put it next to nsurfm,nstrm,nper,kpar on line 13) so the value read from the file is visible to later stages.
2. src/vmec_support/vmecinm_m.f90 — after call nc_open (line 71), read signgs with the optional-variable idiom already used at src/field/boozmn_file.f90:48-50 (if (nf90_inq_varid(ncid, 'signgs', varid) == nf90_noerr)), defaulting to -1 and printing the assumed convention when absent. error stop for any value other than +/-1. Replace line 82 with phi = signgs*phi/(2*pi) and rewrite the comment block at lines 77-81 to state the general rule instead of the signgs = -1 special case.
3. src/vmec_support/spline_vmec_data.f90 — update the comment at lines 67-69 to record sign(torflux) = signgs * sign(phi_edge). No code change: sA_phi (line 73) and aiota (line 675) are already sign-transparent, per the analysis above.
4. src/boozer_chartmap.f90 — next to the torflux attribute (line 286) write an integer signgs global attribute; mirror the read in src/field/boozer_chartmap_io.f90 next to the torflux attribute read (line 72), and add the component to boozer_chartmap_data_t in src/boozer_chartmap_types.f90 (alongside torflux, line 18), defaulting to -1 for files that lack it.
5. src/field/boozmn_reader.f90 — add integer, intent(in), optional :: signgs to the converter; replace line 59 with a signgs-driven expression, default -1; propagate into the chartmap data component from step 4.
6. python/libneo/booz_xform_to_boozer_chartmap.py — same change at line 207 plus a --signgs CLI flag. When unset, use the geometric recovery _derive_psi_prime (defined line 101) to fix the sign: it already returns a signed float(np.median(psi_primes)) from mean(J*B**2)/(G + iota*I) with a signed Jacobian J, so its sign is meaningful. Cross-check its magnitude against abs(phi[-1])/(2*pi) to 1e-6 relative and raise on mismatch. This overlaps Geometric toroidal-flux recovery (derive_psi_prime) in the booz_xform chartmap converter is untested; CI only exercises the phi_b-present branch #307 (the _derive_psi_prime branch is untested, and CI only exercises the phi_b-present path); do both in one PR.
7. python/libneo/vmec.py:115,140 — replace the float(self.phi[-1]) > 0.0 guard with abs(float(self.phi[-1])) > 0.0 and normalise by self.phi[-1], so a negative phiedge no longer silently degrades to linear index mapping.
8. New test test/source/test_vmec_signgs_invariance.f90 implementing scenarios 2 and 3; register in test/CMakeLists.txt next to test_chartmap_matches_vmec (lines 546-550), with target_link_libraries(... neo) and DEPENDS setup_vmec_wout.
9. New script test/scripts/setup_vmec_wout_mirrored.py producing wout_rh.nc from build/test/wout.nc exactly as in the MRE; register as a ctest fixture with DEPENDS setup_vmec_wout, following the setup_vmec_chartmap_python pattern at test/CMakeLists.txt:372-380.
10. New Python test python/tests/test_boozmn_signgs.py extending the synthetic-boozmn construction used by test/source/test_boozmn_reader.f90 (which fills b%phi = phi_full at line 144) to cover both signs.
11. doc/ExtraDocuments/Conventions.md — add a subsection under "Left/right handed systems - The COCOS convention" (heading at line 104) stating torflux = signgs * phi_edge/(2*pi) and that chartmap files now record signgs.
12. Downstream: SIMPLE reads the chartmap through libneo, so no SIMPLE source change is required, but SIMPLE:tools/booz_xform_to_boozer_chartmap.py lines 211-218 (cited in the issue) carries the same hardcoded negation and should be deleted in favour of the libneo copy, or fixed in lockstep. Land the libneo PR first: .github/workflows/downstream-gate.yml dispatches the downstreams listed in ci/downstreams — SIMPLE, NEO-2, MEPHIT, KAMEL and rabe are blocking, NEO-RT is report-only.

Acceptance criteria

• make && cd build && ctest --output-on-failure passes with zero regressions.
• ctest -R test_boozer_converter_vs_simple still matches its hard-coded SIMPLE reference values (test/CMakeLists.txt:558-562). This is the guard that the signgs = -1 path is untouched.
• ctest -R test_vmec_signgs_invariance passes: with wout.nc and mirrored wout_rh.nc, sqrt(g) and Bctrvr_vartheta are antisymmetric and Bctrvr_varphi is invariant, each to 1e-10 relative, at every sample point of scenario 2.
• ncdump -h on a chartmap built from wout_rh.nc shows :signgs = 1 and a torflux of opposite sign to the one from wout.nc; ncdump -h on the wout.nc chartmap shows :signgs = -1 and a torflux numerically unchanged from main.
• grep -rn -- "-phi/(2\*pi)\|-b%phi(\|-float(d\[\"phi\"\])" src/ python/libneo/ returns no unconditional negation outside a signgs-driven expression.
• A wout.nc with signgs deleted still runs and prints the assumed-convention note; one with signgs = 0 stops with an error naming signgs.
• ctest -R test_chartmap_matches_vmec and ctest -R validate_vmec_rho_lcfs_auto unchanged.

References

• S. P. Hirshman and J. C. Whitson, Steepest-descent moment method for three-dimensional magnetohydrodynamic equilibria, Phys. Fluids 26, 3553 (1983) — VMEC and the (s, theta, zeta) chart whose Jacobian sign is signgs.
• A. H. Boozer, Plasma equilibrium with rational magnetic surfaces, Phys. Fluids 24, 1999 (1981) — the identity sqrt(g) B^2 = psi'(G + iota*I) that _derive_psi_prime cites in its own docstring (booz_xform_to_boozer_chartmap.py:104-106).
• O. Sauter and S. Yu. Medvedev, Tokamak coordinate conventions: COCOS, Comput. Phys. Commun. 184, 293 (2013) — the canonical treatment of this class of silent sign error. libneo already implements COCOS classification for g-files in src/magfie/geqdsk_tools.f90 (cocos_t, line 16; classification from line 268), and this issue is the VMEC-side analogue.
• Prior art: STELLOPT's LIBSTELL/Sources/Modules/read_wout_mod.f90 reads the wout signgs variable (NetCDF name constant vn_sgs = 'signgs', long name 'Sign jacobian') into the module variable isigng, i.e. it keeps the value rather than assuming it. That file is already linked from this repo's README.md:103 (URL at README.md:107).

Not asserted here: how any particular downstream code (ASCOT5, DESC, SIMSOPT) handles this internally. python/libneo/nemec.py:16 records that DESC's VMEC reader/writer also uses signgs = -1, which is a statement about the NEMEC writer's assumptions, not evidence that DESC generalises the sign. test/scripts/test_ascot5_compare.py cross-checks the Cerfon-Freidberg analytic solver against ASCOT5's B_GS and touches neither VMEC nor Boozer coordinates, so it says nothing about this problem.

Scope and risk

Single PR, roughly 80-150 lines of source plus two tests. The risk is concentrated in step 2, where a sign slip regresses every VMEC consumer at once; test_boozer_converter_vs_simple is the correct tripwire and must be kept unmodified. Steps 6-7 can be split into a Python-only follow-up if the Fortran change needs to land ahead of the downstream gate.

Related open issues (titles verified): #307 — geometric toroidal-flux recovery (derive_psi_prime) is untested, CI only exercises the phi_b-present branch; fix together with step 6. #171 — chartmap zeta sign support / lifting the vmec = cyl assumption; the same "record the convention in the file" problem on the toroidal angle, so agree the attribute naming with it. #302 — g-file convention checks only on the geqdsk_tools path; the tokamak analogue of this defect. itpplasma/SIMPLE#430 — .bc retirement, verify chartmap and booz_xform conventions. itpplasma/GLISS#10 — read precomputed BOOZ_XFORM files and support asymmetric VMEC equilibria; the lasym side of convention-complete ingestion.

[krystophny]

Two clarifications from an independent native audit:

1. The old vmecinm_m.f90 comment quoted in this issue is demonstrably wrong when it says phi_v=-cylindrical angle. Standard VMEC zeta is the physical geometric cylindrical azimuth. The usual negative Jacobian belongs to the complete (s,theta,zeta) chart. A comment-only correction is pushed at 46c57f5 on agent/eqdsk-booz-xform; it makes no runtime sign change.
2. A fresh standard signgs=-1 VMEC -> booz_xform -> Cartesian gate closes the ordinary path: positions agree within 2.8e-9 m and Cartesian B has relative L2 error 2.6e-5. This validates current behavior for that fixture, but it does not independently validate the proposed signgs=+1 formula.

I therefore agree that the unconditional omission of signgs is a real unsupported-input hazard, while keeping the positive-Jacobian behavior evidence-gated. Before changing production signs, the acceptance suite should run a self-consistent signgs=+1 native fixture, verify the sign of the actual coordinate determinant against the stored metadata, and compare Cartesian B against the corresponding physical reference. A mere edit of the signgs scalar or a reflected geometry without verifying all coupled phi/iota/lambda/B quantities is not sufficient.

[krystophny]

The positive-signgs defect is now independently reproduced with a pure coordinate relabel, avoiding the physical-mirror ambiguity.

Manufactured map

s' = s
theta' = -theta
zeta' = zeta

This leaves Cartesian space and the physical magnetic field unchanged. The exact required transformations are

J' = -J
iota' = -iota
A_theta' = -A_theta
A_zeta' = A_zeta
B_theta' = -B_theta
B_zeta' = B_zeta
B_cartesian' = B_cartesian.

The relabeled wout transforms xn/xn_nyq, Fourier sine/cosine parity, angular lambda, iota/chi, the signed Jacobian harmonics, and u/v field components as tensors; it sets signgs'=+1. It is not a reflected physical equilibrium.

Native execution

Both files were run through SIMPLE's compiled export_boozer_chartmap_tool.x, which uses libneo's native VMEC reader and Boozer transform (SIMPLE 10ac7e7, libneo 46c57f5). Comparing the relabeled chart at theta'=-theta gives:

maximum Cartesian geometry error       9.8429e-9 cm
relative L2 |B| error                   9.3583e-12
torflux original                       -1.334629991641011e6 G cm^2
torflux relabeled reader               -1.334629991641011e6 G cm^2
required torflux relabeled             +1.334629991641011e6 G cm^2
relative Cartesian-B error vs same B    2.0000000000000195
relative Cartesian-B error vs -B        2.6196e-13

The covariants diagnose the same failure:

reader: A_zeta'=-A_zeta, B_theta'=B_theta, B_zeta'=-B_zeta
needed: A_zeta'= A_zeta, B_theta'=-B_theta, B_zeta'= B_zeta

So the reader reverses the complete physical B vector while preserving geometry and magnitude. This is not metadata-only and cannot be hidden by a |B| or orbit-shape test.

The reproduction and machine-readable evidence are in the ITER TC24 audit as:

manufactured_equivalence/manufacture_vmec_theta_relabel.py
manufactured_equivalence/verify_vmec_signgs_native.py
manufactured_equivalence/evidence/vmec_signgs_relabel_defect.json

This closes the evidence caveat in my previous comment. The issue's proposed general rule

torflux = signgs * phi_edge / (2*pi)

has the correct tensor behavior for this relabel: it preserves the current signgs=-1 result and flips A_theta for signgs=+1, restoring the same Cartesian B. The acceptance test should use this coordinate-relabel fixture (or an equivalent one) because it provides both a required invariant and a strong wrong-sign negative control.