From 35e4ae5939a32af1c2effe9a5aee39ceeb2ff983 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Wed, 19 Aug 2026 08:35:35 -0400 Subject: [PATCH 1/2] FFS - IMPROVEMENT - Seed the certified kinetic grid with the resonance-panel physics skeleton Replace the arbitrary constant-size decimation seed with the same structure the NTV psi quadrature panels at: rational windows, located kinetic-resonance surfaces (Omega_l = 0 via KineticForces.kinetic_resonance_psi_nodes -- one source of truth with the torque quadrature), and four interior points per inter-rational span (log-spaced in the axis span to match the core Frobenius structure). Warn loudly when certification exhausts max_rounds instead of returning silently. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01LzbLFQKyuRE5DYZmLokKmk --- src/ForceFreeStates/Kinetic.jl | 66 ++++++++++++++++++-------- src/GeneralizedPerturbedEquilibrium.jl | 14 +++++- 2 files changed, 59 insertions(+), 21 deletions(-) diff --git a/src/ForceFreeStates/Kinetic.jl b/src/ForceFreeStates/Kinetic.jl index 7577e867c..344c966c4 100644 --- a/src/ForceFreeStates/Kinetic.jl +++ b/src/ForceFreeStates/Kinetic.jl @@ -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] diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 39a0d0efa..e25e7fc30 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/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; From bef73a6bad6f3bcd0c7c65d1e05f8655dabf0d4b Mon Sep 17 00:00:00 2001 From: logan-nc Date: Wed, 19 Aug 2026 09:44:40 -0400 Subject: [PATCH 2/2] FFS - MINOR - Raise certification round cap and fix stale seed docstring DIII-D auto-grid runs exhaust 6 rounds in the core spans (warned as designed); 10 rounds lets the certificate reach the spacing floors. Docstring now describes the physics-skeleton seed instead of the removed decimation seed. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01LzbLFQKyuRE5DYZmLokKmk --- src/ForceFreeStates/ForceFreeStatesStructs.jl | 30 +------------------ src/ForceFreeStates/Kinetic.jl | 2 +- 2 files changed, 2 insertions(+), 30 deletions(-) diff --git a/src/ForceFreeStates/ForceFreeStatesStructs.jl b/src/ForceFreeStates/ForceFreeStatesStructs.jl index b493ceab0..9e8a552ee 100644 --- a/src/ForceFreeStates/ForceFreeStatesStructs.jl +++ b/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 diff --git a/src/ForceFreeStates/Kinetic.jl b/src/ForceFreeStates/Kinetic.jl index 344c966c4..599e61419 100644 --- a/src/ForceFreeStates/Kinetic.jl +++ b/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)