EQUIL - NEW FEATURE - Record the resistive-layer domain scan in gpec.h5

The scan now runs whenever kinetic profiles are reachable, with only the
enforcement gated by psihigh_from_layer_overlap, so the file always records what
the layer physics wanted for the domain.

  Equilibrium/psihigh_resolved            domain actually reconstructed
  Equilibrium/psihigh_resistive_layers    where adjacent layers start to overlap
  Equilibrium/ResistiveLayers/            per-surface widths for every scored surface
    rational_psi, rational_m, rational_n, r_s,
    delta_s_abs [m], delta_s_psi [1], extrapolated

Comparing the two scalars answers "did the domain honour the layer physics", which
previously required reading the run log.

Both width forms are kept deliberately: delta_s_abs is the physical thickness and
delta_s_psi is the same width in normalized flux, which is what the overlap
criterion compares against surface spacing. delta_s_abs matches the leaf name
already used in Tearing/LayerWidths per the one-name-per-quantity rule.

psihigh_resolved matches the name PR #367 (immutable input structs) introduces.
It is written explicitly here only because that PR has not landed; once it does,
EquilibriumParameters carries the field and the EQUIL_H5_NAMES loop writes it, so
the explicit line must be deleted or the dataset is created twice. The write site
says so.

Also threads resistivity_model and lnLambda_form from SLAYERControl into the scan;
it had hardcoded SauterNeoModel/:nrl, so a deck selecting another closure would get
a different eta -- and hence a different delta_s -- from the SLAYER stage.

h5 schema test passes (the metadata contract covers the new datasets).

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

Author: logan-nc

src/GeneralizedPerturbedEquilibrium.jl

@@ -109,10 +109,15 @@ function _layer_overlap_profiles(inputs, dir_path::AbstractString)
                       "($(sprint(showerror, err))); the domain scan is skipped."
                 nothing
             end
+            # Every knob that feeds the layer width comes from the same control the SLAYER stage
+            # uses, so delta_s at a shared surface is identical in Equilibrium/ResistiveLayers/
+            # and Tearing/LayerWidths/ rather than silently differing by the resistivity closure.
             loaded !== nothing && return (profiles=loaded.profiles,
                 chi_perp=loaded.chi_perp === nothing ? sl_ctrl.chi_perp : loaded.chi_perp,
                 chi_tor=loaded.chi_tor === nothing ? sl_ctrl.chi_tor : loaded.chi_tor,
-                mu_i=sl_ctrl.mu_i, zeff=sl_ctrl.zeff)
+                mu_i=sl_ctrl.mu_i, zeff=sl_ctrl.zeff,
+                resistivity_model=Runner._build_resistivity_model(sl_ctrl.resistivity_model),
+                lnLambda_form=sl_ctrl.lnLambda_form)
         end
     end
     if haskey(inputs, "KineticForces")
@@ -126,7 +131,9 @@ function _layer_overlap_profiles(inputs, dir_path::AbstractString)
                     T_e=data.T_e, T_i=data.T_i,
                     omega=(data.omega_E === nothing ? zeros(npsi) : data.omega_E),
                     omega_e=zeros(npsi), omega_i=zeros(npsi))
-                return (profiles=profiles, chi_perp=1.0, chi_tor=1.0, mu_i=2.0, zeff=1.0)
+                return (profiles=profiles, chi_perp=1.0, chi_tor=1.0, mu_i=2.0, zeff=1.0,
+                    resistivity_model=Tearing.InnerLayer.SLAYER.SauterNeoModel(),
+                    lnLambda_form=:nrl)
             catch err
                 err isa InterruptException && rethrow()
                 @warn "psihigh_from_layer_overlap: could not read [KineticForces] kinetic_file " *
@@ -162,29 +169,34 @@ end
 # superset, so its measured knot density is still usable), or the accepted trial for an
 # outward move. `eq_config.psihigh` is left untouched when the scan cannot run.
 function _apply_layer_overlap_psihigh!(eq_config, equil, rerun_input, inputs, intr, ctrl;
-    max_trials::Int=4)
+    enforce::Bool=true, max_trials::Int=4)
+    _rec(; kw...) = merge((scan=nothing, recommended=NaN, psihigh_delta_s=NaN, psihigh_visco=NaN,
+            floor=NaN, applied=eq_config.psihigh, enforced=enforce, reason="not_run"), kw)
     if intr.nlow != intr.nhigh
-        @warn "psihigh_from_layer_overlap is ignored for multi-n runs (nlow=$(intr.nlow), nhigh=$(intr.nhigh)); " *
+        enforce && @warn "psihigh_from_layer_overlap is ignored for multi-n runs (nlow=$(intr.nlow), nhigh=$(intr.nhigh)); " *
               "the overlap criterion depends on which n sets the rational spacing. Keeping psihigh=$(eq_config.psihigh)."
-        return equil
+        return equil, _rec(; reason="multi_n_unsupported")
     end
     kin = _layer_overlap_profiles(inputs, intr.dir_path)
     if kin === nothing
-        @info "psihigh_from_layer_overlap: no kinetic profiles configured ([SLAYER] profile_file or " *
+        enforce && @info "psihigh_from_layer_overlap: no kinetic profiles configured ([SLAYER] profile_file or " *
               "[KineticForces] kinetic_file); keeping psihigh = $(@sprintf("%.6f", eq_config.psihigh))."
-        return equil
+        return equil, _rec(; reason="no_kinetic_profiles")
     end
 
     psihigh_safe = eq_config.psihigh
     scan = Tearing.resistive_layer_overlap(equil, kin.profiles;
         n_tor=intr.nlow, psihigh_safe=psihigh_safe,
         chi_perp=kin.chi_perp, chi_tor=kin.chi_tor,
-        mu_i=kin.mu_i, zeff=kin.zeff)
+        mu_i=kin.mu_i, zeff=kin.zeff,
+        resistivity_model=kin.resistivity_model, lnLambda_form=kin.lnLambda_form)
+    _rec2(; kw...) = _rec(; scan=scan, psihigh_delta_s=something(scan.psihigh_delta_s, NaN),
+        psihigh_visco=something(scan.psihigh_visco, NaN), kw...)
     target = scan.psihigh
     if target === nothing
-        @info "psihigh_from_layer_overlap: no resistive-layer overlap within the scanned range; " *
-              "keeping psihigh = $(@sprintf("%.6f", psihigh_safe))."
-        return equil
+        @info "resistive-layer overlap: no overlap within the scanned range; " *
+              "psihigh = $(@sprintf("%.6f", psihigh_safe)) is unconstrained by layer physics."
+        return equil, _rec2(; reason="no_overlap")
     end
     k = scan.first_overlap
     why = "layers of the q=$(scan.m[k])/$(scan.n[k]) and q=$(scan.m[k-1])/$(scan.n[k-1]) surfaces overlap"
@@ -194,38 +206,49 @@ function _apply_layer_overlap_psihigh!(eq_config, equil, rerun_input, inputs, in
     # more likely wrong than the plasma is that pathological.
     n_kept = count(<=(target), scan.psi)
     if n_kept < 2
-        @warn "psihigh_from_layer_overlap: the criterion asks for psihigh = $(@sprintf("%.6f", target)), " *
+        @warn "resistive-layer overlap: the criterion asks for psihigh = $(@sprintf("%.6f", target)), " *
               "which would keep only $n_kept rational surface(s) ($why). Refusing to truncate — " *
               "check the layer widths. Keeping psihigh = $(@sprintf("%.6f", psihigh_safe))."
-        return equil
+        return equil, _rec2(; recommended=target, reason="too_few_surfaces")
     end
 
     # Lower gate: never cut inside the rational window psiedge sits in.
     fl = _psiedge_floor(scan.psi, ctrl.psiedge, psihigh_safe)
+    floor_val = fl === nothing ? NaN : fl.floor
     if fl !== nothing && target < fl.floor
-        eq_config.psihigh = fl.floor
-        @warn "psihigh_from_layer_overlap: the layer criterion asks for psihigh = " *
+        enforce && (eq_config.psihigh = fl.floor)
+        @warn "resistive-layer overlap: the layer criterion asks for psihigh = " *
               "$(@sprintf("%.6f", target)) ($why), but psiedge = $(@sprintf("%.6f", ctrl.psiedge)) " *
               "sits inside the q=$(scan.m[fl.index])/$(scan.n[fl.index]) window, so the domain is " *
               "held at $(@sprintf("%.6f", fl.floor)) to keep the dW edge scan band open. The layer " *
               "criterion is overridden here — the retained surfaces are set by psiedge, not by physics."
-        return equil
+        return equil, _rec2(; recommended=target, floor=floor_val,
+            applied=enforce ? fl.floor : eq_config.psihigh, reason="psiedge_floor")
     end
 
     # Inward: the pass-1 equilibrium already covers the domain, so nothing to verify.
     if target <= psihigh_safe
-        eq_config.psihigh = target
+        enforce && (eq_config.psihigh = target)
         a_top = Tearing.InnerLayer.SLAYER.surface_minor_radius(equil, psihigh_safe)
         a_cut = Tearing.InnerLayer.SLAYER.surface_minor_radius(equil, target)
         @info "psihigh_from_layer_overlap: $(@sprintf("%.6f", psihigh_safe)) -> $(@sprintf("%.6f", target)) " *
               "($why); removes $(@sprintf("%.2f", 1e3 * (a_top - a_cut))) mm of minor radius, " *
-              "$(@sprintf("%.2f", 100 * (a_top - a_cut) / a_top))% of a, keeping $n_kept of $(length(scan.psi)) surfaces"
-        return equil
+              "$(@sprintf("%.2f", 100 * (a_top - a_cut) / a_top))% of a, keeping $n_kept of $(length(scan.psi)) surfaces" *
+              (enforce ? "" : "  [DIAGNOSTIC ONLY -- psihigh_from_layer_overlap is off, psihigh unchanged]")
+        return equil, _rec2(; recommended=target, floor=floor_val,
+            applied=enforce ? target : eq_config.psihigh,
+            reason=enforce ? "applied" : "not_enforced")
     end
 
     # Outward: the domain has to be formed and checked before it can be trusted. The raw psi(R,Z)
     # map lives on DirectRunInput, which plain `efit` uses too -- pass it so the round-trip
     # tracing residual runs on every direct/EFIT deck, not only on efit_by_inversion.
+    if !enforce
+        @info "resistive-layer overlap: layers stay separated beyond the first-pass bound; the " *
+              "criterion would raise psihigh to $(@sprintf("%.6f", target))  [DIAGNOSTIC ONLY -- " *
+              "psihigh_from_layer_overlap is off, psihigh unchanged]"
+        return equil, _rec2(; recommended=target, floor=floor_val, reason="not_enforced")
+    end
     psi_map = rerun_input isa Equilibrium.DirectRunInput ? rerun_input.psi_in : nothing
     psio_ref = rerun_input isa Equilibrium.DirectRunInput ? rerun_input.psio : nothing
     accepted, accepted_equil = psihigh_safe, equil
@@ -268,7 +291,8 @@ function _apply_layer_overlap_psihigh!(eq_config, equil, rerun_input, inputs, in
         @info "psihigh_from_layer_overlap: $(@sprintf("%.6f", psihigh_safe)) -> $(@sprintf("%.6f", accepted)) " *
               "(layers still separated beyond the first-pass bound; $why)"
     end
-    return accepted_equil
+    return accepted_equil, _rec2(; recommended=target, floor=floor_val, applied=accepted,
+        reason=accepted <= psihigh_safe ? "formation_gate" : "applied")
 end
 
 # Re-run the overlap scan on the final equilibrium and report disagreement with the
@@ -280,7 +304,8 @@ function _verify_layer_overlap_psihigh(eq_config, equil, inputs, intr; rtol::Rea
     scan = Tearing.resistive_layer_overlap(equil, kin.profiles;
         n_tor=intr.nlow, psihigh_safe=eq_config.psihigh,
         chi_perp=kin.chi_perp, chi_tor=kin.chi_tor,
-        mu_i=kin.mu_i, zeff=kin.zeff)
+        mu_i=kin.mu_i, zeff=kin.zeff,
+        resistivity_model=kin.resistivity_model, lnLambda_form=kin.lnLambda_form)
     target = scan.psihigh
     target === nothing && return nothing
     rel = abs(target - eq_config.psihigh) / max(eq_config.psihigh, eps())
@@ -471,10 +496,13 @@ function main_from_inputs(
     # formed there. This runs BEFORE grid construction: setup_equilibrium re-clamps psihigh
     # to the closed-flux region internally, and _validate_psi_nodes hard-errors if psihigh
     # moves after the nodes are built.
+    # The scan ALWAYS runs when kinetic profiles are reachable, so gpec.h5 always records what
+    # the resistive-layer logic would have chosen. Only the enforcement is gated by the flag --
+    # which is what makes a discrepancy between the psihigh used and the criterion (a psiedge
+    # floor, a formation back-off, or simply the flag being off) visible after the fact.
     psihigh_pass1 = eq_config.psihigh
-    if eq_config.psihigh_from_layer_overlap
-        equil = _apply_layer_overlap_psihigh!(eq_config, equil, rerun_input, inputs, intr, ctrl)
-    end
+    equil, layer_overlap = _apply_layer_overlap_psihigh!(eq_config, equil, rerun_input, inputs,
+        intr, ctrl; enforce=eq_config.psihigh_from_layer_overlap)
     psihigh_moved = eq_config.psihigh != psihigh_pass1
 
     # Two-pass auto grid: measure the pass-1 equilibrium's curvature (profiles, geometry,
@@ -747,7 +775,8 @@ function main_from_inputs(
             forcing_modes_snapshot,
             gal_data;
             locstab=locstab,
-            ballooning_boundary=ballooning_boundary
+            ballooning_boundary=ballooning_boundary,
+            layer_overlap=layer_overlap
         )
         @info "Results written to $(ctrl.HDF5_filename)"
     end
@@ -957,7 +986,8 @@ function write_outputs_to_HDF5(
     forcing_modes::Union{Nothing,Vector{ForcingTerms.ForcingMode}}=nothing,
     gal_data::Union{GalerkinResult,Nothing}=nothing;
     locstab::Union{FastInterpolations.CubicSeriesInterpolant,Nothing}=nothing,
-    ballooning_boundary=(psi=Float64[], alpha=Float64[], alpha_critical=Float64[])
+    ballooning_boundary=(psi=Float64[], alpha=Float64[], alpha_critical=Float64[]),
+    layer_overlap=nothing
 )
 
     # Idempotent: already done if a PerturbedEquilibrium stage ran. Leaves the stores empty
@@ -1018,6 +1048,31 @@ function write_outputs_to_HDF5(
             val === nothing && continue
             out_h5["Equilibrium/$(get(EQUIL_H5_NAMES, f, String(f)))"] = val
         end
+        # Domain actually reconstructed, and what the resistive-layer physics asked for. The two
+        # differ whenever the psiedge floor, the formation gate, or the enforcement flag
+        # intervened, so comparing them is how a user checks the domain against the layer logic.
+        #
+        # NOTE: `psihigh_resolved` is written explicitly here only until PR #367 (immutable input
+        # structs) lands -- that PR puts `psihigh_resolved` on EquilibriumParameters, at which
+        # point the EQUIL_H5_NAMES loop above writes it and THIS LINE MUST BE DELETED or the
+        # dataset is created twice.
+        out_h5["Equilibrium/psihigh_resolved"] = equil.config.psihigh
+        out_h5["Equilibrium/psihigh_resistive_layers"] =
+            layer_overlap === nothing ? NaN : Float64(layer_overlap.recommended)
+
+        # Per-surface resistive layer widths from the domain scan: every rational surface it
+        # scored, including any located outside the grid. Written whenever the scan could run,
+        # independently of whether the SLAYER stage ran or the truncation was enforced.
+        if layer_overlap !== nothing && layer_overlap.scan !== nothing
+            sc = layer_overlap.scan
+            out_h5["Equilibrium/ResistiveLayers/rational_psi"] = sc.psi
+            out_h5["Equilibrium/ResistiveLayers/rational_m"] = sc.m
+            out_h5["Equilibrium/ResistiveLayers/rational_n"] = sc.n
+            out_h5["Equilibrium/ResistiveLayers/r_s"] = sc.rs
+            out_h5["Equilibrium/ResistiveLayers/delta_s_abs"] = sc.delta_s_m
+            out_h5["Equilibrium/ResistiveLayers/delta_s_psi"] = sc.width_delta_s
+            out_h5["Equilibrium/ResistiveLayers/extrapolated"] = Int.(sc.extrapolated)
+        end
         out_h5["Equilibrium/psi_total"] = equil.psio
         out_h5["Equilibrium/R_axis"] = equil.ro
         out_h5["Equilibrium/Z_axis"] = equil.zo

src/HDF5Schema.jl

@@ -64,6 +64,37 @@ const MAIN_H5_ANNOTATIONS = [
     # --- Equilibrium/ scalars (written per-field when set; superset listed) ---
     "Equilibrium/R_axis" => (; long_name="R-coordinate of the magnetic axis", units="m"),
     "Equilibrium/Z_axis" => (; long_name="Z-coordinate of the magnetic axis", units="m"),
+    "Equilibrium/psihigh_resolved" =>
+        (;
+            long_name="normalized poloidal flux at the outer edge of the reconstructed equilibrium domain; may differ from the requested [Equilibrium] psihigh via the separatrix clamp or the resistive-layer truncation",
+            units="1"
+        ),
+    "Equilibrium/psihigh_resistive_layers" =>
+        (;
+            long_name="normalized poloidal flux at which adjacent rational surfaces' resistive layers begin to overlap; NaN when not computed (no kinetic profiles, multi-n run, or no overlap in range). Compare with psihigh_resolved to check the domain against the layer physics",
+            units="1"
+        ),
+    "Equilibrium/ResistiveLayers/rational_psi" =>
+        (; long_name="normalized poloidal flux of each rational surface scored by the resistive-layer scan", units="1", dims=("rational",)),
+    "Equilibrium/ResistiveLayers/rational_m" =>
+        (; long_name="poloidal mode number m of each scored rational surface", dims=("rational",)),
+    "Equilibrium/ResistiveLayers/rational_n" =>
+        (; long_name="toroidal mode number n of each scored rational surface", dims=("rational",)),
+    "Equilibrium/ResistiveLayers/r_s" =>
+        (; long_name="minor radius of each scored rational surface", units="m", dims=("rational",)),
+    "Equilibrium/ResistiveLayers/delta_s_abs" =>
+        (; long_name="physical resistive layer thickness |δ_s| at each scored rational surface", units="m", dims=("rational",)),
+    "Equilibrium/ResistiveLayers/delta_s_psi" =>
+        (;
+            long_name="the same layer thickness in normalized flux, |δ_s|/|da/dψ|, so it is directly comparable with the spacing between rational surfaces",
+            units="1",
+            dims=("rational",)
+        ),
+    "Equilibrium/ResistiveLayers/extrapolated" =>
+        (;
+            long_name="flag: this surface was located outside the equilibrium grid (edge q-law for a diverted plasma, or the real q profile for an inverse equilibrium) rather than found within it",
+            dims=("rational",)
+        ),
     "Equilibrium/psi_total" => (; long_name="total poloidal flux difference |ψ_axis − ψ_boundary|", units="Wb/rad"),
     "Equilibrium/R_midplane" => (; long_name="R of the boundary at the inboard and outboard midplane crossings", units="m"),
     "Equilibrium/R_extremum" => (; long_name="R of the boundary at its upper and lower Z extrema", units="m"),
@@ -238,7 +269,7 @@ const MAIN_H5_ANNOTATIONS = [
     "SurfaceGeometries/Plasma/z" => (; long_name="Cartesian z of plasma-surface point cloud", units="m"),
     "SurfaceGeometries/Wall/x" => (; long_name="Cartesian x of wall point cloud", units="m"),
     "SurfaceGeometries/Wall/y" => (; long_name="Cartesian y of wall point cloud", units="m"),
-    "SurfaceGeometries/Wall/z" => (; long_name="Cartesian z of wall point cloud", units="m"),
+    "SurfaceGeometries/Wall/z" => (; long_name="Cartesian z of wall point cloud", units="m")
 ]
 
 # Euler-Lagrange operator matrices: same wording per letter, Ideal/ and Kinetic/ variants.
@@ -251,7 +282,7 @@ const _ELM_IDEAL_LETTERS = [
     ("H", "Euler-Lagrange primitive coefficient matrix H"),
     ("F", "Euler-Lagrange derived coefficient matrix F"),
     ("K", "Euler-Lagrange derived coefficient matrix K"),
-    ("G", "Euler-Lagrange derived coefficient matrix G"),
+    ("G", "Euler-Lagrange derived coefficient matrix G")
 ]
 # The kinetic branch overwrites only A, B, C, K, G and adds f0; D, E, H, F are shared
 # unchanged from the ideal set and are not re-emitted.
@@ -261,14 +292,19 @@ const _ELM_KINETIC_LETTERS = [
     ("C", "Euler-Lagrange primitive coefficient matrix C"),
     ("K", "Euler-Lagrange derived coefficient matrix K"),
     ("G", "Euler-Lagrange derived coefficient matrix G"),
-    ("f0", "raw kinetic component matrix f0"),
+    ("f0", "raw kinetic component matrix f0")
 ]
 const ELM_H5_ANNOTATIONS = vcat(
     ["ForceFreeStates/EulerLagrangeMatrices/psi" => (; long_name="normalized poloidal flux ψ_N grid of the operator matrices", scale="psi")],
-    ["ForceFreeStates/EulerLagrangeMatrices/Ideal/$l" =>
-        (; long_name="ideal " * d, dims=("psi", "mode_row", "mode_col"), attach=(1 => "ForceFreeStates/EulerLagrangeMatrices/psi",)) for (l, d) in _ELM_IDEAL_LETTERS],
-    ["ForceFreeStates/EulerLagrangeMatrices/Kinetic/$l" =>
-        (; long_name="kinetic-modified " * d, dims=("psi", "mode_row", "mode_col"), attach=(1 => "ForceFreeStates/EulerLagrangeMatrices/psi",)) for (l, d) in _ELM_KINETIC_LETTERS]
+    [
+        "ForceFreeStates/EulerLagrangeMatrices/Ideal/$l" =>
+            (; long_name="ideal " * d, dims=("psi", "mode_row", "mode_col"), attach=(1 => "ForceFreeStates/EulerLagrangeMatrices/psi",)) for (l, d) in _ELM_IDEAL_LETTERS
+    ],
+    [
+        "ForceFreeStates/EulerLagrangeMatrices/Kinetic/$l" =>
+            (; long_name="kinetic-modified " * d, dims=("psi", "mode_row", "mode_col"), attach=(1 => "ForceFreeStates/EulerLagrangeMatrices/psi",)) for
+        (l, d) in _ELM_KINETIC_LETTERS
+    ]
 )
 
 """