ForceFreeStates - BUGFIX! - Use the fixed start for integrations that begin away from the axis

[logan-nc]

Found while benchmarking the original OMFIT project's error-field case against its Fortran GPEC run. That case starts the integration well inside the plasma rather than at the axis, and with both #457 and #458 applied its innermost resonant-coupling row was still 1.34× Fortran's and the total energies 3 % apart, insensitive to every resolution and tolerance knob in either code.

Cause. initialize_el_at_axis! applies the Frobenius axis start [Glasser 2016 Eq. 51] wherever the integration begins. That start is a power series about the magnetic axis; at an interior surface it is neither the axis-regular solution nor Fortran DCON's start, which is always U₁ = 0, U₂ = I (zero displacement at the first surface). The DIII-D ideal example, which matches Fortran to four digits at psilow = 1e-4, shows the drift directly:

psilow	innermost coupling row, Julia/Fortran norm	lowest W_t, Julia vs Fortran
1e-4	1.000	1e-4
0.01	1.010	0.03 %
0.03	1.042	0.1 %
0.10	1.243	2.3 %
0.30	2.140	7 %
0.30 with the fixed start (frobenius_psi_max = 0)	0.995	0.4 %

On that case the fixed start takes the coupling rows to 1.0000 correlation with Fortran and every per-coil overlap to within 1 % of OMFIT.

Change. One control, frobenius_psi_max (default 0.01), bounds where the Frobenius start is used. Above it the fixed start is used and a warning says so; 0 selects the fixed start everywhere, silently. It replaces the undocumented fixed_axis flag (removed; no example, test or case file set it). Every example in the repo starts at or below 0.01, so nothing tracked moves; the ! is for the interface and for interior-start runs, whose results change.

Verification. test/runtests_eulerlagrange.jl gains a Solovev testset: axis start keeps a nonzero Frobenius U₁ with U₂ = I; a start at ψ_N = 0.3 switches to U₁ = 0, U₂ = I and logs the warning; raising frobenius_psi_max keeps the Frobenius start; frobenius_psi_max = 0 selects the fixed start silently even at the axis. ODE Tests 99/99.

Not touched: which start is physically right for an interior surface. The fixed start is what DCON does and what the published numbers used; the switch reproduces that and stops an axis expansion from being applied silently far from the axis.

cc @matt-pharr

Release note

• Audience: users
• Numerical impact: interior-start runs (psilow or qlow placing the start above ψ_N = 0.01) now begin from the fixed start and change; axis-started runs are unchanged (harness @ e93df8f: diiid_n1 develop vs branch: 47 unchanged; the 5 N/A rows are error-field quantities the main checkout's case TOML tracks that neither ref writes)
• Migration: fixed_axis = true becomes frobenius_psi_max = 0; set frobenius_psi_max above the start to keep the Frobenius start for an interior start

ForceFreeStates uses the Frobenius axis start only when the integration begins at or below the new frobenius_psi_max (default 0.01); an interior start uses Fortran DCON's fixed start (U₁ = 0, U₂ = I) with a warning, and frobenius_psi_max = 0 selects it everywhere. The fixed_axis flag is removed.

Regression report

================================================================
Case: diiid_n1 — DIII-D-like equilibrium, n=1, ideal + perturbed equilibrium
================================================================
Regression Report: diiid_n1
======================================================================================================================================
Ref 1: develop  @ 349a0c26 (2026-09-04)
       env: julia 1.11.7, x86_64-linux-gnu, manifest 06e27666 (pinned), 16 threads/16 BLAS
Ref 2: bugfix/forcefreestates-interior-start-initialization  @ e93df8f0 (2026-09-12)
       env: julia 1.11.7, x86_64-linux-gnu, manifest 06e27666 (pinned), 16 threads/16 BLAS
--------------------------------------------------------------------------------------------------------------------------------------
Quantity                                      develop          bugfix/forcefreestates-interior-start-initialization  Diff       Status
--------------------------------------------------------------------------------------------------------------------------------------
total energy Re(et[1])                        8.012318e-01     8.012318e-01                                          0.0e+00    OK    
total energy Im(et[1])                        4.142529e-05     4.142529e-05                                          0.0e+00    OK    
plasma energy Re(ep[1])                       -1.348486e+00    -1.348486e+00                                         0.0e+00    OK    
vacuum energy Re(ev[1])                       2.149718e+00     2.149718e+00                                          0.0e+00    OK    
vacuum matrix min eigenvalue                  1.873976e-01     1.873976e-01                                          0.0e+00    OK    
plasma energy (all)                           [35 elem]        [35 elem]                                             0.0e+00    OK    
vacuum energy (all)                           [35 elem]        [35 elem]                                             0.0e+00    OK    
total energy (all)                            [35 elem]        [35 elem]                                             0.0e+00    OK    
ODE steps (saved)                             2576             2576                                                  0.0e+00    OK    
ODE steps (total)                             4572             4572                                                  0.0e+00    OK    
q0                                            1.204212e+00     1.204212e+00                                          0.0e+00    OK    
q95                                           4.781723e+00     4.781723e+00                                          0.0e+00    OK    
beta_t                                        1.327024e-02     1.327024e-02                                          0.0e+00    OK    
beta_n                                        1.372511e+00     1.372511e+00                                          0.0e+00    OK    
internal inductance li1                       8.842392e-01     8.842392e-01                                          0.0e+00    OK    
internal inductance li2                       7.080847e-01     7.080847e-01                                          0.0e+00    OK    
internal inductance li3                       7.304433e-01     7.304433e-01                                          0.0e+00    OK    
poloidal beta betap1                          6.680744e-01     6.680744e-01                                          0.0e+00    OK    
poloidal beta betap2                          5.349834e-01     5.349834e-01                                          0.0e+00    OK    
poloidal beta betap3                          5.518761e-01     5.518761e-01                                          0.0e+00    OK    
# singular surfaces                           5                5                                                     0.0e+00    OK    
singular psi locations                        [5 elem]         [5 elem]                                              0.0e+00    OK    
singular q values                             [5 elem]         [5 elem]                                              0.0e+00    OK    
current beta betaj                            4.236478e-01     4.236478e-01                                          0.0e+00    OK    
plasma volume                                 1.829472e+01     1.829472e+01                                          0.0e+00    OK    
plasma current                                1.152130e+00     1.152130e+00                                          0.0e+00    OK    
mpert                                         35               35                                                    0.0e+00    OK    
npert                                         1                1                                                     0.0e+00    OK    
toroidal field bt0                            2.006573e+00     2.006573e+00                                          0.0e+00    OK    
wall field bwall                              3.880145e-01     3.880145e-01                                          0.0e+00    OK    
aspect ratio                                  2.845746e+00     2.845746e+00                                          0.0e+00    OK    
elongation kappa                              1.708350e+00     1.708350e+00                                          0.0e+00    OK    
q profile (checksum)                          0cd285cea88d...  0cd285cea88d...                                       identical  OK    
pressure profile (checksum)                   a1c48b266622...  a1c48b266622...                                       identical  OK    
Mercier D_I profile (checksum)                eeb06744e795...  eeb06744e795...                                       identical  OK    
resistive interchange D_R profile (checksum)  fa37296851f6...  fa37296851f6...                                       identical  OK    
ballooning Delta' profile (checksum)          bb713caf14eb...  bb713caf14eb...                                       identical  OK    
island half-widths                            [5 elem]         [5 elem]                                              0.0e+00    OK    
Chirikov parameter                            [5 elem]         [5 elem]                                              0.0e+00    OK    
||resonant area-weighted field||              5.189179e-04     5.189179e-04                                          0.0e+00    OK    
dominant-coupling singular values             N/A              N/A                                                   N/A        N/A   
|forcing overlap with dominant mode|          N/A              N/A                                                   N/A        N/A   
|delta_nominal| of coil set 1                 N/A              N/A                                                   N/A        N/A   
||ddelta/d(shift)|| over coil sets            N/A              N/A                                                   N/A        N/A   
||ddelta/d(tilt)|| over coil sets             N/A              N/A                                                   N/A        N/A   
PE plasma energy                              3.422677e+00     3.422677e+00                                          0.0e+00    OK    
PE vacuum energy                              3.174510e+00     3.174510e+00                                          0.0e+00    OK    
PE surface energy                             5.826684e+00     5.826684e+00                                          0.0e+00    OK    

Notes for reviewers

• Independent of the error-field stack and of ForcingTerms - BUGFIX! - Take the coil-grid helicity from the equilibrium file's current sign #457/Equilibrium - BUGFIX! - Build the root-area weight with the inverse Fourier transform #458; branched from develop.
• JuliaFormatter normalized EulerLagrange.jl and Riccati/Driver.jl on first touch; the substantive diff is the block in initialize_el_at_axis!, the two docstrings and the control fields.

🤖 Generated with Claude Code

https://claude.ai/code/session_01BzKDmtKYDotJWrcxWf1oJu

[logan-nc]

Amended per review direction: fixed_axis is removed and frobenius_psi_max = 0 selects the fixed start everywhere, so there is one control. Harness rerun at e93df8f: 47 unchanged. Test file 99/99.

[logan-nc]

Stack summary — OMFIT → Julia error-field tolerance migration (updated 2026-09-17)

Nothing in this set merges without a third-party human review, the three small bugfix PRs included.

Full page with the diagram, per-PR table, review order and links to every review package (a private claude.ai artifact until the author shares it; DIII-D-like and synthetic data only): https://claude.ai/code/artifact/408ddd36-e52c-42d8-8f5c-52f680f41da9

Group	PRs	Base
Bugfixes the stack's numbers depend on — review first	#458 root-area weight (careful review: Fourier convention inside one routine), #457 helicity from the current sign, #460 interior-start initialization	develop
Prerequisites	#446 resonant-coupling SVD API, #447 per-coil-set forcing modes	develop
ErrorFields stack, each based on the previous branch	#448 sensitivities → #449 tolerance TOML → #450 sampling → #451 Monte Carlo → #452 risk + scan → #453 plots + phasing → #455 NTV limits → #462 NTV torque against rotation	#448 on develop (carries merges of #446/#447); rebase the chain after the bugfixes and #446/#447 land

Benchmark, qualitatively. The original OMFIT project's case, rebuilt from that project's own run inputs, is reproduced once #458, #457 and #460 are in: dominant-mode singular values agree to 1e-4, singular-coupling rows match Fortran GPEC at 1.0000 correlation and within 0.5 % in norm, and every coil set's dominant-mode overlap agrees within 1 % once both codes use the same converged toroidal coil grid. The DIII-D-like example agrees with Fortran to 1e-4 throughout. Any comparison shown on these PRs uses the DIII-D-like examples only.

Design issue for what comes after the stack: #461 (momentum-balance module fed by tabulated torque surfaces; stationary solvers only).

[d-burg]

Independent check while triaging #440 (which this supersedes — I'm closing it): the frobenius_psi_max = 0.01 default leaves a case on the Frobenius path that needs the fixed start.

#440's fixture is a synthetic diverted DIII-D-like ramp-up (q0 = 3.41, n = 1, no wall, psilow = 0.01, ldp grid mpsi 128, riccati). 0.01 > 0.01 is false, so it keeps the Frobenius start. On this branch at e93df8f:

setting	et[1]	max eig W_p	negative W_p eigs
default (frobenius_psi_max = 0.01)	−1.29e5	4.70e3	1
frobenius_psi_max = 0	+1.504	1.6090e4	0

Fortran DCON on the same deck (same mlow/mhigh/psilim/dmlim): stiff W_p +1.6092e4, no negative eigenvalue. So your fixed start is right and reproduces DCON; the default reports a spurious −1e5 instability.

I understand why the boundary is > at 0.01 — it is what keeps the LAR and DIII-D IMAS examples, which sit at exactly 0.01, unchanged. But your own table already shows 1 % coupling error at ψ_N = 0.01, and on #440's clean neighbouring slice the Frobenius stiff eigenvalue was +2.76e5 against DCON's +1.61e4 (17×) — so the eigenvalue looks wrong everywhere and the threshold mitigates rather than repairs.

Not asking you to widen this PR. Three ways to close it, your call:

1. lower the default;
2. make the comparison >=;
3. keep 0.01 and document that a start at exactly the threshold is on the axis expansion.

If it is the third, say so in the frobenius_psi_max docstring and I will carry the fixture and a @test_broken in a follow-up after this merges.

Filed #464 for the underlying compute_axis_init repair, which neither PR attempts.

Two caveats on my numbers: the fixture is on grid_type = "ldp", not the current auto default, so they want re-confirming on auto; and #440's own regression table is from 2026-09-02 and stale against current develop — the table above was measured today on your head.

[logan-nc]

@d-burg I'm lost in your AI verbosity.

1. are you saying the expansion now causes a spurious instability in the standard DIII-D example? if so, we need to track down when that happened cuz it all benchmarked fine when the feature was added.

2. are you just asking for the default to be lower? I'm fine with that if it still triggers in at least one example and any appropriate tests.

As a side note, is think physically that allowing finite displacement at psilow would be all the more important at higher psilow... but yes the expansion we use is for near zero. correct initialization at higher values will take more thought by future us.

[d-burg]

@logan-nc sorry, I insufficiently edited that one – this PR is still not fully within my understanding and I've tried to learn too many new PRs today...

1. No – not the DIII-D example. That was ForceFreeStates - BUGFIX! - Default to the DCON fixed-axis initial condition #440's synthetic diverted ramp-up (q_0 = 3.41, psilow = 0.01), a case that didn't exist when the Frobenius start was added, and that's the branch we agreed to close since you already covered the needed fix here. The DIII-D example at 1e-4 is unchanged. Claude-confirmed on the auto grid too: "default gives et[1] = -8.7e4 with a negative W_p eigenvalue, frobenius_psi_max = 0 gives +1.5, stiff W_p 1.6026e4 vs DCON's 1.6092e4." I think it all looks good to me

2. Yes. >= (or a default of 0.005) is enough. The Frobenius start still runs on every 1e-4 deck and a10_kinetic at 1e-3. Five of the cases start at exactly 0.01 (IMAS + the four LAR) and would move to the fixed start with the warning, however none are a harness case and no test pins their numbers so I'm assuming that's fine. I agree the expansion is for near-zero, am tracking that in ForceFreeStates - Repair or retire the Frobenius axis start #464. If you and/or your claude finds that shift beneficial, maybe ping me and I'll double check I'm understanding the results

[logan-nc]

gotcha. just make it >=0.005 then! that's good with me.

other than that, is this one good to merge @d-burg ? (or did you skim the diff as a human as well as sicking claude on it?)

[ebursch]

This looks ready to merge.