ForceFreeStates/PerturbedEquilibrium - Automate the poloidal mode range with convergence diagnostics #518

Consolidates the discussion on #193 (draft, v2.0 milestone, to be closed in favour of this issue) so that future work starts from the insights there rather than from that branch. Post-GPEC 2.0.

Problem

The poloidal mode range of ForceFreeStates (FFS) is still the Fortran rule

mlow  = min(n·qmin, 0) − 4 − delta_mlow
mhigh = n·qmax + 2·delta_mhigh        # delta_mhigh is doubled "for consistency with Fortran"

(src/GeneralizedPerturbedEquilibrium.jl), mthvac is a fixed 480 (ForceFreeStates/CoreTypes.jl), and the only convergence test is the manual one from the DCON readme: re-run until the mlow/mhigh components of the eigenfunction fall below ~1e-3, else raise delta_m*. Nothing in gpec.h5 reports whether a run was converged in m.

The most common failure mode in practice is on the PerturbedEquilibrium (PE) side: the physics sets a compact FFS range, but boxy coil sets have a broad poloidal spectrum. ForcingTerms.compute_coil_forcing_modes! projects the coil field onto ffs.mlow:ffs.mhigh and convert_forcing_normalization! skips any mode outside that range silently (ForcingTerms/CoilFourier.jl, 1 <= idx <= mpert || continue). Users discover this as aliasing or missing forcing and re-run with larger delta_m*, paying for a bigger FFS solve they do not need.

Where #193 got to

• Replaced the deltas with _estimate_m_bandwidth: Fourier spectrum of r² at psilim, highest harmonic above m_accuracy × DC. r² rather than the Jacobian because in Hamada coordinates the Jacobian is constant.
• Replaced manual mthvac with _estimate_arc_nonuniformity: mean(ds)/min(ds) on the SFL boundary, mthvac = round_up(max(4·mpert, max(96, 6·mhigh)·arc_factor), 32).
• mlow/mhigh as explicit overrides, m_accuracy as the single knob.
• Validated on Solov'ev only (bandwidth 7, arc_factor 1.64, mthvac 160); DIII-D gave bandwidth 4, arc_factor 4.01 but was blocked on an equilibrium-grid issue in the parent branch. Never tested on a shaping or q95 scan. Consensus in the thread: the estimate needs real investigative effort and is not 2.0 work.

What is already in the code and the literature

A Fourier-convergence test at every rational surface exists on paper and both halves exist in the code. Glasser, Phys. Plasmas 23, 072505 (2016), §III: the eigenvalues of the Frobenius matrix M₀ at a singular surface converge to ±√(−D_I) as the number of retained harmonics grows, and "agreement between these two ways of computing D_I can thus serve as a test for adequate Fourier convergence at a resonant surface" (doi:10.1063/1.4958328). In the repo:

• α from the resonant 2×2 block of m0mat in ForceFreeStates/Surfaces/Asymptotics.jl (currently via sqrt(-det), with a comment that this equals the eigenvalue only when the block is traceless);
• the quadrature D_I from flux-surface averages in LocalStability (behind local_stability_flag) and the GGJ E + F + H − 1/4 in Surfaces/ResistEval.jl / InnerLayer.GGJ.mercier_di.

They are never compared. Jake has been testing this independently in 3D DCON and finds it behaves better on resonance-banded spectra (Fortran-style delta_mband, which does not exist in the Julia code) than on coil-driven "what spectrum do my coils produce" ranges:

(J. Halpern, #193 comment, 2026-10-09.)

Note that in 3D the quadrature comparator does not exist: Glasser, Phys. Plasmas 30, 052502 (2023) (doi:10.1063/5.0144211) states the generalized Mercier criterion can only be evaluated on the rational surface with a given set of resonant harmonics, so there M-convergence is the only check. For tokamaks we have both.

Precedent for a spectral-tail rule: PIES adds harmonics when the amplitude of the highest retained poloidal mode exceeds the precision target (Reiman & Greenside, J. Comput. Phys. 87, 349 (1990), doi:10.1016/0021-9991(90)90257-2). A title-level sweep by code name (MARS-F/K, CASTOR, AEGIS, ELITE, PEST, KINX, NOVA, CAS3D, TERPSICHORE, GATO, SPEC) found no published automated m-range selection in a linear stability code; everyone scans by hand. Not exhaustive.

Proposed directions

1. Diagnostics first, no behaviour change

Always record in gpec.h5, on every run:

• per singular surface: alpha from m0mat and the quadrature D_I, and the mismatch |α² + D_I| (requires running the Mercier quadrature at the singular surfaces regardless of local_stability_flag; it is cheap);
• the FFS edge spectral-tail ratio, max(|ξ_mlow|, |ξ_mhigh|) / max_m |ξ_m| over the least-stable / each requested eigenmode, i.e. the DCON readme rule as a number;
• on the PE side, the fraction of applied-field power outside the FFS range, Σ_{m∉[mlow,mhigh]} |Φ_x,m|² / Σ_m |Φ_x,m|², with a @warn above a threshold instead of the silent continue.

These three numbers are what any auto-logic would have to be calibrated against, and they are useful on their own for database runs.

2. Auto-range for FFS

Keep the q-range core of the Fortran rule and replace the deltas with a tolerance on the diagnostics in (1): extend mlow/mhigh (and re-run FFS, or integrate once with a generous range and report where the tail criterion is met) until the tail ratio and the Mercier mismatch are below m_accuracy. The #193 r²-bandwidth estimate is a reasonable initial guess for the extension, not the criterion. mthvac from the arc-nonuniformity factor can stay as in #193 but needs the shaping/q95 scans it never got. mlow/mhigh remain expert overrides (needed to match a Fortran run exactly).

3. Decouple the PE basis from the FFS basis

This is the part that would remove the number-one re-run cause, and it needs the most thought. The idea: FFS integrates on a physics-motivated range M_F; PE is assembled on a wider forcing-driven range M_P ⊇ M_F; harmonics in M_P \ M_F carry vacuum response only ("behave like vacuum, not dropped"). Some directions:

• Compute the forcing spectrum before FFS setup. ForcingTerms only needs the equilibrium and psilim, so M_P can be known (coil spectrum + a Nyquist margin) when mlow/mhigh are chosen. For file-based forcing M_P is just the span of the supplied modes (today the file path already uses extrema(m_vals) rather than the FFS range).
• Compute the vacuum matrices once on M_P and slice. Green's-function matrix elements do not depend on truncation, so wv and the surface inductance L on M_F are sub-blocks of the M_P matrices. This replaces the two separate VacuumInput(... mlow:mhigh ...) calls in ForceFreeStates/Free.jl and PerturbedEquilibrium/Response.jl with one call plus a view, and is a saving even without the split.
• State the closure for M_P \ M_F precisely, in the inductance language of ResponseMatrices.jl. Λ = inv(2·t₁·(wp + wv)·t₂) is only defined on M_F. The natural "vacuum-like" closure is the rigid-boundary one: ξ_m = 0 for m ∉ M_F. The external harmonics in M_P \ M_F are then neither screened nor amplified by a displacement of their own, but they still couple: external flux in M_P \ M_F drives ξ on M_F through the surface inductance L (which couples harmonics), and ξ on M_F sources flux in all of M_P. On the extended basis the permeability is P = [P_FF P_FX; P_XF I] with both off-diagonal blocks nonzero and only the X–X block trivial. This is different from the energy closure "wp = 0 outside M_F", which would make the plasma infinitely compliant in those modes and is wrong. Write both down in the Λ, L, P notation of Park et al. 2007 and pick one; the rigid-boundary form is the candidate.
• Make it a convergence diagnostic rather than a modelling choice. Running with M_F = M_P (today's behaviour) and with the split, the two must agree as M_F grows. The size of P_FX and the plasma-response fraction at the FFS edge harmonics, |Φ_tot,m − Φ_x,m| / |Φ_x,m| at m = mlow, mhigh, tell whether the vacuum closure beyond M_F is justified for a given case — this is the PE analogue of the readme tail rule.
• Anti-aliasing on the θ grids. Field-space quantities (b·n on the control surface, R-Z reconstruction, resonant fields) are formed on θ grids with products of ξ and equilibrium geometry that spread harmonics; FieldReconstruction.jl already guards with mtheta_rz = max(2·(|mlow| + mpert), 512). With a split basis that rule must be keyed on M_P, and the coil projection (fourier_decompose_bn on mtheta_coil = 480) must satisfy Nyquist for the widest harmonic in M_P, or the "extra" vacuum harmonics are themselves aliased. This is the subtlety flagged in the thread.
• Singular surfaces. Resonant surfaces whose m lies outside M_F are dropped in Surfaces/Finding.jl and SingularCoupling.jl skips m_res ∉ [mlow, mhigh]. A wider M_P adds no resonant surfaces by construction (there is no plasma response there), which is consistent, but the resonant-coupling outputs should say which surfaces were excluded by M_F.

4. Validation

Shaping scan (Solov'ev), q95 and aspect-ratio scans (TJ-based benchmark scripts), DIII-D n = 1–3 with the real coil sets, against the Fortran range and the regression harness; report runtime as well as convergence. AIs over-optimise local conditions, so the acceptance test is a scan, not an example.

References

• Glasser, The direct criterion of Newcomb for the ideal MHD stability of an axisymmetric toroidal plasma, Phys. Plasmas 23, 072505 (2016) — §III, Eq. (46)–(49) and the convergence remark.
• Glasser, Generalized Mercier stability criterion for stellarators, Phys. Plasmas 30, 052502 (2023).
• Reiman & Greenside, Computation of zero β three-dimensional equilibria with magnetic islands, J. Comput. Phys. 87, 349 (1990).
• DCON namelist readme (delta_mlow/delta_mhigh/delta_mband, convergence rule).
• ForceFreeStates: physics-driven auto-computation of mlow/mhigh and mthvac #193 and its thread (Pharr, Halpern, Logan).

Also tracked as a stub in the GPEC Feature Add Ideas doc.

Generated by Claude Code