Equilibrium - BUGFIX! - Place read_imas profiles on the stored psi grid

[calin1989]

Release note

• Audience: users
• Numerical impact: IMAS equilibria whose profiles_1d grid is not uniform in ψ_N, or does not start at psi_axis; none for the regression suite (harness @ 0d056d4)
• Migration: Re-run equilibria read from IMAS files whose profiles_1d grid is not uniform in ψ_N.

read_imas now places F, p and q on the ψ grid stored in profiles_1d.psi. It previously assumed a uniform ψ_N grid, so on non-uniform grids the pressure, and with it β and stability, was evaluated on the wrong flux surfaces.

Fixes #482

Regression report

Regression Report: diiid_n1
===================================================================================================================================================
Ref 1: 11971072cf4f2fa3be41a42bedab506608f5495e  @ 11971072c (2026-10-14)
       env: julia 1.12.7, x86_64-linux-gnu, manifest aa94b997 (pinned), 4 threads/4 BLAS, jit build
Ref 2: 0d056d4d37c9d61e539e76b1cb16357bb084cff5  @ 0d056d4d3 (2026-10-14)
       env: julia 1.12.7, x86_64-linux-gnu, manifest aa94b997 (pinned), 4 threads/4 BLAS, jit build
---------------------------------------------------------------------------------------------------------------------------------------------------
Quantity                                      11971072cf4f2fa3be41a42bedab506608f5495e  0d056d4d37c9d61e539e76b1cb16357bb084cff5  Diff       Status
---------------------------------------------------------------------------------------------------------------------------------------------------
total energy Re(et[1])                        8.008331e-01                              8.008331e-01                              0.0e+00    OK    
total energy Im(et[1])                        7.926479e-05                              7.926479e-05                              0.0e+00    OK    
plasma energy Re(ep[1])                       -1.371895e+00                             -1.371895e+00                             0.0e+00    OK    
vacuum energy Re(ev[1])                       2.172728e+00                              2.172728e+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)                             1415                                      1415                                      0.0e+00    OK    
ODE steps (total)                             1410                                      1410                                      0.0e+00    OK    
q0                                            1.204212e+00                              1.204212e+00                              0.0e+00    OK    
q95                                           4.781723e+00                              4.781723e+00                              0.0e+00    OK    
beta_t                                        1.327026e-02                              1.327026e-02                              0.0e+00    OK    
beta_n                                        1.372510e+00                              1.372510e+00                              0.0e+00    OK    
internal inductance li1                       8.842273e-01                              8.842273e-01                              0.0e+00    OK    
internal inductance li2                       7.080758e-01                              7.080758e-01                              0.0e+00    OK    
internal inductance li3                       7.304338e-01                              7.304338e-01                              0.0e+00    OK    
poloidal beta betap1                          6.680739e-01                              6.680739e-01                              0.0e+00    OK    
poloidal beta betap2                          5.349835e-01                              5.349835e-01                              0.0e+00    OK    
poloidal beta betap3                          5.518760e-01                              5.518760e-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.236468e-01                              4.236468e-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.006572e+00                              2.006572e+00                              0.0e+00    OK    
wall field bwall                              3.880150e-01                              3.880150e-01                              0.0e+00    OK    
aspect ratio                                  2.845751e+00                              2.845751e+00                              0.0e+00    OK    
elongation kappa                              1.708325e+00                              1.708325e+00                              0.0e+00    OK    
q profile (checksum)                          0446dee197d6...                           0446dee197d6...                           identical  OK    
pressure profile (checksum)                   a1e4037665c3...                           a1e4037665c3...                           identical  OK    
Mercier D_I profile (checksum)                558bac8c1f1e...                           558bac8c1f1e...                           identical  OK    
resistive interchange D_R profile (checksum)  781547f9a331...                           781547f9a331...                           identical  OK    
ballooning Delta' profile (checksum)          dddc351e6b37...                           dddc351e6b37...                           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.287814e-04                              5.287814e-04                              0.0e+00    OK    
dominant-coupling singular values             [3 elem]                                  [3 elem]                                  0.0e+00    OK    
|forcing overlap with dominant mode|          1.411600e-04                              1.411600e-04                              0.0e+00    OK    
|delta_nominal| of coil set 1                 7.034885e-05                              7.034885e-05                              0.0e+00    OK    
||ddelta/d(shift)|| over coil sets            1.245468e-04                              1.245468e-04                              0.0e+00    OK    
||ddelta/d(tilt)|| over coil sets             4.176658e-07                              4.176658e-07                              0.0e+00    OK    
PE plasma energy                              3.320581e+00                              3.320581e+00                              0.0e+00    OK    
PE vacuum energy                              3.174510e+00                              3.174510e+00                              0.0e+00    OK    
PE surface energy                             5.507877e+00                              5.507877e+00                              0.0e+00    OK    
PE toroidal torque                            1.119008e-04                              1.119008e-04                              0.0e+00    OK    
NTV torque FGAR [N·m]                         5.499613e-01                              5.499613e-01                              0.0e+00    OK    
NTV kinetic energy dW FGAR [J]                8.231914e-02                              8.231914e-02                              0.0e+00    OK    
Runtime (s)                                   397.9s                                    393.5s                                               --    
||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

Notes for reviewers

Two assertions gate the input: profiles_1d.psi must increase strictly and lie between psi_axis and psi_boundary. Grids that break these assumptions stop with a clear message instead of producing a wrong spline.

🤖 Generated with Claude Code

[logan-nc]

@claude review this

[claude]

Review of #499

The fix is correct and minimal: normalizing profiles_1d.psi by the axis/boundary values is what the spline abscissa should have been. The sign handling is fine, since (psi_1d - psi_axis)/(psi_boundary - psi_axis) is orientation-independent. The regression harness result is included and the PR has an assignee and reviewers.

Points to consider:

1. Clamp can create duplicate knots. If profiles_1d.psi extends past the boundary or before the axis (some IMAS writers include a padded or slightly overshooting last point), the clamp maps several points to exactly 0 or 1. cubic_interp then receives non-strictly-increasing sq_xs, which will error or give garbage. The comment says the clamp is only meant for roundoff. Either clamp only the endpoints, or add an explicit check, e.g. issorted(psi_norm_grid; lt=<=) || error("read_imas: profiles_1d.psi must be strictly monotonic between psi_axis and psi_boundary"). The same check would also catch descending grids that don't match the axis→boundary orientation. Currently nothing validates monotonicity.

2. Test coverage is narrower than the release note. The note says the fix also covers grids that "do not start at psi_axis", but the test grid starts exactly at 0. Add a case with psi_n = range(0.05, 1.0; ...) to exercise the ExtendExtrap path near the axis and the sqrt column. A test for the duplicate-knot or monotonicity case from point 1 would also be useful.

3. Test tolerance. rtol=1e-3 on a linear pressure profile is fine, but only buf[2] is checked. A one-line check on buf[4] (≈ sqrt(0.5)) would cover the sqrt(psi_norm_grid) column that the clamp comment refers to.

Minor: the existing comment "Extract 1D profiles, converting psi..." is still accurate, but a short note that the grid is normalized by the axis/boundary psi would help readers, since the previous code looked equally plausible.

Per CLAUDE.md: this PR must not be merged into develop without approval from a human reviewer (jhalpern30 / logan-nc are requested).

[logan-nc]

@jhalpern30 do you think this really needs a test?

[jhalpern30]

@jhalpern30 do you think this really needs a test?

Yeah I don't think this one needs a unit test. I do however like the suggestion of adding informative assert statements around this change. Its hard to be completely agnostic for all input types, and I'd rather it error out in a useful way and we can then modify it if a new equilibrium type or input comes along. For example, adding in the basic assumptions as assertions here (monotonic, between 0 to 1, anything else) is probably the way to go.

So @calin1989, my suggestions here are then to

1. remove the unit tests
2. add in the assert statements to gate the inputs, and consider if the clamp statement is really needed
3. investigate Nik's comment

Then you can ping me again and it should be good to go

[jhalpern30]

@logan-nc can you confirm if this change fixes the bug you noted in #482? it looks like that was a local SPARC eq so I can't reproduce it

[jhalpern30]

@logan-nc assuming this fixes the issue you had with the SPARC equilibria and will merge this later today

Everything from original review been addressed, and PR has my approval already - only thing holding up the merge

[logan-nc]

TY @calin1989 @jhalpern30