EQUIL - BUGFIX - ldp_i reader: accept F < 0 and an optional psi bounds record

[StuartBenjamin]

GPEC PR report: OFT_interface → develop

This PR was cut down to one issue, the ldp_i reader for OpenFUSIONToolkit/TokaMaker i-files. #297 and #296 have since merged, and develop is merged into this branch (31cb34b1).

The RDCON fixes that were stacked here moved to their own PRs:

• Galerkin factorization failures (d58c2622, dfd771ca) → RDCON - BUGFIX - Galerkin factorization failures: Cholesky falls back to LU, a failed LU stops #301.
• Δ-matrix output format at ≥ 25 surfaces (cc7f996a) → RDCON - BUGFIX - Delta-matrix output formats sized with i0 #302.
• The vacuum zheev real-eigenvalue fix (c9642696) is dropped: GPEC, DCON, RDCON, STRIDE - BUGFIX - Fix LAPACK and ZVODE type and size errors #299 made the same change.

Purpose

Read OpenFUSIONToolkit/TokaMaker i-files correctly with eq_type = "ldp_i". The work is spread over:

• OpenFUSIONToolkit #364 (replaces the closed #357): points lie exactly on their flux surfaces; save_ifile can add FF′, p′ records (profile_derivs) and the psi bounds record (psi_bounds).
• EQUIL - FEATURE - Optional profile_source = "integrate": build F and P in the sq spline from FF′ and p′ (values stays the default) #297 (spline_improvements, merged): can build F and p from the FF′, p′ records with the opt-in profile_source = "integrate". The default is values. For i-files (double precision) the two agree to ≤ 0.2 % on Δ′(2/1) at mpsi ≤ 512, so either works.
• This PR: the reader fixes below.

Changes

1. read_eq_ldp_i takes |F| (8a578482), as read_eq_efit does.
   • inverse_run computes q in proportion to F. An i-file with F < 0 (TokaMaker with B_t < 0) stopped with "Invalid extrapolation near axis".
   • With the fix, flipping the sign of F gives bit-identical output.
   • The i-file format is documented in the reader header.
2. Optional psi bounds record [psi_axis, psi_lcfs] (608e0942).
   • TokaMaker pads the last i-file surface slightly inside the LCFS, so psi_N = 1 used to be the padded surface: psio too small, every psi_N (rational surfaces, psihigh) stretched by 1/(1 − pad).
   • With the record, psi_N and psio use the true axis and LCFS flux. A psihigh above the last file surface is clamped to it, with a note. If the file surfaces fall outside [0, 1] in the new psi_N (psi_axis outside the first surface, or the last surface outside psi_lcfs), the run stops and prints psi_axis, psi_lcfs, psi(0) and psi(mx) (8dc2d0da). psi_axis slightly inside the first surface is not checked. qa is now extrapolated to the true LCFS.
   • Files without the record give bit-identical results. The reader skips the optional FF′, p′ records and tells the records apart by length: a length-mx read of the 2-value bounds record fails, and the reader then rewinds and skips the records it has read. It does not BACKSPACE after the failed read, because the file position after an I/O error is processor dependent: with ifx/ifort, BACKSPACE there misses the bounds record (checked with a stand-alone test). Older builds ignore the new record.

Evidence

A D3D-like TokaMaker i-file (129 × 257) in variants with and without the FF′, p′ and bounds records (bounds 1 % beyond the last surface), run through RDCON on develop and on this branch:

check	result
files without the bounds record, this PR vs develop	bit-identical
files with the record, on develop	bit-identical to those without, so develop ignores it
psio with the record	×1.01 exactly; rational psi_N 2/1, 3/1, 4/1 ÷1.01 exactly
psihigh above the last surface	clamped to 0.990099, with a note
bounds that put the surfaces outside [0, 1] (last surface beyond psi_lcfs, or psi_axis beyond the first surface)	stops with "ldp_i psi bounds disagree with psi" and prints psi_axis, psi_lcfs, psi(0), psi(mx)
every record layout (none; FF′+p′; FF′+p′+bounds; bounds only; FF′+bounds), stand-alone reader test, native and big-endian	bounds found exactly when present, with gfortran 14, ifx 2023.2 and ifort 2023.2 (the earlier BACKSPACE version missed them with ifx/ifort)
F < 0	develop stops; this PR gives output bit-identical to F > 0
DCON DIIID_ideal_example (g-file)	bit-identical to develop
after merging develop (#297): DCON on a TokaMaker i-file (33 × 65) in every record layout, big-endian, F < 0, profile_source = values and integrate (gfortran 16)	bounds found exactly when present; rational psi_N ÷1.01 exactly; F < 0 and big-endian identical to the native F > 0 file

Earlier i-file vs g-file comparison (on the previous version of this PR, with #297): the i-files give Δ′(2/1) = 8.50–8.62 at every mpsi from 128 to 512, while 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.

Status

• make rdcon dcon builds (gfortran 14, Perlmutter)
• Tests above (scripts: ldp_i_reader_scripts.zip)
• Merged develop after EQUIL - FEATURE - Optional profile_source = "integrate": build F and P in the sq spline from FF′ and p′ (values stays the default) #297: read_eq_ldp_i keeps EQUIL - FEATURE - Optional profile_source = "integrate": build F and P in the sq spline from FF′ and p′ (values stays the default) #297's FF′, p′ reads and this PR's |F| and bounds read; tested as above.

How to reproduce

ldp_i_reader_scripts.zip

PY_NC=python3 ./run_variants.sh /path/develop/bin/rdcon /path/pr/bin/rdcon out
PY_NC=python3 ./run_fneg.sh     /path/develop/bin/rdcon /path/pr/bin/rdcon out

🤖 Generated with Claude Code

[logan-nc]

@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

[matt-pharr]

This checks out, I updated the PR and added one minor safety check to your new code. Otherwise good and I am merging.