Equilibrium - FEATURE - Read inverse i-files (ldp_i) such as TokaMaker save_ifile output by StuartBenjamin · Pull Request #511 · OpenFUSIONToolkit/GPEC

StuartBenjamin · 2026-10-07T09:44:09Z

Release note

Audience: users
Numerical impact: none against the code before Equilibrium - BUGFIX! - Build F and P of g-file and IMAS equilibria from the tabulated FF′ and p′ #506: with the "values" default, diiid_n1 and solovev_n1 are identical (harness @ 88f6ecc)
Migration: none

GPEC can now read inverse i-files (eq_type = "ldp_i" or "ifile"), such as OpenFUSIONToolkit TokaMaker's save_ifile output. This is a port of the Fortran read_eq_ldp_i. profile_source (added by #506) now defaults to "values", the tabulated F and P; "derivatives" stays an option.

Regression report

regress --cases diiid_n1,solovev_n1 --refs 0e68a0553,88f6ecc96 (the develop commit #506 branches from vs this head; julia 1.11.7, pinned manifest, 2 threads): diiid_n1 53 unchanged, solovev_n1 22 unchanged. Against #506 itself (b3f412f), diiid_n1 changes (44 of 53 quantities), because #506 reads the g-file with "derivatives"; this also makes the STRIDE Δ′ check in runtests_parallel_integration.jl that #506 fails pass again (94/94). 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: 0e68a0553  @ 0e68a0553 (2026-10-05)
       env: julia 1.11.7, x86_64-linux-gnu, manifest 6e22aba4 (pinned), 2 threads/2 BLAS
Ref 2: 88f6ecc96  @ 88f6ecc96 (2026-10-13)
       env: julia 1.11.7, x86_64-linux-gnu, manifest 6e22aba4 (pinned), 2 threads/2 BLAS
-------------------------------------------------------------------------------------------------
Quantity                                      0e68a0553        88f6ecc96        Diff       Status
-------------------------------------------------------------------------------------------------
total energy Re(et[1])                        8.013468e-01     8.013468e-01     0.0e+00    OK    
total energy Im(et[1])                        4.143487e-05     4.143487e-05     0.0e+00    OK    
plasma energy Re(ep[1])                       -1.348157e+00    -1.348157e+00    0.0e+00    OK    
vacuum energy Re(ev[1])                       2.149503e+00     2.149503e+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)                             1431             1431             0.0e+00    OK    
ODE steps (total)                             1426             1426             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.842212e-01     8.842212e-01     0.0e+00    OK    
internal inductance li2                       7.080713e-01     7.080713e-01     0.0e+00    OK    
internal inductance li3                       7.304295e-01     7.304295e-01     0.0e+00    OK    
poloidal beta betap1                          6.680737e-01     6.680737e-01     0.0e+00    OK    
poloidal beta betap2                          5.349836e-01     5.349836e-01     0.0e+00    OK    
poloidal beta betap3                          5.518763e-01     5.518763e-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)                          2cdacd6807fb...  2cdacd6807fb...  identical  OK    
pressure profile (checksum)                   3d101fa873b1...  3d101fa873b1...  identical  OK    
Mercier D_I profile (checksum)                25ada2b687ca...  25ada2b687ca...  identical  OK    
resistive interchange D_R profile (checksum)  7dc243ce58a9...  7dc243ce58a9...  identical  OK    
ballooning Delta' profile (checksum)          9140d3dfe986...  9140d3dfe986...  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.187848e-04     5.187848e-04     0.0e+00    OK    
dominant-coupling singular values             [3 elem]         [3 elem]         0.0e+00    OK    
|forcing overlap with dominant mode|          1.415721e-04     1.415721e-04     0.0e+00    OK    
|delta_nominal| of coil set 1                 7.055416e-05     7.055416e-05     0.0e+00    OK    
||ddelta/d(shift)|| over coil sets            1.246245e-04     1.246245e-04     0.0e+00    OK    
||ddelta/d(tilt)|| over coil sets             4.180789e-07     4.180789e-07     0.0e+00    OK    
PE plasma energy                              3.422634e+00     3.422634e+00     0.0e+00    OK    
PE vacuum energy                              3.174510e+00     3.174510e+00     0.0e+00    OK    
PE surface energy                             5.825956e+00     5.825956e+00     0.0e+00    OK    
PE toroidal torque                            5.085625e-02     5.085625e-02     0.0e+00    OK    
NTV torque FGAR [N·m]                         5.294859e-01     5.294859e-01     0.0e+00    OK    
NTV kinetic energy dW FGAR [J]                7.920700e-02     7.920700e-02     0.0e+00    OK    
Runtime (s)                                   345.1s           335.5s                      --    
||forcing b~|| (root-area-weighted)           4.683329e-04     4.683329e-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: 0e68a0553  @ 0e68a0553 (2026-10-05)
       env: julia 1.11.7, x86_64-linux-gnu, manifest 6e22aba4 (pinned), 2 threads/2 BLAS
Ref 2: 88f6ecc96  @ 88f6ecc96 (2026-10-13)
       env: julia 1.11.7, x86_64-linux-gnu, manifest 6e22aba4 (pinned), 2 threads/2 BLAS
---------------------------------------------------------------------------------
Quantity                      0e68a0553        88f6ecc96        Diff       Status
---------------------------------------------------------------------------------
total energy Re(et[1])        3.428086e-01     3.428086e-01     0.0e+00    OK    
total energy Im(et[1])        -1.055018e-03    -1.055018e-03    0.0e+00    OK    
plasma energy Re(ep[1])       -1.005939e+01    -1.005939e+01    0.0e+00    OK    
vacuum energy Re(ev[1])       1.040220e+01     1.040220e+01     0.0e+00    OK    
vacuum matrix min eigenvalue  2.171292e+00     2.171292e+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)             405              405              0.0e+00    OK    
ODE steps (total)             403              403              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.621020e-02     4.621020e-02     0.0e+00    OK    
beta_n                        3.215791e+00     3.215791e+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)          90ffc338a738...  90ffc338a738...  identical  OK    
pressure profile (checksum)   5040338aec07...  5040338aec07...  identical  OK    
ca_left (checksum)            0de08e7b7c34...  0de08e7b7c34...  identical  OK    
Runtime (s)                   208.8s           206.9s                      --    
=================================================================================
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_ldp_i (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).
Optional trailing records, told apart by length: FF′ and p′ feed profile_source through file_profiles, as in the g-file reader; [ψ_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 = "ldp", psihigh = 0.995, Hamada)

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 329/329. Scripts: read_ifiles_scripts.zip. The δW scan requested in review is in this comment.

Dependencies: stacked on #506 (this PR's base), which adds profile_source and file_profiles; merge after it. This PR changes #506's default from "derivatives" to "values".

🤖 Generated with Claude Code

…r save_ifile output Port of Fortran read_eq_ldp_i. eq_type = "ldp_i" (alias "ifile") reads the sequential unformatted i-file (real*8 or real*4) into an InverseRunInput. The optional trailing FF′ and p′ records feed profile_source = "derivatives" via file_profiles, as for g-files; without them the tabulated F and P are used. equilibrium_gse! now returns its residual arrays and the θ-integrated residual per surface, so the Grad-Shafranov error can be checked programmatically. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

Checks i-file vs g-file agreement from one TokaMaker solve, GS residual, real*4 vs real*8, the "ifile" alias, and the profile_source fallback. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

matt-pharr · 2026-10-08T14:21:21Z

@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

…heir field-line integral The inverse solver kept the i-file's tabulated q, splined in psi, while the SFL angle nu and the metric are built from the field-line q of the interpolated geometry. Where the input surfaces are coarse near the edge (q ~ log(1-psi)) the two disagree (q up to +0.7 at psi_N 0.9985 on a uniform 129-surface file), which moves psilim and spoils Delta' (2/1: 6.39 vs 8.28 for RDCON). ldp_i/ifile now use the field-line q, as Fortran inverse_run does; LAR, CHEASE and efit_by_inversion keep their input q. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

…ker i-files read_ldp_i tells the optional records apart by length: FF' and p' (npsi each), then [psi_axis, psi_LCFS] (2). With the bounds, psi_N and psio are normalised to the true LCFS, so the padded last surface sits below 1; psihigh is clamped to the outermost input surface for inverse inputs (resolve_psihigh!), on the reader and rerun paths alike. Files without the record read unchanged. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

…q and profile docs LDP_I_KINDS replaces the two hard-coded ("ldp_i", "ifile") lists, so the inverse solver does not repeat the reader's kinds. One comment at q_from_geometry explains both q sources. The profile docstrings no longer say direct-only or x = 1 at the boundary, since the i-file reader calls them with its last node at 1 - pad. psihigh_resolved's docstring says the clamp applies to i-files and reruns. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

… it to the g-file's The ratio to the g-file residual held only because FastInterpolations' periodic spline is wrong on non-uniform θ grids (FastInterpolations.jl#220), which inflates the direct refit's residual. With the fix the g-file residual is 6e-6 on this data, below the i-file's 1e-4, so the ratio check would fail once the fix is released. The absolute bound stays. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

StuartBenjamin · 2026-10-13T06:36:16Z

@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.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

…s" stays an option Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

StuartBenjamin and others added 2 commits October 7, 2026 02:10

Equilibrium - TEST - Add TokaMaker i-file (ldp_i) reader test

6959627

Checks i-file vs g-file agreement from one TokaMaker solve, GS residual, real*4 vs real*8, the "ifile" alias, and the profile_source fallback. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

github-actions Bot added the feature New capability label Oct 7, 2026

matt-pharr self-requested a review October 8, 2026 13:44

matt-pharr assigned StuartBenjamin Oct 8, 2026

StuartBenjamin added this pull request to stack #516 October 8, 2026 22:08

StuartBenjamin and others added 4 commits October 12, 2026 15:27

This was referenced Oct 13, 2026

Equilibrium - g-file reader differentiates FPOL and PRES and discards FFPRIM and PPRIME, putting rounding noise on Δ′ #504

Open

Equilibrium - BUGFIX - Support original and rho grids for inverse equilibria #521

Closed

Equilibrium - MINOR - Tighten the i-file reader docstrings and tests

14445e3

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

StuartBenjamin marked this pull request as draft October 13, 2026 09:37

Equilibrium - MINOR - Default profile_source to "values"; "derivative…

88f6ecc

…s" stays an option Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

StuartBenjamin mentioned this pull request Oct 13, 2026

Equilibrium - BUGFIX - Support original and rho grids for inverse equilibria #523

Merged

StuartBenjamin marked this pull request as ready for review October 13, 2026 11:27

matt-pharr added this to the GPEC v2.0.0 milestone Oct 16, 2026

Navigation Menu

Repository navigation

Equilibrium - FEATURE - Read inverse i-files (ldp_i) such as TokaMaker save_ifile output - #511

Equilibrium - FEATURE - Read inverse i-files (ldp_i) such as TokaMaker save_ifile output#511
StuartBenjamin wants to merge 8 commits into
bugfix/efit-profiles-from-file-derivativesfrom
feature/read_ifiles_on_bugfix

StuartBenjamin commented Oct 7, 2026 •

edited

Loading

Uh oh!

matt-pharr commented Oct 8, 2026 •

edited

Loading

Uh oh!

StuartBenjamin commented Oct 13, 2026 •

edited

Loading

Uh oh!

This branch has not been deployed

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

Conversation

StuartBenjamin commented Oct 7, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Release note

Regression report

Notes for reviewers

Uh oh!

matt-pharr commented Oct 8, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

StuartBenjamin commented Oct 13, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

This branch has not been deployed

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

StuartBenjamin commented Oct 7, 2026 •

edited

Loading

matt-pharr commented Oct 8, 2026 •

edited

Loading

StuartBenjamin commented Oct 13, 2026 •

edited

Loading