Skip to content

Oft interface improvements using ifiles - #298

Open
StuartBenjamin wants to merge 31 commits into
developfrom
OFT_interface
Open

StuartBenjamin wants to merge 31 commits into
developfrom
OFT_interface

Conversation

@StuartBenjamin

@StuartBenjamin StuartBenjamin commented Oct 6, 2026 •

Copy link
Copy Markdown
Collaborator

GPEC PR report: OFT_interface → develop

This PR adds 1 commit on top of PR #297.

Purpose

Getting OpenFUSIONToolkit/TokaMaker i-files accurately read with eq_type = "ldp_i". The work that went into this is spre spread over three PRs:

Changes

  1. read_eq_ldp_i takes |F| (19a8099b), as read_eq_efit does.
    • inverse_run computes q in proportion to F. An i-file with F < 0 (TokaMaker with B_t < 0) therefore stopped with "Invalid extrapolation near axis".
    • With the fix, flipping the sign of F in a file gives the same Δ′ to every digit: 8.50824 for both.
  2. The i-file format is documented in the reader header.

Results: i-file vs g-file

One TokaMaker reference equilibrium (DIII-D-like, q0 = 1.25, q95 = 4.49), written by save_ifile with FF′ and p′ records and by save_eqdsk, run through STRIDE with profile_source = integrate. The last surface is at the same true ψ_N = 0.985 in every run. Integrated GSE is the θ-integrated residual over the θ-integrated source, median over 0.05 < ψ_N < 0.95.

file size GSE local / integrated (mpsi = 128) Δ′(2/1), mpsi = 128 256 512
g-file 129×129 (efit) 0.29 MB 1.7e-4 / 7.1e-5 8.60 8.28 8.61
g-file 257×257 1.1 MB 1.8e-4 / 7.4e-5 8.57 9.24 8.80
g-file 513×513 4.3 MB 1.8e-4 / 1.3e-4 8.52 8.93 8.26
i-file 65×129 (ldp_i) 0.14 MB 2.4e-4 / 4.5e-5 8.58 8.57 8.60
i-file 129×257 0.54 MB 1.5e-4 / 2.4e-5 8.60 8.60 8.60
i-file 257×513 2.1 MB 1.7e-4 / 3.0e-5 8.59 8.50 8.62
ifile_vs_gfile
  • The i-files give Δ′(2/1) = 8.50–8.62 at every mpsi; the g-files scatter over 8.26–9.24 at mpsi ≥ 256 (single-precision ψ(R,Z)).
  • The i-files' integrated GSE is 1.6–5× lower than the g-files'.
  • A 129×257 i-file (0.54 MB) beats the 513×513 g-file (4.3 MB) on both counts.

Status

  • Built; the negative-F test above passes
  • CI Solovev examples. Compared with Zeff_profile_support 99378252, which also has the vacuum fix, on the full chain (OFT_interface = spline_improvements + this commit):
    • The DCON, STRIDE, RDCON and GPEC outputs of the ideal and resistive Solovev examples are bit-for-bit identical.
    • The kinetic example differs by up to 8e-5, and PENTRC by up to 1e-6 in n and T, because PENTRC's kinetic table now uses pchip. The exceptions are the last 2 of 17 PENTRC output points next to the edge (ψ_N > 0.99), where T → 0: there logLambda/nu were already unphysical (negative) with the cubic spline and change value.
    • The example pass/fail pattern is unchanged; the gpec/pentrc/rmatch "failures" in the DIIID examples happen the same way on the base build.

How to reproduce

The scripts and their inputs ship with this PR as OFT_interface_scripts.zip: its README.md has the commands.

StuartBenjamin and others added 30 commits June 21, 2026 21:12
…explitcit indexing control to match Zeff spline assignment
…r frame

VACUUM runs on a theta-reversed boundary, so wv comes back complex-conjugated
relative to wp; conjugate it after mscvac (free_run and free_wvmats). The
eigenvector norm used jmat(jpert-ipert), the transpose of the plasma-matrix
index order; use jmat(ipert-jpert). Up-down symmetric cases are unchanged;
asymmetric ones no longer depend on the theta origin.
…er frame

Same change as DCON: conjugate wv after mscvac (free_run and
free_get_wvac, used by the Galerkin solve) and transpose the norm.
…ier frame

STRIDE called mscvac with complex_flag false, keeping only Re(wv), which is
itself origin-dependent. Keep the full matrix, conjugate it as in DCON, and
transpose the norm.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…ler-supplied knot derivatives

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…' and p' (efit, ldp_i)

New equil.in variable profile_source:
  values    - tabulated f, p; spline-fitted derivatives (default, unchanged)
  hermite   - tabulated f, p with the tabulated derivatives as knot slopes
  integrate - f**2/2 and p integrated inward from the boundary values
              (as in Julia GPEC PR #506), tabulated derivatives as slopes
read_eq_efit no longer discards FFPRIM and PPRIME, and fixes an
out-of-bounds index when the g-file has no q profile. read_eq_ldp_i
reads optional trailing FF' and p' records. Falls back to values with a
warning when the derivatives are non-finite or disagree in sign.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…nverse sq construction

sq_in and sq are fit with spline_fit_hermite so the profile slopes
from profile_source are not re-derived by spline differentiation.
newq0 rescaling updates the f slope consistently (f' -> f'/ffac).
With profile_source = values, output is bit-for-bit unchanged.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Against a TokaMaker equilibrium with known profiles, integrate gives F'' errors
of 0.1-0.8% (0.1 < psi_n < 0.85) that fall with resolution, where values and
hermite give 2-47% that grow with mpsi (single-precision f amplified by
differentiation). STRIDE Delta' is correspondingly more consistent across
g-file resolutions. Set profile_source="values" to recover the previous
behaviour.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
… end-slope order

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…dcon-vacuum-theta-frame)

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
A second assignment of 4 titles reallocated sq%title, so the dump record
was 6 bytes short and read_eq_dump failed (I/O past end of record).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
equil_out_dump appends the supplied-derivative flags with sq%fs1(:,1:2), then
eqfun; read_eq_dump reads them when present (older dumps read as before) and
no longer fails when the equilibrium is re-read (STRIDE psilim reform).
DCON efit -> dump -> dcon reproduces f, mu0p, q, D_I, D_R exactly.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…bulated data without derivatives

Fritsch-Carlson interior slopes with scipy PchipInterpolator's end slopes, so
Python front ends (bouquet, TPS) can reproduce GPEC's interpolant exactly.
Unit test checks values and derivatives against scipy.

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

The tabulated Zeff (RDCON MRE terms) and the PENTRC kinetic input table carry
pedestal-scale gradients; a cubic spline rings there, pchip does not.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Derivative sign taken from f and p together, so a flat p (or f) no longer
forces the values fallback; the dump writes eqfun only when it exists;
shorter equil.in entry; unit test drops redundant wrappers.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Tabulated values paired with tabulated derivatives put the values' rounding
noise into f'' (the worst option in the comparison). spline_fit_hermite stays,
used by integrate to keep the derivatives.

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

With f < 0 (e.g. a TokaMaker equilibrium with B_t < 0) inverse_run built a
negative q and stopped ('Invalid extrapolation near axis'). Flipping the sign of
f in an i-file now gives identical results.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
ev is COMPLEX but zheev returns REAL eigenvalues, so W_v_eigenvalue in the
netcdf output had the eigenvalues packed pairwise into the real and imaginary
parts of the first half of ev, with stale values in the rest. zheev also got
lwork=2*mpert+1 for a work array of 2*mpert-1. Use a REAL buffer and the
correct lwork. Only W_v_eigenvalue changes; all other outputs (Delta', W_t,
eigenvectors) are bit-identical on the DIII-D resistive example.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…rf info)

gal_solve ignored info from zgbtrf/zpbtrf, so an indefinite matrix with solver='cholesky'
gave a garbage Delta' with no warning (g147131 n=1: 2/1 = 2.3e7). Now program_stop with the
info value. LU results unchanged (bit-identical Delta' on g147131).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…Galerkin matrix is not positive definite

zpbtrf info was ignored, so solver='cholesky' on an indefinite matrix gave a garbage Delta'
(g147131 n=1: 2/1 = 2.3e7). Now the lower band is kept, and on zpbtrf failure the same
Hermitian matrix is solved by LU (2/1 = 7.3787; solver='LU' gives 7.3641, its boundary rows
differ). zgbtrf failure stops the run. solver='LU' results are bit-identical.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…flowed to '**' at 25+ singular surfaces and crashed rdcon writing delta_gw.out)

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

logan-nc commented Oct 9, 2026

Copy link
Copy Markdown
Contributor

@StuartBenjamin we are trying to wrap up fortran dev... fill out the reviewer and assignee and labels metadata asap to get this in before 1.6. I suggest Matt or Daniel. Same for your other PRs.

Proceed without me

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

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants