Solve free boundary once and for all

[jurasic-pf]

NESTOR has it's problems, but is to date the best free boundary solver we have for VMEC. The B^2 pressure field it produces inherently doesn't assume any sheet currents by construction, and qualitatively "leads to smooth convergence dynamics". Fundamentally better solvers for the virtual casing / exterior field contribution exist, e.g. BIEST but have poor convergence dynamics when coupled to VMEC that have to be understood first.

BIEST is an interesting candidate, if we can make it converge robustly and penalize sheet currents.

See the verification calculations done in

• https://arxiv.org/abs/2412.05680
• https://pubs.aip.org/aip/pop/article-abstract/32/4/043906/3344075/Verification-of-the-Shafranov-shift-in-free?redirectedFrom=fulltext

In particular a high pressure equilibrium with large shafranov shift shpuld be a challenging test case, because the virtual casing contribution is then particularly large

Goal: Get a robust and accurate free boundary solver for VMEC

[CharlesCNorton]

on it 👍

[jurasic-pf]

One important verification should be that actual surface currents only appear when the pressure jump at the boundary is non zero. (not to be confused with the virtual surface currents from virtual casing)

[CharlesCNorton]

These two observations come from validating the AEGIS virtual-casing solver (#629) against BIEST, the reference recommended in that thread. AEGIS is superseded by the native BIEST and Vac2 solvers (#637, #638) and #629 is closed; the observations below are general and bear on how free-boundary solver quality is judged for this issue.

Field tangency, not fixed-boundary magnitude accuracy, governs in-loop quality. The free-boundary force couples the vacuum field to the boundary only through |B_ext|²/2, so it enforces pressure balance rather than n·B_ext = 0. A solver's in-loop behaviour therefore depends on whether its field is tangent by construction. NESTOR, Vac2, and BIEST are each tangent by construction, and on cth_like up to β ≈ 2% their magnetic axes agree to within ~2 mm. AEGIS evaluates the virtual-casing field without solving for tangency, so its field is not tangent: on a 3D NCSX case it matched BIEST's magnitude to 0.02% (direction cosine 0.999991) while its boundary normal-field residual was n·B/|B| = 2e-2 against BIEST's 1e-13, and the coupled equilibrium's axis was displaced by 4% (R_axis 1.591 vs 1.530).

delbsq does not rank solvers by correctness. On that NCSX case AEGIS converged to a lower delbsq than NESTOR (2.4e-3 vs 3.4e-3), consistent with the delbsq improvement reported for AEGIS earlier in #629, but at the displaced, less-tangent boundary above. delbsq measures a solver's consistency with its own vacuum field, so convergence to an incorrect boundary can still yield a small value. The normal-field residual and axis position measured against BIEST separate the correct boundary from the self-consistent one.

[unalmis]

FYI, I had written a differentiable, BIEST style gmres solver in links 1 and 2 that is tangent by construction and found it converged well in my testing.