Skip to content

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

Open
StuartBenjamin wants to merge 2 commits into
bugfix/efit-profiles-from-file-derivativesfrom
feature/read_ifiles_on_bugfix
Open

StuartBenjamin wants to merge 2 commits into
bugfix/efit-profiles-from-file-derivativesfrom
feature/read_ifiles_on_bugfix

Conversation

@StuartBenjamin

Copy link
Copy Markdown
Collaborator

Release note

  • Audience: users
  • Numerical impact: none (harness @ 6959627)
  • Migration: none

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

Regression report

regress --cases diiid_n1 --refs bugfix/efit-profiles-from-file-derivatives,feature/read_ifiles_on_bugfix: all 53 tracked quantities identical to the base branch.

Regression Report: diiid_n1
==========================================================================================================================================
Ref 1: bugfix/efit-profiles-from-file-derivatives  @ b3f412f35 (2026-10-06)
       env: julia 1.11.7, x86_64-linux-gnu, manifest f1e0eb38 (pinned), 2 threads/1 BLAS
Ref 2: feature/read_ifiles_on_bugfix  @ 6959627de (2026-10-06)
       env: julia 1.11.7, x86_64-linux-gnu, manifest f1e0eb38 (pinned), 2 threads/1 BLAS
------------------------------------------------------------------------------------------------------------------------------------------
Quantity                                      bugfix/efit-profiles-from-file-derivatives  feature/read_ifiles_on_bugfix  Diff       Status
------------------------------------------------------------------------------------------------------------------------------------------
total energy Re(et[1])                        7.959956e-01                                7.959956e-01                   0.0e+00    OK    
total energy Im(et[1])                        1.244153e-04                                1.244153e-04                   0.0e+00    OK    
plasma energy Re(ep[1])                       -1.361080e+00                               -1.361080e+00                  0.0e+00    OK    
vacuum energy Re(ev[1])                       2.157076e+00                                2.157076e+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)                             1419                                        1419                           0.0e+00    OK    
ODE steps (total)                             1414                                        1414                           0.0e+00    OK    
q0                                            1.204207e+00                                1.204207e+00                   0.0e+00    OK    
q95                                           4.781704e+00                                4.781704e+00                   0.0e+00    OK    
beta_t                                        1.328123e-02                                1.328123e-02                   0.0e+00    OK    
beta_n                                        1.373645e+00                                1.373645e+00                   0.0e+00    OK    
internal inductance li1                       8.842320e-01                                8.842320e-01                   0.0e+00    OK    
internal inductance li2                       7.080798e-01                                7.080798e-01                   0.0e+00    OK    
internal inductance li3                       7.304383e-01                                7.304383e-01                   0.0e+00    OK    
poloidal beta betap1                          6.686244e-01                                6.686244e-01                   0.0e+00    OK    
poloidal beta betap2                          5.354245e-01                                5.354245e-01                   0.0e+00    OK    
poloidal beta betap3                          5.523312e-01                                5.523312e-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.238951e-01                                4.238951e-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.006569e+00                                2.006569e+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)                          566f5ea79f6e...                             566f5ea79f6e...                identical  OK    
pressure profile (checksum)                   7342e0453497...                             7342e0453497...                identical  OK    
Mercier D_I profile (checksum)                a74297148c20...                             a74297148c20...                identical  OK    
resistive interchange D_R profile (checksum)  547606911e8d...                             547606911e8d...                identical  OK    
ballooning Delta' profile (checksum)          8272768c3101...                             8272768c3101...                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.245521e-04                                5.245521e-04                   0.0e+00    OK    
dominant-coupling singular values             [3 elem]                                    [3 elem]                       0.0e+00    OK    
|forcing overlap with dominant mode|          1.421101e-04                                1.421101e-04                   0.0e+00    OK    
|delta_nominal| of coil set 1                 7.082240e-05                                7.082240e-05                   0.0e+00    OK    
||ddelta/d(shift)|| over coil sets            1.245837e-04                                1.245837e-04                   0.0e+00    OK    
||ddelta/d(tilt)|| over coil sets             4.176884e-07                                4.176884e-07                   0.0e+00    OK    
PE plasma energy                              3.425341e+00                                3.425341e+00                   0.0e+00    OK    
PE vacuum energy                              3.174511e+00                                3.174511e+00                   0.0e+00    OK    
PE surface energy                             5.875246e+00                                5.875246e+00                   0.0e+00    OK    
PE toroidal torque                            -5.064459e-02                               -5.064459e-02                  0.0e+00    OK    
NTV torque FGAR [N·m]                         5.446590e-01                                5.446590e-01                   0.0e+00    OK    
NTV kinetic energy dW FGAR [J]                8.179981e-02                                8.179981e-02                   0.0e+00    OK    
Runtime (s)                                   396.5s                                      376.2s                                    --    
||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

Changes

  • read_ldp_i (src/Equilibrium/ReadEquilibrium.jl):
    • It reads the i-file into an InverseRunInput, the same input the CHEASE readers build.
    • Reals may be real8 or real4; the reader tells them apart by record length. The Fortran reader only takes real*8.
    • The optional trailing FF′ and p′ records, which TokaMaker writes, feed profile_source = "derivatives" through file_profiles, as the g-file reader does. Without them the reader warns and uses the tabulated F and P.
  • Rerun from gpec.h5: the reader stores the same InverseIngest as the CHEASE readers, so this works with no extra code.
  • equilibrium_gse! now returns its residual arrays: (; xs, ys, flux_x, flux_y, source, total, error, errori). Before, it discarded them, so the Grad-Shafranov residual could not be checked from code. Its HDF5 output is unchanged.
  • Unit test: "Load TokaMaker i-file (ldp_i)" in test/runtests_equil.jl, with 118 KB of data in test/test_data/TokaMaker_ifile/. That data is g65.geqdsk and i33x65.ifile (real8 and real4), all from one TokaMaker solve. The test checks:
    • the i-file agrees with the g-file
    • its GS residual is small, and five times below the g-file's
    • real4 agrees with real8
    • the alias "ifile" works
    • the derivatives and values routes for profile_source, and the fallback when the records are absent
  • Docs: the eq_type list, the profile_source section and the InverseIngest docstring now mention i-files.

Notes for reviewers

Stacked on bugfix/efit-profiles-from-file-derivatives. This PR targets that branch, not develop. It uses profile_source and file_profiles, which the bugfix adds. It should merge after the bugfix, and it reaches develop together with it. The harness compares against the bugfix branch so that only this PR's changes show; a run against develop moves 44 quantities, all from the bugfix.

Setup of the comparisons below. One TokaMaker reference equilibrium (DIII-D-like H-mode, q0 = 1.25, q95 = 4.49) was written both by save_ifile and by save_eqdsk. GPEC settings: grid_type = "ldp", psihigh = 0.995, Hamada coordinates. "GSE" is equilibrium_gse!'s local residual, normalized by the largest term.

i-file and g-file give the same equilibrium. The axis agrees to the last printed digit (1e-5 m), as does psio (0.259291). q agrees to about 1e-5 at ψ = 0.1–0.95. F and P match the exact TokaMaker profiles to 1e-7–1e-6 and about 1e-5.

GS residual (stock TokaMaker DIII-D mesh, 4 cm plasma cells), mpsi = 256.

file typical surface (median over θ-maxima, 0.05 < ψ < 0.95) worst surface
g-files 129 / 257 / 513 1.2e-3 / 4.3e-4 / 2.8e-4 6–9e-3
i-files 65×129 / 129×257 / 257×513 1.8e-5 / 1.7e-5 / 2.4e-5 0.9–2e-2

The worst-surface residual comes from TokaMaker's solution, not from the i-file.

  • It sits at ψ ≈ 0.78–0.99, at a few fixed poloidal angles (θ ≈ 0.045, 0.895, 0.97). This is where the pedestal p′ steepens.
  • The g-file shows the same jump at the same ψ and the same angles.
  • It hardly changes with i-file packing (pack_lcfs), 257 surfaces, or 1025 angles. Finer GPEC grids (larger mpsi or mtheta) make it slightly worse, not better.
  • Refining TokaMaker's plasma mesh removes it:
TokaMaker plasma cell size i-file 129×257, ψ 0.5–0.75 ψ 0.8–0.9 ψ 0.9–0.99 g-file 257, ψ 0.5–0.9
4 cm (stock mesh) 3–7e-5 1.5–4e-3 1.2–2.8e-2 4e-4 to 1.6e-3
2 cm (34,706 cells) 7e-6 to 3e-5 5e-5 to 5e-4 2.7e-3 to 2.3e-2 6e-4 to 1.4e-3
1 cm (101,560 cells) 6e-6 to 2.5e-5 5e-5 to 1e-4 4e-4 to 4.4e-3 1e-3 to 1.6e-3
  • The i-file's residual follows the solve's accuracy, while the g-file stays at about 1e-3, limited by its R,Z grid. On the 1 cm mesh the i-file beats the g-file at every ψ.
  • Recommendation for TokaMaker users: 2 cm plasma cells for routine work, 1 cm when the pedestal matters.
  • save_ifile's tracer tolerance is not a factor. Every point is Newton-snapped onto the exact ψ contour (ifile_snap).

profile_source = "derivatives" gives no gain on these files. "values" gives an equal or smaller F′ error: for the i129x257 file, 5.2e-4 against 1.2e-3. TokaMaker now writes both formats in double precision, so the tabulated F has no rounding for the derivative route to fix. The Julia reader also uses FF′ only to rebuild F at the grid points, then takes F′ from a spline. The Fortran reader instead passes FF′/F through as the spline slopes. This reading of the cause is a hypothesis; it bears on the bugfix PR, which makes "derivatives" the default.

Reproduce. The scripts ship as
read_ifiles_scripts.zip; its README.md has the commands.


Assignee and a named human reviewer required before merge; otherwise open as draft.

StuartBenjamin and others added 2 commits October 7, 2026 02:10
…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>
@github-actions github-actions Bot added the feature New capability label Oct 7, 2026
@matt-pharr
matt-pharr self-requested a review October 8, 2026 13:44
@matt-pharr

matt-pharr commented Oct 8, 2026 •

Copy link
Copy Markdown
Collaborator

@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
StuartBenjamin added this pull request to stack #516 October 8, 2026 22:08

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

feature New capability

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants