ForceFreeStates - FEATURE! - Cap the integration domain where resistive layers overlap

[amlakbekaluachule]

Release note

• Audience: users
• Numerical impact: one shipped deck moves, diiid_slayer_n1, from two input changes: Zeff(ψ) is now derived from the kinetic file's n_i (Zeff 1 → about 1.9 at its surfaces), and the shipped DIII-D-like kinetic file's χ⊥ is doubled. Growth rates go 2/1 194.1 → 224.8 s⁻¹ (+16%), 3/1 −247.5 → −283.1 s⁻¹ (+14%), 4/1 −1000.8 → −1116.7 s⁻¹ (+12%); Lundquist S falls 22–31%, P_perp rises 38–57%, D_norm, tauk, Q_e, Q_i and Q_root move 8–12%, ω ≤ 0.01%. The other 227 tracked quantities across the 16 cases shared with develop are unchanged. Results also change, on no shipped deck, for runs with set_psilim_via_dmlim = true and a finite qhigh, where dmlim could previously step past qhigh (qhigh = 4.0 gave qlim = 4.2, now 3.2) (harness @ 484a95a)
• Migration: SLAYER results change for any deck whose kinetic file carries n_i: Zeff is now derived from it (new [SLAYER] zimp, default 6) instead of the scalar [SLAYER] zeff, which remains only as the fallback when the file has neither n_i nor a zeff dataset. To prescribe Zeff, add a zeff dataset to the kinetic file. A deck with set_psilim_via_dmlim = true and a finite qhigh truncates one rational lower whenever qhigh sits on a rational or less than dmlim/n above one; raise qhigh above (rational + dmlim)/n to keep the old domain. LayerWidths gains delta_norm, delta_visco and delta_dr, so code constructing it positionally must supply them. build_slayer_inputs and build_ggj_inputs accept zeff as a scalar or a callable of ψ. The shipped DIII-D-like kinetic file's chi_e is doubled.

GPEC can now locate the point past which adjacent rational surfaces' resistive layers overlap — the outer limit of the matched-asymptotic treatment, after Fitzpatrick, Nucl. Fusion 2025, Sect. 5.9 — and, when psilim_from_layer_overlap = true, use it as an upper bound on the stability integration domain. psihigh stays the upper limit: the cap can only narrow the domain. The SLAYER layer and the overlap scan now take Zeff(ψ) from the kinetic profiles.

What is wired

piece	where
opt-in control	[ForceFreeStates] psilim_from_layer_overlap (default false, off on every shipped deck except the truncation example)
the scan	Tearing.resistive_layer_overlap, run whenever [SLAYER] names a readable profile_file, once per n; the innermost overlap point binds
the criterion	surface k overlaps k−1 when ψ_k − w_k/2 < ψ_{k−1} + w_{k−1}/2, with w the closed-form Eq. (100) width (no Riccati solve); the cut is at ψ_{k−1}
the cap	sing_lim!(...; psilim_cap=...): an upper bound on qlim, applied before dmlim, which cannot step past it
the record	ForceFreeStates/LayerOverlap/, written whenever the scan ran
the log	a warning whenever psilim lands past the overlap point (flag on or off), and the scan's own caveats and notes

The scan uses the run's own [SLAYER] plasma: profiles, χ⊥(ψ), Zeff(ψ), bt, mu_i, resistivity model and Coulomb-logarithm form, through SLAYER's own loader.

New since the last review

1. Surfaces beyond the grid are placed on the X-point separatrix law, not a fit (1ec489c6b, !).
The fitted edge law assumed q diverges at ψ_N = 1. A reconstruction may put its boundary slightly inside the X-point flux — ψ_N = 1.0001 on the shipped DIII-D-like file — so extrapolated surfaces were misplaced, and one could be dropped when the fit and the traced edge q disagreed across an integer. Equilibrium.separatrix_q_law takes the law from the equilibrium's own 2D flux:

q(ψ) = q_anchor + Σ_i A_i ln[(ψ_Xi − ψ_anchor) / (ψ_Xi − ψ)],    A_i = F / (2π R_Xi sqrt(−det H_i))

• Anchored on the traced q at the outermost grid surface; nothing is fitted.
• Saddles are found by following each poloidal-field minimum on the outermost traced flux surface to the critical point it leads to, with one term per saddle reached. The nearest one sets where q diverges. On single-null EFIT reconstructions this also picks up the far secondary X-point (ψ_N 1.07–1.13), and on one of five checked a null inboard of the plasma (ψ_N 1.57). These extra terms change the law's shear by at most 5% and do not move the divergence.
• Rejected when its shear differs from the traced one by more than ×1.5 over the outer four knots; not built for inverse equilibria (they carry their own q) or when no saddle is found. A limited plasma with an X-point nearby is not rejected: on the one limited EFIT reconstruction checked, the law was accepted with the X-point at ψ_N = 1.093 (shear ratio 0.86). The scan then places one rational past the boundary, where no closed surface exists, and warns that an overlap could not be established; the domain is unchanged. Telling a limited plasma from a padded file needs the limiter contour, which GPEC does not use. The docstring and the scan's message still say "limited plasma" for the no-law case and will be corrected with the next source change.
• The fitted edge_q_law remains only as the grid-refinement edge-density gate.

2. psihigh is the upper limit, with a graceful fallback. At most one rational past ψ_N = 1 is placed, only to judge whether the outermost surface inside is clean. The recommendation never exceeds ψ_N = 1, and when no overlap can be established the domain keeps its requested limit and the run says so.

3. Zeff(ψ) from the kinetic profiles (75457a1b3, !). The scan had used Zeff = 1, deuterium and default closures whatever [SLAYER] said; SLAYER itself used a scalar Zeff. Both now use one profile, in this order: a zeff dataset in the kinetic file; else zimp − (n_i/n_e)(zimp − 1) with the new [SLAYER] zimp; else the scalar [SLAYER] zeff. The SLAYER stage logs which. This is what moves diiid_slayer_n1 (Zeff 1 → about 1.9 on the shipped profiles). It is independent of truncation and wants a SLAYER reviewer's eye.

4. The scan's messages are now shown. They were collected but never logged. Messages that qualify the verdict (caveats: an overlap that could not be established, a deciding surface placed on the law, a backstop reached) are logged as warnings; per-surface detail (notes) as info.

Verified end to end

examples/DIIID-like_truncation_example (DIII-D-like equilibrium, psihigh = 0.9999, cap on; about two minutes):

Layer-overlap scan (n = 1): q = 9/1 lies beyond the equilibrium's own boundary (ψ = 1), inside its X-point flux ψ = 1.000100; it is placed only to judge whether the outermost surface inside is clean, and needs kinetic profiles past ψ = 1 to be scored
Resistive-layer overlap caps the domain: qlim 8.352 -> 7.000 (psi 0.999900 -> 0.998508)
Setting psilim via dmlim: initial qlim = 7.000, dmlim = 0.200
Filtered 2 singular surface(s) outside integration domain: [([7], 7.0), ([8], 8.0)]
qlim = 6.200, psilim = 0.995

The q = 7 and q = 8 layers are the first adjacent pair to overlap (sum of half-widths over separation = 1.24; the 6|7 pair is at 0.56), so the cap is the q = 7 surface and dmlim then stops 0.2 past q = 6. The resulting psilim = 0.99498 reproduces the ideal deck's hand-chosen psihigh = 0.995: Re(ep[1]) −1.37193 against −1.37187, Re(et[1]) 0.80082 against 0.80084. The harness case diiid_truncation_n1 pins this run, including the LayerOverlap outputs.

Re-run at this head with lower limits: at psihigh = 0.999 the cap still binds (qlim 7.247 → 7.000, then dmlim to 6.2); at 0.998 and 0.995 it is inert and qlim is 6.2 from dmlim alone.

Tests (runtests_layer_overlap_fitzpatrick.jl, 9 + 42 + 29 passing):

• Fitzpatrick's JET example: overlap at Ψ = 0.9985 for n = 1 and 0.9952 for n = 4.
• The X-point flux (ψ_N = 1.0001), coefficient and shear ratio on the DIII-D-like file; a law built on a grid stopping at ψ_N = 0.995 lands within 5% of where a longer grid traces q = 7 and q = 8; the scan's overlap point agrees to 5e-5 whichever grid is used.
• Zeff from each of the three sources; [SLAYER] settings reaching the scan.
• A cap on a rational, and 1000 ulps either side, gives the same qlim.

Limits

• Kinetic profiles end at ψ_N = 1. A rational between the boundary and the X-point flux cannot be scored, so an overlap between it and the outermost surface inside cannot be established; the domain then keeps psihigh and the run warns. This arises only for g-files written with boundary padding, which puts ψ_N = 1 inside the X-point flux: TokaMaker's save_eqdsk with lcfs_pad and the default truncate_eq=True cuts the equilibrium at 1 − pad and renormalises, so the X-point lands at ψ_N = 1/(1 − pad) (the shipped example is padded by 1e-4). On real DIII-D cases the first overlapping pair sits 2e-5 to 2e-4 inside the separatrix, so a truncating pad of 1e-3 hides it. EFIT-written files are not affected: 34 of 34 diverted ones checked have the X-point at ψ_N = 1 to within 6e-8. Writing with truncate_eq=False avoids it: the flux map, F and p then reach the separatrix, and only the last entry of the file's q column is a stand-in (q at 1 − pad), which GPEC does not read because it integrates q from the flux map.
• Untested: disconnected double null (no file to hand). Not modelled: snowflake and x-point-target divertors; the shear check is the backstop.
• χ⊥ is a proxy. Eq. (100) wants the one-fluid energy diffusivity; the kinetic file's chi_e is the electron power-balance value. The shipped file's χ⊥ is doubled so the truncation deck sits clear of threshold (a test-enabling choice, recorded in the file's provenance attribute).
• The NTV path does not read the new zeff dataset, and the X-point flux and per-surface Zeff are not written to gpec.h5. Both are left for follow-up issues.
• TOML entry point only: the flag is inert, with a warning, on the programmatic solve(prob, alg) path.
• The two width channels (Riccati delta_s and Eq. 100) disagree on where overlap starts. Both are written; only Eq. (100) sets the domain, following the paper.

Regression report

regress --refs 44e0afcc4,484a95ae2 --force — develop as merged into the branch against this head (the merge commit), all 18 cases; julia 1.11.6, CI-pinned manifest, 8 threads:

diiid_error_field              Summary: 10 unchanged
diiid_multi_n                  Summary: 19 unchanged
diiid_n1                       Summary: 53 unchanged
diiid_n1_riccati               Summary: 17 unchanged
diiid_n1_riccati_precompiled   Summary: 17 unchanged
diiid_slayer_n1                Summary: 9 changed, 8 unchanged
                               SLAYER Lundquist S           3.367e+07 (31.10%)  ** CHANGED **
                               SLAYER D_norm                2.465e-01 (10.56%)  ** CHANGED **
                               SLAYER P_perp                2.214e+01 (44.02%)  ** CHANGED **
                               SLAYER tauk                  7.006e-06 (10.56%)  ** CHANGED **
                               SLAYER Q_e                   1.431e-01 (7.79%)   ** CHANGED **
                               SLAYER Q_i                   1.120e-01 (8.62%)   ** CHANGED **
                               SLAYER Q_root [2/1,3/1,4/1]  7.380e-02 (10.57%)  ** CHANGED **
                               SLAYER ω_Hz [2/1,3/1,4/1]    1.326e+00 (0.01%)   ** CHANGED **
                               SLAYER γ_Hz [2/1,3/1,4/1]    1.159e+02 (11.58%)  ** CHANGED **
diiid_truncation_n1            new on this branch; no develop baseline (see below)
efit_fixedbdy_separatrix       Summary: 5 unchanged
gal_resistive_diiid            Summary: 10 unchanged
gal_resistive_pe               Summary: 8 missing/N/A   (on both refs)
ggj_ray_q500i                  Summary: 4 unchanged
ggj_reference                  Summary: 4 unchanged
solovev_kinetic_calculated     Summary: 17 unchanged
solovev_kinetic_multiion       Summary: 6 unchanged
solovev_kinetic_ntv            Summary: 6 unchanged
solovev_kinetic_nuzero         Summary: 14 unchanged
solovev_multi_n                Summary: 15 unchanged
solovev_n1                     Summary: 22 unchanged

236 quantities compared, 227 unchanged. The nine that move are all diiid_slayer_n1. Unchanged there: the surfaces, r_s, the shear and iota_e, which do not depend on Zeff or χ⊥. The two causes pull the growth rates in opposite directions: measured separately on this branch before the develop sync, the doubled χ⊥ alone lowered the 2/1 rate by 11% and the Zeff profile alone raised it by 30%.

Every other case is unchanged, including diiid_n1, whose NTV reads neither. gal_resistive_pe extracts nothing on either ref, develop included, so that predates this PR. Growth rates are e-folding rates in s⁻¹; the harness label γ_Hz is a misnomer. The harness exits non-zero only because diiid_truncation_n1's example directory does not exist on develop.

diiid_truncation_n1 is new, so it has values at this head only:

layer-overlap cap applied              1
layer-overlap psilim                   9.985081e-01
first overlapping surface index        7
diffusive-resistive layer widths (Δψ)  [7 elem]
qlim                                   6.200000e+00
psilim                                 9.949767e-01
mpert                                  37
# singular surfaces                    5
singular q values                      [5 elem]
singular psi locations                 [5 elem]
plasma energy Re(ep[1])               -1.371927e+00
total energy Re(et[1])                 8.008153e-01

The width array has 7 elements (q = 2–8) where the previous head had 15: surfaces past the first rational beyond ψ_N = 1 are no longer placed. Runtimes are not quoted; the run was not made under controlled load.

Notes for reviewers

• The scan is driven at rs_method = :flux, not the SLAYER default. Eq. (100) is not covariant, so it is valid only in the toroidal-flux label of the paper's Eq. (30).
• Acceptance thresholds introduced here, for explicit sign-off: the ×1.5 shear bound on the separatrix law (measured ratios are 0.93–1.08 on five EFIT reconstructions and 0.96–1.01 on the 35 other equilibria checked; a law missing one of two equal X-points sits at 0.5), and the test tolerances of 5% on surface placement and 5e-5 on the overlap point.
• delta_s_abs appears in a third HDF5 group, by the one-name-per-quantity rule.

Earlier history of this PR (bugs found, review responses) — unchanged since the last review

• dmlim could step past qhigh and the cap (14d60b6b9, !): qhigh = 4.0 gave qlim = 4.2; it now steps back against min(qmax, qhigh, cap).
• The scan used the wrong χ (f5c60e779, !): every width was computed with χ = 1 while SLAYER used the file's χ(ψ).
• The criterion no longer depends on the Riccati solve (cd6dce306, !): a failed delta_s solve used to drop a surface and hide an overlap.
• A single failing surface aborted the scan; the empty-scan path threw; a catch-all hid code errors (6781b639f, c02a3c43d, 8cfbbee7e).
• Multi-n runs cap at the innermost overlap point across nn_low:nn_high.
• The domain-limit convention is ψ[k−1], the paper's own 1 − ε_c; ψ[k−1] + w/2 exceeded ψ = 1 on 4 of 10 equilibria checked.
• Edge-density gate for the auto grid (cf85db8c5, !): shear consistency within ×2 replaces a tuned r² floor.
• Critical-Δ citations corrected (9c1937bda): Connor, Ham, Hastie & Liu, PPCF 57, 065001 (2015).
• Jacobian independence (asked by @logan-nc): r_s invariant to 1e-16 and the Δψ widths to ~7e-6 under pest / hamada / equal_arc.
• Raising the ideal and Riccati decks to psihigh = 0.9995 was declined: it adds the q = 7 surface and moves their results substantially (pe_toroidal_torque changes sign), and with the flag on the cap would pull them back to the domain they already have.

Validation

• Develop merged into this branch at 484a95ae2, no conflicts
• Full test suite at this head passes
• Docs build locally
• Regression harness at this head (above)

Important

Do not merge without a fresh third-party human review. Both earlier approvals predate the changes listed under "New since the last review".

🤖 Generated with Claude Code

[amlakbekaluachule]

@logan-nc this PR does two things:

1. Fixes a bug in slayer_layer_thickness (LayerThickness.jl)
2. Adds set_resistive_width_based_psihigh (new file,SetResistiveWidthPsihigh.jl)

Let me know if anything looks off.

[d-burg]

I synced this branch with feature/tearing-growthrates (merge 8a771555) and retargeted the PR base from develop to feature/tearing-growthrates. The diff now shows your 2 files instead of 65.

Worth knowing why the retarget was needed: LayerThickness.jl doesn't exist on develop at all — the whole SLAYER module is tearing-branch work, and develop only carries a Slayer.jl stub. So the PR couldn't have merged there regardless of how clean it was.

The merge hit two conflicts, both caused by src/Tearing/InnerLayer/SLAYER/ being relocated to src/InnerLayer/SLAYER/ on the parent branch (the InnerLayer module had to move above ForceFreeStates to break a dependency cycle). Nothing to do with your changes — resolved keeping your δ_FKR / δ_visco constructor arguments. The package precompiles and the SLAYER tests pass on the merged branch.

Two things worth a look before review:

1. SetResistiveWidthPsihigh.jl is never loaded

SLAYER.jl includes LayerParameters.jl, Riccati.jl, LayerThickness.jl and LayerInputs.jl (lines 51–54) but not SetResistiveWidthPsihigh.jl, so the helper isn't reachable:

julia> using GeneralizedPerturbedEquilibrium
julia> isdefined(GeneralizedPerturbedEquilibrium.InnerLayer.SLAYER, :set_resistive_width_based_psihigh)
false

Adding include("SetResistiveWidthPsihigh.jl") alongside the others should wire it in, plus an export if it's meant to be public API.

2. δ_FKR / δ_visco aren't in the LayerWidths docstring

The # Fields list in LayerThickness.jl (lines 133–139) still stops at d_beta. The repo convention is to document struct fields in the docstring rather than as inline comments, so these want a line each saying what they are and their units.

Minor: I also stripped some trailing whitespace from both files while resolving the merge — the pre-commit hygiene hook rejects it.

[logan-nc]

@d-burg whats the status of this? Should it be reviewed now? Or should it still be marked as draft?

[amlakbekaluachule]

Fixed both flagged issues — SetResistiveWidthPsihigh.jl is now included/exported, and δ_FKR/δ_visco are documented in the LayerWidths docstring.

[logan-nc]

@claude review this

[claude]

Review

Reminder (per repo policy): this PR must not be merged without third-party human review, regardless of any automated approval. Also — please run the regression harness (regress --cases <relevant_case> --refs develop,local) before merge; this PR touches LayerWidths/slayer_layer_thickness, which is exactly the kind of numerical-output change the harness is meant to catch.

Correctness / design concerns

1. Hardcoded example data file inside an exported library function. set_resistive_width_based_psihigh unconditionally loads examples/DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode_kinetic.h5 (SetResistiveWidthPsihigh.jl:39) regardless of the equilibrium/kinetic data actually carried by the caller's r. Since the function is exported from both SLAYER and InnerLayer, any user who calls it on a different equilibrium silently gets DIII-D-like kinetic profiles instead of their own, and the relative path only resolves when the working directory happens to be the repo root. build_slayer_inputs already accepts a profiles::KineticProfiles argument (LayerInputs.jl) — this function should take the same, rather than reading a fixed file. If this really is meant to be a one-off sanity-check script for that specific example (as the PR title suggests), it probably shouldn't be a permanent exported module API at all.

2. Reimplements existing rational-surface-finding instead of reusing it. find_psi_for_q brute-forces a 200,000-point grid of q_of_psi_real per call (called once per m in q_scan, i.e. up to ~39× → ~7.8M spline evaluations) followed by a hand-rolled 60-iteration bisection. src/ForceFreeStates/Sing.jl already has _find_rational_surfaces/rational_psi_nodes, which locate q=m/n surfaces via Roots.jl's find_zero(..., Roots.Brent()) directly on profiles.q_spline, segmented between q-extrema (so it's also reverse-shear-safe, which the new bisection is not — it just finds the first sign change on a dense scan). Roots.jl is already a project dependency. Per CLAUDE.md's minimal-change-discipline section, this should reuse rational_psi_nodes/_find_rational_surfaces instead of re-deriving an independent, less robust q(ψ) inversion.

3. omega_star_at (SetResistiveWidthPsihigh.jl:59-67) duplicates _omega_star_at in LayerInputs.jl almost verbatim. Same formula, same sign convention, same structure — this is a second copy of logic that already exists and is already exercised by build_slayer_inputs.

4. Extrapolation edge case: A = (q2 - q1) / (log(1-ψ1) - log(1-ψ2)) (SetResistiveWidthPsihigh.jl:9) divides by a log-difference that goes to -Inf if x_knots[end] == 1.0 exactly, making A = 0; the target_q > q_last branch then computes exp((q_last-target_q)/A) → division by zero. In practice psihigh < 1 is typical, but there's no guard or comment noting this assumption, and this is precisely the kind of degenerate input this function is meant to help users avoid — worth an explicit check/error message.

Physics traceability

5. δ_FKR = r_s · S^(-1/3) is labeled "FKR resistive-inertial length scale" (LayerThickness.jl docstring), but the classic Furth-Killeen-Rosenbluth (1963) resistive-inertial ("constant-ψ") tearing-layer width scaling is usually quoted as δ ∝ r_s S^{-2/5}, not S^{-1/3} (the S^{-1/3} scaling is closer to the Δ'→∞ / resistive-internal-kink limit). No equation citation is given, and neither docs/development/references.md nor docs/resources/ currently lists an FKR reference — only the 2026 Burgess two-fluid-layer paper is present. Given CLAUDE.md's explicit encouragement to cite the source equation for physics kernels, please add the citation (equation number) this formula comes from, and double check the exponent — this looks worth a pass from the fortran-physics-reviewer agent before merge.

Style / consistency

6. New struct fields δ_FKR, δ_visco (LayerThickness.jl:155-156) use unicode Greek + mixed casing, while every other field in LayerWidths/SLAYERParameters uses ASCII delta_*/d_* snake_case (delta_s, delta_s_m, d_beta). This is inconsistent with the rest of the struct and forces awkward unicode access like lw.δ_visco in the new call site — recommend delta_fkr/delta_visco to match convention.
7. Stray blank line before end in the LayerWidths struct definition (between the new fields and end) — minor formatting nit.
8. q1_val (SetResistiveWidthPsihigh.jl:63) is bound to dqdpsi_hybrid(ψ_m) (a derivative), reusing the name shortly after q1, q2 were bound to actual q-values at ψ1, ψ2 a few lines above — easy to misread; consider dqdpsi_val or similar.
9. Parameter name r for the response/model object (set_resistive_width_based_psihigh(r; ...)) is easy to confuse with the various rs/rs_val (minor radius) variables used throughout the same function.

Documentation / tests

10. Both new exported functions (set_resistive_width_based_psihigh, print_layer_width_table) have no docstrings. Per CLAUDE.md, Documenter CI fails on missing_docs for any exported symbol without a docstring — this should fail the docs build as-is.
11. No test coverage was added for the two new LayerWidths fields or for the new 158-line function (diff touches only src/, no files under test/).

None of the above touch the pre-existing slayer_layer_thickness/LayerWidths core computation beyond the two new fields, so the blast radius of a fix should be contained to the new file plus the docstring/struct nits.

[logan-nc]

Status

Rewritten from the original 172-line contribution, which had zero call sites and compared layer widths in metres against positions in normalized flux. That file is gone.

Feature — [Equilibrium] psihigh_from_layer_overlap, opt-in, default off. Forms a first-pass equilibrium, finds where adjacent rational surfaces' resistive layers overlap, checks the reconstruction can actually deliver the resulting domain, and re-forms. Floored by the rational-q window containing psiedge so the dW scan band stays open. The scan runs whenever kinetic profiles are readable regardless of the flag, so gpec.h5 always records Equilibrium/psihigh_resistive_layers beside psihigh_resolved, plus per-surface widths under Equilibrium/ResistiveLayers/.

Default is off because enabling it crosses a rational surface on a diverted plasma: on DIIID-like_ideal_example the domain goes 0.995 → 0.993281, 5 surfaces → 4, mpert 35 → 34, moving et[1] 76% and NTV torque ~1800%. That is the criterion working, not a defect, but not a change to make on every user's behalf silently.

Nine defects fixed along the way, all found by trying to make the feature trustworthy:

defect	consequence
field-line ODE retcode never checked	truncated solves consumed as closed flux surfaces — silently wrong equilibria
round-trip check sampled only at ψ knots	structurally blind to inter-knot ringing; passed while et[1] reached −47707
surface_da_dpsi stencil clamped at 1−1e-4	the real limit on near-boundary accuracy, not psi_cap
run_slayer passed b0exp as B_T	Lundquist number off by ~95% — split out as #399
trial formations unguarded	crash instead of back-off, psihigh left corrupted
round-trip unavailable on the efit path	the path every shipped deck uses
gpec.h5 replay crash	stored relative profile paths don't resolve in the output dir
edge q-law duplicated, applicable to limited plasmas	one shared helper, with a guard that actually discriminates
inverse equilibria extrapolated	fabricated 3 non-existent surfaces on the CHEASE fixture

Verification — 122 + 67 + 59 + 32 + 14 tests; docs clean; the metres→flux conversion triple-checked three independent ways (the binding pair is 1.821 mm apart with half-widths summing to 2.338 mm, argued in pure metres); the back-off gate validated against ForceFreeStates energies as ground truth, and it flips one step before they degrade.

Still open

• Final frozen six-case regression. diiid_slayer_n1 needs re-pinning — for the InnerLayer.SLAYER - BUGFIX! - Use the physical toroidal field, not the b0exp normalization #399 bt fix, not for this feature.
• Full suite against the current tree; the IMAS deck has never been run with the flag on (left off, with a comment explaining why).

[logan-nc]

@d-burg — ready for review, but please take #399 first.

#399 is the one-line bt fix (SLAYER was running at b0exp = 1.0 T instead of the physical ~1.95 T, so the Lundquist number was off by ~95%). It is cherry-picked here as 4304caa, so once #399 merges this branch's copy resolves cleanly and the diff shrinks to this PR's own work.

That ordering also separates the regression movement. diiid_slayer_n1 needs re-pinning, and with #399 landed first it is unambiguous which part is the bt correction and which is this feature — measured separately, #399 alone accounts for S +94.5%, D_norm 35%, γ 4.4%.

Worth your attention specifically:

• The psiedge floor governs every DIII-D deck, so the layer criterion never actually sets the domain in the shipped examples — each asks for ~0.9892 and is held at the q=6 window top. Working as specified, but you may want a deck that exercises the physics path.
• The bt fix changes the layer widths themselves, so any Tearing/LayerWidths/ numbers you have from before InnerLayer.SLAYER - BUGFIX! - Use the physical toroidal field, not the b0exp normalization #399 were computed at half the true field.
• delta_s_abs now appears in two h5 groups. They agree to ~1e-6 and the difference in scope is documented in hdf5-conventions.md; flag it if that still reads as ambiguous.

Everything except the opt-in flag is a bug fix that stands on its own — the field-line ODE retcode was never checked, and the round-trip gate was structurally blind to spline ringing (it passed while et[1] hit −47707).

[logan-nc]

@d-burg, as claude says this has hit a point where it needs your expert review. Workflows are in place to decide the truncation from the resonant layer widths. The issue is that it is reporting big widths and wanting to truncate way in (~0.9892 for DIII-D). I find this suspect, since Fitzpatrick had been finding things like to be very far out past 0.9999. In the process of adding the truncation logic, claude did convince me the units of the layer widths had been wrong in the original overlap-checking-loop (comparing layer widths of m to rational spacing in psin) and a conversion from m to psin was added... so maybe double check that?

Do you think we need to try and reproduce the exact case from Fitzpatricks paper on this as a sanity check?

[logan-nc]

@d-burg I was excited to see the commits roll in! Whats the status now? Still working through any todos?

[d-burg]

Update at 642caffee — what changed, where the cut fires, and how far to trust it

Both approvals are at 6a75d0f50. Results-changing commits have landed since, so this needs a fresh review, not a re-stamp. The PR description has been rewritten for this head; this note is the short version.

Changed since the last update

Change	Commit	Why
Surfaces beyond the grid sit on the X-point separatrix law, taken from the 2D flux	1ec489c6b !	The fitted law assumed q diverges at ψ_N = 1. The X-point is at ψ_N = 1.0001 on the shipped file and further out on some reconstructions, which misplaced surfaces and could drop one.
psihigh is always the upper limit	1ec489c6b	Only one rational past ψ_N = 1 is placed, to judge the last surface inside. If no overlap can be established the domain keeps its limit and the run warns.
Zeff(ψ) from the kinetic file, for SLAYER and the scan	75457a1b3 !	The scan used Zeff = 1 and default closures whatever [SLAYER] said. New [SLAYER] zimp; a zeff dataset in the file overrides. Moves diiid_slayer_n1; independent of truncation.
Scan caveats and notes are logged	1ec489c6b	They were collected and never shown.

Earlier changes (the dmlim bound, the closed-form criterion, the χ⊥ fix, the truncation deck) are unchanged and summarised in the description.

Is it firing?

Yes, on examples/DIIID-like_truncation_example (psihigh = 0.9999, flag on):

Layer-overlap scan (n = 1): q = 9/1 lies beyond the equilibrium's own boundary (ψ = 1), inside its X-point flux ψ = 1.000100; it is placed only to judge whether the outermost surface inside is clean, and needs kinetic profiles past ψ = 1 to be scored
Resistive-layer overlap caps the domain: qlim 8.352 -> 7.000 (psi 0.999900 -> 0.998508)
Setting psilim via dmlim: initial qlim = 7.000, dmlim = 0.200
Filtered 2 singular surface(s) outside integration domain: [([7], 7.0), ([8], 8.0)]
qlim = 6.200, psilim = 0.995

Where. The q = 7 and q = 8 layers are the first adjacent pair to overlap, so the cap is the q = 7 surface, and dmlim then steps inside it to qlim = 6.2. Five rational surfaces remain (q = 2–6) of seven.

How. Two adjacent layers overlap when the sum of their half-widths exceeds their separation in ψ. That ratio is 0.56 for the 6|7 pair and 1.24 for 7|8.

On real equilibria

Run outside the repository on DIII-D H-mode discharges, n = 1, χ⊥ = 4 m²/s, Zeff(ψ) from the measured profiles. Two sets of equilibria for the same discharges, which differ in where the boundary sits relative to the X-point.

Set 1: boundary within 2e-5 of the X-point flux (three discharges, two times each). The separatrix law applies on all six (shear ratio 0.987–1.001).

case	first overlapping pair	ratio	overlap point ψ_N	qlim at psihigh 0.998	at 0.999	at 0.9999
A	5|6	1.14	0.99847	4.2 (unchanged)	5.2 → 4.2	6.2 → 4.2
B	5|6	1.30	0.99876	4.2 (unchanged)	4.2 (unchanged)	5.2 → 4.2
C	5|6	1.10	0.99862	4.2 (unchanged)	4.2 (unchanged)	5.2 → 4.2
D	5|6	1.25	0.99880	4.2 (unchanged)	4.2 (unchanged)	5.2 → 4.2
E	6|7	2.59	0.99969	4.2 (unchanged)	5.2 (unchanged)	6.2 → 5.2
F	6|7	2.56	0.99970	4.2 (unchanged)	5.2 (unchanged)	6.2 → 5.2

• The overlap point is the same to five digits whichever psihigh the grid stops at.
• At psihigh = 0.998 the cap never binds: dmlim already stops inside it. It binds on one case at 0.999 and on all six at 0.9999.

Set 2: boundary at 99.9% of the X-point flux (the same three discharges and one more, eight slices). q reaches only 4.7–5.5 at ψ_N = 1, and no pair inside overlaps for any χ⊥ from 1 to 64 m²/s (closest ratio 0.48 at χ⊥ = 4). The next rational lies between the boundary and the X-point, where there are no profiles, so it cannot be scored: the scan reports that no overlap was established and the domain keeps psihigh. This is the fallback path, exercised by real data.

So the shipped example fires, and Set 2 does not, because of where each equilibrium's boundary sits, not because of the plasma.

Sensitivity to χ⊥ and Zeff (Set 1):

• Widths go as χ⊥^(1/4), so the result is quantised. The 5|6 pair starts to overlap at χ⊥ ≈ 1.4–2.7 m²/s on A–D and 4.9 on E, F; every case reaches qlim = 4.2 by χ⊥ = 16. A factor of two in χ⊥ moves the cut by at most one rational.
• Zeff = 1 instead of the measured profile lowers the ratios by about 20% and changes the answer on two of six cases (A: 5|6 → 6|7; C: no overlap established). The measured Zeff matters.
• Multiplying χ_e by two as a stand-in for χ_e + χ_i is not recommended: the layer width wants the pressure-gradient-weighted mean of the two channels, not their sum.

How confident are we?

Solid

• Off by default, and safe when off: the scan then only writes diagnostics and warnings.
• Can only narrow the domain: psihigh, qhigh and dmlim still apply inside the cap.
• Matches the paper: Fitzpatrick's JET overlap points, ψ = 0.9985 (n = 1) and 0.9952 (n = 4).
• Independent of the grid: the overlap point does not move with psihigh, in the repository's test and on the six real cases.
• Fails towards "no cap", visibly: unreadable profiles, an unscorable surface, a limited edge or a law that fails its shear check each leave the domain alone and say so in the log.

Open

• Rationals beyond the reconstruction boundary cannot be scored without profiles there (Set 2).
• χ⊥ is a proxy (electron power-balance χ_e for the one-fluid χ_E), and the shipped file's is doubled so the test deck sits clear of threshold.
• Untested: disconnected double null. Not modelled: snowflake and x-point-target divertors.
• The two width channels disagree on where overlap starts; only Eq. (100) sets the domain.
• The NTV path does not yet read the new zeff dataset.

Bottom line: ready as an opt-in diagnostic and optional cap. Not ready to be on by default.

🤖 Generated with Claude Code

[logan-nc]

@d-burg great figures. Good Zeff catch... whats the status in human words? close???

[d-burg]

@logan-nc yes of course, just wanted to lay everything out for the record. The "Set 2: boundary at 99.9% of the X-point flux" point is the most concerning to me and I finally just got to the bottom of it. As you can see, in the default example TokaMaker g-file, every rational above the boundary q (8.79) lies at psi_N > 1 (these are closed surfaces inside the separatrix that the file drops)...

... and another view here:

The origin of this is the truncate_eq=True flag in the Tokamaker save_eqdsk, which I've always just left on by default. Looking into the mechanism, it cuts the equilibrium at 1 − lcfs_pad (1e-4 here) and renormalises psi_N to that surface, so the X-point sits at psi_N = 1/(1 − pad) and everything between is dropped from the file.

psi_pad is always active on the last node of the geqdsk saved q profile out of Tokamaker, but GPEC only looks at the other non-q profiles, so it seems like truncate_eq=False will give us more accurate bookkeeping. Anyone who wants this feature to fire or run more reliably on a Tokamaker equilibrium will need to turn that off. EFIT geqdsks are not affected since they aren't truncated/padded in this way and their psi_N=1.0 is consistent.

Should I regenerate the D3D-like example without the profile padding, or just file an issue so this merges asap? With my chi_e boost + Zeff fix + small non-default pad of 1e-4 (vs the default 1e-2), the feature does fire on this truncate_eq=True equilibrium even though it doesn't fire on most other Tokamaker equilibria I've tried.

Other than this issue:

• Evan's InnerLayer.SLAYER - FEATURE! - Add perpendicular Prandtl number options and an ion heat diffusivity profile #501 will be helpful for using the correct perpendicular diffusivity D_perp instead of the existing proxy electron thermal chi_perp_e, though for that we do need a chi_i. Also, D_perp is a weighted mean that usually lands lower than chi_perp_e, which means that this overlap criterion will shift even farther out. We could either stack this one on InnerLayer.SLAYER - FEATURE! - Add perpendicular Prandtl number options and an ion heat diffusivity profile #501 or extend InnerLayer.SLAYER - FEATURE! - Add perpendicular Prandtl number options and an ion heat diffusivity profile #501 to correct the overlap calculation once this merges
• A limited plasma with a nearby X-point gets a spurious "overlap not established" warning, may need future work

[logan-nc]

Wow. What a tricky nuance to find in the TokaMaker equilibrium. Good sleuthing @d-burg.

I'd say that since this currently fires on the equilibrium we have, lets merge it now. This PR is getting old and a bit unwieldy, so lets avoid more creep. Similarly, I vote we do not stack this on #501 - it's better to get this foundation of truncation logic in, then naturally update the calculations as the layer models improve. If the #501 updates push the truncation too far out, then we may need to think about updating the example equilibrium to have truncate_eq=False - but lets not invent problems for ourselves just yet.

You are the official assignee, so I'll let you do the honors. Thanks for getting this over the finish line. It's a really great feature!

[logan-nc]

Wooooohooo!!! 🎉

Thanks @d-burg!