Pressure lost due to a saturated island

[logan-nc]

This PR adds a simple estimate of the pressure lost when an RMP driven island saturates to the saturated width predicted by Callen's model.

This can be used to compare to EPED to predict if an RMP island is expected to decrease the pedestal pressure enough to stabilize ELMs. It can also be benchmarked by comparing to TM1 results.

Example output when I boost the DIII-D ideal example to 1e6 amp I-coils:

...
Reading table from file: 
   g147131.02300_DIIID_KEFIT.kin
 Computing total resonant fields
           psi            q      singflx singlfx_crit     chirikov     w_island          w_v     w_v_crit        w_sat        w_min         dP/P
     5.936E-01        2.000    1.809E-01    4.729E-04        7.825    1.760E+00    3.008E-01    1.223E-01    1.261E-01    3.311E-03    2.378E-01
     8.186E-01        3.000    1.180E-01    2.396E-04       10.102    1.107E+00    2.285E-01    1.024E-01    1.028E-01    3.632E-03    2.761E-01
     9.282E-01        4.000    2.597E-02    5.624E-04        6.817    4.108E-01    1.766E-01    6.900E-02    9.333E-02    2.168E-03    1.068E+00
     9.884E-01        5.000    2.483E-01    9.979E-04       10.039    6.050E-01    1.112E-01    3.471E-02    8.126E-02    9.568E-04    8.048E+00
   split cpu time =            3  seconds
   total cpu time =           25  seconds
 Computing Clebsch displacements
...

In other words, that huge of an RMP would be expected to clobber the pedestal completely.
I like printing dP/P because it allows the user to eyeball if say ~10% of the pressure is eaten by an island at the pedestal top.
The netcdf saves both dP and P separately.

A known issue with the implementation is that each island is currently treated as perturbing the original pressure profile. When islands overlap, this is an over estimate. We should really work our way from the edge inwards and modify the pressure profile successively.

[logan-nc]

@priyanshlunia there is some bug in that last commit but it should be clear what I am trying to do. Can you fix it and plot the modified pressure profile from saturated islands? Post a screenshot when you do! TY

[logan-nc]

any progress @priyanshlunia?

[logan-nc]

@priyanshlunia 😞

[priyanshlunia]

Went through the WIP commit (GPEC - WIP - Attempting to get a fully modified pressure profile output) and found three distinct bugs in the p_mod (modified pressure profile) code in gpout_singfld (gpec/gpout.f), plus a stray debug edit. All are now fixed. Here's the breakdown.

Background

The goal of the WIP commit is a full radial pressure profile p_mod(psi_n) that flattens the pressure gradient (dP/dpsi -> 0) inside each saturated island, then writes it to gpec_profile_output_n*.nc for plotting. The relevant target file is initialized in gpout_init_netcdf: its full-radial axis is the dimension psi_n, of size mstep, with coordinate values psifac(1:mstep) (the DCON integration grid). It also carries psi_n_rational (size msing) for the per-rational-surface quantities.

Bug 1 — grid/size mismatch (the crash)

p_mod was declared DIMENSION(0:mpsi) and evaluated on the equilibrium grid sq%xs(0:mpsi), but then written to the psi_n NetCDF variable, whose length is mstep on the psifac grid. These are different grids and different lengths — in the DIII-D example mpsi+1 = 129 vs mstep = 1548 — so nf90_put_var fails (and even if the sizes coincided, the coordinate mapping would be wrong).

Fix: declare p_mod as DIMENSION(mstep) and evaluate the interpolating spline on psifac(1:mstep) so it lands on the file's actual psi_n axis (both the Callen branch and the no-island ELSE branch):

!  before
REAL(r8), DIMENSION(0:mpsi) :: p_mod
...
DO ipsi=0,mpsi
   CALL spline_eval(spl,sq%xs(ipsi),0)
   p_mod(ipsi) = spl%f(1)
ENDDO

!  after
REAL(r8), DIMENSION(mstep) :: p_mod
...
DO ipsi=1,mstep
   CALL spline_eval(spl,psifac(ipsi),0)
   p_mod(ipsi) = spl%f(1)
ENDDO

Bug 2 — astat variable collision (silently dropped area_rational)

astat is set from the area_rational existence check and later consumed to decide whether to define and write area_rational:

astat = nf90_inq_varid(fncid, "area_rational", a_id)   ! ~line 2065
...
IF(astat/=nf90_noerr) THEN ... define area_rational ... ENDIF   ! ~line 2110
...
IF(astat/=nf90_noerr) THEN CALL check(nf90_put_var(fncid, a_id, area)) ENDIF   ! ~line 2152

The WIP commit inserted astat = nf90_inq_dimid(fncid, "psi_n", ps_id) between the define and the write, clobbering astat. Because psi_n always exists in this file, astat became nf90_noerr, so the IF(astat/=nf90_noerr) guard on the area write was skipped and area_rational was left as NetCDF fill values.

Fix: don't reuse astat. psi_n is defined unconditionally in gpout_init_netcdf, so no existence guard is needed — inquire, define, and write p_mod unconditionally, leaving the area_rational status untouched:

CALL check( nf90_inq_dimid(fncid, "psi_n", ps_id) )
CALL check( nf90_def_var(fncid, "p_mod", nf90_double, (/ps_id/), pm_id) )
CALL check( nf90_put_att(fncid, pm_id, "units", "Pa") )
CALL check( nf90_put_att(fncid, pm_id, "long_name",
$      "Modified pressure with saturated island correction") )
...
CALL check( nf90_put_var(fncid, pm_id, p_mod) )   ! unconditional

Bug 3 — sign error in the trapezoidal integration

The profile is reconstructed by integrating dP/dpsi inward from a fixed edge value. The recurrence added the (negative) gradient:

p_tmp(ipsi) = p_tmp(ipsi+1) + 0.5_r8 * (dp_tmp(ipsi+1) + dp_tmp(ipsi)) / REAL(npsi_tmp,r8)

Since dP/dpsi < 0, this makes the pressure decrease toward the core (and go negative) — backwards. It must subtract:

p_tmp(ipsi) = p_tmp(ipsi+1) - 0.5_r8 * (dp_tmp(ipsi+1) + dp_tmp(ipsi)) / REAL(npsi_tmp,r8)

Also — reverted a stray debug edit

The WIP commit had bumped the il/iu correction-coil currents in docs/examples/DIIID_ideal_example/coil.in (1e3 -> 1e5/1e6). That was scratch used to drive large islands; it's reverted to the baseline 1 kA so the shipped example is unchanged.

Verification

Built (gfortran + Homebrew NetCDF + Accelerate) and ran the DIII-D ideal example end-to-end (dcon then gpec). Automated checks on gpec_profile_output_n1.nc:

• len(p_mod) == len(psi_n) == 1548 (= mstep) — Bug 1
• p_mod finite, positive, monotonically non-increasing on the baseline (no-island) case — Bug 3
• area_rational = [34.8, 44.0, 49.2, 52.7] populated (not fill values) — Bug 2

With the coil currents raised enough to saturate islands, p_mod correctly flattens across each [psi_rational +/- hw_sat] region and drops below the equilibrium profile by the cumulative pressure lost, and dP_res/P_res per surface are consistent with that offset.

Note: for a saturated edge island (e.g. q=5) dP/P can exceed 1 under very strong drive — an artifact of the flat-top estimate near the plasma edge where P -> 0, not a code issue.

A run with the DIII-D ideal example is shown below:

[logan-nc]

@claude review this PR

[claude]

Claude finished @logan-nc's task in 8m 4s —— View job

---

Review: Pressure lost due to a saturated island

Reviewed gpec/gpout.f diff against origin/develop (base 7c08271, 5 commits, +112/-16). I also re-verified the three bugs @priyanshlunia previously found and fixed (grid/size mismatch, astat collision, integration sign) — all three fixes are correctly in place, and the debug coil.in edit is not present in this diff.

Bug: ASCII output format descriptor mismatch (gpec/gpout.f:1984)

WRITE(out_unit,'(1x,f6.3,16(es17.8e3))')
     $           singtype(ising)%q,singtype(ising)%psifac,
     $           ... (13 more items) ...
     $           hw_sat(ising),hw_min(ising),
     $           dP_r(ising),P_r(ising)

The item list now has 18 values (q through P_r(ising)), but the format only provides 17 slots: f6.3 (1) + 16(es17.8e3) (16) = 17. Before this PR the format matched exactly (16 items ↔ 15(es17.8e3)); adding dP_r/P_r needed the repeat count bumped by 2 (15→17), but it was only bumped by 1 (15→16).

Since the I/O list has more items than the format has descriptors, Fortran performs format reversion: it reverts to the last parenthesized repeat group (16(es17.8e3)) and starts a new record to output the remaining item. So P_r(ising) will be written on its own line, separate from the rest of the row, corrupting the column alignment of gpec_singfld_n*.out for every rational surface.

Fix: change 16(es17.8e3) → 17(es17.8e3).

(The corresponding header format '(1x,a6,18(1x,a16))' and the verbose stdout format '(1x,es13.3,f13.3,2es13.3,f13.3,6es13.3)' were both incremented correctly — only the ascii data-row format has the off-by-one.)

Fix this →

Minor / style

• gpec/gpout.f:1844,1846 — MAX(0.0, respsi - hw_sat(ising)) / MIN(1.0, respsi + hw_sat(ising)) use bare default-kind literals 0.0/1.0 mixed with REAL(r8) arguments. This is functionally harmless (0 and 1 are exact in any FP kind), but it's inconsistent with the identical clamp pattern a few lines later (gpec/gpout.f:1920,1922, MAX(0.0_r8, ...) / MIN(1.0_r8, ...)) and with the project's CLAUDE.md convention of always suffixing REAL(r8) literals with _r8.
• gpec/gpout.f:1849 — comment typo: "latetr" → "later".
• gpec/gpout.f:1844,1846 — these spline_eval(sq, ..., 1) calls request first-derivative evaluation (order 1) but only sq%f(2) (the value) is used; order 0 would suffice. Trivial cost, not blocking.

Notes (not bugs)

• The new full-profile p_mod output (written to gpec_profile_output_n*.nc) actually already handles overlapping saturated islands correctly — dp_tmp is zeroed across the union of [psifac-hw_sat, psifac+hw_sat] for every rational surface before the cumulative integration, so overlapping regions aren't double-flattened. The overlap caveat noted in the PR description ("each island is currently treated as perturbing the original pressure profile... an over estimate") still applies to the per-surface dP_r/P_r diagnostics (ascii + dP_res/P_res in the netCDF), since each is computed independently against the unmodified sq pressure spline — worth a one-line doc/comment clarifying that only the per-surface diagnostics carry that limitation, not p_mod.
• npsi_tmp = 1000 gives a fixed grid spacing of 0.001 in psi_n for the flattening/integration grid. For a saturated island narrower than that (possible for small hw_sat near a densely-packed set of surfaces), no grid point may fall inside [psileft, psiright], silently skipping the flattening for that surface. Not an issue for the DIII-D example shown, but worth keeping in mind for future edge cases.

I was not able to build/run the code in this sandbox (network and script-execution access is restricted here), so this review is based on static reading of the diff and surrounding code only — please compile and run the DIII-D example to confirm the ascii-format fix before merging.

[logan-nc]

@priyanshlunia please address the review comments above, including the minor ones which all seem like good changes to me.

So now there is a npsi_tmp and a mstep grid? Should we just use mstep throughout?