EQUIL/EXAMPLES - Default psihigh_from_layer_overlap on across every deck

Flips EquilibriumConfig.psihigh_from_layer_overlap to true and sets it
explicitly, annotated, in all 17 decks. Kinetic profiles are wired via a
[SLAYER] profile_file entry using each deck's OWN family file -- never a
[KineticForces] block, which would flip needs_kinetic_profiles and change the
auto-grid knot density, an unrelated results channel. A [SLAYER] section
carrying only profile_file is inert: slayer_ctrl.enabled defaults false and the
stage returns before running or writing anything.

Deck sweep, equilibrium stage, all 17:
  MOVED  DIIID-like_SLAYER          0.9995 -> 0.993281
  MOVED  DIIID-like_ideal           0.995  -> 0.993281
  no-op  DIIID gal_resistive, _pe   0.993 (already at the floor)
  no-op  Solovev x5, LAR x4, a10
  n/a    DIIID-like_ideal_IMAS      needs a runtime dd; flag left off anyway

The analytic decks no-op exactly as predicted: limited plasmas with finite edge
q, so the edge q-law fit is rejected and the two well-separated surfaces never
overlap. LAR has no legitimate kinetic file, so those four no-op on the missing
profiles instead.

NOTE FOR REVIEW: on all four DIII-D decks the psiedge floor OVERRIDES the layer
criterion. Each asks for ~0.9892 and each is held at the q=6 window top. The
feature as configured therefore reduces to "clamp psihigh to the first rational
surface outboard of psiedge"; the layer physics sets the domain nowhere in the
shipped examples.

IMAS deck keeps the flag off with a comment: its rerun path re-runs read_imas and
must resolve the same psihigh on both passes (_validate_psi_nodes errors
otherwise), which is what resolving a new domain changes; "imas" is also outside
EFIT_KINDS so there is no separatrix clamp and no raw psi map.

benchmarks/equil_psihigh_scan.jl now records the structured equilibrium_edge_quality
verdict per point instead of a hand-rolled round-trip, captures et/ep/ev, and adds
a gate-validation pass that fails loudly if any psihigh exists where the gate
passes while the energies are already unusable. It forces the flag off so the
ladder tests the psihigh it sets.

docs: corrects a naming bug I introduced -- workflow.md called the variable
dW_edge_scan_start, but that rename is an unmerged draft (PR #382) and the code
says psiedge. Documents the floor, the cross-section dependency on
[ForceFreeStates] psiedge, and that truncations are reported in mm as well as psi.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01PSrf6JCViFfVzqzkQ66o6b

Author: logan-nc

benchmarks/equil_psihigh_scan.jl

@@ -6,9 +6,14 @@ methods near the separatrix.
 
 Key outputs:
   - Last successful psihigh for each method
-  - Max roundtrip error in outer 10% (ψ > 0.90) vs. psihigh
+  - The `equilibrium_edge_quality` verdict (jacobian / roundtrip / ringing) vs. psihigh
   - q monotonicity violations in outer 10%
-  - et[1] (free-boundary MHD stability eigenvalue) vs. psihigh
+  - et[1], ep[1], ev[1] (free-boundary MHD energies) vs. psihigh
+  - GATE VALIDATION: does the edge-quality verdict flip at or before the energies go bad?
+
+This doubles as the validation of the psihigh back-off gate used by
+`[Equilibrium] psihigh_from_layer_overlap`. That flag is forced OFF here so the
+ladder actually tests each requested psihigh instead of the feature resolving its own.
 
 Physics note: In a diverted tokamak (DIIID), q → ∞ as ψ → 1. The EFIT q-profile
 extrapolates to a finite qa, which is physically wrong near the separatrix. This
@@ -29,7 +34,7 @@ using GeneralizedPerturbedEquilibrium.Equilibrium
 using TOML, Printf, Statistics, HDF5
 
 example_path = length(ARGS) > 0 ? ARGS[1] : joinpath(@__DIR__, "../examples/DIIID-like_ideal_example")
-config_path  = joinpath(example_path, "gpec.toml")
+config_path = joinpath(example_path, "gpec.toml")
 
 psihigh_values = [0.980, 0.985, 0.990, 0.993, 0.995, 0.996, 0.997, 0.998, 0.999, 0.9995, 0.9999, 1.0]
 methods = ["efit", "efit_arclength", "efit_by_inversion"]
@@ -38,20 +43,22 @@ function make_config(path::String, eq_type::String, psihigh::Float64)
     raw = TOML.parsefile(path)
     raw["Equilibrium"]["eq_type"] = eq_type
     raw["Equilibrium"]["psihigh"] = psihigh
+    raw["Equilibrium"]["psihigh_from_layer_overlap"] = false  # the ladder sets psihigh; don't let the feature move it
     base = dirname(path)
     return Equilibrium.EquilibriumConfig(raw["Equilibrium"], base)
 end
 
 """
 Run the full ForceFreeStates pipeline for a given eq_type and psihigh, returning
-real(et[1]). Uses a temporary directory so the example directory is not modified.
-Returns NaN on failure.
+(et[1], ep[1], ev[1]) as reals. Uses a temporary directory so the example directory
+is not modified. Returns all-NaN on failure.
 """
-function run_ffs_et1(config_path::String, eq_type::String, psihigh::Float64)::Float64
+function run_ffs_energies(config_path::String, eq_type::String, psihigh::Float64)
     example_dir = dirname(config_path)
     raw = TOML.parsefile(config_path)
     raw["Equilibrium"]["eq_type"] = eq_type
     raw["Equilibrium"]["psihigh"] = psihigh
+    raw["Equilibrium"]["psihigh_from_layer_overlap"] = false  # the ladder sets psihigh; don't let the feature move it
     raw["ForceFreeStates"]["force_termination"] = true   # write gpec.h5 after FFS, skip PE
     raw["ForceFreeStates"]["vac_flag"] = true
     raw["ForceFreeStates"]["verbose"] = false
@@ -69,22 +76,24 @@ function run_ffs_et1(config_path::String, eq_type::String, psihigh::Float64)::Fl
         try
             GeneralizedPerturbedEquilibrium.main([tmpdir])
             h5open(joinpath(tmpdir, "gpec.h5"), "r") do h5
-                et = read(h5["ForceFreeStates/FreeBoundaryStability/eigenmode_energies"])
-                return real(et[1])
+                g = h5["ForceFreeStates/FreeBoundaryStability"]
+                return (real(read(g["eigenmode_energies"])[1]),
+                    real(read(g["eigenmode_plasma_energies"])[1]),
+                    real(read(g["eigenmode_vacuum_energies"])[1]))
             end
         catch
-            return NaN
+            return (NaN, NaN, NaN)
         end
     end
 end
 
 # JIT warmup: run equilibrium + FFS once before timing
 println("JIT warmup (equilibrium + FFS)...")
-run_ffs_et1(config_path, "efit", 0.990)
+run_ffs_energies(config_path, "efit", 0.990)
 
-println("=" ^ 70)
+println("="^70)
 println("psihigh robustness scan (equilibrium + stability)")
-println("=" ^ 70)
+println("="^70)
 
 rows = []
 
@@ -109,31 +118,24 @@ for method in methods
         q_at_psihigh = NaN
         q_edge_slope = NaN
         et1 = NaN
+        ep1 = NaN
+        ev1 = NaN
+        gate_ok = false
+        gate_failed = "formation-threw"
+        jac_ok = false
+        ringing = NaN
 
         if success
             raw_profile = Equilibrium.read_efit(cfg)
             psi_xs = pe.rzphi_xs
-            mtheta = length(pe.rzphi_ys) - 1
-            psio = pe.psio
-
-            # Roundtrip error in outer 10%
-            errors = Float64[]
-            for ipsi in 1:length(psi_xs)
-                psi_xs[ipsi] < 0.90 && continue
-                ψ = psi_xs[ipsi]
-                for itheta in 1:8:mtheta+1
-                    θ = pe.rzphi_ys[itheta]
-                    r2  = pe.rzphi_rsquared((ψ, θ))
-                    off = pe.rzphi_offset((ψ, θ))
-                    rfac = sqrt(max(r2, 0.0))
-                    η = 2π * (θ + off)
-                    R = pe.ro + rfac * cos(η)
-                    Z = pe.zo + rfac * sin(η)
-                    ψ_recon = 1.0 - raw_profile.psi_in((R, Z)) / psio
-                    push!(errors, abs(ψ_recon - ψ))
-                end
-            end
-            roundtrip_max_edge = isempty(errors) ? NaN : maximum(errors)
+
+            # One implementation of the edge verdict, shared with the psihigh back-off gate.
+            eq = Equilibrium.equilibrium_edge_quality(pe; psi_map=raw_profile.psi_in, psio=pe.psio)
+            gate_ok = eq.ok
+            gate_failed = isempty(eq.failed) ? "-" : join(eq.failed, "+")
+            jac_ok = eq.jacobian_ok
+            ringing = eq.ringing
+            roundtrip_max_edge = eq.roundtrip_err
 
             # q monotonicity in outer 10%
             q_profile = [pe.profiles.q_spline(ψ) for ψ in psi_xs]
@@ -149,25 +151,27 @@ for method in methods
                                (psi_xs[end] - psi_xs[end-1])
             end
 
-            # et[1]: full FFS run
-            et1 = run_ffs_et1(config_path, method, psihigh)
+            # Downstream sanity: full FFS run
+            et1, ep1, ev1 = run_ffs_energies(config_path, method, psihigh)
         end
 
         status = success ? "OK" : "FAIL"
         et1_str = isnan(et1) ? "   FAIL" : @sprintf("%+.4f", et1)
-        @printf("  psihigh=%.4f  %-4s  t=%.2fs  rt_edge=%.2e  q_mono_viol=%d  q(end)=%.2f  et[1]=%s\n",
-            psihigh, status, runtime, roundtrip_max_edge,
-            Int(isnan(q_mono_violations_edge) ? -1 : q_mono_violations_edge),
-            isnan(q_at_psihigh) ? -1.0 : q_at_psihigh,
-            et1_str)
-
-        push!(rows, (
-            method=method, psihigh=psihigh, success=success, runtime=runtime,
-            roundtrip_max_edge=roundtrip_max_edge,
-            q_mono_violations_edge=q_mono_violations_edge,
-            q_at_psihigh=q_at_psihigh, q_edge_slope=q_edge_slope,
-            et1=et1, error_msg=err_msg
-        ))
+        @printf("  psihigh=%.4f  %-4s  gate=%-4s(%-16s) rt=%.2e  ring=%-8.3g q_viol=%d  et[1]=%s\n",
+            psihigh, status, gate_ok ? "PASS" : "FAIL", gate_failed, roundtrip_max_edge, ringing,
+            Int(isnan(q_mono_violations_edge) ? -1 : q_mono_violations_edge), et1_str)
+
+        push!(
+            rows,
+            (
+                method=method, psihigh=psihigh, success=success, runtime=runtime,
+                gate_ok=gate_ok, gate_failed=gate_failed, jacobian_ok=jac_ok, ringing=ringing,
+                roundtrip_max_edge=roundtrip_max_edge,
+                q_mono_violations_edge=q_mono_violations_edge,
+                q_at_psihigh=q_at_psihigh, q_edge_slope=q_edge_slope,
+                et1=et1, ep1=ep1, ev1=ev1, error_msg=err_msg
+            )
+        )
     end
 end
 
@@ -181,7 +185,7 @@ end
 println("\nResults saved to: $output_csv")
 
 # Summary
-println("\n" * "=" ^ 70)
+println("\n" * "="^70)
 println("Last successful psihigh per method:")
 for method in methods
     method_rows = filter(r -> r.method == method && r.success, rows)
@@ -193,14 +197,69 @@ for method in methods
     end
 end
 
-println("\n" * "=" ^ 70)
+# ---------------------------------------------------------------------------
+# GATE VALIDATION -- does the edge-quality verdict track downstream sanity?
+#
+# The gate is validated if, walking psihigh outward, it never says PASS while the
+# ForceFreeStates energies are already garbage. A PASS-with-garbage band means the
+# metric is blind to the failure it exists to catch, and the psihigh back-off gate
+# must not be trusted (nor the feature defaulted on) until that is fixed.
+# ---------------------------------------------------------------------------
+println("\n" * "="^70)
+println("GATE VALIDATION: edge-quality verdict vs. ForceFreeStates sanity")
+println("="^70)
+et_finite(r) = isfinite(r.et1) && isfinite(r.ep1) && isfinite(r.ev1)
+verdict_bad = NamedTuple[]
+for method in methods
+    mrows = sort(filter(r -> r.method == method, rows); by=r -> r.psihigh)
+    # Reference energy from the innermost psihigh that produced finite energies.
+    ref = findfirst(r -> r.success && et_finite(r), mrows)
+    et_ref = ref === nothing ? NaN : mrows[ref].et1
+    println("\n--- $method  (reference et[1] = $(isnan(et_ref) ? "n/a" : @sprintf("%+.4f", et_ref))) ---")
+    @printf("  %-9s %-7s %-6s %-18s %-11s %-11s %s\n",
+        "psihigh", "formed", "gate", "gate failed on", "et[1]", "|et/et_ref|", "assessment")
+    for r in mrows
+        sane = r.success && et_finite(r) &&
+               (isnan(et_ref) || abs(r.et1) < 100 * max(abs(et_ref), eps()))
+        ratio = (isnan(et_ref) || !isfinite(r.et1)) ? NaN : abs(r.et1 / et_ref)
+        assess =
+            !r.success ? "formation threw" :
+            (sane && r.gate_ok) ? "consistent (both good)" :
+            (!sane && !r.gate_ok) ? "consistent (both bad)" :
+            (!sane && r.gate_ok) ? "*** GATE BLIND: passes but energies bad ***" :
+            "conservative (gate fails, energies still ok)"
+        (!sane && r.gate_ok && r.success) && push!(verdict_bad, r)
+        @printf("  %-9.4f %-7s %-6s %-18s %-11s %-11s %s\n", r.psihigh,
+            r.success ? "yes" : "THREW", r.gate_ok ? "PASS" : "FAIL", r.gate_failed,
+            isfinite(r.et1) ? @sprintf("%+.4f", r.et1) : "nonfinite",
+            isnan(ratio) ? "-" : @sprintf("%.3g", ratio), assess)
+    end
+end
+
+println("\n" * "="^70)
+if isempty(verdict_bad)
+    println("GATE VALIDATED: no psihigh where the gate passes but the energies are bad.")
+else
+    println("GATE NOT VALIDATED -- $(length(verdict_bad)) point(s) where the gate PASSES")
+    println("but the ForceFreeStates energies are already unusable:")
+    for r in verdict_bad
+        @printf("   %-20s psihigh=%.4f  et[1]=%s\n", r.method, r.psihigh,
+            isfinite(r.et1) ? @sprintf("%+.4f", r.et1) : "nonfinite")
+    end
+    println("\nFix the edge-quality metric before relying on the back-off gate.")
+end
+println("="^70)
+
+println("\n" * "="^70)
 println("et[1] vs psihigh (free-boundary stability eigenvalue):")
 println("  (positive = stable, negative = unstable)")
 @printf("  %-22s  %s\n", "psihigh", join([@sprintf("%-22s", m) for m in methods]))
 for psihigh in psihigh_values
-    vals = [let r = findfirst(x -> x.method == m && x.psihigh == psihigh, rows)
-                r === nothing || isnan(rows[r].et1) ? "      -" : @sprintf("%+.4f", rows[r].et1)
-            end for m in methods]
+    vals = [
+        let r = findfirst(x -> x.method == m && x.psihigh == psihigh, rows)
+            r === nothing || isnan(rows[r].et1) ? "      -" : @sprintf("%+.4f", rows[r].et1)
+        end for m in methods
+    ]
     @printf("  %-22.4f  %s\n", psihigh, join([@sprintf("%-22s", v) for v in vals]))
 end
 println()

docs/src/workflow.md

@@ -52,7 +52,7 @@ to different stages and are computed at different times:
 |:--- |:--- |:--- |:--- |
 | `psihigh` | `[Equilibrium]` | The **equilibrium spline domain** — where the ψ grid ends and every spline stops being defined | The deck, then optionally the separatrix clamp and the layer-overlap scan below |
 | `psilim` | internal (`ForceFreeStatesInternal`) | The **stability integration domain** — where the Euler-Lagrange / Δ′ integration stops. Always ≤ `psihigh` | `sing_lim!`, from `qhigh`, `dmlim`, and the equilibrium's `qmax` |
-| `dW_edge_scan_start` | `[ForceFreeStates]` | Start of a **diagnostic** dW(ψ) scan. Moves nothing | The deck |
+| `psiedge` | `[ForceFreeStates]` | Start of the **diagnostic** dW(ψ) scan over `[psiedge, psilim]`. Moves nothing itself, but floors the layer-overlap truncation below | The deck |
 
 `psihigh` is chosen before the equilibrium is built; `psilim` is derived from the finished
 equilibrium. Shrinking `psihigh` therefore also pulls `psilim` in, but not the reverse.
@@ -71,6 +71,11 @@ by default:
 `truncate_at_dW_peak` is **not** in this list. It moves `psilim`, not `psihigh`, is documented
 as making Δ′ and δW unreliable, and is retained only for experimental work.
 
+!!! note "Pending rename"
+    PR #382 proposes renaming `psiedge` to `dW_edge_scan_start`, precisely because it reads like
+    the plasma edge when it is only a diagnostic scan start. That PR is still a draft, so the
+    name in the code and in every deck today is `psiedge`.
+
 ### Choosing `psihigh` from resistive-layer physics
 
 With `psihigh_from_layer_overlap = true`, the deck's `psihigh` becomes a *first-pass* bound —
@@ -84,7 +89,13 @@ field-line-tracing noise — and the final domain is computed:
    Widths are converted from meters to normalized flux before being compared against surface
    separations.
 3. Where two neighbouring layers overlap, neither surface has a well-separated inner region,
-   so the domain is cut at the inner edge of the inner one.
+   so the domain is cut at the inner edge of the inner one — **floored** at
+   `min(psihigh, top of the rational-q window containing psiedge)`, i.e. the first rational
+   surface outboard of `psiedge`. The dW edge scan runs over `[psiedge, psilim]`, so the domain
+   has to stay far enough out for its maximum to be findable. When the floor binds, the run
+   warns that the retained surfaces are set by `psiedge` rather than by the layer physics.
+   Setting `psiedge >= psihigh` (the documented "scan disabled" sentinel) removes the floor.
+   The scan also refuses to truncate at all if fewer than two rational surfaces would survive.
 4. If that asks for a larger domain than the first pass, form it and check it with
    [`equilibrium_edge_quality`](@ref Equilibrium.equilibrium_edge_quality); on failure, back off toward the first-pass bound and
    warn, naming the check that failed.
@@ -93,7 +104,12 @@ field-line-tracing noise — and the final domain is computed:
 
 The scan needs kinetic profiles, taken from `[SLAYER] profile_file` when present and
 `[KineticForces] kinetic_file` otherwise; with neither configured it is skipped with a message.
-Single-`n` runs only. The domain is chosen once, from the first-pass equilibrium, and the final
+It also reads `psiedge` from `[ForceFreeStates]` for the floor — so although the switch lives in
+`[Equilibrium]`, the behaviour depends on two settings in other sections. Single-`n` runs only.
+
+Truncations are reported in millimetres of minor radius as well as in ψ. A cut that reads as
+large in ψ is usually small in physical terms, because ψ is heavily compressed near the edge:
+on the DIII-D-like deck, 0.9995 → 0.98919 removes 3.6 mm, 0.66% of the minor radius. The domain is chosen once, from the first-pass equilibrium, and the final
 equilibrium is re-scanned as a check — it is deliberately not iterated to a fixed point, since
 the auto grid it sits on top of does not itself reach one.
 

examples/DIIID-like_SLAYER_example/gpec.toml

@@ -11,7 +11,7 @@ jac_type = "hamada"            # Coordinate system (hamada, pest, boozer, equal_
 grid_type = "log_asymptotic"   # Radial grid packing type
 psilow = 1e-4                  # Lower limit of normalized flux coordinate
 psihigh = 0.9995               # Upper limit of normalized flux coordinate (bounds the equilibrium splines; the integration bound psilim is derived from qhigh/dmlim)
-psihigh_from_layer_overlap = false  # Choose the final psihigh from resistive-layer overlap instead of taking it from this file; when true the psihigh above is the first-pass (spline-safe) bound
+psihigh_from_layer_overlap = true   # Resolve the final psihigh from resistive-layer overlap; the psihigh above is then the first-pass (spline-safe) bound, floored by the rational-q window psiedge sits in
 mpsi = 0                       # Number of radial grid points (0 = auto-compute from psi_accuracy)
 psi_accuracy = 0.001           # Target absolute error in q for auto-mpsi
 mtheta = 256                   # Number of poloidal grid points

examples/DIIID-like_gal_resistive_example/gpec.toml

@@ -11,6 +11,7 @@ jac_type = "hamada"            # Coordinate system (hamada, pest, boozer, equal_
 grid_type = "ldp"              # Radial grid packing type (ldp = linear-derivative packing toward rationals)
 psilow = 1e-4                  # Lower limit of normalized flux coordinate
 psihigh = 0.993                # Upper limit of normalized flux coordinate (0.993 stays clear of the separatrix; truncating at ≳0.998 is numerically unreasonable here)
+psihigh_from_layer_overlap = true   # Resolve the final psihigh from resistive-layer overlap; the psihigh above is then the first-pass (spline-safe) bound, floored by the rational-q window psiedge sits in
 mpsi = 128                     # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 psi_accuracy = 0.001           # Target relative accuracy of splined profile derivatives for the auto grid
 mtheta = 256                   # Number of poloidal grid points
@@ -83,3 +84,9 @@ gal_eta                = [8e-8, 8e-8, 8e-8, 8e-8]          # Per-surface resisti
 gal_rho                = [3.3e-7, 3.3e-7, 3.3e-7, 3.3e-7]  # Per-surface mass density ρ [kg/m³] (core→edge)
 gal_rotation           = [1.0, 1.0, 1.0, 1.0]             # Per-surface rotation frequency f [Hz] (core→edge); forced eigenvalue γ_s = 2πi·n·f
 gal_gamma              = 1.6666666666666667               # Ratio of specific heats Γ for the resistive-layer coefficients (5/3)
+
+[SLAYER]
+# Kinetic profiles for the resistive-layer overlap scan that resolves psihigh
+# ([Equilibrium] psihigh_from_layer_overlap). The tearing analysis itself stays off
+# (enabled defaults false), so only profile_file is read here.
+profile_file            = "../DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode_kinetic.h5"  # Kinetic profile file (HDF5 GPEC kinetic schema or ASCII table)

examples/DIIID-like_gal_resistive_pe_example/gpec.toml

@@ -10,6 +10,7 @@ jac_type = "hamada"            # Coordinate system (hamada, pest, boozer, equal_
 grid_type = "ldp"              # Radial grid packing type (ldp = linear-derivative packing toward rationals)
 psilow = 1e-4                  # Lower limit of normalized flux coordinate
 psihigh = 0.993                # Upper limit of normalized flux coordinate (0.993 stays clear of the separatrix; truncating at ≳0.998 is numerically unreasonable here)
+psihigh_from_layer_overlap = true   # Resolve the final psihigh from resistive-layer overlap; the psihigh above is then the first-pass (spline-safe) bound, floored by the rational-q window psiedge sits in
 mpsi = 128                     # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 psi_accuracy = 0.001           # Target relative accuracy of splined profile derivatives for the auto grid
 mtheta = 256                   # Number of poloidal grid points
@@ -101,3 +102,9 @@ compute_singular_coupling = true        # Compute singular layer coupling metric
 verbose = true                          # Enable verbose logging
 write_outputs_to_HDF5 = true            # Write perturbed equilibrium outputs to HDF5
 reg_spot = 0.05                         # Regularization width for singular surfaces (0 = disabled)
+
+[SLAYER]
+# Kinetic profiles for the resistive-layer overlap scan that resolves psihigh
+# ([Equilibrium] psihigh_from_layer_overlap). The tearing analysis itself stays off
+# (enabled defaults false), so only profile_file is read here.
+profile_file            = "../DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode_kinetic.h5"  # Kinetic profile file (HDF5 GPEC kinetic schema or ASCII table)