Equilibrium - BUGFIX - Converge PEST Δ′ in the radial grid #524

Problem

On a direct (traced) equilibrium, Δ′ in PEST coordinates does not converge as mpsi grows. The same equilibrium in hamada coordinates does. Δ′ at Mercier index α is sensitive to coefficient error at the surface, amplified by roughly x^(−2α), so a small geometry error that is invisible in the energies shows up in Δ′.

Solovev, q0 = 1.2 (one rational at q = 2, α = 0.886), wall a = 0.05, etol 1e-10, eulerlagrange_tolerance 1e-12, Riccati, ldp grid. Run on #398's head (81006c3a6):

• Hamada moves 0.6% from 256 to 1024. PEST moves 15%, in the wrong direction.
• et₁ agrees across all six runs to 2e-6, so this is a Δ′-only defect.
• The exact-coefficient reference for this case is −150.73 (local analytic patch, settings differ slightly).

Likely sources

1. J from a spline θ-derivative. PEST forms the Jacobian as (1 + ff_deriv_val[4]) * y_out[end, 2] * 2π * psio (DirectEquilibrium.jl, equilibrium_solver), i.e. by differentiating a periodic cubic spline of ∫dl/B_p. Hamada does not need that derivative. Before Equilibrium - BUGFIX! - Sample flux surfaces on common straight-fieldline angles, cap core EL knots, and suppress near-axis kinetic terms #398, ⟨∂J/∂ψ⟩ from this formula was off by 88% rms against the exact Solovev value.
2. θ resolution of the traced surface. After Equilibrium - BUGFIX! - Sample flux surfaces on common straight-fieldline angles, cap core EL knots, and suppress near-axis kinetic terms #398 each surface is sampled only at the mtheta + 1 straight-fieldline nodes. A prototype that fixed this case to ≤1e-4 across mpsi 128–2048 in both PEST and hamada used 2048 field-line points per surface plus direct J. On Equilibrium - BUGFIX! - Sample flux surfaces on common straight-fieldline angles, cap core EL knots, and suppress near-axis kinetic terms #398, direct J alone recovers about a third of the gap at mpsi 1024.

Suggested approach

• Compute J directly at each θ node from the field there: J = 2π ψ₀ ∮(w dl/B_p) / w(θ), where w is the Jacobian weight (power_bp, power_b, power_r, power_rc). Snippet below; it is a drop-in replacement for the one line above.
• Then find the remaining error: rerun the table with mtheta doubled and quadrupled, and with a tighter dense-output tolerance on the field-line solve. Whichever closes the PEST/hamada gap at mpsi 1024 is the second half of the fix.
• Check arclength_fieldline_int (DirectEquilibriumArcLength.jl) gets the same J.

rfac_j = sqrt(max(ff_val[1], 0.0)); eta_j = 2π * (theta + ff_val[2])
r_j = ro + rfac_j * cos(eta_j); z_j = zo + rfac_j * sin(eta_j)
direct_get_bfield!(bf_j, r_j, z_j, raw_profile.psi_in, raw_profile.sq_in, sqd_j, psio; derivs=1)
bp_j = sqrt(bf_j.br^2 + bf_j.bz^2); b_j = sqrt(bp_j^2 + (bf_j.f / r_j)^2)
w_j = bp_j^equil_params.power_bp * b_j^equil_params.power_b / (r_j^equil_params.power_r * rfac_j^equil_params.power_rc)
rzphi_nodes[ipsi, itheta, 4] = y_out[end, 5] / w_j * 2π * psio

(bf_j = DirectBField() and sqd_j = deriv1(raw_profile.sq_in) are allocated once before the ψ loop.)

Acceptance

• Solovev α 0.89: PEST and hamada Δ′ agree to ≤1e-3 relative at mpsi 256 and 1024, and each is flat in mpsi to the same level.
• DIII-D-like n = 1 Riccati: PEST and hamada Δ′ agree on all five surfaces once Equilibrium - BUGFIX! - Build F and P of g-file and IMAS equilibria from the tabulated FF′ and p′ #506's profiles are in.
• Every direct-equilibrium PEST result moves, so the PR takes ! and re-pins the harness.

Dependencies

• Stacks on Equilibrium - BUGFIX! - Sample flux surfaces on common straight-fieldline angles, cap core EL knots, and suppress near-axis kinetic terms #398 (common straight-fieldline sampling; approved, needs its conflicts resolved). Equilibrium - BUGFIX! - Sample flux surfaces on common straight-fieldline angles, cap core EL knots, and suppress near-axis kinetic terms #398 removes the sparse-sampling noise that dominated before; this issue is what is left.
• Independent of Equilibrium - BUGFIX! - Build F and P of g-file and IMAS equilibria from the tabulated FF′ and p′ #506 (Equilibrium - g-file reader differentiates FPOL and PRES and discards FFPRIM and PPRIME, putting rounding noise on Δ′ #504), ForceFreeStates - BUGFIX! - Use trace-aware Frobenius exponents at resonant surfaces #517 and the crossing work (ForceFreeStates - BUGFIX - Decouple the ideal-crossing jump width from singfac_min #481).