PE/EQUIL - BUGFIX! - Withhold unvalidated kinetic tearing diagnostics and free the grid they constrained

src/Equilibrium/GridRefinement.jl

@@ -418,7 +418,7 @@ Build the refined pass-2 ψ grid from a formed pass-1 equilibrium: measured-curv
 density (`_knot_density`), equidistribution, a global minimum-spacing floor
 (`enforce_min_spacing`), and rational-surface bracketing (`bracket_mandatory_nodes`). `tau` is
 the target interpolation accuracy (`psi_accuracy`); `kin` optionally supplies kinetic profiles
-whose pedestal gradients attract knots; `mandatory` lists rational-surface ψ values to bracket; `pinned` lists ψ values inserted as plain knots without a cleared zone (kinetic-resonance surfaces);
+whose pedestal gradients attract knots; `mandatory` lists rational-surface ψ values to bracket (`bracket_mandatory = false` skips the bracketing and keeps the fine knots inside instead); `pinned` lists ψ values inserted as plain knots without a cleared zone (kinetic-resonance surfaces);
 `singfac_min` and `n_min` (smallest |n| in the run) set each surface's matching half-stencil
 `dpsi = singfac_min/(n_min·|q′|)`, and the bracket half-width is `bracket_coef·dpsi` (floored at
 `min_spacing`). Rational surfaces are bracketed, not pinned: a knot on the surface would make the
@@ -429,6 +429,7 @@ function refined_psi_grid(equil::PlasmaEquilibrium;
     kin::Union{Nothing,KineticProfileSplines}=nothing,
     mandatory::Vector{Float64}=Float64[],
     pinned::Vector{Float64}=Float64[],
+    bracket_mandatory::Bool=true,
     singfac_min::Float64=1e-4,
     n_min::Int=1,
     bracket_coef::Float64=BRACKET_COEF,
@@ -455,7 +456,13 @@ function refined_psi_grid(equil::PlasmaEquilibrium;
     # pinned node inside a rational's bracket zone is cleared by it (the Δ′ clean-interval
     # requirement wins locally; the rational's own dense floor resolves that neighbourhood).
     grid = isempty(pinned) ? grid : merge_mandatory_nodes(grid, pinned)
-    isempty(mandatory) && return grid
+    # `bracket_mandatory = false` keeps the locally-uniform fine patch around each rational
+    # instead of clearing it: the cleared zone exists solely so the ideal Δ′ stencil sees a
+    # consistent 3rd derivative across the surface, so it is only worth its cost when that Δ′ is
+    # actually delivered (kinetic runs withhold it — see PerturbedEquilibriumControl
+    # `unvalidated_kinetic_tearing`). The knots inside are what resolve the kinetic matrices'
+    # near-rational structure.
+    (isempty(mandatory) || !bracket_mandatory) && return grid
     min_half_widths = [max(bracket_coef * singfac_min / (n_min * abs(equil.profiles.q_deriv(m))), min_spacing) for m in mandatory]
     return bracket_mandatory_nodes(grid, mandatory, min_half_widths, min_spacing)
 end

src/GeneralizedPerturbedEquilibrium.jl

@@ -273,9 +273,19 @@ function main_from_inputs(
             isempty(pinned) ||
                 @info "Pinning $(length(pinned)) kinetic-resonance surfaces into the ψ grid: $(round.(sort(pinned); digits=3))"
         end
+        # Rational bracketing serves the ideal Δ′ stencil. Kinetic runs withhold that Δ′ (see
+        # PerturbedEquilibriumControl `unvalidated_kinetic_tearing`), so unless the user opts back
+        # in, skip the bracketing and let the fine near-rational knots resolve the kinetic
+        # matrices instead. Read from the raw inputs: the PE control struct is built later.
+        pe_raw = get(inputs, "PerturbedEquilibrium", Dict{String,Any}())
+        wants_kinetic_tearing = Bool(get(pe_raw, "unvalidated_kinetic_tearing", false))
+        bracket_mandatory = !(ctrl.kinetic_factor > 0) || wants_kinetic_tearing
+        bracket_mandatory ||
+            @info "Kinetic run: rational-surface bracketing disabled (the ideal Δ′ it serves is withheld); " *
+                  "near-rational knots retained for the kinetic matrices"
         psi_nodes = Equilibrium.refined_psi_grid(equil;
             tau=eq_config.psi_accuracy, kin=kinetic_profiles, mandatory=mandatory, pinned=pinned,
-            singfac_min=ctrl.singfac_min, n_min=n_min)
+            bracket_mandatory=bracket_mandatory, singfac_min=ctrl.singfac_min, n_min=n_min)
         rerun_input = if additional_input !== nothing
             # Analytic *Config, IMAS dd, or prebuilt RunInput — all re-formable. The IMAS
             # path re-runs read_imas, which must resolve the same psihigh both passes;

src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl

@@ -26,8 +26,10 @@ Medium Priority (defer for MWE):
   - `singular_point_method::String` - Method for singular point treatment (default: "standard")
 
 Regularization:
-    # High Priority (MWE)
+# High Priority (MWE)
+
   - `reg_spot::Float64` - Regularization width for singular surface smoothing (default: 0.05). Set to 0 to disable. Must be ≥ 0.
+  - `unvalidated_kinetic_tearing::Bool` - Emit the jump-derived tearing diagnostics (Δ′, resonant current, island widths, Chirikov) in **kinetic** runs. Default `false`: those quantities are withheld because their extraction assumes ideal outer-region asymptotic matching at an ideal rational surface, which kinetic terms invalidate — the surfaces are kinetically displaced (`kinsing`) and the layer response is not the ideal one, so the numbers would be delivered with no validated meaning. Setting `true` restores them for research use, at the user's risk, with a warning. Ideal runs are unaffected. Enabling it also restores rational-surface bracketing of the ψ grid (see `Equilibrium.refined_psi_grid`), which exists to keep the ideal Δ′ stencil consistent.
 """
 @kwdef struct PerturbedEquilibriumControl
     # High Priority (MWE)
@@ -48,6 +50,7 @@ Regularization:
     # Regularization width for singular surface smoothing (matches Fortran gpec.f reg_spot).
     # Set to 0 to disable regularization. Must be non-negative.
     reg_spot::Float64 = 5e-2
+    unvalidated_kinetic_tearing::Bool = false
 end
 
 """
@@ -121,6 +124,7 @@ Metadata [n_rational] — identifies each (surface, n) row:
 
 Control-surface forcing/response spectra [numpert_total], in the three Pharr (2026) field
 representations (all tesla; no flux/weber is stored):
+
   - `forcing_b`/`response_b` - bare normal field b (Σ⁻¹·b̃)
   - `forcing_b_rootarea`/`response_b_rootarea` - root-area-weighted field b̃ (coordinate-invariant)
   - `forcing_b_area`/`response_b_area` - area-weighted field b̄ (= S·b̃; flux is Φ = A·b̄)
@@ -185,18 +189,18 @@ well-conditioned flux-space inductances L, Λ:
     rational_surface_idx::Vector{Int} = Int[]
 
     # Control-surface forcing/response spectra in the three weightings of field representations [numpert_total], tesla
-    forcing_b::Vector{ComplexF64}           = ComplexF64[]  # bare normal field b (forcing Φ_x)
-    forcing_b_rootarea::Vector{ComplexF64}  = ComplexF64[]  # root-area-weighted field b̃ (coordinate-invariant)
-    forcing_b_area::Vector{ComplexF64}      = ComplexF64[]  # area-weighted field b̄
-    response_b::Vector{ComplexF64}          = ComplexF64[]  # bare normal field b (response Φ_tot = P·Φ_x)
+    forcing_b::Vector{ComplexF64} = ComplexF64[]  # bare normal field b (forcing Φ_x)
+    forcing_b_rootarea::Vector{ComplexF64} = ComplexF64[]  # root-area-weighted field b̃ (coordinate-invariant)
+    forcing_b_area::Vector{ComplexF64} = ComplexF64[]  # area-weighted field b̄
+    response_b::Vector{ComplexF64} = ComplexF64[]  # bare normal field b (response Φ_tot = P·Φ_x)
     response_b_rootarea::Vector{ComplexF64} = ComplexF64[]  # root-area-weighted field b̃
-    response_b_area::Vector{ComplexF64}     = ComplexF64[]  # area-weighted field b̄
+    response_b_area::Vector{ComplexF64} = ComplexF64[]  # area-weighted field b̄
 
     # Control surface matrices [numpert_total × numpert_total], root-area-weighted field (b̃) space
-    plasma_inductance::Matrix{ComplexF64}  = zeros(ComplexF64, 0, 0)  # Λ̃ (field space)
+    plasma_inductance::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0)  # Λ̃ (field space)
     surface_inductance::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0)  # L̃ (field space)
-    permeability::Matrix{ComplexF64}       = zeros(ComplexF64, 0, 0)  # P̃ = R⁻¹·Λ·L⁻¹·R
-    reluctance::Matrix{ComplexF64}         = zeros(ComplexF64, 0, 0)  # ϱ̃ = R†·L⁻¹·(Λ−L)·L⁻¹·R
+    permeability::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0)  # P̃ = R⁻¹·Λ·L⁻¹·R
+    reluctance::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0)  # ϱ̃ = R†·L⁻¹·(Λ−L)·L⁻¹·R
     rootarea_to_area_weight::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0)  # S = Σ/√A at psilim: b̃→b̄ recovery operator
     surface_area::Float64 = 0.0  # scalar control-surface area A = ∫J|∇ψ|dθ (flux: Φ = A·b̄; conform R = S·A)
 

src/PerturbedEquilibrium/SingularCoupling.jl

@@ -134,7 +134,7 @@ function _solution_at(
 
     # Same-side candidate nodes around the bracket, trimmed to the 4 nearest psi.
     side = sign(psi - psi_surf)
-    idxs = [j for j in max(1, il-3):min(nstep, ir+3) if sign(odet.psi_store[j] - psi_surf) == side]
+    idxs = [j for j in max(1, il - 3):min(nstep, ir + 3) if sign(odet.psi_store[j] - psi_surf) == side]
     while length(idxs) > 4
         abs(odet.psi_store[idxs[1]] - psi) > abs(odet.psi_store[idxs[end]] - psi) ? popfirst!(idxs) : pop!(idxs)
     end
@@ -307,7 +307,7 @@ function compute_singular_coupling_metrics!(
         state.C_penetrated_area_weighted_field = zeros(ComplexF64, n_rational, numpert_total)
     else
         state.C_penetrated_area_weighted_field = zeros(ComplexF64, 0, 0)
-        @warn "No inner-layer B_pen supplied; penetrated field not computed." maxlog=1
+        @warn "No inner-layer B_pen supplied; penetrated field not computed." maxlog = 1
     end
     state.C_delta_prime = zeros(ComplexF64, n_rational, numpert_total)
     state.rational_psi = zeros(Float64, n_rational)
@@ -493,6 +493,34 @@ function compute_singular_coupling_metrics!(
     # Phase 5: Island diagnostics from applied resonant vectors
     compute_island_diagnostics!(state, n_rational)
 
+    # Kinetic guard: every quantity above descends from the outer-solution jump across the
+    # rational surface, whose reading as a tearing index assumes ideal asymptotic matching at an
+    # ideal rational surface. Kinetic terms displace the surfaces (`kinsing`) and change the
+    # layer response, so that reading is unvalidated — withhold rather than deliver numbers with
+    # no established meaning. Surface metadata (ψ, q, m, n, area) and the inner-layer penetrated
+    # field are kept; they do not come from the jump.
+    if ffit.kinetic_populated && !ctrl.unvalidated_kinetic_tearing
+        state.C_delta_prime = zeros(ComplexF64, 0, 0)
+        state.delta_prime = ComplexF64[]
+        state.C_resonant_current = zeros(ComplexF64, 0, 0)
+        state.resonant_current = ComplexF64[]
+        state.C_resonant_area_weighted_field = zeros(ComplexF64, 0, 0)
+        state.resonant_area_weighted_field = ComplexF64[]
+        state.C_island_width_sq = zeros(ComplexF64, 0, 0)
+        state.island_width_sq = ComplexF64[]
+        state.island_half_width = Float64[]
+        state.chirikov_parameter = Float64[]
+        @warn "Kinetic run: withholding the jump-derived tearing diagnostics (Δ′, resonant current, " *
+              "resonant area-weighted field, island widths, Chirikov). Their extraction assumes ideal " *
+              "outer-region matching at an ideal rational surface, which kinetic terms invalidate; no " *
+              "validated kinetic formulation exists in this code. Set [PerturbedEquilibrium] " *
+              "unvalidated_kinetic_tearing = true to emit them anyway (research use, at your own risk)."
+    elseif ffit.kinetic_populated
+        @warn "Kinetic run with unvalidated_kinetic_tearing = true: Δ′ and the island/current diagnostics " *
+              "derived from the resonant-surface jump have NO validated kinetic formulation. Do not " *
+              "publish these numbers without independent verification."
+    end
+
     if ctrl.verbose
         max_island = isempty(state.island_half_width) ? 0.0 : maximum(state.island_half_width)
         max_dp = isempty(state.delta_prime) ? 0.0 : maximum(abs.(real.(state.delta_prime)))