EQUIL - FEATURE - Optional profile_source = "integrate": build F and P in the sq spline from FF′ and p′ (values stays the default)

[StuartBenjamin]

Builds on #296. Adds an opt-in profile_source = "integrate" that builds F and p from the g-file's FF′ and p′ columns, the method of Julia GPEC #506. The default, values, is unchanged.

Problem

read_eq_efit reads FFPRIM and PPRIME but never uses them: it splines F and p and differentiates them instead. When a g-file holds single-precision values, as TokaMaker's do, and F varies by only ~2 % across the plasma, differentiating twice turns the rounding into tens of percent of noise in F″, which Δ′ depends on.

Changes

• spline_fit_hermite(spl, endmode, keep) (equil/spline.f): the usual fit, but keeps caller-supplied knot slopes for the columns flagged in keep.
• profile_source in equil.in, for efit and for ldp_i files with FF′/p′ records:
  • values (default): as before.
  • integrate: F²/2 and p integrated inward from the boundary, with the file's derivatives kept as knot slopes. Falls back to values with a warning if the derivatives are non-finite, F = 0 somewhere, the integrated F² ≤ 0, or only one of FF′, p′ has the opposite sign to its values.
• The slopes are carried through direct_run/inverse_run, newq0, sq_out and the dump file.
• read_eq_ldp_i reads optional FF′/p′ records. read_eq_efit no longer writes out of bounds when the g-file has no q column.
• eq_type = "dump" failed for direct equilibria ("I/O past end of record": a second, 4-entry sq%title assignment in direct_run shortened the record). Fixed; dumps now also store the slopes and eqfun, and older dumps still read.

Results

Default (values): bit-identical to #296 on DIIID_ideal (DCON), DIIID_resistive (RDCON) and the CI Solovev chains.

integrate: rms error of F″ (0.1 < ψ_N < 0.85) against an exact TokaMaker equilibrium written as g-files (Python replica of GPEC's sq fit, in the zip):

g-file	mpsi	values	integrate
257	128	1.8 %	0.8 %
257	512	12 %	~0.2 %
513	512	29 %	~0.1 %

values gets worse as the grid is refined; integrate converges. On STRIDE Δ′(2/1) it cuts the spread across three g-file resolutions from 1.43 to 0.96 (mpsi 256) and from 0.80 to 0.54 (mpsi 512). On double-precision i-files the two settings give STRIDE Δ′ within 0.6 % (2/1 within 0.2 %). When the values are already accurate there is no noise to remove, and integrate can shift pedestal surfaces, which is why values stays the default.

Tests

• make -C regression/spline_tests run (manual, not in CI): spline_fit_hermite keeps the supplied slopes exactly, and unflagged columns equal spline_fit.
• Comparisons above: spline_improvements_scripts.zip (commands in its README).

Dependencies

Contains #296 except its last commit (6eda483, RDCON only; no overlap). Conflicts with #300 only where both append code (equil/spline.f: keep both, renumber spline_fit_pchip to subprogram 25; the spline-test makefile and .gitignore: build both tests).

---

🤖 Generated with Claude Code

[d-burg]

A flag before this is reviewed: the Julia PR this ports (OpenFUSIONToolkit/GPEC#506) is now parked, and I would not make integrate the default here. Your OFT change (OpenFUSIONToolkit/OpenFUSIONToolkit#358) turns out to be the cleaner fix. Details are in OpenFUSIONToolkit/GPEC#504.

I ran this branch (a95a365b, STRIDE, psilow 0.01, psihigh 0.998, mtheta 512) on bouquet's public D3D-like example: four equilibria, each written by the single-precision writer and by the #358 writer from the same solver state.

1. With double-precision files, values is already clean. Same state, two g-files whose inputs differ only in the last digits; rms gap in Δ′(2/1) over three states:

mpsi	single, values	single, integrate	double, values	double, integrate
256	0.18	0.026	0.023	0.034
512	1.10	0.018	0.067	0.022

2. integrate changes files that have no noise to remove. integrate minus values on the double-precision files, mpsi 256:

state	2/1	3/1	4/1	5/1	6/1
A	+0.04	+0.04	0.00	+7.4 (of 182)	+19.1 (of −63)
B	−0.04	−0.04	+0.02	−4.6 (of 896)
draw 1	−0.01	+0.03	+0.14	−0.7 (of 245)
draw 2	+0.54	−0.01	+0.61	−2.9 (of 636)

At mpsi 512 draw 2 gives +0.15 on 2/1 and +1.5 on 4/1.

The cause is the quadrature. Integrating FF′ and p′ needs an assumed shape between nodes (a cubic spline here and in #506). Near steep or kinked pedestal profiles that differs from the file's own F and P. EFIT-family files tabulate FPOL and PRES as the trapezoid integral of FFPRIM and PPRIME, so no single rule matches every writer. In the Julia tests the cubic, trapezoid and Akima rules agree on Δ′(2/1) to 0.1 but differ by several units on surfaces in the pedestal foot.

Suggestion. Keep profile_source = "values" as the default, and keep integrate as an opt-in for g-files already written in single precision, with a note that outer surfaces can move. The rest of the PR (the Hermite fit, the dump fix, pchip for tabulated inputs, the reader bug fix) is not affected by this.

I can share the eight g-files and the run script if you want to reproduce it.

[StuartBenjamin]

Thanks, agreed. values is the default again (a47fc020), and integrate stays as an opt-in for single-precision g-files; the equil.in entry says that outer surfaces can move. With the default, the DIIID examples and the CI Solovev chains are bit-identical to the base branch. I also moved the pchip fit for tabulated Zeff and kinetic profiles to its own PR (#300), since it changes default results for runs that read a kinetic_file and has nothing to do with profile_source. 78977ab9 adds small review cleanups (bit-identical). The branch was force-pushed for the split.