FFS - NEW FEATURE - Certified adaptive psi grid for the calculated kinetic matrices

src/ForceFreeStates/ForceFreeStatesStructs.jl

@@ -233,7 +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_tol::Float64` - Relative tolerance for the certified adaptive ψ grid of the calculated kinetic matrices. When > 0, the expensive kinetic kernel is evaluated on a seed grid (the ideal coefficient-spline knots) and intervals are refined until the spline-predicted **total** (ideal + kinetic) matrices match fresh evaluations within `kinetic_grid_tol · max|T|` for every element of every consumed family, including the non-Hermitian adjoint combination. `0` (default) evaluates the kernel on every equilibrium knot, today's behaviour. Applies only to `kinetic_source = "calculated"`
+  - `kinetic_grid_tol::Float64` - Relative tolerance for the certified adaptive ψ grid of the calculated kinetic matrices. When > 0, the expensive kinetic kernel is evaluated on a physics-skeleton seed (rational windows, located kinetic-resonance surfaces, and a few interior points per inter-rational span) and intervals are refined until the spline-predicted **total** (ideal + kinetic) matrices match fresh evaluations within `kinetic_grid_tol · max|T|` for every element of every consumed family, including the non-Hermitian adjoint combination. `0` (default) evaluates the kernel on every equilibrium knot, today's behaviour. Applies only to `kinetic_source = "calculated"`
   - `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).
@@ -484,71 +484,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)
@@ -557,11 +532,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

@@ -22,7 +22,7 @@ always contains the seed, so structure resolved by the ideal grid is retained.
 """
 function certified_kinetic_grid(seed::Vector{Float64}, evaluate::Function,
     ideal_scales::NTuple{6,Float64}, tol::Float64, rationals::Vector{Float64};
-    max_rounds::Int=6, verbose::Bool=true)
+    max_rounds::Int=10, verbose::Bool=true)
 
     xs = copy(seed)
     kw, kt = evaluate(xs)
@@ -125,6 +125,10 @@ function certified_kinetic_grid(seed::Vector{Float64}, evaluate::Function,
         kt = newkt[1:row, :, :]
         uncertified = BitVector(newunc)
     end
+    nunc = count(uncertified)
+    nunc > 0 && @warn "Certified kinetic grid: $nunc interval(s) still uncertified after $max_rounds rounds " *
+          "(first near ψ=$(round(xs[findfirst(uncertified)]; digits=4))) -- results may be under-resolved there; " *
+          "raise max_rounds or loosen kinetic_grid_tol"
     @info "Certified kinetic grid: $(length(seed)) seed -> $(length(xs)) knots " *
           "($nevals kernel evaluations, $added refined, tol=$tol)"
     return xs, kw, kt
@@ -149,6 +153,11 @@ Dispatches on `ctrl.kinetic_source`:
     is loaded before KineticForces, so a direct import would invert the
     dependency order.
 
+When `ctrl.kinetic_grid_tol > 0` on the "calculated" path, the kernel is driven over a
+certified adaptive grid instead of the full equilibrium grid; `resonance_psis` supplies the
+located kinetic-resonance surfaces (the same Ω_ℓ = 0 nodes the NTV ψ quadrature panels at)
+as mandatory seed points.
+
 Both paths apply `ctrl.kinetic_factor` as a global scale before the FKG Schur
 reduction.
 """
@@ -158,7 +167,8 @@ function make_kinetic_matrix(
     ffit::FourFitVars,
     intr::ForceFreeStatesInternal,
     metric::MetricData;
-    calculated_source::Union{Nothing,Function}=nothing
+    calculated_source::Union{Nothing,Function}=nothing,
+    resonance_psis::Vector{Float64}=Float64[]
 )
     xs = metric.xs
     mpsi = length(xs)
@@ -174,28 +184,44 @@ function make_kinetic_matrix(
             "`calculated_source=KineticForces.compute_calculated_kinetic_matrices` explicitly."
         )
         if ctrl.kinetic_grid_tol > 0
-            # Certified adaptive grid: seed with the ideal coefficient-spline knots and refine
-            # until the total matrices are spline-predictable to kinetic_grid_tol everywhere.
-            # Seed with a coarse skeleton of the ideal coefficient-spline knots: the kernel is the
-            # expensive part, so the seed must not scale with the equilibrium grid. Every ~seed
-            # gap the certificate cannot vouch for is refined, so coarse seeding trades cheap
-            # certificates for expensive blanket evaluation. Endpoints and rational-window knots
-            # are always retained.
+            # Certified adaptive grid: seed with the physics skeleton and refine until the
+            # total matrices are spline-predictable to kinetic_grid_tol everywhere. The seed is
+            # the same structure the NTV ψ quadrature panels at: the rational windows (retained
+            # from the ideal coefficient grid), the located kinetic-resonance surfaces, and a few
+            # interior points per inter-rational span (log-spaced in the axis span to match the
+            # core's Frobenius power-law structure). The kernel is the expensive part, so the
+            # seed must not scale with the equilibrium grid; every gap the certificate cannot
+            # vouch for is refined, so coarse seeding trades cheap certificates for expensive
+            # blanket evaluation. Seed adequacy is not correctness-critical (certification is) --
+            # the seed's job is to pre-place knots where sharp structure is expected so it cannot
+            # alias between midpoints.
             base = isempty(ffit.matrix_xs) ? metric.xs : ffit.matrix_xs
-            seed_max = 49
-            if length(base) > seed_max
-                rats = [sng.psifac for sng in intr.sing]
-                keep_idx = falses(length(base))
-                keep_idx[1] = keep_idx[end] = true
-                stride = max(1, (length(base) - 1) ÷ (seed_max - 1))
-                keep_idx[1:stride:end] .= true
-                for (i, x) in pairs(base)
-                    any(abs(x - r) <= Equilibrium.RATIONAL_RES_RADIUS for r in rats) && (keep_idx[i] = true)
+            rats = [sng.psifac for sng in intr.sing]
+            lo, hi = base[1], base[end]
+            seed = [lo, hi]
+            for x in base   # rational-window knots: the anti-aliasing floor where layers live
+                any(abs(x - r) <= Equilibrium.RATIONAL_RES_RADIUS for r in rats) && push!(seed, x)
+            end
+            append!(seed, filter(p -> lo < p < hi, resonance_psis))
+            anchors = sort!(filter(r -> lo < r < hi, unique(rats)))
+            spans = [lo; anchors; hi]
+            nbetween = 4   # interior points per span: cubic-spline support everywhere before round one
+            for i in 1:length(spans)-1
+                a, b = spans[i], spans[i+1]
+                if i == 1 && b > 4 * a
+                    append!(seed, exp.(range(log(a), log(b); length=nbetween + 2))[2:end-1])
+                else
+                    append!(seed, range(a, b; length=nbetween + 2)[2:end-1])
                 end
-                seed = base[keep_idx]
-            else
-                seed = copy(base)
             end
+            sort!(seed)
+            # Merge near-coincident points (e.g. a resonance on a rational) but never the endpoints.
+            merged = [seed[1]]
+            for x in seed[2:end]
+                x - merged[end] > Equilibrium.RATIONAL_RES_SPACING / 8 && push!(merged, x)
+            end
+            merged[end] = hi
+            seed = merged
             hint = Ref(1)
             ideal_scales = ntuple(ic -> begin
                     sp = (ffit.amats, ffit.bmats, ffit.cmats, ffit.dmats_prim, ffit.emats_prim, ffit.hmats)[ic]

src/GeneralizedPerturbedEquilibrium.jl

@@ -423,8 +423,20 @@ function main_from_inputs(
             KineticForces.compute_calculated_kinetic_matrices(
                 c, e, i, m, f;
                 kf_ctrl=kf_ctrl, kinetic_profiles=kinetic_profiles, psis=psis)
+        # Locate the kinetic-resonance surfaces (Ω_ℓ = 0) that seed the certified matrix
+        # grid -- the same nodes the NTV ψ quadrature panels at (one source of truth).
+        resonance_psis = Float64[]
+        if ctrl.kinetic_source == "calculated" && ctrl.kinetic_grid_tol > 0
+            for n_res in intr.nlow:intr.nhigh
+                n_res == 0 && continue
+                append!(resonance_psis, KineticForces.kinetic_resonance_psi_nodes(
+                    kinetic_profiles, equil;
+                    n=n_res, nl=kf_ctrl.nl, zi=kf_ctrl.zi, mi=kf_ctrl.mi,
+                    electron=kf_ctrl.electron, wdfac=kf_ctrl.wdfac))
+            end
+        end
         make_kinetic_matrix(ctrl, equil, ffit, intr, metric;
-            calculated_source=calculated_cb)
+            calculated_source=calculated_cb, resonance_psis=resonance_psis)
 
         # Find kinetically-displaced singular surfaces (zeros of det(F̄)) for ODE crossings.
         # Matches Fortran ksing_find (sing.f:1486-1616). singfac_min > 0 gates crossings;