Fix the sign of the cross terms in NESTOR's analytic add-back

[CharlesCNorton]

NESTOR subtracts a tangent-plane kernel numerically (RegularizedIntegrals) and adds its Fourier coefficients back analytically (SingularIntegrals). The analytic side used the opposite sign of the metric and curvature cross terms guv, auv, so mode (m, n) received the coefficient of (m, −n). The error vanishes on axisymmetric surfaces, is first order in the helical shaping, and does not decrease with resolution; VMEC 8.52 has the same sign in analyt.f90.

On an exact exterior Neumann problem (interior dipole line, helical torus) the residual drops from 4.9e-3 to 7e-6 and now converges with resolution; on the cth_like β=0.2% equilibrium the NESTOR boundary moves from 0.43 mm rms to 0.08 mm rms from a 10-digit BIEST solve of the same problem, and a W7-X standard-configuration equilibrium at β=2.2% moves by 1.7 mm on the axis and 1.3 mm rms on the boundary. educational_VMEC with the same two-line change reproduces the corrected results to 1e-11.

New tests: the assembled analytic coefficients against Gauss–Legendre quadrature of the subtracted kernels (singular_integrals_test, fails with the old sign) and the interior-dipole exact solution (nestor_test). The golden scalars in LasymFreeBoundaryMatchesEducationalVmec come from the patched educational_VMEC, the multigrid iteration guard is 328, and wout_cth_like_free_bdy.nc is regenerated (it matches the patched educational_VMEC at the reference test's tolerances).

Fixes the NESTOR part of #628.

[chatgpt-codex-connector]

Codex usage limits have been reached for code reviews. Please check with the admins of this repo to increase the limits by adding credits.
Credits must be used to enable repository wide code reviews.

[jurasic-pf]

Wow, well done on spotting this issue - you mention that you did a convergence check against the analytical solution, could you share some plots for it, since you already generated the data? In particular I'm curious how the B.n decays qualitatively (is/was there clear structures visible in the error? How did this go unnoticed?), and rms decay in log-log

[CharlesCNorton]

Plots from the exact test (a line of dipoles on the axis of a helically deformed torus; the exact vacuum field on the surface is zero, so the total field NESTOR reconstructs is the error). NESTOR imposes B·n itself, so the residual is the tangential field.

The old-sign residual is 6% of |B_ext| at every resolution and lives entirely in the low modes of the deformation, as (m, +1)/(m, −1) pairs of equal weight, which is the signature of the (m, n) ↔ (m, −n) swap. With the fix it drops to 0.1–0.4% and decays close to h⁻² up to nzeta 48, then plateaus around 1e-3 in |B| (1e-6 in |B|²) with an n = 0 structure; that floor is not the mgrid interpolation (a 601² grid gives the same numbers) and is three orders below the defect, so I have not chased it. rms |B_total| / rms |B_ext| for nzeta 18, 24, 36, 48, 72 (mpol = nzeta/3, ntor = nzeta/6): old 6.2e-2 throughout; fixed 3.9e-3, 2.4e-3, 1.2e-3, 7.6e-4, 1.1e-3.

It went unnoticed because the error vanishes on axisymmetric surfaces, is smooth and shares the symmetry of the geometry, moves ordinary equilibria by fractions of a millimetre, and Fortran VMEC and VMEC++ agree with each other since both carry it; only a solver-independent reference exposes it.

[jurasic-pf]

Haha, so fast! :)

[jurasic-pf]

Really great work catching that. I guess the remaining m=0 structure is primarily due to the expansion error of the singularity subtraction, what do you think?
And what's your intuition on why the error starts increasing around nzeta >50? I would suspect the derivatives of the geometric terms start to blow up around there so truncation error becomes dominant.

[CharlesCNorton]

@jurasic-pf Yes. The plateau is the auto grid tying ntheta to 2·mpol+6: at mpol 8, ntor 4, nzeta 24 the residual is 2.4e-3 / 7.2e-4 / 4.2e-4 for ntheta 22 / 44 / 88 (nzeta 24 → 96 only gives 4.2e-4 → 3.3e-4), and at ntheta 88 it is 4.2e-4 / 1.7e-4 / 1.5e-4 for mpol 8 / 12 / 16, about first order in ntheta as expected for a remainder that is only C⁰ at the singular point.

[CharlesCNorton]

@jurasic-pf The rise is mpol, not nzeta (the ladder ties mpol to nzeta/3): at fixed ntheta the residual is 1.5e-4 for mpol 16 and 5.7e-4 for mpol 24, mpol 32 diverges for both signs and also on the axisymmetric 1/R test, while nzeta 24 → 96 at fixed mpol improves it. That points at the analytic side at high l — the T_l recurrences and the alternating cmn sums with kL = mpol+1+ntor — rather than the geometry, whose derivatives are exact here (three Fourier modes).

[jurasic-pf]

The PR convincingly demonstrates that there's a bug that needs fixing, but I wonder where it actually is - NESTOR and VMEC use different handedness in the coordinates, so there's a few terms that are lacking a signgs flip. I wonder if the mixed sincos terms here

vmecpp/src/vmecpp/cpp/vmecpp/free_boundary/surface_geometry/surface_geometry.cc

Line 421 in 7ecfce8

snv[o] = signOfJacobian * (rub[d] * zvb[d] - zub[d] * rvb[d]);

or here

vmecpp/src/vmecpp/cpp/vmecpp/free_boundary/surface_geometry/surface_geometry.cc

Line 426 in 7ecfce8

guv[o] = 2.0 * (rub[d] * rvb[d] + zub[d] * zvb[d]) / s_.nfp;

have the sign convention backwards and we are just patching the consequence here

[CharlesCNorton]

@jons-pf Yes, it is the same issue. Negating guv and auv on the way into the analytic part exchanges ap ↔ am and Ap ↔ Am, i.e. T^± and S^±, which is exactly what the call-site swap does in the m ≥ 1, n ≥ 1 branch (the m = 0 / n = 0 branch uses T^+ + T^- and is unaffected); the explicit un-swap, implemented as a check, is bit-identical to this change. Intentional or not is settled by the exact solution: with the swap the assembled coefficients of mode (m, n) are those of the mirror kernel (m, −n), and NESTOR's residual on an exact exterior solution is 4.9e-3 at every resolution; without it 7e-6 and converging. If you want the form that mirrors the Fortran call, it can be written as the explicit swap instead.

[CharlesCNorton]

@jurasic-pf Checked: the stored guv equals 2 X_θ·X_ζ/nfp evaluated analytically on the helical torus to 3e-17, so greenf subtracts the physical kernel; the mismatch is on the analytic side, and the analysum2 argument swap jons found is its origin.

[CharlesCNorton]

fixing conflicts