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

Activity

1. d-burg

   added

   bugSomething isn't workingSomething isn't working

   on Oct 5, 2026
2. d-burg commented on Oct 5, 2026

   

   d-burg

   on Oct 5, 2026

   CollaboratorAuthor

   More actions

   Correction to the first caveat above. The large 4/1 changes (13 as written, 35 or 50 rebuilt) came from rebuilding F alone, which leaves it inconsistent with the as-written P near the edge. With F and P both rebuilt from FFPRIM and PPRIME:

   • the same-state gap stays closed on all three surfaces (rms B − A: 2/1 0.16, 3/1 0.17, 4/1 0.22; as written 3.70, 0.67, 1.20);
   • 4/1 returns to or near its as-written value in three of five pairs (for example −16.7 → −16.6) and is 4–7 away in the others.

   So a fix should take F and P together from FFPRIM and PPRIME, never one alone. The end-node features remain an open but smaller question. Rebuilding P alone moves Δ′(2/1) by ≤ 0.27 on eight of the ten files; the kinetic EFIT 257 file moves by 2.1, cause not yet identified.
3. d-burg commented on Oct 9, 2026

   

   d-burg

   on Oct 9, 2026

   CollaboratorAuthor

   More actions

   Update: the writer fix is enough for new files; the GPEC-side change (#506) is parked

   OpenFUSIONToolkit/OpenFUSIONToolkit#358 (merged 2026-10-07) makes TokaMaker's gs_save_eqdsk write doubles instead of single precision. I tested it with GPEC unchanged.

   Test. bouquet's public D3D-like H-mode example: four equilibria, each written by the old and the new writer from the same solver state (one OFT build with only that change added). develop 0e68a0553, n = 1, Riccati, settings as in the issue.

   Same state, two g-files whose inputs differ only in the last digits. rms gap in Δ′ over three states:

   grid	writer	2/1	3/1	4/1
   ldp 256	single	0.29	0.83	1.08
   ldp 256	double	0.014	0.021	0.088
   ldp 512	single	1.28	1.98	1.26
   ldp 512	double	0.07	0.08	0.13
   auto	single	1.76	2.25	1.59
   auto	double	0.42	0.14	0.41

   • On fixed grids the gap falls by about a factor of 20.
   • The auto-grid residual is the grid's own response to small input changes (one state at −0.72 on 2/1, two near 0.1). That is a separate matter from this issue.
   • develop on double-precision files agrees with Equilibrium - BUGFIX! - Build F and P of g-file and IMAS equilibria from the tabulated FF′ and p′ #506 to 0.04 on 2/1, 3/1 and 4/1 (ldp 256).

   Why #506 is parked. Integrating FF′ and p′ needs an assumed shape between nodes, and near profile kinks that assumption moves Δ′ by more than the noise it removes.

   • Public example, 6/1 surface, ldp 256: −20.5 (single) and −20.3 (double) on develop, −11.2 with Equilibrium - BUGFIX! - Build F and P of g-file and IMAS equilibria from the tabulated FF′ and p′ #506.
   • On the five same-state pairs in the issue, Δ′(4/1) differs from the value given by the file's own F and P by 5.1 rms with the cubic rule of Equilibrium - BUGFIX! - Build F and P of g-file and IMAS equilibria from the tabulated FF′ and p′ #506, 4.9 with trapezoid and 1.5 with Akima. Δ′(2/1) is the same for all three rules to within 0.1.
   • No rule suits every writer. EFIT and kinetic-EFIT files tabulate FPOL and PRES as the trapezoid integral of FFPRIM and PPRIME (trapezoid reproduces them to rounding); TokaMaker files with piecewise-linear profiles do not follow that rule.
   • This corrects a caveat in the issue: the end nodes of FFPRIM are not what moves the edge surface. Changing the axis node has no effect on Δ′, and the boundary node moves 4/1 by at most 0.4. The kinks in the pedestal foot are the cause.

   Still open.

   • g-files already written in single precision (older TokaMaker output, some EFIT output) keep the noise. A reader warning when FPOL is exactly single precision would at least make it visible.
   • PSIRZ matters too: single-precision ψ(R,Z) alone puts about 0.03 on Δ′(2/1) per file on a fixed grid. Only the writer can fix that, and VAC - REFACTOR - Give the Vacuum module its own output struct #358 does.
   • Nine printed digits leave rounding of about 1e-9. I have not tested g-files finer than 257 nodes, where the amplification is larger.
   • Fortran GPEC behaves the same way; see the comment on EQUIL - FEATURE - Optional profile_source = "integrate": build F and P in the sq spline from FF′ and p′ (values stays the default) PrincetonUniversity/GPEC#297.
4. StuartBenjamin commented on Oct 13, 2026

   

   StuartBenjamin

   on Oct 13, 2026

   Last edited by StuartBenjamin

   Collaborator

   More actions

   @d-burg thanks for the update. I tried to reproduce the "integrate moves the pedestal surfaces" result with files where I know the answer: one TokaMaker state written both as double-precision g-files and as i-files. An i-file can be written with as many surfaces as you like, so I could follow each rule as the file gets finer.

   Setup

   • D3D-like H-mode, high-res mesh, FE order 4. FF′ and P′ are piecewise linear on 257 nodes, as in bouquet's run grid, so there are kinks at every node and in the pedestal foot.
   • Files: g-files with 257 and 513 nodes; i-files with 129–1025 surfaces.
   • Codes: Julia GPEC Galerkin (with the i-file reader of Equilibrium - FEATURE - Read inverse i-files (ldp_i) such as TokaMaker save_ifile output #511) and Fortran RDCON; mpsi 512, mtheta 512, psihigh 0.999 with dmlim truncation.

   Your effect reproduces on the 257-node g-file:

   	values	integrate (#506 cubic rule)
   Δ′(2/1)	8.242	8.181
   Δ′(4/1)	−17.45	−17.63
   Δ′(5/1)	110.7	116.5

   But both rules converge to the same limit as the file gets finer:

   Julia GPEC, values / integrate	2/1	4/1	5/1
   g-file 513	8.163 / 8.157	−16.90 / −17.48	106.6 / 104.0
   i-file 513	8.157 / 8.137	−17.72 / −17.54	100.44 / 100.55
   i-file 1025	8.153 / 8.151	−17.64 / −17.61	100.15 / 100.18

   What this shows:

   • On the 257-node g-file, 5/1 is off by +10 (values) and +16 (integrate). At 513 nodes it is +6 and +4. Neither rule is closer at both resolutions.
   • So on these files "integrate − values" measures the resolution error of the g-file at the pedestal surfaces. It is not an integrate defect, and the file's own F and P are no better a reference.
   • A check without GPEC agrees. I rebuilt F′, F″, P′ and P″ from the file nodes with each rule and compared them with TokaMaker's own profiles over ψ_N 0.85–0.995:
     • values, cubic, trapezoid and Hermite-integrate are all wrong by the same order (within ×1.5);
     • Equilibrium - BUGFIX! - Build F and P of g-file and IMAS equilibria from the tabulated FF′ and p′ #506's cubic rule is the worst of the four, by a little.

   The Fortran integrate rule has a weakness. Fortran's integrate (#297) uses integrated values together with the file's FF′/F as Hermite knot slopes. That makes the values and slopes slightly inconsistent, and it is the most erratic rule at 4/1 in RDCON (−17.1 … −18.6 across files). RDCON's 4/1 and 5/1 also scatter by ±0.4 and ±10 with smooth profiles, so part of that is the solver, not the profiles.

   What I'm doing with it:

   • Fortran PERTURBED EQUILIBRIUM - BUG FIX - Use Λ† in reluctance to match Fortran gpresp_reluct #297 now defaults to values, as you suggested, with integrate opt-in.
   • For TokaMaker i-files (binary real*8, 513 or more surfaces), values and integrate agree through 5/1 in Julia GPEC, so either is fine there.
   • For 9-digit g-files the pedestal Δ′ is resolution-limited at 257–513 nodes whatever the rule. Tabulated values with FF′ slopes ("hermite") is the best rule on i-files but the worst on 9-digit g-files, because print rounding gets into F″.

   Plots and scripts are available if useful.