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

[amlakbekaluachule]

Release note

• Audience: users
• Numerical impact: one shipped deck moves (harness @ af9a8ba): the DIII-D-like kinetic profile's χ⊥ is doubled, which moves diiid_slayer_n1 — P_perp ×2 exactly, the 2/1 growth rate 194.6 → 173.6 Hz (−11%), the 3/1 and 4/1 rates by −9% and −8%, Q_root ≤ 0.3%, ω ≤ 0.1%. The other 216 tracked quantities across the 15 cases shared with develop are unchanged. Results also change, on no shipped deck, for (1) 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), and (2) real reconstructions on grid_type = "auto" whose edge the old q-law gate wrongly rejected.
• Migration: a deck with dmlim on and a finite qhigh now 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. sing_lim! gains an optional psilim_cap keyword and is otherwise unchanged. The layer-overlap scan reads its profiles only from [SLAYER] profile_file, so a deck without one no longer writes ForceFreeStates/LayerOverlap/. The shipped DIII-D-like kinetic file's chi_e is doubled, so SLAYER results computed from it differ from earlier versions.

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 psilim. Its layer widths use the same profiles and χ⊥(ψ) as the SLAYER analysis. The overlap point is recorded whenever the scan runs, and GPEC warns whenever psilim lands past it. The truncation logic is documented in the new "Integration domain truncation" section of the stability docs.

What is actually wired

piece	where
opt-in control	[ForceFreeStates] psilim_from_layer_overlap (default false)
the scan	Tearing.resistive_layer_overlap, run in prepare_force_free_states! whenever the run's [SLAYER] section names a readable profile_file, once per n in nn_low:nn_high; the innermost overlap point binds
the cap	sing_lim!(...; psilim_cap=...) — an upper bound on qlim, applied before dmlim, which can no longer step past it
the record	ForceFreeStates/LayerOverlap/ — always written when the scan ran
the warning	fires whenever psilim lands past the overlap point, flag on or off

The scan runs regardless of the flag; only applying the cap is gated, so a run always shows where layer physics would have cut the domain. Profiles, and χ⊥(ψ)/χ_φ(ψ), come from the [SLAYER] profile_file through SLAYER's own loader, so the scan and the SLAYER analysis see identical inputs.

Why it stays off by default, and why the DIII-D decks keep their psihigh

With the shipped DIII-D-like profiles the Eq. (100) overlap point is ψ ≈ 0.9985, at the q = 7 surface. That is outside the ideal and Riccati decks' psihigh = 0.995, so the cap could not bind there even if enabled, and inside the SLAYER deck's 0.9995, where enabling it would exclude q = 7. The flag is off on every shipped deck, so none of them changes; the SLAYER deck does now log the "past the overlap point" warning, because its dmlim domain ends near ψ = 0.9989.

Raising the ideal and Riccati decks to psihigh = 0.9995, as suggested, would move results substantially, because the domain gains the q = 7 surface:

diiid_n1, 0.995 → 0.9995
msing / mpert	5 → 6 / 35 → 36
pe_toroidal_torque	−0.0506 → +0.0492 (sign flip)
et_imag	+1.23e-4 → −1.71e-6 (sign flip)
ep_real, ev_real	−37 %, −24 %
NTV torque, resonant field	−26 %, −21 %
q0, q95	unchanged to 8 figures

It would also now put the overlap point inside their domain, so with the flag on the cap would pull them straight back to q ≤ 6 — the domain they already have. The change is therefore not made here; where the cap binds is covered by a dedicated deck (below).

Verified end to end

examples/DIIID-like_truncation_example runs the DIII-D-like equilibrium at psihigh = 0.9999 with the cap on — ideal free-boundary stability only, about two minutes:

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, so the cap sits at the q = 7 surface itself, and dmlim then stops 0.2 of a rational spacing past the last surface inside it, q = 6. The resulting domain, psilim = 0.99498, reproduces the ideal deck's hand-chosen psihigh = 0.995: ep_real −1.348513 vs −1.34833, et_real 0.8010979 vs 0.801353.

The new harness case diiid_truncation_n1 tracks that run, including the LayerOverlap outputs (applied, psilim_overlap, first_overlap_index, width_dr) that no case tracked before.

Correction to an earlier version of this description, which said dmlim's trunc "backs off to the rational below" the cap. It did here, but only because q(psilim_cap) happened to evaluate to 7.99999…: the cap lands exactly on a rational by construction, and a value rounding the other way gave qlim = 8.2, i.e. past the cap. Fixed in 14d60b6b9 (below); the capped domain no longer depends on that rounding.

The domain-limit convention

The scan cuts at ψ[k−1], the position of the last surface with a well-separated inner region. This is the paper's own 1 − ε_c, and it is what runtests_layer_overlap_fitzpatrick.jl reports in order to reproduce the paper's Ψ = 0.9985 / 0.9952.

Settled against ten equilibria (the shipped deck plus nine IDA DIII-D reconstructions) rather than by preference — the two half-width alternatives both fail:

convention	result
ψ[k−1] − w/2 (the original)	drops the last surface the criterion found clean; half a layer width more conservative than the paper
ψ[k−1]	always inside the plasma, always inside the first contaminated surface
ψ[k−1] + w/2	exceeds ψ = 1 on 4 of 10, and on those same 4 overshoots surface k itself

Full per-case table in the thread.

Bugs found while making the feature trustworthy

1. The edge q-law gate rejected 8 of 9 real DIII-D reconstructions (d4ce058ff, carries !). EDGE_FIT_MIN_R2 = 0.99 was calibrated against the synthetic TokaMaker deck (r²_log 0.9957, 0.9989). Real EFIT/CAKE reconstructions carry noise in the edge q and fit at 0.929–0.977, so all were refused: no surfaces placed past the grid, and the scan reported "no overlap" — a silent fail-open on exactly the case the criterion exists for. No absolute threshold can separate the classes (the limited a10 deck sits at 0.972, inside the diverted band); the r²_log > r²_linear comparison does, and is unchanged. Lowered to 0.90 at the time; since replaced by a physical shear-consistency check with a 0.5 backstop (see "Edge q-law gate" below). The ! is because the same gate controls the auto-grid edge density floor.

2. Reverse shear crashed the layer build — split out as #431, being pre-existing on develop. #431 has since merged; the cherry-picked copy resolved away in the develop sync merge (7e9e66508).

3. A single failing surface aborted the whole scan (6781b639f). The per-surface guard caught only ArgumentError.

4. The no-surfaces path threw instead of returning an empty scan (c02a3c43d, found in adversarial review). The early return built LayerOverlapScan with 13 arguments for a 15-field struct — a MethodError the pipeline's blanket catch downgraded to a spurious "scan failed" warning. Fixed, with a regression test forcing the path (m_max=0, extrapolation off), and the catch now logs the backtrace so a programming error is distinguishable from a physics refusal.

5. dmlim could step past qhigh and the cap (14d60b6b9, carries !). sing_lim! stepped the (last rational + dmlim)/n offset back only while it exceeded qmax, so qhigh = 4.0 gave qlim = 4.2, and a cap on a rational could be exceeded depending on float rounding. It now steps back against min(qmax, qhigh, cap): qhigh = 4.0 and 4.1 both give 3.2, and a cap on the exact rational and 1000 ulps either side all give the same qlim. 8 of the 12 new assertions fail against the old step-back. The existing cap-wiring tests all ran with set_psilim_via_dmlim = false, which is why this was not caught.

6. The scan's catch-all hid code errors (8cfbbee7e). resistive_layer_overlap already records expected per-surface failures in its notes and rethrows anything else; the outer catch downgraded those to a warning (bug 4 was found through exactly that). Only the kinetic-file read is guarded now.

Perpendicular diffusivity

The scan used the wrong χ (f5c60e779, carries !). It re-read the kinetic file with its own loader and passed no diffusivity, so every layer width was computed with chi_perp = chi_tor = 1.0, while the SLAYER analysis of the same surfaces used the file's χ(ψ). On the shipped profiles χ⊥ is 1.25–3.29 m²/s, and since delta_dr ∝ χ⊥^(1/4) the scan's widths were 6–35% narrower than SLAYER's own. The scan now takes its profiles, HDF5 group and χ(ψ) from [SLAYER] through SLAYER's own loader.

The shipped χ⊥ profile is doubled (ce29f7d99, carries !), from 1.25–3.29 to 2.50–6.59 m²/s; chi_phi is unchanged and the file's provenance attribute records the scaling. On this equilibrium the overlap cut is marginal — it first moves one surface inward at χ⊥ ×1.26 and next at ×~30 — so ×2 sits clear of both, centres the profile on TJ's default of 4 m²/s, and makes the truncation deck exercise a deeper cut. It is a deliberate change to input data, justified by how uncertain this proxy is (next paragraph), not by new physics.

Which χ each formula wants — recorded, and deferred. The layer width goes through τ_E = r²/χ_E, the one-fluid perpendicular energy diffusivity (every Fitzpatrick paper defines it that way). The critical-Δ island width W_d ∝ (χ_⊥/χ_∥)^(1/4) describes electron-temperature flattening against an electron χ_∥, so its χ_⊥ is the electron channel. GPEC feeds both from one profile, the kinetic file's chi_e, which is autoONETWO's electron power-balance CHI_E: the right channel for the critical Δ and a proxy for χ_E in the layer width. This matches TJ, which uses one Chip for both and assumes T_e = T_i. Separating the two needs an ion chi_i in the kinetic schema, and the power-balance ion χ available upstream is not yet usable (negative across 0.3 ≤ ρ < 0.95 on the case checked), so it is deferred and recorded in the stability docs as a to-do.

Critical-Δ citations corrected (9c1937bda). Cited in eight places as "Connor-Hastie-Helander 2015"; it is Connor, Ham, Hastie & Liu, PPCF 57, 065001 (2015), Eqs. (59) and (61), after Lütjens, Luciani & Garbet, Phys. Plasmas 8, 4267 (2001). The W_d loop and its short/long-mean-free-path χ_∥ blend follow Fitzpatrick, Phys. Plasmas 32, 062509 (2025).

Edge q-law gate

EDGE_FIT_MIN_R2 = 0.90 sat just under the worst real reconstruction then measured (0.929), so it was fitted to that sample, and the relative r2_log > r2_linear test was already doing the discriminating. It is replaced as the pathology guard by a physical check (cf85db8c5, carries !): the fitted A in q = q̄ + A(ln(1−ψ) − ū) is the local logarithmic shear, so it must match the profile's own −(1−ψ)q′ over the fitted knots, within a factor of 2. Across 14 real DIII-D reconstructions from IDA runs the ratio is 0.884–1.160 (all accepted), the synthetic DIII-D-like deck 0.948, and the limited Solovev edge 0.666. EDGE_FIT_MIN_R2 = 0.5 stays as an explicit backstop against a fit that describes nothing, documented as such.

Jacobian independence

@logan-nc asked for confirmation that delta_s gives a Jacobian-independent ψ width. Measured under pest / hamada / equal_arc: r_s invariant to 1e-16, da/dψ to 3e-14, |δ_s| and both Δψ widths to ~7e-6 (the equilibrium solve's own noise — the inputs are invariant to 1e-14). Same under :flux. Confirmed clean.

The check did find that :fsa and :halfwidth are coordinate-dependent (up to 7% in r_s); those labels live in #403, are not defaults, and are not used here.

Response to the review-package questions (@logan-nc, 2026-09-15)

question	outcome
convention vs dmlim	kept as agreed; dmlim now provably stays inside the cap (bug 5)
two width channels	Eq. (100) kept as the criterion; the
shipped decks stay a no-op	yes — flag off everywhere; the psihigh = 0.9995 harmonisation is declined for the reasons above
Solovev auto-grid classification	added to runtests_grid_refinement.jl (9f2010cfd): Solovev is rejected as limited on both passes, including with the absolute r² floor disabled, and the DIII-D-like equilibrium is accepted as a positive control
always-on cost	kept; the scan is ~15 layer solves per n

Notes for reviewers

• The scan is driven at rs_method = :flux, not the SLAYER default. Eq. (100) is not covariant — δ ∝ (q/|dq/dr|)^(1/2) — so it is only valid in the coordinate it was derived in, the toroidal-flux label of the paper's Eq. (30).
• delta_s_abs appears in a third HDF5 group. The one-name-per-quantity rule says this is right, but @logan-nc asked to be told if it stayed ambiguous.
• Still open, not gated on this PR: the two width channels (|δ_s| Riccati vs Eq. 100) disagree on where the layers first overlap. Both are reported.
• Deferred: separate χ_E (layer width) and electron χ⊥ (critical Δ) inputs; see Perpendicular diffusivity.
• Seen in passing, not addressed: the Newcomb zero-crossing check (evaluate_stability_criterion!) reports a non-Hermitian W-inverse at ~40 integration steps on the DIII-D-like equilibrium at any psihigh. It runs after every forward and Riccati integration whatever vac_flag is. Correction: an earlier version of this description, and the message of 6de1b8f03, said running the truncation deck with vac_flag = true avoided it; it does not.

Regression report

regress --cases <all 17 cases> --refs cb18e386e,af9a8badd --force — develop as merged into the branch vs this head (the merge commit); julia 1.11.6, pinned manifest:

diiid_error_field           Summary: 8 unchanged
diiid_multi_n               Summary: 19 unchanged
diiid_n1                    Summary: 53 unchanged
diiid_n1_riccati            Summary: 17 unchanged
diiid_slayer_n1             Summary: 4 changed, 13 unchanged
                            SLAYER P_perp                5.030e+01 (100.00%)  ** CHANGED **
                            SLAYER Q_root [2/1,3/1,4/1]  2.251e-03 (0.32%)    ** CHANGED **
                            SLAYER ω_Hz [2/1,3/1,4/1]    9.325e+00 (0.08%)    ** CHANGED **
                            SLAYER γ_Hz [2/1,3/1,4/1]    8.353e+01 (8.28%)    ** 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; see below)
ggj_ray_q500i               Summary: 4 unchanged
ggj_reference               Summary: 4 unchanged
solovev_kinetic_calculated  Summary: 14 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

220 quantities compared, 216 unchanged. The four changes are the doubled χ⊥ on diiid_slayer_n1: P_perp is exactly ×2 on every surface, and the growth rates move 2/1 194.6 → 173.6 Hz, 3/1 −249.2 → −226.8 Hz, 4/1 −1008.8 → −925.3 Hz — far above the root-search jitter. Every SLAYER input that does not depend on χ⊥ (S, D_norm, tauk, iota_e, r_s) is unchanged, and so is every other case, including diiid_n1, whose NTV does not read χ. So the SLAYER-only profile path and the edge-gate change move nothing on any other deck.

gal_resistive_pe extracts nothing on either ref, develop included, so it predates this PR; it is being looked at separately.

Runtimes are not quoted: the run was not made under controlled load.

diiid_truncation_n1 is new, so it runs at this head only:

Ref 1: af9a8badd
layer-overlap cap applied              1
layer-overlap psilim                   9.985081e-01
first overlapping surface index        7
diffusive-resistive layer widths (Δψ)  [15 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.348513e+00
total energy Re(et[1])                 8.010979e-01

Changes from adversarial review

• The criterion no longer depends on the Riccati solve (cd6dce306, carries !). The scan
  called slayer_layer_thickness, which always runs the stiff riccati_del_s solve, and dropped
  any surface whose delta_s came back non-finite -- even though the criterion uses only the
  closed-form delta_dr. A deleted surface left _first_overlap_limit comparing its neighbours
  across a gap twice the physical one, hiding an overlap and reporting too wide a domain, silently.
  delta_s is now best-effort (NaN plus a note, its own channel only); a surface that cannot be
  scored on the DR channel is skipped with an explicit "not conclusive" note, matching the m_max
  and max_layer_solves backstops. Measured unchanged: psihigh = 0.9997535788 over 14 surfaces,
  no notes, 3.9 s per n.
• The truncation case pins physics again (6de1b8f03). It ran vac_flag = false, which writes
  no energies, so only bookkeeping was tracked on the one deck where the cap binds. The deck now
  runs free-boundary and the case pins ep_real and et_real.
• refined_psi_grid's psihigh keyword and _truncate_density are removed (9848fa332) --
  left from the design that re-formed the equilibrium, unreachable and untested since the qlim
  cap replaced it. (The Migration note below no longer mentions them.)
• The flag is inert on the programmatic solve(prob, alg) path, which resolves no kinetic
  file; documented in the control docstring and the stability page (a3e7c51fa).
• LayerOverlapScan gains a keyword constructor -- the positional form takes 15 arguments in a
  fixed order, which is how the earlier empty-scan return came to pass 13.

Validation

• Develop merged into this branch at 620e4c531, no conflicts
• runtests_layer_overlap_fitzpatrick.jl at this head: JET reproduction 9/9 (Ψ = 0.9985111 vs paper 0.9985 for n = 1, 0.9949738 vs 0.9952 for n = 4, plus the no-ODE width path asserted equal to the full one), cap wiring 28/28
• runtests_grid_refinement.jl at this head: 66/66, including a test that the shear check is live
• runtests_slayer_params.jl 71/71, runtests_kinetic.jl 277/277, runtests_slayer_inputs.jl 62/62
• Edge-fit shear check on 14 real DIII-D reconstructions from IDA runs: ratio 0.884–1.160, all accepted; Solovev 0.666, rejected
• Multi-n end to end on the DIII-D-like equilibrium at psihigh = 0.9999: n = 1 overlaps at 0.99975, n = 2 at 0.99935, nn = 1:2 binds at n = 2
• Docs build locally with the new cross-references resolving
• Full suite: left to CI on this push

Important

Do not merge without human 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?

[jhalpern30]

Please ping me again when this is ready for review - it seems to still be being developed. Turning off auto-notifications on this for now to avoid my inbox being blown up

[logan-nc]

Same @d-burg. It says

d-burg marked this pull request as ready for review 15 hours ago

above but I see no comment explaining how/why it got promoted from draft. The SLAYER stack is big rn so I am a bit lost - hold our hands through this please 🙏

[d-burg]

I put together a review package for this PR — the overlap criterion and cut conventions as one diagram, the always-recorded / opt-in-applied wiring, the ten-equilibrium convention evidence, and the bugs found en route.

https://claude.ai/code/artifact/71b5b47a-7b76-45e8-9714-bd43041bcde4

[logan-nc]

@d-burg lets merge this one way or another this week. Responses to the review package artifact questions:

Does the convention choice survive contact with dmlim?

I am fine with this for now. The result is perhaps a very small fraction of "ideal-ish" plasma getting truncated out - it is still much better than what we had before and pushing further out into the known-stochastic region feels unjustified as well as numerically harder.

The two width channels still disagree by about two surfaces. |δ_s| (Riccati, the dispersion-relation equivalent sheet thickness) cuts at 0.98978; Eq. (100) (the spatial decay scale) cuts at 0.99939.

I believe @d-burg declared Eq. (100) to be the more consistent one to the paper and thus agree this is where we should start. Besides, 0.98978 is pretty deep and starts to cut into the pedestal... so more thought would be needed about what that means → future task to be considered if/when this has been used a bit more so we have more experience.

Should the shipped decks stop being a no-op?

No. All DIII-D like examples should get psihigh = 0.9995 like already done in the SLAYER and IMAS examples. It is an op when this is used.

The edge_q_law gate change is exercised by diiid_n1 and diiid_slayer_n1 (both auto grids, both unchanged) but by no limited deck on an auto grid — every Solovev/LAR case uses grid_type = "ldp" with fixed mpsi. The limited-plasma branch is argued from measured r² values, not covered by a harness case. Should a Solovev auto-grid case be added?

Nah. Too much of a change. Lets just add a super quick test that does solovev with auto-grid and confirms it is correctly classified as limited - no full example needed.

Cost lands on runs that never asked for the feature.

This is fast and it's always valuable to have the basic outputs - so when users come to developers with "my run is funny" we can quickly look it up and see if the truncation they used was "reasonable". This is the number one gotcha for new users.

[logan-nc]

Approved once the task list associated with the response to the review package questions is addressed.

[jhalpern30]

I gave a very cursory review and did not read the entire chain. These comments are all optional and I'll leave it up to you if you decide any of them are worth it, its just things I saw when skimming the code; sorry if anything was repeated from the above. Don't want to add too much work on a PR that's been sitting for a while

The implementation specifics (i.e. dmlim logic, opting in to get truncation, warning the user if they requested a psihigh past truncation, no-ops on most shipped decks) all seem ok with me

[jhalpern30]

This should probably just be in the docstring (which needs to be updated with the extra argument anyway)

[d-burg]

Moved into the docstring, which now documents psilim_cap and the dmlim bound; the inline comment is down to one line (14d60b6).

Doing this turned up a real bug. dmlim stepped back only against qmax, so it could land past the cap, and past qhigh too: qhigh = 4.0 gave qlim = 4.2. And because the cap sits exactly on a rational, whether it was exceeded came down to how q(psilim_cap) rounded. dmlim now steps back until it is inside min(qmax, qhigh, cap), so qhigh = 4.0 and 4.1 both give 3.2; the tests put a cap on the exact rational and 1000 ulps either side. No shipped deck moves.

[jhalpern30]

Why not just take the most conservative truncation limit? Like min(layer_overlap(n) for n in nlow:nhigh)?

[logan-nc]

Agreed

[d-burg]

Done in 8cfbbee. The scan now runs for each n in nn_low:nn_high and the innermost overlap point binds. On the DIII-D-like equilibrium at psihigh = 0.9999, n=1 overlaps at 0.99975 and n=2 at 0.99935, so nn = 1:2 binds at n=2. The recorded LayerOverlap group is the binding scan, and its n dataset says which.

[jhalpern30]

Is there an actual reason for the try catch here or is this an AI artifact?

[logan-nc]

skimming these comments - but I think the resistive_layer_overlap itself might throw errors? I agree this setup is opaque and could/should be cleaner

[d-burg]

Fair question: it was hiding real errors (the empty-scan MethodError was downgraded to a "scan failed" warning). resistive_layer_overlap already records expected per-surface failures in its notes and rethrows anything else, so the outer catch-all is gone (8cfbbee). The only guard left is on reading the kinetic file, since that is external input to an optional diagnostic and a bad file should skip the scan rather than abort someone's stability run.

[jhalpern30]

Is false the right default here? If this is a true physical limit, do we want people to have to know that they need to turn this on or they might accidentally set truncation out too far?

[d-burg]

Keeping it false, now with numbers. On the shipped DIII-D-like decks the overlap point is at ψ ≈ 0.9997, beyond every deck's psihigh, so the cap cannot bind there and turning it on changes nothing (measured: identical results flag on vs off, even at psihigh = 0.9995).

On setting the truncation out too far: the scan always runs and is recorded under ForceFreeStates/LayerOverlap/, and GPEC now warns whenever psilim lands past the overlap point, flag on or off (14d60b6). DIIID-like_truncation_example (psihigh = 0.9999) is the case where the cap binds, with a harness case tracking it.

[jhalpern30]

If this is a property of the equilibrium, should this be done during setup_equilibrium? And along the same lines, would it make sense for the layer_overlap to be a member of the PlasmaEquilibrium struct? Completely up to you, just a thought. I guess its a balance of scope/struct bloat vs. carrying along layer_overlap everywhere

[d-burg]

I'd keep it out of the equilibrium. The overlap point depends on n (now on the whole nn range) and on the kinetic profiles, neither of which PlasmaEquilibrium carries, so it belongs to the run rather than the equilibrium. It is computed once in prepare_force_free_states! and only goes to sing_lim! and the HDF5 writer, so it isn't carried around anywhere else.

[logan-nc]

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