Restore ASCOT5 adapters on top of canonical split

[krystophny]

Summary

• reintroduce the Python ASCOT5 helpers and pytest coverage on top of the canonical-only base branch
• wire the public package export so downstream users see
• revive the detailed ASCOT5 TODO for follow-up work

Testing

• source activate.sh && cmake --build build -j
• source activate.sh && pip install -e .
• source activate.sh && cd build && ctest --output-on-failure
• source activate.sh && pytest -q

[slopqueue]

ORCHESTRATE_REVIEW_STATUS: request_changes

All verification complete. Here's the review verdict.

---

PR #130 Review: ASCOT5 field conversion helpers

Verdict: Request changes — two critical bugs in write_b3ds_hdf5

Critical: psi not transposed in HDF5 output

File: python/libneo/ascot5/__init__.py:560

a5py's B_3DS.write_hdf5 (tag 5.6.1) transposes psi from (nr, nz) → (nz, nr) before writing:

# a5py reference
psi = np.transpose(psi)  # (nr, nz) -> (nz, nr)
g.create_dataset("psi", (psi_nz, psi_nr), data=psi, dtype="f8")

libneo writes psi directly without transposing:

# libneo (line 560)
g.create_dataset("psi", data=data["psi"], dtype="f8")  # writes (nr, nz)

When nr != nz, ASCOT5 reads the array as (nz, nr) and transposes it to (nr, nz) — getting a transposed psi grid. The test test_field_from_vmec_generates_b3ds uses nr=nz=24 (square), so the bug is masked. The mgrid test uses nr=6, nz=5 but never reads the HDF5 back.

Fix: g.create_dataset("psi", data=np.transpose(data["psi"]), dtype="f8")

Critical: psi grid metadata not written

File: python/libneo/ascot5/__init__.py:538-567

a5py writes six psi-grid datasets that ASCOT5's C reader expects:

# a5py reference — all six written
g.create_dataset("psi_rmin", ...)  # psi_rmax, psi_nr, psi_zmin, psi_zmax, psi_nz

libneo writes none of them. ASCOT5 will either fail or silently use defaults.

Fix: Add the six datasets, defaulting to the B-grid values (same as a5py does when psi_rmin=None):

g.create_dataset("psi_rmin", data=(data["b_rmin"],), dtype="f8")
g.create_dataset("psi_rmax", data=(data["b_rmax"],), dtype="f8")
g.create_dataset("psi_nr", data=(data["b_nr"],), dtype="i4")
g.create_dataset("psi_zmin", data=(data["b_zmin"],), dtype="f8")
g.create_dataset("psi_zmax", data=(data["b_zmax"],), dtype="f8")
g.create_dataset("psi_nz", data=(data["b_nz"],), dtype="i4")

Minor: unused Iterable import

File: python/libneo/ascot5/__init__.py:7

Iterable imported but never used. Dict and Tuple are used; Iterable is not.

Minor: zaxis_cs likely wrong in _load_vmec_metadata

File: python/libneo/ascot5/__init__.py:75-78

if "zaxis_cs" in ds.variables:
    axis_z = float(ds.variables["zaxis_cs"][0])

In VMEC, zaxis_cs[0] is the n=0 sine coefficient (always 0). The axis Z coordinate is zaxis_cc[0] (the constant term). This branch should be removed — always use zaxis_cc[0].

Test gap: no round-trip verification

Neither test reads the HDF5 file back with a5py's B_3DS.read() to verify the layout. A round-trip test with nr != nz would catch the psi transpose bug. The existing tests only check that the file exists and the group name starts with B_3DS.

Verified correct

• br/bphi/bz HDF5 layout: np.transpose(field.br, (2, 1, 0)) matches a5py exactly: (nr, nphi, nz) → (nz, nphi, nr)
• Gauss→Tesla: applied at export only; internal fields stay in Gauss
• b_phimax formula: phi_grid_deg[-1] + delta_phi = period boundary. Operator precedence is correct (ternary binds looser than +)
• field_from_vmec psi from iphi==0: valid — psi is a flux surface quantity, independent of phi
• field_from_mgrid transpose: (nphi, nz, nr) → (nr, nphi, nz) matches B3DSField internal layout
• _fill_nan_by_local_mean: correct iterative fill with increasing window
• Newton/Picard/coarse-search solver: correct clamping, wrapping, and fallback chain
• _initial_guess: reasonable radial projection using outer surface edge

Recommendation

Block merge until the two critical bugs are fixed. The psi transpose and missing grid metadata will produce incorrect ASCOT5 input files whenever nr != nz. Add a round-trip test with non-square grids to prevent regression.

[slopqueue]

Review verdict: Request changes

Summary: The PR adds Python ASCOT5 field export helpers (field_from_vmec, field_from_mgrid, write_b3ds_hdf5) with a frozen B3DSField dataclass, a multi-strategy Newton/Picard/coarse-search flux-coordinate solver, doc, and pytest coverage. The VMEC path is well-tested with value-level assertions; the Fortran wrapper interface at src/f2py_interfaces/f2py_vmec_wrappers.f90:8-19 correctly matches the Python unpacking (17 intent-out args, 3 inputs). The optional import in python/libneo/__init__.py:33-36 correctly tolerates missing h5py/_magfie. However, the required gate CI is red (downstream NEO-2 and rabe builds fail), and the mgrid path has a unit inconsistency that would produce a 10000x too-small field.

Findings:

1. [blocker] CI: gate — The gate required check failed. NEO-2 CI: Fatal Error: Cannot open module file 'hdf5_tools.mod' for reading at (1): No such file or directory in h5merge_multispec.f90:14. rabe CI: Fatal Error: Cannot open module file 'boozmn_reader.mod' in boozer_field.f90:7. The hdf5_tools module exists in libneo (src/hdf5_tools/hdf5_tools.F90:1) but downstream builds can't find its .mod. These failures stem from the canonical-split refactoring on the feature branch (c4c01b8 Decouple VMEC wrappers from canonical helpers), not from this PR's Python-only diff, but the gate tests this PR's ref and must be green to merge.
2. [major] python/libneo/ascot5/__init__.py:567-569 — mgrid B-field unit bug. write_b3ds_hdf5 multiplies all B components by GAUSS_TO_TESLA (1e-4). For the VMEC path (Gauss from Fortran) this is correct. But field_from_mgrid (lines 474-476) reads mgrid data in Tesla (SI): src/field/mgrid.f90:4 confirms mgrid stores "physical cylindrical components" and rmin/rmax are in metres. Multiplying Tesla by 1e-4 yields a field 10000x too small. The doc (doc/ascot5.md:18-19) says B3DSField stores Gauss, but field_from_mgrid never converts. Fix: convert mgrid data to Gauss (multiply by 1e4) before storing in B3DSField, or skip the conversion in write_b3ds_hdf5 for mgrid-origin fields.
3. [minor] python/libneo/ascot5/__init__.py:148 — _solve_flux_coordinates early exit math.hypot(r_target, z_target) < 1.0e-6 checks against (R=0, Z=0) instead of the magnetic axis (axis_r, axis_z), which the function doesn't receive. For typical stellarator grids (R ~ 0.5–2 m) this never triggers. The caller's axis_dist <= 0.03 fallback at line 407 mitigates the practical impact.
4. [minor] python/libneo/ascot5/__init__.py:58 — b_phimax in as_dict uses phi[1] - phi[0] for the last cell width, which is only correct for uniformly spaced phi grids. For mgrid files with a phi variable (test_ascot5_mgrid.py:48 creates one), non-uniform spacing would give an incorrect b_phimax.
5. [minor] test/python/test_ascot5_mgrid.py:57-89 — No value-level test for mgrid B-field. The test checks shapes and metadata only, never verifying output field values or units. A value assertion (e.g., comparing g["br"] against the known synthetic input) would catch the unit bug in finding 2.

Verdict: request_changes — the gate CI must be green and the mgrid unit inconsistency must be resolved before merge; the VMEC path and test structure are otherwise sound.