Migrate compute_nabla_s to the direct VMEC metric and regenerate golden

[krystophny]

Context

PR #100 cut RABE over to the merged libneo boozer_sub. To keep the golden record valid with no logic change, libneo restored the splined sqrt_g_ss output gated on use_B_r (itpplasma/libneo#335), and compute_nabla_s reads that splined value, as before.

The more accurate alternative

nabla_s = |grad s| can be computed directly from the VMEC metric at the Boozer->VMEC back-mapped angles, bypassing the field3d spline interpolation of sqrt_g_ss:

call boozer_to_vmec(s, theta_B, phi_B, theta_V, phi_V)
call splint_vmec_data(s, theta_V, phi_V, ...)        ! R, dR/d*, dZ/d*
call metric_tensor_vmec(...)                          ! g_vmec, sqrt_g_vmec
nabla_s = sqrt(g(2,2)*g(3,3) - g(2,3)**2)/abs(sqrt_g_vmec)

All three calls are public in libneo boozer_sub / spline_vmec_sub. This drops a spline-interpolation layer, so it is strictly closer to the true metric.

Difference and tradeoff

• At the golden grid (grid_refinement=6, multharm=6) the splined and direct values differ by ~1.5e-6 — the exact uniform drift seen on Use libneo boozer_sub; drop rabe local converter #100, since every output carries the dr_dAtheta proportional to 1/nabla_s prefactor.
• At coarse grids the gap is larger (~1e-2 at multharm=3).
• Cost: a boozer_to_vmec Newton solve per nabla_s evaluation.

Action when adopted

1. Replace compute_nabla_s with the direct-metric path (already prototyped on Use libneo boozer_sub; drop rabe local converter #100 history).
2. Regenerate the golden record (make golden_update) on the pinned nix toolchain — outputs shift by ~1.5e-6.
3. Revisit boozer_sub: restore sqrt_g_ss output, gated on use_B_r libneo#335: if RABE is the only sqrt_g_ss consumer, the gated restore can be dropped from boozer_sub.

[krystophny]

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

Still open and still actionable. compute_nabla_s at src/vmec/boozer_field.f90:349 reads the splined sqrt_g_ss exactly as the issue describes, and its body is unchanged since #100 (e1142da). Two corrections to the issue text, and one new finding:

1. The splined path carries three interpolation stages, not one.
2. #108 (4f70aa4) has since added booz_xform and chartmap field types, and those paths appear to leave nabla_s unassigned today. That turns this into a two-part change.

Paths below are prefixed [rabe] or [libneo]. libneo line numbers are quoted at the pinned revision 54e372f2ef15cf6cf0f36d734500c64d10cd5ada ([rabe] cmake/libneo.cmake:29), which is what RABE actually builds against.

Purpose

Every RABE output is proportional to dr_dAtheta, and dr_dAtheta is a pure function of the surface-averaged nabla_s:

! [rabe] src/offset/coefficients.f90:114
dr_dAtheta = real(sign_sqrtg, dp)/(average%nabla_s*psi_tor_edge)

So any interpolation error in nabla_s shows up as a uniform multiplicative drift on Lambda_A, Lambda_B, Lambda_S and nu_star_crit (the four 1-D output variables, app/main.f90:146-164) — it never averages out, and it is indistinguishable from a physics error in a NEO benchmark. sqrt(g^{ss}) is available in closed form from the VMEC metric at any point, so paying an interpolation floor for it is a pure loss of accuracy. This is the same class of fix as #92 ("spline floor of field-line |B| maxima; evaluate maxima spectrally"): move a quantity that has an exact evaluator off the interpolation chain, so the error floor is set by the Fourier/VMEC representation rather than by RABE's grid.

Current state

RABE side. [rabe] src/vmec/boozer_field.f90:349-369:

    subroutine compute_nabla_s(self, theta, phi, nabla_s)
        class(boozer_field_t), intent(in) :: self
        real(dp), intent(in) :: theta, phi
        real(dp), intent(out) :: nabla_s

        real(dp) :: dummy(32), sqrt_g_ss

        if (.not. self%fixed_to_surface) &
            error stop "compute_nabla_s: call fix_to_surface first"

        call splint_boozer_coord(self%fixed_stor, theta, phi, 0, &
                                 dummy(1), dummy(2), dummy(3), &
                                 ... &
                                 sqrt_g_ss)
        nabla_s = sqrt_g_ss/cm2m     ! :367

use_B_r = .true. is set in exactly one place, [rabe] src/vmec/boozer_field.f90:71 inside init_from_vmec (imported at :63). init_from_boozmn (:85-93) and init_from_chartmap (:100-108) do not set it.

libneo side. The exact evaluator already exists — it is what builds the grid in the first place, [libneo] src/boozer_converter.F90:1006-1031:

    function get_sqrt_g_ss_contravariant(s, theta, varphi) result(sqrt_g_ss)
        use spline_vmec_sub, only: splint_vmec_data
        use spline_vmec_sub, only: metric_tensor_vmec
        ...
        sqrt_g_ss = sqrt(g_vmec(2, 2)*g_vmec(3, 3) - g_vmec(2, 3)**2) &
                    /abs(sqrt_g_vmec)
    end function get_sqrt_g_ss_contravariant

It is private to boozer_sub (the public list, [libneo] src/boozer_converter.F90:16-25, does not contain it) — so RABE must reimplement those lines locally, or libneo must export it. boozer_to_vmec is public (:20), and splint_vmec_data / metric_tensor_vmec live in spline_vmec_sub ([libneo] src/vmec_support/spline_vmec_data.f90:654, :1193), which has no private statement, so both are accessible.

How many stages the splined path actually has. get_sqrt_g_ss_contravariant is evaluated on the VMEC-angle grid ([libneo] src/boozer_converter.F90:828-829), then Lagrange re-interpolated onto the Boozer poloidal grid (:903, perqua_t(7,...)), then onto the Boozer toroidal grid, then stored (:953) and finally packed into the 3D batch spline (:1244) that splint_boozer_coord evaluates. The issue text says "drops a spline-interpolation layer"; it actually drops two Lagrange re-interpolations and a 3D spline evaluation. That strengthens the case and explains why the coarse-grid gap (~1e-2 at multharm=3) is so much larger than the fine-grid one.

Build wiring is already in place. [rabe] src/vmec/CMakeLists.txt:6-8 already links LIBNEO::vmec_support, and that target compiles spline_vmec_data.f90 ([libneo] src/vmec_support/CMakeLists.txt:1-7), so use spline_vmec_sub needs no CMake change.

Latent bug found while checking scope. splint_boozer_coord assigns its optional sqrt_g_ss_B only under a gate, at [libneo] src/boozer_converter.F90:296 and :346:

            if (present(sqrt_g_ss_B) .and. boozer_state%use_B_r) &
                sqrt_g_ss_B = y_eval(2)

For the chartmap path this is not merely unset by omission — build_boozer_from_chartmap explicitly sets use_B_r = .false. ([libneo] src/boozer_converter.F90:1367) and builds a single-quantity spline, field3d_num_quantities = 1 (:1410), so y_eval(2) is not even populated. [libneo] src/field/boozmn_reader.f90 contains no use_B_r reference at all and so leaves the module default .false. ([libneo] src/vmec_support/canonical_coordinates_mod.f90:103).

Therefore for field_type = 'booz_xform' or 'chartmap', RABE's local sqrt_g_ss in compute_nabla_s is never written — it is an uninitialized local. A full run should then abort at [rabe] src/fieldline/surface_average.f90:67-69 with "error: average sqrt_g11 is NaN.", or produce garbage. No test catches this: [rabe] test/integration/vmec/test_chartmap.f90 checks only Bmod at three grid nodes, and the only three compute_nabla_s call sites under test/ (test_surface_area.f90:55, test_against_neo_qa.f90:204, test_against_neo_qh.f90:202) all use init_from_vmec. Reproduce this before writing code — it changes the shape of the fix.

The prototype referenced by the issue is not in this checkout. git log --all -S'boozer_to_vmec' returns nothing, so the "#100 history" version is not in reachable history. Write it fresh.

Reproduction / evidence

[rabe] test/integration/vmec/test_surface_area.f90 is already the right oracle: it integrates abs(sqrtg*nabla_s) over a surface and compares to simsopt references (:38-42), at grid_refinement=3 (:32) — exactly the coarse regime where the issue claims a ~1e-2 gap. Its tolerance is reltol = 1e-3 (:9), loose enough to hide it.

cd /path/to/rabe
make CONFIG=Release
ctest -R TestSurfaceArea --test-dir build/test --output-on-failure   # passes today at 1e-3

Quantify the gap before changing anything. Add a scratch program that, for s = 0.5 on test/integration/vmec/input/wout_LandremanPaul2021_QA_reactorScale_lowres_reference.nc, loops over a 32x32 (theta_B, phi_B) grid and prints, per point, the splined value from compute_nabla_s next to the direct value from boozer_to_vmec + splint_vmec_data + metric_tensor_vmec. Expected, per the issue: max relative difference O(1e-2) at grid_refinement=3, O(1.5e-6) at grid_refinement=6. Attach that table to the PR — it is the evidence that justifies regenerating the golden record.

Reproduce the chartmap NaN:

cd build/test/integration/vmec        # `input/` is copied here by CMakeLists.txt:26
printf '&rabe_config\n field_type = "chartmap",\n field_file = "input/circ_chartmap.nc",\n M_pol = -1.0d0, N_tor = 4.0d0, s_tor = 0.5,\n sign_sqrtg = -1.0d0, max_n_fieldlines = 20, n_eta = 50\n/\n' > rabe.in
<path-to>/rabe.x     # expect: "error: average sqrt_g11 is NaN."

Golden run and update ([rabe] Makefile:42-49, fixture test/golden/input/rabe.in, comparison at rtol=1e-10, atol=0.0 in test/golden/compare_golden.py:37-38):

make golden           # ctest -L golden; expect GoldenRecordCompare to FAIL by ~1.5e-6 after the change
make golden_update    # cp build/test/golden/rabe.nc test/golden/expected/

Run make golden_update only inside nix develop on the pinned toolchain (flake.nix:4, gfortran 11 / glibc 2.35). #90 traced ~4.6e-7 cross-platform drift in Lambda_A to adaptive quadrature amplifying libm differences; regenerating on a host toolchain bakes that drift into the committed record.

Behaviour (BDD)

Scenario: direct metric agrees with the spline on a well-resolved grid
  Given wout_LandremanPaul2021_QA_reactorScale_lowres_reference.nc loaded via
        init_from_vmec(..., radial_spline_order=5, angular_spline_order=5,
                       grid_refinement=6)
  And   the field fixed to s = 0.5
  When  compute_nabla_s is evaluated on a 32x32 (theta_B, phi_B) grid
  Then  every value is strictly positive
  And   every value agrees with the pre-change splined value to within 1e-5 relative

Scenario: direct metric beats the spline on a coarse grid
  Given the same wout file loaded via init_from_vmec(..., grid_refinement=3)
  When  test/integration/vmec/test_surface_area.f90 runs with reltol tightened to 1e-5
  Then  the surface area at all five s values (0.1, 0.3, 0.5, 0.7, 0.9) matches the
        simsopt references in test_surface_area.f90:38-42
  And   the same test with the old splined compute_nabla_s fails at that tolerance

Scenario: golden record shifts by the expected uniform factor and no more
  Given test/golden/input/rabe.in unchanged (QH wout, five s_tor values)
  When  make golden is run on the nix-pinned toolchain before regeneration
  Then  GoldenRecordCompare fails at rtol=1e-10 on Lambda_A, Lambda_B, Lambda_S
        and nu_star_crit
  And   every reported relative difference is below 1e-5
  And   the relative shift has the same sign for all four variables on a given
        surface, and the magnitudes agree within a factor 2

Scenario: golden record passes after regeneration
  Given test/golden/expected/rabe.nc regenerated with make golden_update inside
        nix develop
  When  make golden is run again on the same toolchain
  Then  GoldenRecordCompare passes at the unchanged rtol=1e-10, atol=0.0

Scenario: non-VMEC field types fail loudly instead of silently
  Given rabe.in with field_type = "chartmap" and field_file = "input/circ_chartmap.nc"
  When  rabe.x runs
  Then  it aborts with the message
        "compute_nabla_s: nabla_s requires field_type = 'vmec_nc'"
  And   it does not abort with "error: average sqrt_g11 is NaN."
  And   it does not produce finite-looking but uninitialized output

Implementation plan

1. [rabe] src/vmec/boozer_field.f90 — add a private module procedure
   nabla_s_direct(stor, theta_B, phi_B) result(nabla_s). Extend the existing
   use boozer_sub, only: list at :5 with boozer_to_vmec, and add
   use spline_vmec_sub, only: splint_vmec_data, metric_tensor_vmec at procedure
   scope (matching the file's existing style at :63). Body — the argument
   lists are exactly those of the libneo private original, so copy them verbatim:

       real(dp) :: theta_V, varphi_V
       real(dp) :: dummy(10)
       real(dp) :: R, dR_ds, dR_dtheta, dR_dphi
       real(dp) :: dZ_ds, dZ_dtheta, dZ_dphi
       real(dp) :: g_vmec(3, 3), sqrt_g_vmec
   
       call boozer_to_vmec(stor, theta_B, phi_B, theta_V, varphi_V)
       call splint_vmec_data(stor, theta_V, varphi_V, &
                             dummy(1), dummy(2), dummy(3), dummy(4), dummy(5), &
                             R, &
                             dummy(6), dummy(7), &
                             dR_ds, dR_dtheta, dR_dphi, &
                             dZ_ds, dZ_dtheta, dZ_dphi, &
                             dummy(8), dummy(9), dummy(10))
       call metric_tensor_vmec(R, dR_ds, dR_dtheta, dR_dphi, &
                               dZ_ds, dZ_dtheta, dZ_dphi, g_vmec, sqrt_g_vmec)
       nabla_s = sqrt(g_vmec(2, 2)*g_vmec(3, 3) - g_vmec(2, 3)**2) &
                 /abs(sqrt_g_vmec)/cm2m

   Two things that are easy to get wrong:

   • The first argument of boozer_to_vmec and splint_vmec_data is the flux
     label s, not rho. Both routines take rho_tor = sqrt(r) internally
     ([libneo] src/boozer_converter.F90:433), so pass self%fixed_stor
     directly, as get_sqrt_g_ss_contravariant does at :828-829.
   • Keep the /cm2m. libneo VMEC geometry is in cm, so the raw quantity is
     1/cm and /cm2m (i.e. *100) converts to 1/m, exactly as the current
     line :367 does.

2. [rabe] src/vmec/boozer_field.f90 — add a logical :: use_direct_metric = .false. component to boozer_field_t (components at :28-33); set it .true. in init_from_vmec only. compute_nabla_s (:349) branches on it: nabla_s_direct when true, and in the .false. branch
   error stop "compute_nabla_s: nabla_s requires field_type = 'vmec_nc'".
   This replaces today's silent uninitialized read with a loud failure and keeps the change reviewable.

3. File the follow-up for the non-VMEC fallback. Do not fold it into this PR. The clean fix lives in libneo: have build_boozer_from_chartmap ([libneo] src/boozer_converter.F90:1319-1418) populate sqrt_g_ss_grid from the chartmap geometry and set use_B_r, or add a chartmap-native sqrt(g^{ss}) evaluator; the boozmn path needs the equivalent. Open one RABE issue and one libneo issue, and link them here.

4. [rabe] test/integration/vmec/test_nabla_s_direct.f90 (new) — standalone program in the repo idiom (program, use utils, only: not_same, error stop on failure; see test_chartmap.f90:44-60). Loads the QA reference wout at grid_refinement=6, fixes s = 0.5, sweeps a 16x16 (theta_B, phi_B) grid, and asserts (a) nabla_s > 0 everywhere, (b) agreement with hard-coded reference values generated once from the direct path at reltol = 1e-12, (c) round-trip consistency: boozer_to_vmec followed by vmec_to_boozer (both public, [libneo] src/boozer_converter.F90:20) returns the input angles to 1e-12.

5. [rabe] test/integration/vmec/CMakeLists.txt — register it in the existing three-line idiom:

   add_executable(test_nabla_s_direct.x test_nabla_s_direct.f90)
   target_link_libraries(test_nabla_s_direct.x vmec_lib utils_lib LIBNEO::boozer)
   add_test(TestNablaSDirect test_nabla_s_direct.x)
   set_tests_properties(TestNablaSDirect PROPERTIES LABELS "quick")

   (LIBNEO::boozer is needed only if the test calls boozer_to_vmec directly — compare the test_booz_xform_signs.x entry, which does.)

6. [rabe] test/integration/vmec/test_surface_area.f90:9 — tighten reltol from 1e-3 to whatever the measured direct-path residual against the simsopt references supports, with a factor ~3 of headroom (target 1e-5). Record the measured value in a comment. This is the regression lock that prevents a silent revert to the splined path.

7. make golden_update inside nix develop, then commit test/golden/expected/rabe.nc in the same PR as the source change, with the before/after difference table from the "Reproduction" section in the commit message.

8. [rabe] CLAUDE.md — the "Keeping Documentation in Sync" section triggers on new derived types, modules, make targets and test labels, so this change strictly requires no edit. Update the boozer_field.f90 bullet under "Libraries (src/)" anyway to record that compute_nabla_s uses the direct VMEC metric for vmec_nc and aborts for the other two loaders — that is exactly the kind of surprise a future reader needs.

9. libneo cleanup (separate PR, after this one lands). The final bullet of the issue — dropping the use_B_r-gated sqrt_g_ss restore from boozer_sub (boozer_sub: restore sqrt_g_ss output, gated on use_B_r libneo#335) — is only safe once step 3 is settled, because the chartmap/boozmn path may end up needing that same grid. Grep for other consumers before removing: [libneo] src/boozer_converter.F90:296, :346, :753, :775, :828, :903, :953, :972, :1205, :1244. Order: RABE PR -> confirm no consumer -> libneo removal PR -> bump the RABE pin at [rabe] cmake/libneo.cmake:29.

Acceptance criteria

• make test_all passes with no new failures.
• ctest -R TestNablaSDirect --test-dir build/test passes (new test from step 4).
• ctest -R TestSurfaceArea --test-dir build/test passes at the tightened reltol (step 6), and demonstrably fails at that tolerance with the old splined path — show both runs in the PR.
• make golden passes against the regenerated test/golden/expected/rabe.nc at the unchanged rtol=1e-10 (test/golden/compare_golden.py:37).
• The golden shift is documented: a table of old vs new Lambda_A, Lambda_B, Lambda_S, nu_star_crit for all five surfaces in test/golden/input/rabe.in, with relative differences at the ~1.5e-6 level and no sign flips.
• make golden_update was run inside nix develop; state the toolchain in the PR body.
• field_type = "chartmap" and field_type = "booz_xform" abort with the explicit error stop from step 2, and the step 3 follow-up issues are linked.
• test/integration/vmec/test_against_neo_qa.f90 and test_against_neo_qh.f90 still pass. Note these are weak evidence, not a discriminator: their reltol is 1.4e-2 (test_against_neo_qa.f90:11), four orders above the expected shift, and the live nabla_s_ref values are the ones at :99-104 (the :95-98 assignment is dead, immediately overwritten). Their header states the references come from "the NEO code" via a converted .bc file. Treat them as a non-regression guard only. If you want a real independent check, re-run the internal NEO-2 comparison suite (test/internal, label internal, git submodule) and report whether the direct path moves nabla_s toward those references.

Scope and risk

Steps 1-2 and 4-8 are a single reviewable PR (roughly 150 lines of source plus the regenerated golden .nc); step 3's libneo work and step 9's sqrt_g_ss removal are separate PRs with a strict ordering dependency.

Main risks:

• Cost. The Newton solve in boozer_to_vmec ([libneo] src/boozer_converter.F90:635-680, niter = 100, epserr = 1.0e-14) now runs once per nabla_s evaluation, including every QUADPACK quadrature node reached through nabla_s_over_B_squared ([rabe] src/fieldline/fieldline_integrands.f90:63-73). Watch .github/workflows/benchmark.yml. If runtime regresses badly, cache the back-mapped VMEC angles per (theta, phi) rather than reverting.
• Golden regeneration environment. Must happen on the pinned toolchain or Cross-platform golden-record drift in Lambda_A traced to adaptive quadrature + libm (env, not code) #90 resurfaces as a spurious cross-platform failure.
• Scope creep into libneo. Resist it; step 3 is a separate issue by design.

Related open issues: #92 (same "replace a splined value with exact evaluation" pattern, applied to field-line |B| maxima — the two together should move the whole error floor), #30 (slow convergence of <sqrt(g11)> — the same quantity, but a quadrature question rather than an interpolation one; this change will alter the numbers #30 measures), #90 (golden-record drift; read before regenerating), and itpplasma/libneo#335 (the gated sqrt_g_ss restore this issue proposes to retire).

References

• Grassler, G. S., Albert, C. G., Kasilov, S. V. & Kernbichler, W. (2026), Asymptotic modeling of the bootstrap and Ware pinch effect in near omnigenous stellarators, Phys. Plasmas 33(7), 070703, doi:10.1063/5.0332431 — the paper RABE implements; defines Lambda_A and Lambda_B (README.md:331).
• Albert, C. G. et al. (2025), On the convergence of bootstrap current to the Shaing-Callen limit in stellarators, J. Plasma Phys. 91(3), E77, doi:10.1017/S0022377825000200 — the finite-boundary-layer correction Lambda_S computed at src/offset/coefficients.f90:23-72 (README.md:333).
• Boozer, A. H. (1981), Plasma equilibrium with rational magnetic surfaces, Phys. Fluids 24, 1999 — the coordinate system and the Jacobian/metric relations sqrt(g^{ss}) is derived from.
• Hirshman, S. P. & Whitson, J. C. (1983), Steepest-descent moment method for three-dimensional magnetohydrodynamic equilibria, Phys. Fluids 26, 3553 — VMEC; the R(s,theta,phi), Z(s,theta,phi) representation that splint_vmec_data splines and metric_tensor_vmec differentiates.
• Landreman, M. & Paul, E. (2022), Magnetic fields with precise quasisymmetry for plasma confinement, Phys. Rev. Lett. 128, 035001 — the QA/QH equilibria used by the golden record and the surface-area test.

Note on other codes: libneo's own source comment ([libneo] src/boozer_converter.F90:1230-1232) records that SIMPLE leaves use_B_r off and consumes the lean single-quantity Bmod spline, and that sqrt_g_ss and B_r exist in the batch spline only for use_B_r callers such as RABE's nabla_s. That asymmetry is what this issue removes on the RABE side. No claim is made here about how booz_xform, DESC, SIMSOPT or any other external code handles the same quantity — that was not verified.