DCON, RDCON, STRIDE - BUGFIX - Trace-aware Frobenius exponents at resonant surfaces

[StuartBenjamin]

Fortran port of Julia GPEC #517.

Near a rational surface the solutions go as |ψ − ψ_s|^p, where p are the eigenvalues of the 2×2 resonant block of M0: p = tr/2 ∓ √(tr²/4 − det). DCON, RDCON and STRIDE use ∓√(−det), which assumes tr = 0. It isn't: on the DIII-D example tr is imaginary, 1e-7 to 8e-4 (largest at q = 5 and 6). Δ′ is read off a coefficient |dψ|^(2α) below the big solution, so this small error becomes a large spurious imaginary part on the Δ′ diagonal: STRIDE gives Δ′(q=5) = −2412 − 47180i.

Changes

• sing_vmat/sing1_vmat store the powers tr/2 ∓ √(tr²/4 − det) in singp%power. The zeroth-order eigenvectors, sing_get_ua, sing_get_dua, sing_get_ua_cut and the sing1_* versions use them.
• RDCON's sing_order_ceiling uses 2 Re √(tr²/4 − det), as ForceFreeStates - BUGFIX! - Use trace-aware Frobenius exponents at resonant surfaces OpenFUSIONToolkit/GPEC#517 does.
• The left side reuses the right-side powers: flipping sig gives M0_L = −D M0_R D with D = diag(1, −1), so the left powers are −p.
• RDCON appends the powers to each surface's header record in vmat.bin, and RMATCH reads them. Old RMATCH reads new files. New RMATCH on an old file rereads the header without the powers (by rewind and skip, since the file position after a failed read is processor dependent) and uses ±α, giving bitwise the old ua.
• singp%alpha = SQRT(-CMPLX(di)) was single precision (no KIND); now KIND=r8 in DCON, RDCON and STRIDE. α no longer sets the powers, so this only changes the α written to vmat.bin.

Results

DIII-D examples (STRIDE: DIIID_ideal; RDCON: DIIID_resistive, Cholesky), Δ′ diagonal at q = 2…6:

	develop	this PR
STRIDE	7.394+0.098i, −4.619−0.28i, −17.15+5.5i, −2412−47180i, −113.6−438i	7.391, −4.618, −17.11, −2412.4−0.048i, −113.6 (Im ≤ 2e-5 relative)
RDCON, Cholesky	3.854+0.007i, −10.02−0.083i, −23.50−0.62i, −3346−13773i, −988−102i	3.853+0.008i, −10.02+0.006i, −23.50+0.11i, −3395−21.3i, −988.0+1.7i

• STRIDE: the diagonal imaginary parts drop to round-off; the real parts move ≤ 0.2 %.
• RDCON Galerkin: the imaginary part falls 650× at q = 5 but not to round-off, as for the Julia Galerkin solver in ForceFreeStates - BUGFIX! - Use trace-aware Frobenius exponents at resonant surfaces OpenFUSIONToolkit/GPEC#517. RDCON's LU Δ′ on this example is sensitive to round-off (rewriting x^p as 1/x^(−p) moves it by up to 12 % of the largest entry), so the table uses Cholesky.
• W eigenvalues (DCON, STRIDE, RDCON) move ≤ 4e-7 relative.
• The off-diagonal Δ′ entries are non-Hermitian in both builds, for another reason.

Not changed: the sing1 degenerate branch (sing1_flag) still assumes ±1/2, and RMATCH match.f keeps −0.5 ± √(−di). As with OpenFUSIONToolkit/GPEC#517, a physics check against Glasser 2016 would be welcome. The RDCON numbers above still include the issues fixed in #303 and #304.

Scripts: frobenius_trace_scripts.zip: scripts/run.sh RUNSET <case> <code> (edit its paths), then scripts/collect.sh TESTROOT > tests.txt and scripts/plot_imdp.py. vmat.bin compatibility: scripts/rmatch_read/run_read_tests.sh NEW_WT DEVELOP_WT TESTROOT.

---

🤖 Generated with Claude Code

[matt-pharr]

This checks out and I verified that it is the same as Julia #517 OpenFUSIONToolkit/GPEC#517.