TestAgainstNEO_QH zero-component abstol is toolchain-sensitive

[krystophny]

Goal

TestAgainstNEO_QH asserts an analytically-zero spline derivative component against abstol = 1e-10, while the sibling symmetry-zero check uses abstol_for_zero = 1e-3. On gfortran 16.1.1 (Arch/CachyOS) the component evaluates to 1.577e-10 and the test fails; Ubuntu CI (gfortran 13/14) stays below the bound. The failure reproduces bit-identically on a clean build of current main with the pinned libneo, independent of any open PR, so the bound encodes toolchain-specific roundoff rather than physics.

Decide the intended tolerance semantics for exact-zero reference components in this test and apply them consistently (e.g. route the zero components through the same near-zero floor the dB_dx check uses, or document the 1e-10 bound as CI-toolchain-specific).

Reproduction

cmake -S . -B build -G Ninja && cmake --build build -j
ctest --test-dir build/test -R TestAgainstNEO_QH --output-on-failure

On gfortran 16.1.1: dB/dphi = 1.5772720057194170E-010 vs reference 0.0, abstol = 1e-10 -> FAIL.

Non-goals

• No blanket tolerance widening; only the exact-zero reference components are affected.

[krystophny]

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

The specific failure reported here — TestAgainstNEO_QH tripping on dB/dphi = 1.577e-10 against abstol = 1e-10 — is already fixed on main, by 4f70aa4 (2026-07-03, "Support booz_xform and chartmap field input in addition to .bc (#108)"), one day after this issue was filed. That commit added

! Symmetry-zero dB_dx components (exact 0 in the NEO reference) come out
! as spline roundoff, measured <= 1.6e-10 with an optimized libneo build.
real(dp), parameter :: abstol_roundoff = 1e-9

(test/integration/vmec/test_against_neo_qh.f90:14-16) and switched the dB_dx assertion at test/integration/vmec/test_against_neo_qh.f90:191-192 from abstol_in=abstol to abstol_in=abstol_roundoff.

I am not closing it, because the second half of the stated goal — "decide the intended tolerance semantics for exact-zero reference components and apply them consistently" — is demonstrably not done, and the sibling QA test still carries the identical latent trip-wire. The spec below is scoped to that residual only.

Purpose

Both NEO comparison tests validate boozer_field_t against NEO reference field evaluations on quasi-symmetric configurations. At the stellarator-symmetry point (theta = 0, phi = 0) the angular derivatives of |B| vanish exactly, so the reference values are hard zeros and the code returns pure interpolation roundoff. A hard absolute bound on such a component tests the compiler and libm, not the field representation — exactly the environment-versus-code distinction already established in #90. Right now the repo holds three different, undocumented policies for the same situation across two neighbouring files, so the next toolchain bump will re-open this as a new "mystery test failure" rather than being caught by an explicit, documented floor.

Current state

The three inconsistent policies:

1. QH, B_theta (analytically zero, refs ~1e-16): reltol_in=0.0_dp, abstol_in=abstol_for_zero with abstol_for_zero = 1e-3 (test_against_neo_qh.f90:13, :144-145). A 1e-3 floor on a quantity whose true value is 1e-16 is thirteen orders of slack — it would not catch a real regression. QA does the same with abstol_for_zero = 1e-4 (test_against_neo_qa.f90:12, :145-146).
2. QH, dB_dx: blanket abstol_roundoff = 1e-9 applied to the whole 3-vector (test_against_neo_qh.f90:191-192), which widens the floor for the non-zero components too. This is harmless in practice (the smallest non-zero QH reference component is 3.85688857474143776e-03, test_against_neo_qh.f90:90, so reltol = 1.2e-2 dominates) but it contradicts this issue's own non-goal, "no blanket tolerance widening".
3. QA, dB_dx: the third component is sliced out entirely —

! Skip dB/dphi (component 3) — for QA it is at the spline noise floor
if (not_same(dB_dx_b(1:2), dB_dx_ref(case, 1:2), &
             reltol_in=reltol, abstol_in=abstol)) then

(test_against_neo_qa.f90:192-194) — while dB_dx_ref(1, 2) = 0.0 (test_against_neo_qa.f90:79-81) is still compared against abstol = 1e-10 (test_against_neo_qa.f90:11). That is the same assertion that failed here for QH, still live in QA, and it is only untripped because CI never runs a modern compiler: .github/workflows/test.yml:25-34 runs both the build and make test_all inside nix develop .#default, and flake.nix:5-7 pins nixpkgs to nixos-22.11 (gfortran 11, glibc 2.35, per the comment on flake.nix:5-6).

not_same_array (src/utils/utils.f90:71-95) is an any() over abs(a1-a2) > reltol*abs(a1) + abstol, so a single component decides the whole check; there is no per-component tolerance vector.

Candidate primitive for the floor, and its caveat

The repo already exposes a deferred binding for field accuracy: field_base declares

!! \brief Relative accuracy of B_mod evaluations for this field representation.
real(dp) function rel_accuracy_B(self)

(src/field/field_base.f90:57-63, deferred at :16), and boozer_field_t implements it as a constant 1e-9_dp (src/vmec/boozer_field.f90:371-375, bound at :46, so bfield%rel_accuracy_B() is callable from the tests). The derivative returned to the tests is dB_dx = dlnB_dx*B_mod (src/vmec/boozer_field.f90:311 and :329).

Be honest about what this is. rel_accuracy_B is documented as the relative accuracy of B_mod, not of its angular derivatives; spline derivative error is generally worse than the value error. So rel_accuracy_B()*B_mod is a stated convention with a single named knob, not a derivation. It is still a strict improvement over three unexplained magic constants, and numerically it is comfortably above the observed roundoff: ~5.58e-9 at QH case 1 (bmod_ref(1) = 5.58389078029500130e+00, test_against_neo_qh.f90:70) and ~5.70e-9 at QA case 1 (bmod_ref(1) = 5.70341902299086545e+00, test_against_neo_qa.f90:66), versus the reported 1.577e-10.

Note the tension with #92, which measured a field-line |B| maxima spline error of 7.6e-8 (QA, cubic, median) — two orders above the hard-coded 1e-9. Do not silently treat 1e-9 as measured truth. Either keep it and say in the comment that it is a convention, or (out of scope here) open a follow-up to re-derive the constant.

Reproduction / evidence

CI runs nix develop .#default --command make test_all (.github/workflows/test.yml:34), i.e. ctest ... -L "quick|slow" (Makefile:33-34); TestAgainstNEO_QA and TestAgainstNEO_QH both carry label quick (test/integration/vmec/CMakeLists.txt:11, :16). Reproduce the residual on a non-pinned toolchain:

# do NOT use `nix develop`; the point is to use the system gfortran
gfortran --version                       # e.g. 16.1.1
cmake -S . -B build -DCMAKE_BUILD_TYPE=Debug && cmake --build build -j
ctest --test-dir build/test -R 'TestAgainstNEO_Q[AH]' --output-on-failure

Fixtures test/integration/vmec/input/wout_LandremanPaul2021_QA_reactorScale_lowres_reference.nc and ..._QH_...nc are copied into the build tree by test/integration/vmec/CMakeLists.txt:1. cmake/libneo.cmake fetches libneo at pin 54e372f2ef15cf6cf0f36d734500c64d10cd5ada unless -DLIBNEO_PATH= is given, so the build needs network access.

Measure before changing anything. Add a temporary print *, case, dB_dx_b immediately before the dB_dx check at test_against_neo_qa.f90:192, run TestAgainstNEO_QA, and record dB_dx_b(2) and dB_dx_b(3) at case 1 (stor = 0.1, theta = 0.0, phi = 0.0, from test_against_neo_qa.f90:44-46). Do the same for B_theta_b at test_against_neo_qa.f90:145 and test_against_neo_qh.f90:144. Paste all four numbers into the PR description; they are the inputs to steps 4 and 5 and the justification for whatever bound you write down.

Behaviour (BDD)

Scenario: symmetry-zero derivative components are held to a documented floor
  Given test/integration/vmec/test_against_neo_qh.f90 with dB_dx_ref(1,2) = dB_dx_ref(1,3) = 0.0
  When TestAgainstNEO_QH is built with a host gfortran >= 15 and run
  Then the dB_dx check at case 1 compares those two components against
       bfield%rel_accuracy_B()*bmod_ref(1) = 5.58e-9, not against 1e-10
  And the test passes

Scenario: the same policy holds in the QA sibling
  Given test/integration/vmec/test_against_neo_qa.f90 currently slices dB_dx_b(1:2)
  When the slice is removed and all three components are checked
  Then dB_dx_b(3) at every case is compared against the same zero-floor rule
  And TestAgainstNEO_QA passes on both gfortran 11 (nix pin) and the host compiler

Scenario: non-zero components keep their original strictness
  Given dB_dx_ref(3,2) = 3.85688857474143776e-03 in test_against_neo_qh.f90
  When that reference value is multiplied by 1.05 in a scratch build
  Then TestAgainstNEO_QH fails, proving the floor did not swallow a real regression

Scenario: pinned CI is unchanged
  Given the nixos-22.11 / gfortran 11 devShell in flake.nix
  When `nix develop .#default --command make test_all` runs
  Then all previously passing tests still pass and the golden record is untouched

Implementation plan

1. src/utils/utils.f90 — add a third specific procedure not_same_array_zero_aware and register it in the existing not_same generic interface (src/utils/utils.f90:6-9):

   function not_same_array_zero_aware(array_1, array_2, reltol_in, zero_floor_in)

   Elementwise over i: if array_2(i) == 0.0_dp, the element fails when abs(array_1(i)) > zero_floor; otherwise it fails when abs(array_1(i) - array_2(i)) > reltol*abs(array_2(i)). Aggregate with any(), matching not_same_array (:92-93). Return .true. when either argument contains a NaN, reusing contains_nan (src/utils/utils.f90:97-111).
   Note the deliberate change from not_same_array, which scales reltol by abs(array_1) — the computed value. Scaling by the reference is the correct convention when the reference is the trusted quantity. Because the new procedure has a distinct dummy-argument name (zero_floor_in vs abstol_in), generic resolution stays unambiguous; call it with keyword arguments.
   utils_lib (src/utils/CMakeLists.txt:1) is already linked by both NEO test targets (test/integration/vmec/CMakeLists.txt:9, :14), so no build wiring is needed. Do not put this in test/helpers/ — helpers_lib (test/helpers/CMakeLists.txt:1) is not linked by the NEO tests, and adding it would drag in field_lib, myplot_lib and friends for a pure numeric helper.

2. test/unit/test_not_same_zero_floor.f90 (new) — unit-test the helper directly:

   • exact-zero reference, computed 1.577e-10, zero_floor = 5.6e-9 → passes;
   • the same with zero_floor = 1e-10 → fails;
   • non-zero reference 3.85688857474143776e-03, computed 4.04973300347850965e-03 (+5 %), reltol = 1.2e-2 → fails (|diff| = 1.93e-4 > 4.63e-5);
   • non-zero reference 3.85688857474143776e-03, computed 3.86e-3, reltol = 1.2e-2 → passes;
   • NaN in either argument → fails.
     Register in test/unit/CMakeLists.txt with add_executable / target_link_libraries(... utils_lib) / add_test / set_tests_properties(TestNotSameZeroFloor PROPERTIES LABELS "quick"), following the pattern of the neighbouring entries (e.g. test/unit/CMakeLists.txt:1-4).

3. test/integration/vmec/test_against_neo_qh.f90 — delete abstol_roundoff (:16) and its comment (:14-15); declare real(dp) :: zero_floor and set it inside the case loop as zero_floor = bfield%rel_accuracy_B()*bmod_ref(case); replace the call at :191-192 with the new helper (reltol_in=reltol, zero_floor_in=zero_floor). Add a comment stating that rel_accuracy_B is the B_mod accuracy and is being used as the derivative floor by convention.

4. test/integration/vmec/test_against_neo_qh.f90, B_theta — replace the check at :144-145 with the same helper. B_theta_ref is numerically zero (~1e-16, :58-62), so it goes through the zero branch. Tighten abstol_for_zero (:13) to the measured B_theta_b magnitude times a safety factor of ~10, recorded in a comment; delete the constant if the helper's zero_floor covers it. Do the same in test_against_neo_qa.f90 (:12, :145-146) — note abstol_for_zero = 1e-4 there is live, so do not delete it without replacing the call site.

5. test/integration/vmec/test_against_neo_qa.f90, dB_dx — remove the dB_dx_b(1:2) slice and the "Skip dB/dphi" comment (:192-194), check all three components with the new helper, and update the four print * diagnostic lines at :195-200 to print the full 3-vectors again.

6. .github/workflows/test.yml — add a second job that does not use nix, so this class of regression is visible:

   test-host-gfortran:
     runs-on: ubuntu-latest
     steps:
     - uses: actions/checkout@v4
     - run: sudo apt-get update && sudo apt-get install -y gfortran cmake ninja-build pkg-config libnetcdf-dev libnetcdff-dev
     - run: cmake -S . -B build -DCMAKE_BUILD_TYPE=Debug && cmake --build build -j
     - run: ctest --test-dir build/test -R 'TestAgainstNEO_Q[AH]|TestStellaratorSymmetryQA' --output-on-failure

   The apt list mirrors the nix toolchain set (flake.nix:14-22) minus the Python golden-record deps, which this job does not need. Keep this job out of the golden label path — the golden record must stay on the pinned toolchain per flake.nix:5-7. If it goes red for unrelated reasons on ubuntu-latest, keep the job but do not let it gate merges until it has been green once.

7. CLAUDE.md — add one bullet under ### Testing (CLAUDE.md:17-28) stating the rule: analytically-zero reference components are compared against field%rel_accuracy_B() * B_mod via not_same(..., zero_floor_in=...), never against a bare absolute constant.

No libneo change is required. Do not touch test/golden/expected/rabe.nc or test/golden/compare_golden.py — this issue is strictly about the NEO field-evaluation comparisons.

Acceptance criteria

• nix develop .#default --command make test_all passes with no new failures (pinned gfortran 11 path unchanged).
• ctest --test-dir build/test -R 'TestAgainstNEO_Q[AH]' --output-on-failure passes on a host gfortran >= 15 build.
• TestAgainstNEO_QA checks all three dB_dx components; rg -n 'dB_dx_b\(1:2\)' test/integration/vmec/test_against_neo_qa.f90 returns nothing.
• No bare abstol = 1e-10 is applied to any reference component whose value is literally 0.0 in either NEO test.
• New unit test test/unit/test_not_same_zero_floor.f90 is registered with label quick and passes under make test.
• The four measured magnitudes from the "Reproduction / evidence" section are recorded in the PR description, and every remaining tolerance constant in both NEO tests has a comment giving its measured basis.
• Falsification check recorded in the PR: multiplying dB_dx_ref(3, 2) in test_against_neo_qh.f90 by 1.05 makes TestAgainstNEO_QH fail.
• make golden on the pinned toolchain still reports zero drift against test/golden/expected/rabe.nc.
• Added/changed lines follow the repo style (CLAUDE.md:106-110: 88 columns, 4-space indent).

References

No specific literature reference for the tolerance convention itself; this is an implementation detail internal to the code. The reltol*|reference| + abstol split with an explicit floor for analytically-zero quantities is ordinary numerical-testing practice (numpy.testing.assert_allclose's rtol/atol, Boost.Test's tolerance), not a published method.

The underlying symmetry — odd-parity quantities vanishing identically at the stellarator-symmetry point of a stellarator-symmetric equilibrium — is standard; the formalism is set out in R. L. Dewar and S. R. Hudson, "Stellarator symmetry", Physica D 112, 275 (1998), in the Boozer coordinates of A. H. Boozer, "Plasma equilibrium with rational magnetic surfaces", Phys. Fluids 24, 1999 (1981). Neither paper says anything about numerical tolerances; they are cited only for why the reference values are exact zeros.

Scope and risk

Single small PR (two test files, one helper procedure in src/utils/utils.f90, one unit test, one CI job). Main risk is the new non-pinned CI job being flaky on newer compilers for unrelated reasons — see step 6 for the mitigation. Closely related: #90 (golden-record drift is environment, not code — same root-cause class) and #92 (evaluating |B| and its maxima spectrally would drop the spline floor and change the right magnitude for these floors; if #92 lands first, re-measure before hard-coding anything).