[StuartBenjamin]

Release note

• Audience: users
• Numerical impact: none. diiid_n1 and solovev_n1 are identical to develop (harness @ 15a4b7b)
• Migration: none

GPEC can now read inverse i-files (eq_type = "ifile"), such as OpenFUSIONToolkit TokaMaker's save_ifile output. This is a port of the Fortran read_eq_ldp_i.

Regression report

regress --cases diiid_n1,solovev_n1 --refs 84b7e4800,15a4b7b97 (develop vs this head; julia 1.12.6, pinned manifest, 8 threads): diiid_n1 53 unchanged, solovev_n1 22 unchanged. No harness case reads an i-file; the unit tests cover the reader.

Full report

================================================================
Case: diiid_n1 — DIII-D-like equilibrium, n=1, ideal + perturbed equilibrium
================================================================

Regression Report: diiid_n1
=================================================================================================
Ref 1: 84b7e4800  @ 84b7e4800 (2026-10-20)
       env: julia 1.12.6, x86_64-linux-gnu, manifest 4bdd795f (pinned), 8 threads/8 BLAS, jit build
Ref 2: 15a4b7b97  @ 15a4b7b97 (2026-10-20)
       env: julia 1.12.6, x86_64-linux-gnu, manifest 4bdd795f (pinned), 8 threads/8 BLAS, jit build
-------------------------------------------------------------------------------------------------
Quantity                                      84b7e4800        15a4b7b97        Diff       Status
-------------------------------------------------------------------------------------------------
total energy Re(et[1])                        8.006123e-01     8.006123e-01     0.0e+00    OK    
total energy Im(et[1])                        7.924925e-05     7.924925e-05     0.0e+00    OK    
plasma energy Re(ep[1])                       -1.372261e+00    -1.372261e+00    0.0e+00    OK    
vacuum energy Re(ev[1])                       2.172874e+00     2.172874e+00     0.0e+00    OK    
vacuum matrix min eigenvalue                  1.873975e-01     1.873975e-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)                             800              800              0.0e+00    OK    
ODE steps (total)                             795              795              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.842269e-01     8.842269e-01     0.0e+00    OK    
internal inductance li2                       7.080757e-01     7.080757e-01     0.0e+00    OK    
internal inductance li3                       7.304340e-01     7.304340e-01     0.0e+00    OK    
poloidal beta betap1                          6.680739e-01     6.680739e-01     0.0e+00    OK    
poloidal beta betap2                          5.349836e-01     5.349836e-01     0.0e+00    OK    
poloidal beta betap3                          5.518764e-01     5.518764e-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.236480e-01     4.236480e-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.708322e+00     1.708322e+00     0.0e+00    OK    
q profile (checksum)                          ea980cf7e8cf...  ea980cf7e8cf...  identical  OK    
pressure profile (checksum)                   75be5949b7c1...  75be5949b7c1...  identical  OK    
Mercier D_I profile (checksum)                cafa40526b39...  cafa40526b39...  identical  OK    
resistive interchange D_R profile (checksum)  c70e46071c9d...  c70e46071c9d...  identical  OK    
ballooning Delta' profile (checksum)          029149d57e8c...  029149d57e8c...  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.288864e-04     5.288864e-04     0.0e+00    OK    
dominant-coupling singular values             [3 elem]         [3 elem]         0.0e+00    OK    
|forcing overlap with dominant mode|          1.411634e-04     1.411634e-04     0.0e+00    OK    
|delta_nominal| of coil set 1                 7.035049e-05     7.035049e-05     0.0e+00    OK    
||ddelta/d(shift)|| over coil sets            1.245519e-04     1.245519e-04     0.0e+00    OK    
||ddelta/d(tilt)|| over coil sets             4.176407e-07     4.176407e-07     0.0e+00    OK    
PE plasma energy                              3.320677e+00     3.320677e+00     0.0e+00    OK    
PE vacuum energy                              3.174510e+00     3.174510e+00     0.0e+00    OK    
PE surface energy                             5.509017e+00     5.509017e+00     0.0e+00    OK    
PE toroidal torque                            1.119454e-04     1.119454e-04     0.0e+00    OK    
NTV torque FGAR [N·m]                         5.483495e-01     5.483495e-01     0.0e+00    OK    
NTV kinetic energy dW FGAR [J]                8.163740e-02     8.163740e-02     0.0e+00    OK    
Runtime (s)                                   355.2s           340.1s                      --    
||forcing b~|| (root-area-weighted)           4.683330e-04     4.683330e-04     0.0e+00    OK    
resonant area-weighted field b^r              [5 elem]         [5 elem]         0.0e+00    OK    
=================================================================================================
Summary: 53 unchanged


================================================================
Case: solovev_n1 — Solovev analytical equilibrium, n=1, ideal stability
================================================================

Regression Report: solovev_n1
=================================================================================
Ref 1: 84b7e4800  @ 84b7e4800 (2026-10-20)
       env: julia 1.12.6, x86_64-linux-gnu, manifest 4bdd795f (pinned), 8 threads/8 BLAS, jit build
Ref 2: 15a4b7b97  @ 15a4b7b97 (2026-10-20)
       env: julia 1.12.6, x86_64-linux-gnu, manifest 4bdd795f (pinned), 8 threads/8 BLAS, jit build
---------------------------------------------------------------------------------
Quantity                      84b7e4800        15a4b7b97        Diff       Status
---------------------------------------------------------------------------------
total energy Re(et[1])        1.090560e-02     1.090560e-02     0.0e+00    OK    
total energy Im(et[1])        2.141300e-10     2.141300e-10     0.0e+00    OK    
plasma energy Re(ep[1])       -1.041640e+01    -1.041640e+01    0.0e+00    OK    
vacuum energy Re(ev[1])       1.042730e+01     1.042730e+01     0.0e+00    OK    
vacuum matrix min eigenvalue  2.170703e+00     2.170703e+00     0.0e+00    OK    
plasma energy (all)           [32 elem]        [32 elem]        0.0e+00    OK    
vacuum energy (all)           [32 elem]        [32 elem]        0.0e+00    OK    
total energy (all)            [32 elem]        [32 elem]        0.0e+00    OK    
ODE steps (saved)             287              287              0.0e+00    OK    
ODE steps (total)             285              285              0.0e+00    OK    
q0                            1.900006e+00     1.900006e+00     0.0e+00    OK    
q95                           3.147422e+00     3.147422e+00     0.0e+00    OK    
beta_t                        4.621027e-02     4.621027e-02     0.0e+00    OK    
beta_n                        3.217775e+00     3.217775e+00     0.0e+00    OK    
# singular surfaces           2                2                0.0e+00    OK    
singular psi locations        [2 elem]         [2 elem]         0.0e+00    OK    
singular q values             [2 elem]         [2 elem]         0.0e+00    OK    
mpert                         32               32               0.0e+00    OK    
npert                         1                1                0.0e+00    OK    
q profile (checksum)          af93c6269245...  af93c6269245...  identical  OK    
pressure profile (checksum)   e252b3b40ede...  e252b3b40ede...  identical  OK    
ca_left (checksum)            06afdeb86cdc...  06afdeb86cdc...  identical  OK    
Runtime (s)                   279.8s           267.6s                      --    
=================================================================================
Summary: 22 unchanged

Notes for reviewers

GPEC had no reader for the i-file, the inverse format TokaMaker's save_ifile writes. An i-file stores the flux surfaces themselves, so the equilibrium does not depend on a g-file's R,Z grid.

Changes

• read_ifile (ReadEquilibrium.jl) builds the same InverseRunInput and InverseIngest (for gpec.h5 reruns) as the CHEASE readers. Reals may be real8 or real4 (the Fortran reader takes only real*8). F and P are the tabulated values.
• Optional trailing records, told apart by length: FF′ and p′ are skipped; [ψ_axis, ψ_LCFS] (save_ifile(psi_bounds=True), ALL - IMPROVEMENT - Self-describing HDF5 metadata (issue #226, PR B) #364) normalizes ψ_N to the true LCFS.
• resolve_psihigh! clamps psihigh to the outermost file surface, for i-files and on the inverse rerun path.
• i-file equilibria take q from the field-line integral of their geometry, as Fortran inverse_run does; the tabulated q, splined in ψ, is inaccurate where edge surfaces are coarse. Other inverse inputs keep their input q.
• equilibrium_gse! returns its residual arrays (HDF5 output unchanged).

Results (one TokaMaker DIII-D-like H-mode solve, written by save_ifile and save_eqdsk; grid_type = "core_edge_packed" (then named "ldp"), psihigh = 0.995, Hamada; measured before the rebase onto develop)

• i-file and g-file give the same axis (to 1e-5 m), psio and q (to ~1e-5 at ψ_N 0.1–0.95).
• The i-file's Grad-Shafranov residual is ~2e-5 on typical surfaces at all resolutions tested (FastInterpolations ≤ 0.4.19). Its peaks (1–2e-2 at ψ_N 0.78–0.99) come from the TokaMaker solve: the g-file has them too, and 1 cm instead of 4 cm TokaMaker plasma cells lower them to ≤ 4.4e-3.
• Δ′(2/1), against Fortran RDCON on the same files (with the FastInterpolations periodic-spline fix, fix(cubic): wrong super-diagonal in PeriodicBC system on non-uniform grids ProjectTorreyPines/FastInterpolations.jl#220):

i-file	tabulated q	field-line q	RDCON
uniform ψ, 129 surfaces	6.392	8.285	8.280
65×129	7.696	8.306	8.280
edge-packed 129×257	8.339	8.404	8.318

Tests: two testsets in test/runtests_equil.jl on 137 KB of data from the same solve (test/test_data/TokaMaker_ifile/); runtests_equil.jl 320/320. Scripts: read_ifiles_scripts.zip. The δW scan requested in review is in this comment.

Dependencies: none. This PR was stacked on #506 and fed the i-file's FF′ and p′ records to its profile_source option. #506 is parked, so the PR is rebased onto develop without it; the reader skips those records until profile_source lands.

---

🤖 Generated with Claude Code

[matt-pharr]

@StuartBenjamin I am working on a code review for this, but in the mean time, this is a significant feature addition. In order to make sure everything is working well (and this should be cheap because the code is reasonably fast) I would really love to see some convergence scans showing convergence of n=1 delta W between gfiles and ifiles, or something similar you can use to convince us that this feature works properly.

Furthermore, this ifile_snap thing seems important! Hopefully we can get that merged into Tokamaker soon

[StuartBenjamin]

@matt-pharr here is the g-file vs i-file scan you asked for.

Setup

• One TokaMaker solver state (D3D-like H-mode, high-res DIII-D mesh, FE order 4), written two ways:
  • double-precision g-files (OFT VAC - REFACTOR - Give the Vacuum module its own output struct #358);
  • cut-cell i-files (OFT ALL - IMPROVEMENT - Self-describing HDF5 metadata (issue #226, PR B) #364, with the FF′/P′ and ψ-bounds records).
• n = 1 free-boundary δW (et[1]), no wall, forward integrator, ψ_lim set by dmlim 0.2.
• This branch with ForceFreeStates - BUGFIX! - Use trace-aware Frobenius exponents at resonant surfaces #517 and ForceFreeStates - BUGFIX! - Trigger the Gaussian reduction on the conditioning of U₁ and factor U₁ in the plasma response matrix #520 merged in.
• Two profile shapes:
  • "kinked": FF′ and P′ piecewise linear on 257 nodes, as bouquet's run-grid profiles are;
  • "smooth": the same profiles resampled smooth.

δW, kinked profiles (values; mtheta 512)

file	mpsi 512	1024	2048
g-file, 257 nodes	0.8058	0.8081	0.8097
g-file, 513 nodes	0.8093	0.8114	0.8124
i-file, 513 surfaces	0.8089	0.8112	0.8125
i-file, 1025 surfaces	–	0.8113	0.8123

• The 513-node g-file and the i-files agree to ≤ 0.05 % at every mpsi. The 257-node g-file is 0.3 % low.
• The larger error is GPEC's own radial grid, not the file type:
  • δW still rises 0.15–0.3 % per doubling of mpsi;
  • at mpsi 256 every file gives about 0.758 (6 % low).
• mtheta 1024 instead of 512 changes δW by < 2e-5.
• The smooth profiles show the same picture: g513 0.8127 / 0.8138, i513 0.8129 / 0.8139 at mpsi 1024 / 2048.
• Fortran RDCON's δW on the same files (mpsi 512): i-files 0.8118–0.8129, g513 0.8121.
• Δ′ agrees too. At mpsi 512 (Galerkin, without ForceFreeStates - BUGFIX! - Use trace-aware Frobenius exponents at resonant surfaces #517/ForceFreeStates - BUGFIX! - Trigger the Gaussian reduction on the conditioning of U₁ and factor U₁ in the plasma response matrix #520), Δ′(2/1) is 8.163 (g513), 8.157 (i513) and 8.153 (i1025). RDCON gives 8.182, 8.157 and 8.156.

ifile_snap. No longer needed. OFT #364's cut-cell writer puts every i-file point on the FE flux surface to ~1e-14 of the flux range, with no field-line tracing. OFT #357 (which had the snap) is closed.

What changed in this PR since you looked:

• the q of i-file equilibria now comes from the field-line integral, not the tabulated q. Δ′(2/1) was 6.39 instead of 8.29 on an unpacked i-file;
• the optional ψ-bounds record that OFT ALL - IMPROVEMENT - Self-describing HDF5 metadata (issue #226, PR B) #364 writes is now read.

Scripts and data are available if you want to rerun any of it.

[matt-pharr]

Made necessary changes to get this merged. Thank you for the very useful contribution that will aid tokamaker + GPEC coupling efforts!

[matt-pharr]

Violates the "no people's initials as code names" rule 😡

[matt-pharr]

I recommend you amend the name. You can note in the docstring that it is Don Pearlstein's format, but it should be called something more just like read_ifile