From aefa3cf1ff153c4c5ecc69ae2b6e81a9f25901ca Mon Sep 17 00:00:00 2001 From: logan-nc Date: Tue, 18 Aug 2026 08:44:07 -0400 Subject: [PATCH 1/3] FFS - IMPROVEMENT - Cap the core knot density of the EL coefficient splines The Euler-Lagrange coefficient splines inherited every knot of the equilibrium grid. A cubic spline's third-derivative jumps at knots scale as (node error)/dpsi^3, so the equilibrium's near-axis packing (dpsi ~ 1e-6) amplifies even tolerance-level node error into huge C2 kinks, and the adaptive integrator's step size becomes slaved to the knot spacing -- measured directly: core jump magnitudes grow ~34x per mpsi doubling while a (tol/J)^(1/4) step model reproduces the observed step-count ladder. The coefficients are near-cylindrical in the core and do not need that packing. Build their splines on a subset of the equilibrium grid with core density capped at dpsi >= 0.05*psi below psi = 0.1; node values are unchanged, only knot density. Measured on DIII-D stripped decks (route-a base, mpsi 512/1024): accepted EL steps 2768 -> 2188 and 4403 -> 3034, per-doubling growth 1.59x -> 1.39x, warm run -19% at mpsi=512, with et[1] unchanged to 3e-8 relative and the Riccati BVP Delta-prime diagonal unchanged to 4e-7 on all five surfaces. Interim fixed-cap form; a matrix-curvature knot selection is planned to replace the fixed rule, with this commit as the fallback. Co-Authored-By: Claude Opus 5 (1M context) --- src/ForceFreeStates/Fourfit.jl | 39 +++++++++++++++++++++++----------- 1 file changed, 27 insertions(+), 12 deletions(-) diff --git a/src/ForceFreeStates/Fourfit.jl b/src/ForceFreeStates/Fourfit.jl index 1ce7fde50..3eef94da5 100644 --- a/src/ForceFreeStates/Fourfit.jl +++ b/src/ForceFreeStates/Fourfit.jl @@ -496,23 +496,38 @@ function make_matrix(equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStates # FastInterpolations now natively supports complex values - no need to split real/imag # Create complex series interpolants with per-column extrap BC - ffit.amats = cubic_interp(metric.xs, Series(amats_flat); ffit.itp_opts...) - ffit.bmats = cubic_interp(metric.xs, Series(bmats_flat); ffit.itp_opts...) - ffit.cmats = cubic_interp(metric.xs, Series(cmats_flat); ffit.itp_opts...) - ffit.dmats_prim = cubic_interp(metric.xs, Series(dmats_flat); ffit.itp_opts...) - ffit.emats_prim = cubic_interp(metric.xs, Series(emats_flat); ffit.itp_opts...) - ffit.hmats = cubic_interp(metric.xs, Series(hmats_flat); ffit.itp_opts...) - ffit.fmats_lower = cubic_interp(metric.xs, Series(fmats_lower_flat); ffit.itp_opts...) - ffit.fmats_prim = cubic_interp(metric.xs, Series(fmats_prim_flat); ffit.itp_opts...) - ffit.fmats_gal = cubic_interp(metric.xs, Series(fmats_gal_flat); ffit.itp_opts...) - ffit.gmats = cubic_interp(metric.xs, Series(gmats_flat); ffit.itp_opts...) - ffit.kmats = cubic_interp(metric.xs, Series(kmats_flat); ffit.itp_opts...) + # Decouple the coefficient-spline knots from the equilibrium grid in the packed core: + # a cubic spline's third-derivative jumps scale as (node error)/dpsi^3, so equilibrium-grade + # core packing amplifies tolerance-level node error into jumps that slave the EL step size. + # The coefficients are near-cylindrical there; cap the density at dpsi >= 0.05*psi below psi=0.1. + keep = Int[1] + for i in 2:length(metric.xs)-1 + x = metric.xs[i] + if x >= 0.1 || (x - metric.xs[keep[end]]) >= 0.05 * x + push!(keep, i) + end + end + push!(keep, length(metric.xs)) + mxs = metric.xs[keep] + @info "EL coefficient-spline grid: $(length(metric.xs)) -> $(length(mxs)) knots after core density cap" + + ffit.amats = cubic_interp(mxs, Series(amats_flat[keep, :]); ffit.itp_opts...) + ffit.bmats = cubic_interp(mxs, Series(bmats_flat[keep, :]); ffit.itp_opts...) + ffit.cmats = cubic_interp(mxs, Series(cmats_flat[keep, :]); ffit.itp_opts...) + ffit.dmats_prim = cubic_interp(mxs, Series(dmats_flat[keep, :]); ffit.itp_opts...) + ffit.emats_prim = cubic_interp(mxs, Series(emats_flat[keep, :]); ffit.itp_opts...) + ffit.hmats = cubic_interp(mxs, Series(hmats_flat[keep, :]); ffit.itp_opts...) + ffit.fmats_lower = cubic_interp(mxs, Series(fmats_lower_flat[keep, :]); ffit.itp_opts...) + ffit.fmats_prim = cubic_interp(mxs, Series(fmats_prim_flat[keep, :]); ffit.itp_opts...) + ffit.fmats_gal = cubic_interp(mxs, Series(fmats_gal_flat[keep, :]); ffit.itp_opts...) + ffit.gmats = cubic_interp(mxs, Series(gmats_flat[keep, :]); ffit.itp_opts...) + ffit.kmats = cubic_interp(mxs, Series(kmats_flat[keep, :]); ffit.itp_opts...) # TODO: set powers # Do we need this yet? Only called if power_flag = true # Jacobian Fourier band ψ-spline, used for the power normalization in Free.jl - ffit.jmats = cubic_interp(metric.xs, Series(jmats_flat); ffit.itp_opts...) + ffit.jmats = cubic_interp(mxs, Series(jmats_flat[keep, :]); ffit.itp_opts...) return ffit end From a737b3d9d25218593c143e1ac25eb658f5bb5b61 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Tue, 18 Aug 2026 15:20:14 -0400 Subject: [PATCH 2/3] FFS - MINOR - Log the coefficient-spline grid only when the cap removes knots On the production two-pass auto grid the cap is a no-op (its core spacing already satisfies the density rule), so the unconditional message was noise. Co-Authored-By: Claude Opus 5 (1M context) --- src/ForceFreeStates/Fourfit.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/ForceFreeStates/Fourfit.jl b/src/ForceFreeStates/Fourfit.jl index 3eef94da5..928a25866 100644 --- a/src/ForceFreeStates/Fourfit.jl +++ b/src/ForceFreeStates/Fourfit.jl @@ -509,7 +509,8 @@ function make_matrix(equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStates end push!(keep, length(metric.xs)) mxs = metric.xs[keep] - @info "EL coefficient-spline grid: $(length(metric.xs)) -> $(length(mxs)) knots after core density cap" + length(mxs) < length(metric.xs) && + @info "EL coefficient-spline grid: $(length(metric.xs)) -> $(length(mxs)) knots after core density cap" ffit.amats = cubic_interp(mxs, Series(amats_flat[keep, :]); ffit.itp_opts...) ffit.bmats = cubic_interp(mxs, Series(bmats_flat[keep, :]); ffit.itp_opts...) From 0c17097b5ce9cddec9a83183f7627f14653baded Mon Sep 17 00:00:00 2001 From: logan-nc Date: Tue, 18 Aug 2026 22:16:53 -0400 Subject: [PATCH 3/3] FFS - IMPROVEMENT - Ground the density cap in the Frobenius scale and protect rationals Reframes the cap as log-uniform sampling of the Frobenius region: near the axis every component is a power law in psi, and cubic interpolation of psi^p on a log-uniform grid errs by ~(p*c)^4/384, so c = 0.05 resolves even the steepest spectrum component (p = mmax/2) to ~2e-4 while the physics responds far below that. Two generalization guards, motivated by cross-equilibrium testing: the capped region now ends at the innermost rational surface when that sits inside psi = 0.1, and no knot inside a rational's RATIONAL_RES_RADIUS window is ever removed -- preserving the Delta'-stencil structure for decks (e.g. higher n) whose rationals reach the core. Both guards are no-ops on every current case, verified: DIII-D m512 and Solovev m512 reproduce the previous cap's step counts and knot sets exactly. Cross-equilibrium check of the rule itself: tj_analytic_direct m1024 (analytic, traced) 1470 -> 959 steps at et[1] 1.3e-9; LAR m1024 (inversion path, clean geometry) 951 -> 875 at et[1] identical to 8 digits; Solovev ldp m512/m1024 ~unchanged steps at ~4.5e-7 absolute et[1] shift (a +-10.4 cancellation amplifies this to 3e-5 relative). Co-Authored-By: Claude Opus 5 (1M context) --- src/ForceFreeStates/Fourfit.jl | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/src/ForceFreeStates/Fourfit.jl b/src/ForceFreeStates/Fourfit.jl index 928a25866..c7fcd2988 100644 --- a/src/ForceFreeStates/Fourfit.jl +++ b/src/ForceFreeStates/Fourfit.jl @@ -499,11 +499,22 @@ function make_matrix(equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStates # Decouple the coefficient-spline knots from the equilibrium grid in the packed core: # a cubic spline's third-derivative jumps scale as (node error)/dpsi^3, so equilibrium-grade # core packing amplifies tolerance-level node error into jumps that slave the EL step size. - # The coefficients are near-cylindrical there; cap the density at dpsi >= 0.05*psi below psi=0.1. + # Near the axis every component is a Frobenius power law in psi, and power laws are scale-free, + # so log-uniform sampling (dpsi >= c*psi) resolves them at constant relative accuracy: cubic + # interpolation of psi^p on that grid errs by ~(p*c)^4/384, so c = 0.05 resolves even the + # steepest spectrum component (p = mmax/2) to ~2e-4 while the physics responds far below that + # (the steep components carry vanishing solution amplitude). The capped region ends at the + # innermost rational surface (or psi = 0.1, whichever is smaller) and never removes a knot + # inside a rational's resolution window, preserving the Delta'-stencil structure the + # equilibrium grid encodes (GridRefinement RATIONAL_RES_RADIUS). + cap_edge = 0.1 + rationals = [s.psifac for s in intr.sing] + isempty(rationals) || (cap_edge = min(cap_edge, minimum(rationals) - Equilibrium.RATIONAL_RES_RADIUS)) + in_rational_window(x) = any(abs(x - r) <= Equilibrium.RATIONAL_RES_RADIUS for r in rationals) keep = Int[1] for i in 2:length(metric.xs)-1 x = metric.xs[i] - if x >= 0.1 || (x - metric.xs[keep[end]]) >= 0.05 * x + if x >= cap_edge || in_rational_window(x) || (x - metric.xs[keep[end]]) >= 0.05 * x push!(keep, i) end end