FFS - EXPERIMENT - One-round bisection refinement of the calculated kinetic matrix grid

User-proposed safety net, trial implementation: kernel-evaluate every interval
midpoint above the validity band (one extra threaded kernel pass), compare
against the increment splines with a per-interval LOCAL scale (the lesson from
the parked global-scale certificate), insert missed midpoints as knots reusing
their kernel values, and log the residual map unconditionally. Knob
kinetic_grid_bisect (default 0 = off).

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01LzbLFQKyuRE5DYZmLokKmk

Author: logan-nc

src/ForceFreeStates/ForceFreeStatesStructs.jl

@@ -233,6 +233,7 @@ gpec.toml.
   - `numsteps_init::Int` - Initial array size for ODE data storage
   - `numunorms_init::Int` - Initial array size for solution normalization data
   - `singfac_min::Float64` - Fractional distance from rational q at which ideal jump condition is enforced
+  - `kinetic_grid_bisect::Float64` - EXPERIMENTAL one-round bisection refinement of the calculated kinetic matrix grid. Every interval midpoint above the near-axis validity band is kernel-evaluated; the measured values are compared against the increment splines with a per-interval local scale, and midpoints missed by more than this relative tolerance become knots (their kernel values are reused, so refinement costs exactly one extra kernel pass). A residual map is logged regardless of insertions. `0` (default) disables
   - `set_psilim_via_dmlim::Bool` - Truncate the integration domain at `(last_rational_q + dmlim) / n` rather than at `qhigh` / `psihigh`. Fortran STRIDE found that truncating ~20 % above the outermost rational (`dmlim = 0.2`) avoids a numerical kink instability in δW that appears when the integration ends too close to or just below a rational surface. **For diverted equilibria where q → ∞ at the separatrix** (e.g. DIII-D geqdsks, the bulk of production use) this costs negligible physical domain because rationals get arbitrarily dense near the LCFS — `set_psilim_via_dmlim = true` is the safe and recommended default. **For limited circular / analytical equilibria with finite q at the edge** (Solovev, LAR scans), rationals are sparse and 20 % above the last rational chops off too much edge, so set `set_psilim_via_dmlim = false` and let `qhigh` / `psihigh` control the truncation. Multi-`n` runs are not supported by this truncation (the "outermost rational + dmlim / n" depends on which `n`); when `set_psilim_via_dmlim = true` with `nn_low != nn_high`, `sing_lim!` warns and falls back to `qhigh` / `psihigh`. Default `true`.
   - `dmlim::Float64` - Distance beyond last rational surface (normalised ∈ [0,1) in units of 1/n). Only used when `set_psilim_via_dmlim` is true. Fortran STRIDE convention is 0.2 (truncate 20 % of one rational-surface spacing above the last surface), retained here.
   - `sing_order::Int` - Order of singular layer (Frobenius) expansion at rational surfaces. Default 6 (Fortran STRIDE convention for Δ' calculations; lower values trade accuracy for speed).
@@ -270,6 +271,7 @@ gpec.toml.
     numsteps_init::Int = 4000
     numunorms_init::Int = 100
     singfac_min::Float64 = 1e-4   # Matches Fortran STRIDE; required nonzero for the Riccati path.
+    kinetic_grid_bisect::Float64 = 0.0  # EXPERIMENTAL: one bisection-refinement round of the calculated kinetic matrix grid; midpoints of every interval above the validity band are kernel-evaluated and inserted where the spline misses them by more than this relative tolerance (local scale). 0 = off
     set_psilim_via_dmlim::Bool = true   # Safe default for diverted equilibria (most production use); set false for limited/analytical (LAR, Solovev). Auto-skipped for multi-n. See docstring.
     dmlim::Float64 = 0.2
     sing_order::Int = 6
@@ -481,71 +483,46 @@ and a small set of temporary matrices and factors used to compute singular-layer
   - `numpert_total::Int` - Total number of Fourier mode combinations (m × n) used in the calculation.
 
   - `numunorms_init::Int` - Initial allocation size for the number of normalization operations recorded.
-
   - `msing::Int` - Number of singular surfaces in the equilibrium (used to size asymptotic coefficient arrays).
-
   - `numsteps_init::Int` - Initial allocation size for the number of integration steps to store.
-
   - `step::Int` - Current integration step index (1-based, like `istep` in the original Fortran).
-
   - `psi_store::Vector{Float64}` - Stored psi values at each saved integration step (length `numsteps_init`).
-
   - `q_store::Vector{Float64}` - Stored q values at each saved integration step (length `numsteps_init`).
-
   - `u_store::Array{ComplexF64,4}` - Stored solution arrays at each saved step with shape
     `(numpert_total, numpert_total, 2, numsteps_init)` (complex solution state used by the solver).
-
   - `du_store::Array{ComplexF64,3}` - dΞ_ψ/dψ (the u₁ block only) at each saved step, shape
     `(numpert_total, numpert_total, step)`. Empty until `materialize_derivative_stores!` fills it,
     except on the galerkin-matched path which supplies the analytic derivative at construction.
     du₂/dψ is never stored densely — its only consumer evaluates it on demand at bracket nodes.
-
   - `xi_s_store::Array{ComplexF64,3}` - Clebsch displacement Ξ_s at each saved step, eq. 18 of Glasser 2016,
     shape `(numpert_total, numpert_total, step)`. Empty until materialized, same as `du_store`.
-
   - `u_store_el_basis::Bool` - True when `u_store` holds the Euler-Lagrange state `(u₁, u₂)`, so the
     derivative kernel can be re-applied to it. False on the sparse parallel path, whose stored columns
     are chunk-endpoint Riccati matrices; `materialize_derivative_stores!` refuses to run there.
-
   - `du_store_populated::Bool` - True once `du_store`/`xi_s_store` hold valid data in the final
     (post-transform, post-normalization) basis. Set by `materialize_derivative_stores!` or by the
     galerkin-matched constructor; stays false where the stores cannot be materialized, e.g. the
     sparse parallel path whose solution is in the Riccati basis.
-
   - `crit_store::Vector{Float64}` - Stored crit parameter values (smallest eigenvalue of W⁻ꜝ) (length `numsteps_init`).
-
   - `ca_r::Array{ComplexF64,4}` - Asymptotic coefficients just to the right of each singular surface
     with shape `(numpert_total, numpert_total, 2, msing)`.
-
   - `ca_l::Array{ComplexF64,4}` - Asymptotic coefficients just to the left of each singular surface
     with shape `(numpert_total, numpert_total, 2, msing)`.
-
   - `edge_scan::EdgeScanState` - Edge dW scan state and results. Initialized as a disabled sentinel (N_edge=0) and replaced by `findmax_dW_edge!` when a scan runs.
-
   - `psifac::Float64` - Current normalized flux coordinate for the integrator.
-
   - `q::Float64` - Safety factor value at `psifac` (current q during integration).
-
   - `u::Array{ComplexF64,3}` - Current working solution arrays with shape `(numpert_total, numpert_total, 2)`.
-
   - `ising_start::Int` - Index of the starting singular surface to be crossed during integration.
-
   - `psimax::Float64` - Maximum psi value for which the integrator is allowed to run in next integration region.
-
   - `needs_crossing::Bool` - Flag indicating whether a rational surface needs to be crossed after the current integration region.
-
   - `nzero::Int` - Count of detected zero crossings (used for diagnostics).
-
   - `new::Bool` - Flag indicating whether a new `unorm0` should be computed after a fixup.
 
     # Initialization parameters
-
   - `unorm::Vector{Float64}` - Current norms of the solution vectors (length `numpert_total`).
-
   - `unorm0::Vector{Float64}` - Reference/initial norms of the solution vectors (length `numpert_total`).
 
     # Saved data throughout integration
-
   - `ifix::Int` - Number of normalization operations performed (index into normalization arrays).
 
 # Total ODE solver steps taken (all steps, not just saved ones)
@@ -554,11 +531,8 @@ and a small set of temporary matrices and factors used to compute singular-layer
 
   - `sing_flag::Vector{Bool}` - Boolean flags indicating which stored normalizations correspond to singular solutions    # Edge dW scan state and results (disabled sentinel when psiedge >= psilim, i.e. no edge scan)
     (length `numunorms_init`).
-
   - `zeroed_idx::Vector{Vector{Int}}` - For each ideal rational surface jump, a vector of indices of solutions that were zeroed.    # Data for integrator
-
   - `fixfac::Array{ComplexF64,3}` - Fix-up factors for Gaussian reduction with shape `(numpert_total, numpert_total, numunorms_init)`.
-
   - `fixstep::Vector{Int64}` - Step indices (psi step positions) at which normalization/fixups were performed (length `numunorms_init`).
 """
 @kwdef mutable struct OdeState

src/ForceFreeStates/Kinetic.jl

@@ -1,3 +1,75 @@
+
+"""
+    bisect_refine_kinetic_grid(xs, kw, kt, evaluate, tol, psi_c) → (xs, kw, kt)
+
+One bounded bisection-refinement round for the calculated kinetic matrix grid (experimental).
+Every interval midpoint above `2·psi_c` (outside the near-axis validity band) is kernel-evaluated
+in one threaded batch; for each family the increment spline on the current knots predicts the
+midpoint, and the max-element residual is normalized by a **per-interval local scale** (the max
+increment magnitude over the interval's endpoints and midpoint, floored at 1% of the family's
+global maximum — the local scale is the lesson from the parked global-scale certificate, the
+floor keeps near-zero crossings from flagging on noise). Midpoints whose worst-family residual
+exceeds `tol` become knots; their kernel values are reused, so the round costs exactly one extra
+kernel pass regardless of how many knots it inserts. The residual distribution is logged
+unconditionally — the map of what the ideal-driven grid misses is the diagnostic even when
+nothing is inserted.
+"""
+function bisect_refine_kinetic_grid(xs::Vector{Float64}, kw::Array{ComplexF64,3}, kt::Array{ComplexF64,3},
+    evaluate::Function, tol::Float64, psi_c::Float64)
+    lo = 2 * psi_c
+    cand = [i for i in 1:(length(xs)-1) if xs[i] >= lo]
+    isempty(cand) && return xs, kw, kt
+    mids = [0.5 * (xs[i] + xs[i+1]) for i in cand]
+    kw_m, kt_m = evaluate(mids)
+    np2 = size(kw, 2)
+    worst = zeros(length(mids))
+    buf = Vector{ComplexF64}(undef, np2)
+    for ic in 1:6
+        gscale_w = maximum(abs, @view(kw[:, :, ic]))
+        gscale_t = maximum(abs, @view(kt[:, :, ic]))
+        sp_w = cubic_interp(xs, Series(@view(kw[:, :, ic])))
+        sp_t = cubic_interp(xs, Series(@view(kt[:, :, ic])))
+        for (k, m) in pairs(mids)
+            i = cand[k]
+            for (sp, arr, arrm, gs) in ((sp_w, kw, kw_m, gscale_w), (sp_t, kt, kt_m, gscale_t))
+                gs == 0 && continue
+                sp(buf, m)
+                scale = max(maximum(abs, @view(arr[i, :, ic])), maximum(abs, @view(arr[i+1, :, ic])),
+                    maximum(abs, @view(arrm[k, :, ic])), 0.01 * gs)
+                r = maximum(abs(buf[j] - arrm[k, j, ic]) for j in 1:np2) / scale
+                worst[k] = max(worst[k], r)
+            end
+        end
+    end
+    q = [round(sort(worst)[max(1, ceil(Int, f * length(worst)))]; sigdigits=2) for f in (0.5, 0.9, 0.99, 1.0)]
+    ins = findall(>(tol), worst)
+    @info "Kinetic bisection map: $(length(mids)) midpoints probed, residual quantiles " *
+          "(50/90/99/100%) = $q; $(length(ins)) exceed tol=$tol" *
+          (isempty(ins) ? "" : " at ψ=$(round.(mids[ins]; digits=4))")
+    isempty(ins) && return xs, kw, kt
+    keep = sort(ins)
+    nnew = length(xs) + length(keep)
+    xs2 = Vector{Float64}(undef, nnew)
+    kw2 = Array{ComplexF64,3}(undef, nnew, np2, 6)
+    kt2 = Array{ComplexF64,3}(undef, nnew, np2, 6)
+    row = 0
+    ki = 1
+    for i in eachindex(xs)
+        row += 1
+        xs2[row] = xs[i]
+        kw2[row, :, :] .= kw[i, :, :]
+        kt2[row, :, :] .= kt[i, :, :]
+        while ki <= length(keep) && cand[keep[ki]] == i
+            row += 1
+            xs2[row] = mids[keep[ki]]
+            kw2[row, :, :] .= kw_m[keep[ki], :, :]
+            kt2[row, :, :] .= kt_m[keep[ki], :, :]
+            ki += 1
+        end
+    end
+    return xs2, kw2, kt2
+end
+
 """
     make_kinetic_matrix(ctrl, equil, ffit, intr, metric;
                         calculated_source=nothing)
@@ -57,6 +129,13 @@ function make_kinetic_matrix(
             mpsi = length(xs)
             kw_flat, kt_flat = calculated_source(ctrl, equil, intr, metric, ffit; psis=xs)
         end
+        if ctrl.kinetic_grid_bisect > 0
+            xs, kw_flat, kt_flat = bisect_refine_kinetic_grid(
+                collect(xs), kw_flat, kt_flat,
+                psis -> calculated_source(ctrl, equil, intr, metric, ffit; psis=psis),
+                ctrl.kinetic_grid_bisect, axis_validity_psi_c)
+            mpsi = length(xs)
+        end
         kw_flat .*= ctrl.kinetic_factor
         kt_flat .*= ctrl.kinetic_factor
     else