Equilibrium - BUGFIX! - Build the root-area weight with the inverse Fourier transform

[logan-nc]

⚠️ NEEDS CAREFUL REVIEW: this changes the Fourier convention used inside one routine

Nothing in Utilities/FourierTransforms.jl changes. The forward transform ft(data) = basis·data/N with basis = exp(−i(mθ − nν)) and the inverse inverse(ft, modes) = adjoint(basis)·modes are untouched, and their round-trip test still passes (20/20). What changes is that compute_sqrtamat stops rolling its own backward step and calls that inverse.

History of the routine, so the reviewer can judge whether it was ever right:

commit	date	backward step in compute_sqrtamat	library inverse at that time
02b0f0b7 (first version)	2026-06-11	(1/N) Σ_m c_m exp(−imθ_j) via cslth/snlth	exp(+imθ) (forward was exp(−imθ)/N, Fortran iscdftf)
ece3e442 (complex overhaul)	2026-06-30	(basis · c)/N = the same exp(−imθ)/N	conj(basis)·modes = exp(+imθ)
1c955043 (basis transposed)	2026-06-30	(transpose(basis) · c)/N, mechanically adjusted	adjoint(basis)·modes = exp(+imθ)
this PR		inverse(ft, c)	unchanged

At every step the routine's backward kernel was the forward kernel exp(−imθ) with a 1/N, while the library's inverse was exp(+imθ) without one. The two refactors carried the hand-rolled line across faithfully; the mismatch dates from the first version. The routine was never tested directly: the develop runtests_coordinate_invariant.jl uses "a generic invertible complex operator standing in for the √weight", and the scripts/test_power_norm_invariance.jl the old docstring cited is not in the tree.

Why the early benchmarks could not see it. Reflecting m → −m is unitary and the 1/N is a scalar, so every norm or Parseval-type identity holds for the old operator up to a constant. The difference only shows in the structure. On the DIII-D example surface (ψ_N 0.99, mtheta 256, 35 modes; script in the review package):

check	old operator	this PR
S_old[:, m] == S_new[:, −m] / N over the 25 mirrored modes	max relative residual 0.0
diagonal	max 0.027	mean 6.8129 = ⟨√(J∣∇ψ∣)⟩ = 6.8129
Hermitian residual	0.34	1.5e-16
‖S b‖² / ∫∣b∣² w² dθ	1.1e-5 (≈ 1/N²)	0.996 (truncated basis; 1.0 on the complete basis in the test)
vs Fortran J_surf_2 on the benchmark surface	zero diagonal, norm ~300× small	correlation 1.0000, diagonal within 0.05 %

The first row is the whole bug in one line: the old column for mode m is the correct column for mode −m divided by N, exactly. The Fortran operator is the independent reference for which of the two is right.

---

compute_sqrtamat builds the √(J|∇ψ|) convolution operator Σ that the coordinate-invariant field space rests on: the b̃→b̄ operator S = Σ/√A, the flux conform R = S·A, field_space_response_matrices, the root-area-weighted forcing and response fields written to PerturbedEquilibrium/, and the conformed singular-coupling matrices. It took each unit mode to θ-space with transpose(basis)/N instead of the inverse transform adjoint(basis), so the matrix came out with its rows reflected in m and scaled by 1/mtheta: a zero diagonal and a norm about 300× too small, where Fortran GPEC's J_surf_2 is a near-identity weight (diagonal ≈ 0.92).

Found while benchmarking the error-field assessment of the original OMFIT project against its Fortran GPEC run, where the dominant-mode singular values were 190× too small and every per-coil overlap was off by factors of 10–900. With the fix the operator matches the Fortran J_surf_2 on that surface to the grid difference (correlation 1.0000, diagonal within 0.05 %) and the per-coil overlaps land on the Fortran values (median ratio 0.95; the remaining scatter is a separate helicity bug fixed in its own PR).

What changes

• src/Equilibrium/CoordinateInvariant.jl: the inverse step of compute_sqrtamat uses adjoint(ft.basis) (the documented inverse transform) with no 1/N; the docstring now states the correct identity and the Fortran correspondence.
• test/runtests_coordinate_invariant.jl: a testset on the Solovev surface checks the diagonal equals ⟨√(J|∇ψ|)⟩, Hermiticity, each column as the transform of the weighted unit mode, the weighted Parseval identity on the complete basis, and rootarea_to_area_weight · area_to_rootarea_weight = I. The old code fails all but Hermiticity.
• regression-harness/cases/diiid_n1.toml: tracks ‖forcing b̃‖ (PerturbedEquilibrium/forcing_b_root_area), since no tracked quantity on develop depended on the operator.

What moves: every field-space output (ResponseMatrices/ field-space matrices, forcing_b_root_area, forcing_b_area, response_b_*, SingularCoupling/C_* as stored) and anything read from them. Flux-space quantities and the harness's existing 47 quantities are unchanged.

Review package (the operator before/after on the DIII-D example surface and the reflection check): linked in the comments as a claude.ai artifact.

cc @matt-pharr

Release note

• Audience: users
• Numerical impact: field-space outputs change (harness @ 76808df (message-only rewrite of 244b2e6): 47 unchanged, 1 changed — the new ‖forcing b̃‖ quantity, 30.6 → 4.683e-4; identical to the run at 26db56c)
• Migration: none

The root-area weight operator (Σ/√A, Fortran's J_surf_2) was built with the wrong inverse transform, reflecting it in m and scaling it by 1/mtheta; all coordinate-invariant field-space outputs and the conformed singular-coupling matrices were wrong. Fixed and tested against the weighted Parseval identity and against Fortran GPEC.

Regression report

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/equilibrium-sqrtamat-inverse-fourier  @ 26db56c8 (2026-09-12)
       env: julia 1.11.7, x86_64-linux-gnu, manifest 06e27666 (pinned), 16 threads/16 BLAS
----------------------------------------------------------------------------------------------------------------------------------------------
Quantity                                      develop          bugfix/equilibrium-sqrtamat-inverse-fourier  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           
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           
PE toroidal torque                            5.087465e-02     5.087465e-02                                 0.0e+00              OK           
NTV torque FGAR [N·m]                         5.296762e-01     5.296762e-01                                 0.0e+00              OK           
NTV kinetic energy dW FGAR [J]                7.924971e-02     7.924971e-02                                 0.0e+00              OK           
||forcing b~|| (root-area-weighted)           3.062518e+01     4.683328e-04                                 3.062e+01 (100.00%)  ** CHANGED **
resonant area-weighted field b^r              [5 elem]         [5 elem]                                     0.0e+00              OK           
==============================================================================================================================================
Summary: 1 changed, 47 unchanged

🤖 Generated with Claude Code

https://claude.ai/code/session_01BzKDmtKYDotJWrcxWf1oJu

[logan-nc]

Commit message of the fix reworded (history rewritten, content identical: 244b2e6 → 76808df) to remove case-specific values; the review package artifact now shows the operator on the DIII-D example surface instead. The benchmark case is proprietary; only agreement percentages are quoted.

[ebursch]

Definitely a good add, suprised it wasn't in there already.

[ebursch]

The implemented changes look good to me, I think using the built in fourier library consistently makes sense and the added test is great. I'm a bit worried about the remaining differences noted in my earlier comment, but my intuition is that they stem from a different part of the workflow. That said, maybe have claude look into it and double check it is unrelated to this issue and if not then I think this piece is good to go.

[logan-nc]

Still working through everything, but one initial thought is that there remain some significant differences between the GPEC and the fixed Julia results, considering this is a log-scale plot. Are two of the arrays really an order of magnitude different between versions?

I see a factor of ~2 off for one pair of arrays and perfect agreement for another in that plot. The legend says there is still a helicity bug. I agree we should confirm that the agreement gets better if the same plot is remade with that bug fixed. Maybe this is a case where trying to make independent PRs for simplicity was the wrong call and we should stack - can you confirm the plot improved when #457 is merged into here? If so, I am fine merging as a stack as long as you are @ebursch

(Edited to remove device-specific coil names; the per-coil comparison will be redone on the DIII-D-like example against Fortran GPEC so it can be shown here.)

[ebursch]

Do you think it's a factor of 2 even with the log scale @logan-nc ?

[logan-nc]

For the historical record, @ebursch was right → we are investigating in more detail now

[logan-nc]

Review package for this PR (republished at a new link; DIII-D-like and synthetic data only): https://claude.ai/code/artifact/3697420a-2578-4693-a150-c1890d39a657

[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).

[logan-nc]

On the earlier question about the per-coil plot: that figure compared against a run made before #457 was applied, and it used data that cannot be shown here, so it has been withdrawn. The comparison has been redone where it can be shown: the DIII-D-like error-field example against Fortran GPEC, both codes on matched DCON and coil-grid settings, with this fix, #457 and #460 applied. For the C-coil as built and for single F-coil hoops shifted 1 mm (x and y) or tilted 0.1°, the dominant-mode overlap from Julia's sensitivity table matches Fortran's finite difference of two runs to 0.01–6 % (the 6 % on the weakest hoop), and the dominant singular value to 0.04 % (core window ψ_N ≤ 0.9). Figure in the #448 review package: https://claude.ai/code/artifact/7eb87dcf-6a2d-4532-8bca-f58ceb76968e. On stacking: the three bugfixes are independent of each other and of the stack, but the stack's numbers need all three, so reviewing #458, #457 and #460 first and merging them ahead of the stack is the intent.

[ebursch]

It looks like these new artifacts are not public @logan-nc

[ebursch]

@logan-nc The develop branch plot is now all zeros? It used to be diagonal w/ slope = -1. There haven't been any commits since the first version, so it seems a bit suspicious.

[logan-nc]

@ebursch its just a colorbar normalization issue - the remake normalized the two panels to the same scale and they are off by ~1/M - claude is remaking it now but I believe it still looks the same as before structure wise

[ebursch]

Okay I think this one is ready to go now. @logan-nc