Read Boozer .bc directly via a bc_field_t (no NEO-2/SuiteSparse dependency)

[krystophny]

Motivation

rabe-vs-NEO-2 benchmarks currently compare two different Boozer transforms of the same equilibrium. NEO-2 reads .bc files generated through booz_xform; rabe converts the VMEC wout internally with the libneo converter (adopted in #44). The measured difference between the two transforms sits at the QI offset signal level and orders of magnitude above the QA signal (see the converter-fidelity gate results in #92), so the converter confounds every benchmark. Reading the same .bc file NEO-2 uses removes that confound.

The earlier neo_field_t served this purpose but pulled in NEO-2 and SuiteSparse, so #74 moved it to the internal test submodule for the public release.

Proposed design

A bc_field_t extending field_t that reads a .bc file directly, with no NEO-2 or SuiteSparse code:

• |B| by direct Fourier sum from bmnc (and bmns for non-symmetric files).
• Jacobian from the current columns: sqrt(g) proportional to (G + iota*I)/B^2 in Boozer coordinates.
• Covariant components B_phi, B_theta from the Jpol/nper and Itor columns.
• Radial spline of the Fourier coefficients across surfaces.
• nabla_s mode switch: default uses the correct stream-function sign (phi_cyl = phi_B - v); an opt-in neo2_bc mode reproduces the legacy opposite-sign convention of existing .bc data, for benchmarking against historical NEO-2 results (see the 36 use simples vmec to boozer converter #44 commit message for the mechanism).

The evaluation core could live in libneo (itpplasma/libneo#305) with bc_field_t as a thin adapter here.

References

• 36 use simples vmec to boozer converter #44 (converter adoption and the v-sign convention)
• 58 prepare codebase for public release remove mandatory suitesparse dependency #74 (removal of neo_field_t from the public tree)
• QA benchmark discrepancy: spline floor of field-line |B| maxima; evaluate maxima spectrally #92 (converter-fidelity gate results)
• Standalone Boozer field evaluator from .bc files (lift neo_magfie design, no SuiteSparse) libneo#305 (standalone evaluator)

[krystophny]

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

Still relevant, and closer to done than the issue text suggests: a working bc_field_t already exists on the unmerged branch origin/feat/bc-field, and libneo's .bc reader is already inside the pin main uses. What is missing is the rebase onto the field_type dispatch introduced by 4f70aa4, an API rename the branch predates, and a real fixture-backed test. This spec is written as "land the branch", not "start from scratch".

Purpose

rabe's neoclassical coefficients are evaluated entirely from Boozer-space geometry (src/fieldline/, src/offset/), so a benchmark against NEO-2 currently differences two independent VMEC→Boozer transforms rather than two physics models. Reading the same .bc file NEO-2 reads removes the transform from the comparison entirely.

What rabe#92 actually concluded (quoting its measurement table, not paraphrasing it):

• Spline-interpolation floor on field-line |B| maxima, spline vs direct Fourier sum, s = 0.5, QA (LandremanPaul2021 reactor scale): cubic median 7.6e-8, but quintic median 3.6e-9 — and rabe runs quintic (ns_tp = 5).
• QA signal level is ~1e-8. So at rabe's actual spline order the pure interpolation floor sits below the signal, and QA benchmark discrepancy: spline floor of field-line |B| maxima; evaluate maxima spectrally #92 says so explicitly: "at quintic order (which rabe uses, ns_tp = 5) the pure interpolation floor drops to 3.6e-9, marginally BELOW the QA signal."
• QA benchmark discrepancy: spline floor of field-line |B| maxima; evaluate maxima spectrally #92's closing conclusion is that the converter difference, not the spline order, sets the QA floor: "the libneo-vs-booz_xform converter difference is at or above the QI signal and two to three orders above the QA signal ... the QA benchmark floor is set by the VMEC-to-Boozer converter, not by the spline interpolation this issue originally targeted."

That is the argument for #98, and it is stronger than a spline-floor argument: #92 concludes the converter is the confound, and #98 removes the converter. Do not repeat the claim that the 7.6e-8 cubic number is rabe's operating floor, and do not describe 7.6e-8 vs 1e-8 as "two orders" — #92 calls it "a factor ~8".

Secondary, and much weaker: the branch's bc_field_t declares rel_accuracy_B = 1e-12_dp where boozer_field_t declares 1e-9_dp (src/vmec/boozer_field.f90:374). These are hard-coded constants consumed by the tolerance logic, not measured accuracies. Treat the change as "the .bc path asserts a tighter tolerance", not as evidence of measured improvement.

Current state

main accepts exactly three field formats. src/read_file.f90:142-146:

if (len(trim(field_type)) == 0 .or. &
    (trim(field_type) /= 'vmec_nc' .and. &
     trim(field_type) /= 'booz_xform' .and. &
     trim(field_type) /= 'chartmap')) then
    print *, "field_type must be 'vmec_nc', 'booz_xform', or 'chartmap'"

dispatched in app/main.f90:77-86 onto init_from_vmec / init_from_boozmn / init_from_chartmap (src/vmec/boozer_field.f90:35-37). There is no .bc path, despite the title of 4f70aa4 ("Support booz_xform and chartmap field input in addition to .bc"): that commit added the other two formats and renamed boozer_field_init into the three init_from_* bindings; the ".bc" in its title refers to nothing that exists in the tree.

Four facts that shape the work:

1. The libneo reader is already in the pin. libneo/src/field/bc_file.f90 (bc_data_t, read_bc_file, 212 lines, added in libneo commit 6dbe6a5, which landed as libneo PR #365) parses the Strumberger/NEO-2 .bc format and is compiled into the LIBNEO::boozer target (libneo/CMakeLists.txt:338-358), which sits above the LIBNEO_BUILD_NUMERICS early-return at libneo/CMakeLists.txt:363-364 — so no BLAS/LAPACK, no fortnum, no SuiteSparse. 6dbe6a5 is an ancestor of rabe's current pin 54e372f2 (cmake/libneo.cmake:29), verified with git merge-base --is-ancestor. No libneo bump is needed.

   (Whether this fully or partly closes Standalone Boozer field evaluator from .bc files (lift neo_magfie design, no SuiteSparse) libneo#305 was not verified here — the commit references PR #365, and the relationship to #305 is an inference. Check #305 before claiming it is done.)

2. The existing boozmn detour is not a substitute. libneo/python/libneo/bc_to_booz_xform.py can convert .bc → boozmn, which rabe reads. But libneo/src/field/boozmn_reader.f90 uses only bmnc, buco, bvco, iota, phi — rmnc, zmns, pmns are read and dropped — and then calls build_boozer_from_chartmap, which sets use_B_r = .false. (libneo/src/boozer_converter.F90:1367, inside the subroutine beginning at line 1319). splint_boozer_coord assigns its sqrt_g_ss_B argument only under that flag (libneo/src/boozer_converter.F90:296-297 and 346-347), while init_from_boozmn (src/vmec/boozer_field.f90:85-93) never sets use_B_r = .true. the way init_from_vmec does (src/vmec/boozer_field.f90:71). So compute_nabla_s (src/vmec/boozer_field.f90:349-369) leaves its local sqrt_g_ss unwritten on the boozmn and chartmap paths. Also, interp_modes in boozmn_reader.f90:145-189 is piecewise linear in rho ((1-frac)*bmnc_h(mn,k) + frac*bmnc_h(mn,k+1)), which reintroduces exactly the radial error Read Boozer .bc directly via a bc_field_t (no NEO-2/SuiteSparse dependency) #98 wants splined away.

3. A module-name collision blocks the obvious naming. test/helpers/bc_file.f90 declares module bc_file — the same name as libneo's. The branch already handles this by renaming it to bc_writer. Its only consumer in rabe is test/unit/test_read_namelist.f90.

4. The branch predates a rename and will not compile as-is. test/integration/vmec/test_bc_roundtrip.f90 calls field_nc%boozer_field_init(nc_filename, radial_spline_order=5, angular_spline_order=5, grid_refinement=3). boozer_field_init does not exist anywhere on main (grep -rn boozer_field_init src/ test/ app/ returns nothing); 4f70aa4 split it into init_from_vmec. The call site must be rewritten during the rebase.

The branch origin/feat/bc-field (9 commits, merge-base 8969e9d, +808/−23) contains src/field/bc_field.f90 (370 lines) implementing every bullet of the proposed design: batch radial splines (construct_batch_splines_1d) of bmnc, rmnc, zmns, vmns and of (iota, B_theta, B_phi); sqrtg = -psi_tor_edge*(iota*B_theta + B_phi)/B**2 (line 223); currents mapped by -mu0/(2 pi); and the sign switch at lines 260-263:

v_factor = -two_pi/self%nfp
if (self%neo2_nabla_s) v_factor = -v_factor
p_t = v_factor*v_t
p_p = 1.0_dp + v_factor*v_p

It predates field_type, so it dispatches on the filename extension in a hand-written init_field in app/main.f90 (branch lines 216-247) and converts field to class(field_t), allocatable with two select type shims (fix_field_to_surface, get_field_iota). The extension sniffing and the shims must be discarded; the class(field_t), allocatable conversion must be kept.

Actual rebase conflict set. Measured with git merge-tree --write-tree --name-only main origin/feat/bc-field:

file	result
CLAUDE.md	CONFLICT
app/main.f90	CONFLICT
cmake/libneo.cmake	CONFLICT
test/integration/vmec/CMakeLists.txt	CONFLICT
src/read_file.f90	auto-merges clean
test/helpers/bc_writer.f90	auto-merges clean

Reproduction / evidence

Confirm gap (1), the missing .bc input:

cd /path/to/rabe && make CONFIG=Release
mkdir -p /tmp/bc && cd /tmp/bc
cp /path/to/libneo/python/tests/circ.bc .
printf "&rabe_config\n field_file='circ.bc',\n field_type='bc',\n s_tor=0.5,\n M_pol=-1.0d0, N_tor=4.0d0, sign_sqrtg=-1.0d0,\n max_n_fieldlines=20, n_eta=100\n/\n" > rabe.in
/path/to/rabe/build/rabe.x

Observed: field_type must be 'vmec_nc', 'booz_xform', or 'chartmap' then error stop. Expected after this issue: rabe.nc with finite Lambda_A, Lambda_B, Lambda_S.

Confirm gap (2), the unwritten nabla_s on the existing detour:

python -m libneo.bc_to_booz_xform /path/to/libneo/python/tests/circ.bc boozmn_circ.nc
# rabe.in with field_type='booz_xform', field_file='boozmn_circ.nc'
/path/to/rabe/build-debug/rabe.x     # Debug build: -finit-real=nan, cmake/compiler.cmake:15

Expected observation, an abort, not NaN output: nothing writes sqrt_g_ss, so with -finit-real=nan the surface-averaged nabla_s is NaN, dr_dAtheta is NaN, and execution stops at src/offset/coefficients.f90:116-119 with error: dr_dAtheta is NaN.. Note that the zero guard just above at src/offset/coefficients.f90:107 does not fire (abs(NaN) < machine_eps is false) — the NaN guard at 116 is what catches it. In Release without -finit-real=nan the uninitialised value is whatever is on the stack, so the run may instead produce silent, non-reproducible garbage that passes both guards. That silent-garbage case is the actual hazard.

Verification recipe for the new path: the branch's test/unit/test_bc_field_nabla_s.f90 builds a circular-cross-section torus with R0 = 10.0, a = 1.0, B0 = 5.0, s0 = 0.49, r(s) = a*sqrt(s), B = B0*sqrt(R0/R), and flux = pi*a^2*B0, for which |nabla s| = 2*sqrt(s)/a holds exactly; it asserts B_mod against B0/sqrt(1 + eps*cos(theta)) to reltol 1e-12 and nabla_s to 1e-10 over an 8x4 (theta, phi) grid. That is the acceptance test to keep.

Behaviour (BDD)

Scenario: .bc input is accepted through the field_type dispatch
  Given rabe.in with field_type = 'bc' and field_file = 'circ.bc'
  When build/rabe.x runs in a directory containing libneo's python/tests/circ.bc
  Then rabe.nc is written with finite Lambda_A, Lambda_B, Lambda_S at s_tor = 0.5
   And no "field_type must be" error is printed

Scenario: analytic |nabla s| on a circular torus
  Given the synthetic .bc written by test/unit/test_bc_field_nabla_s.f90
        with a = 1.0, B0 = 5.0, R0 = 10.0, vmns = 0.0
  When bc_field_t%fix_to_surface(0.49) then compute_nabla_s(theta, phi)
       over theta = i*2*pi/8 for i in 0..7 and phi = j*2*pi/4 for j in 0..3
  Then every value matches 2*sqrt(0.49)/1.0 to relative tolerance 1e-10
   And compute_B_mod matches 5.0/sqrt(1 + (sqrt(0.49)/10.0)*cos(theta))
       to relative tolerance 1e-12

Scenario: the v-sign switch is a no-op when vmns vanishes
  Given the vmns = 0 file above
  When one bc_field_t is initialised with neo2_nabla_s = .false.
       and another with neo2_nabla_s = .true.
  Then their compute_nabla_s values agree to relative tolerance 1e-14

Scenario: the v-sign switch is equivalent to flipping vmns
  Given a second .bc file identical except vmns = -0.05
  When it is read with neo2_nabla_s = .false.
   And the vmns = +0.05 file is read with neo2_nabla_s = .true.
  Then the two compute_nabla_s values agree to relative tolerance 1e-14

Scenario: VMEC round trip preserves the field
  Given input/wout_LandremanPaul2021_QH_reactorScale_lowres_reference.nc
        loaded into boozer_field_t via init_from_vmec with grid_refinement = 3
  When it is written to roundtrip.bc over 9 surfaces linspace(0.1, 0.9)
       with m_max = n_max = 16 and reloaded as bc_field_t
  Then at s = 0.5 over a 17x13 angle grid B_mod and sqrtg agree to relative 1e-5
   And iota, R and the covariant components agree to relative 1e-7

Scenario: the boozmn path no longer leaves nabla_s unwritten
  Given rabe.in with field_type = 'booz_xform'
  When build/rabe.x runs
  Then it either produces a metric-derived nabla_s or aborts with an explicit
       "boozmn path does not provide |grad s|" message, never an unwritten value

The third and fourth scenarios are what test_bc_field_nabla_s.f90 actually asserts (branch lines 75-96). The branch does not assert that the two neo2_nabla_s settings differ by any particular amount, and it cannot: the analytic case has vmns = 0, where the switch is provably a no-op. Do not write an acceptance criterion of the form "the default satisfies the analytic check while .true. does not".

Implementation plan

1. Rebase origin/feat/bc-field onto main. Merge-base is 8969e9d; 4f70aa4 and 5f035b7 land afterwards. Resolve the four conflicts listed in the table above. src/read_file.f90 merges cleanly.

2. src/field/bc_field.f90 — take from the branch unchanged except: delete the unused grow subroutine at branch lines 354-368 (it has no call site in the file; the equivalent grow now lives in libneo src/field/bc_file.f90:196-210 alongside first_surface_pass), and re-run fprettify with the repo's .fprettify config for the 88-column rule.

3. src/read_file.f90 — two independent edits:

   • From the branch diff: add logical, public, protected :: neo2_nabla_s to the module variables, add neo2_nabla_s to the namelist /rabe_config/ list, and default it to .false. immediately after unsafe_mode = .false..
   • Written fresh, not from the branch: add 'bc' to the field_type whitelist at lines 142-146 and to the error message. The branch has no whitelist change at all — it predates field_type and dispatches on the file extension instead. Discard the branch's extension sniffing.

4. app/main.f90 — do not import the branch's init_field / fix_field_to_surface / get_field_iota shims. Instead extend the existing select case (trim(field_type)) at line 77 with case ('bc'). field must change from type(boozer_field_t) (line 34) to class(field_t), allocatable because bc_field_t and boozer_field_t are unrelated concrete types; keep psi_tor_edge, R, nfp as local scalars captured at load time, as the branch does.

   fix_to_surface and get_iota are bindings on boozer_field_t but not on field_t (src/field/field_base.f90:12-17 declares exactly six deferred procedures: compute_B_sqrtg_dB_dx, compute_B_and_dB_dx, compute_B_mod, compute_nabla_s, rel_accuracy_B, get_covariant_components). Two options:

   • select type in app/main.f90 — smallest diff, touches one file.
   • Add the two as deferred bindings on field_t — cleaner, but note this requires implementing them on all six current implementors, not two: src/field/fourier_field.f90:9, src/vmec/boozer_field.f90:27, test/helpers/mock_field.f90:7, test/helpers/mock_perturbed_field.f90:8, test/helpers/anti_sigma_field.f90:7, test/helpers/horner_fourier_field.f90:8. Confirm with grep -rn "extends(field_t)" src/ test/ before committing to this route.

   Pick select type unless the reviewer asks otherwise; the six-implementor cost was previously understated.

5. src/field/CMakeLists.txt — add bc_field.f90 to field_lib and add LIBNEO::boozer to target_link_libraries (branch diff, verbatim). Discard the branch's cmake/libneo.cmake change entirely: the branch sets the pin to 6dbe6a5, which is older than main's 54e372f2, so taking it would be a downgrade. 54e372f2 already contains bc_file.f90.

6. test/helpers/bc_file.f90 → test/helpers/bc_writer.f90 (module bc_writer), update test/helpers/CMakeLists.txt and the sole consumer test/unit/test_read_namelist.f90. Required: the current name shadows libneo's bc_file module, which bc_field.f90 imports.

7. test/unit/test_bc_field_nabla_s.f90 — take from the branch as-is; it only depends on field_lib and utils_lib.
   test/integration/vmec/test_bc_roundtrip.f90 — take from the branch and change call field_nc%boozer_field_init(...) to call field_nc%init_from_vmec(...) (same argument list: radial_spline_order=5, angular_spline_order=5, grid_refinement=3). Without this the test does not compile against main.
   Register both with LABELS "quick" per the branch diff to test/unit/CMakeLists.txt and test/integration/vmec/CMakeLists.txt (test names TestBcFieldNablaS, TestBcRoundtrip).

8. New: test/integration/vmec/input/circ.bc — copy libneo/python/tests/circ.bc (109 KB) into that directory. No CMake change is needed for the copy: test/integration/vmec/CMakeLists.txt:1 already does file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/input DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) for the whole directory. Add test/integration/vmec/test_bc_circ.f90 asserting iota, psi_tor_edge, and B_mod at three fixed (s, theta, phi) points against values dumped from libneo.boozer.BoozerFile, and register it as TestBcCirc with LABELS "quick". The branch only has synthetic fixtures; a real file exercises the header/mode-table parser.

9. src/read_file.f90 guard — error stop if neo2_nabla_s = .true. with field_type /= 'bc', so the flag cannot be silently ignored.

10. CLAUDE.md, README.md, test/golden/input/README.md — document field_type = 'bc' and neo2_nabla_s. CLAUDE.md:114-116 mandates this: "update ... This file (CLAUDE.md) — if adding/removing/renaming modules, libraries, derived types, make targets, test labels, or changing I/O formats".

11. Separate PR (do not fold in): the boozmn/chartmap nabla_s hole. Either make init_from_boozmn set use_B_r = .true. and have libneo's boozmn_reader build sqrt_g_ss from the rmnc/zmns/pmns it currently discards (libneo/src/field/boozmn_reader.f90, libneo/src/boozer_converter.F90:1367), or make rabe's compute_nabla_s abort when use_B_r is false. This is a libneo-first change; file it against itpplasma/libneo.

Acceptance criteria

• make test_all (which runs ctest --test-dir build/test -L "quick|slow") passes with TestBcFieldNablaS, TestBcRoundtrip, and TestBcCirc present.
• ctest -R TestBcFieldNablaS --test-dir build/test passes: nabla_s within 1e-10 relative of 2*sqrt(0.49)/1.0, B_mod within 1e-12, over the full 8x4 grid.
• ctest -R TestBcRoundtrip --test-dir build/test passes at tol_B = tol_sqrtg = 1e-5, tol_flux_funcs = 1e-7, tol_dB_ang = 1e-3, tol_dB_ds = 1e-2.
• make golden unchanged: test/golden/expected/rabe.nc still reproduced at the rtol=1e-10 in test/golden/compare_golden.py:37 — the type(boozer_field_t) → class(field_t), allocatable conversion must not perturb the vmec_nc path.
• The MRE in "Reproduction / evidence" runs to completion and writes finite coefficients.
• With vmns = 0, neo2_nabla_s = .true. and .false. agree to 1e-14 relative; with vmns /= 0, reading the -vmns file with .false. equals reading the +vmns file with .true. to 1e-14 relative.
• grep -rn "module bc_file" src test matches nothing in rabe (only libneo owns that name).
• grep -rn "boozer_field_init" src test app matches nothing (the branch's stale API name is gone).
• pre-commit run --all-files clean.

References

• Boozer, A. H., Plasma equilibrium with rational magnetic surfaces, Phys. Fluids 24, 1999 (1981) — the coordinate system and the sqrt(g) ∝ (G + iota*I)/B^2 Jacobian identity the issue's third bullet invokes. This is the one physics citation the change actually rests on.
• Kernbichler, W. et al., Solution of the drift kinetic equation in stellarators and tokamaks, Plasma Phys. Control. Fusion 58, 104001 (2016) — NEO-2, the reference code whose .bc input this issue is about matching.

Everything else here is an implementation detail internal to rabe and libneo and has no literature reference. In particular, no claim is made about how ASCOT5, BEAMS3D, SFINCS, SIMSOPT, DESC or GPEC represent their fields: an earlier draft of this comment asserted that "BEAMS3D and SFINCS consume boozmn/.bc harmonics directly rather than re-splining |B| on a 3D grid" and that booz_xform uses "radially interpolated coefficients". Neither was verified, both are probably wrong (BEAMS3D works on a cylindrical R/phi/Z spline grid; booz_xform transforms surface-by-surface with no radial interpolation), and neither is needed to justify the change. Do not reintroduce them.

Scope and risk

Two PRs: (a) the rebase-and-land of feat/bc-field with the field_type = 'bc' reconciliation, the boozer_field_init → init_from_vmec fixup, and the bc_writer rename — one reviewable PR, the code already exists; (b) the boozmn/chartmap nabla_s fix in libneo, which is independent and should not block (a). The one real regression risk in (a) is the type(boozer_field_t) → class(field_t), allocatable conversion in app/main.f90, which the golden record covers directly.

Related open issues: itpplasma/libneo#305 (the standalone evaluator; its reader half appears to have landed as libneo PR #365 — confirm before closing anything), #92 (the converter-fidelity measurement whose actual conclusion motivates this issue; see the Purpose section for what it does and does not say), #103 (the sibling compute_nabla_s accuracy question on the VMEC path), itpplasma/NEO-RT#44 (the other consumer named in libneo#305).