From de0db4f551b4774df959109e146dfffd6a5532b3 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Sun, 16 Aug 2026 18:28:59 -0400 Subject: [PATCH 1/4] EQUIL - IMPROVEMENT - Sample every flux surface at the same straight-fieldline angles Each surface was traced independently and then splined on that surface's OWN solver-chosen abscissae before being resampled onto the common theta grid, so the resample error was uncorrelated between neighbouring surfaces -- white noise in psi that grid refinement amplifies rather than reduces. Neither reltol nor abstol touched it, because it is remap interpolation error rather than integration error. The trace now returns its dense solution, and equilibrium_solver root-solves (Brent, bracketed by the monotone jac-weighted flux integral) for the angle at which the normalised straight-fieldline angle reaches each target node, evaluating there. Every surface is sampled at identical abscissae and the resample error at the output nodes is zero. The arclength tracer returns nothing for the solution and keeps the previous path. Measured on DIII-D stripped decks at eulerlagrange_tolerance 1e-10, accepted Euler-Lagrange steps fall 2309/3977/7638 -> 1881/2768/4403 for mpsi 256/512/1024, a 42% reduction at mpsi=1024, and the per-doubling growth drops from 1.92x to 1.59x. Surface geometry residuals now converge with refinement instead of sitting on a floor, and every knot-to-knot correlation turns positive. Co-Authored-By: Claude Opus 5 (1M context) --- src/Equilibrium/DirectEquilibrium.jl | 78 +++++++++++++++---- src/Equilibrium/DirectEquilibriumArcLength.jl | 4 +- 2 files changed, 67 insertions(+), 15 deletions(-) diff --git a/src/Equilibrium/DirectEquilibrium.jl b/src/Equilibrium/DirectEquilibrium.jl index 961a0f7f5..586d81001 100644 --- a/src/Equilibrium/DirectEquilibrium.jl +++ b/src/Equilibrium/DirectEquilibrium.jl @@ -227,6 +227,32 @@ function direct_position!(raw_profile::DirectRunInput) return ro, zo, rs1, rs2 end +""" + eta_at_sfl_angle(sol, y_out, x, total_x) -> Float64 + +Integration angle η at which the normalised straight-fieldline angle ∫jac·dl/Bp reaches `x`. + +`y_out[:, 5]` is monotone in η, so it brackets the root to one solver step and Brent converges in +a handful of dense-output evaluations. Used to sample every flux surface at the *same* SFL angles +instead of resampling each surface's own solver steps (issue #376). +""" +function eta_at_sfl_angle(sol, y_out::Matrix{Float64}, x::Float64, total_x::Float64) + x <= 0 && return y_out[1, 1] + x >= 1 && return y_out[end, 1] + target = x * total_x + hi = searchsortedfirst(view(y_out, :, 5), target) + hi = clamp(hi, 2, size(y_out, 1)) + lo = hi - 1 + eta_lo, eta_hi = y_out[lo, 1], y_out[hi, 1] + f(eta) = sol(eta)[4] - target + flo, fhi = f(eta_lo), f(eta_hi) + # Degenerate bracket (repeated η, or the root sitting exactly on a step) needs no solve. + flo == 0 && return eta_lo + fhi == 0 && return eta_hi + (flo * fhi > 0 || eta_hi <= eta_lo) && return eta_lo + (eta_hi - eta_lo) * (target - y_out[lo, 5]) / max(y_out[hi, 5] - y_out[lo, 5], eps()) + return find_zero(f, (eta_lo, eta_hi), Roots.Brent()) +end + """ direct_fieldline_int(psifac, raw_profile, ro, zo, rs2) @@ -250,9 +276,12 @@ from 1:5 rather than 0:4 as in Fortran. - `y_out[:, 4]`: ∫(dl/(R²Bp)) - `y_out[:, 5]`: ∫(jac*dl/Bp) + - `sol`: the dense ODE solution, so callers can evaluate the trace at prescribed SFL angles + rather than resampling this surface's own solver steps (`nothing` for tracers without it). + - `bfield`: A `DirectBField` object with values at the integration start point. """ -function direct_fieldline_int(psifac::Float64, raw_profile::DirectRunInput, ro::Float64, zo::Float64, rs2::Float64)::Tuple{Matrix{Float64},DirectBField} +function direct_fieldline_int(psifac::Float64, raw_profile::DirectRunInput, ro::Float64, zo::Float64, rs2::Float64) # Find the starting point on the flux surface (outboard midplane) psi0_guess = raw_profile.psio * (1.0 - psifac) @@ -291,10 +320,12 @@ function direct_fieldline_int(psifac::Float64, raw_profile::DirectRunInput, ro:: callback = DiscreteCallback((u, t, i) -> true, refine_affect!; save_positions=(true, false)) prob = ODEProblem{true}(direct_fieldline_der!, u0, (0.0, 2π), params) - sol = solve(prob, Vern9(); callback=callback, reltol=equil_config.etol, abstol=1e-8, dt=2π / 200, adaptive=true, dense=false) + # Dense output lets the caller evaluate the trace at the SFL angles it actually wants, instead + # of splining this surface's solver-chosen steps and resampling (issue #376). + sol = solve(prob, Vern9(); callback=callback, reltol=equil_config.etol, abstol=1e-8, dt=2π / 200, adaptive=true, dense=true) sol_matrix = reduce(hcat, sol.u::Vector{Vector{Float64}})' - return hcat(sol.t::Vector{Float64}, sol_matrix), bfield + return hcat(sol.t::Vector{Float64}, sol_matrix), bfield, sol end """ @@ -495,18 +526,39 @@ robustness. ff_deriv_val = zeros!(pool, Float64, 4) for ipsi in (mpsi+1):-1:1 # outermost to innermost - y_out, bfield = fieldline_int(psi_nodes[ipsi], raw_profile, ro, zo, rs2) + y_out, bfield, sol = fieldline_int(psi_nodes[ipsi], raw_profile, ro, zo, rs2) checkpoint!(pool, Float64) - # Fit data into temporary straight fieldline poloidal angle splines - ff_x_nodes = acquire!(pool, Float64, size(y_out, 1)) - @. ff_x_nodes = @view(y_out[:, 5]) / y_out[end, 5] - - ff_fs_nodes = acquire!(pool, Float64, size(y_out, 1), 4) - @. ff_fs_nodes[:, 1] = @view(y_out[:, 3])^2 - @. ff_fs_nodes[:, 2] = @view(y_out[:, 1]) / (2π) - ff_x_nodes - @. ff_fs_nodes[:, 3] = bfield.f * (@view(y_out[:, 4]) - ff_x_nodes * y_out[end, 4]) - @. ff_fs_nodes[:, 4] = @view(y_out[:, 2]) / y_out[end, 2] - ff_x_nodes + # Straight-fieldline angle x = normalised ∫jac·dl/Bp, monotone in the integration angle η. + # + # Sampling x at this surface's own solver steps and resampling onto theta_nodes leaves a + # resample error that is uncorrelated between neighbouring surfaces, i.e. white noise in ψ + # that grid refinement then amplifies (issue #376). With dense output we instead solve for + # the η where x hits each target node and evaluate there, so every surface is sampled at + # the same abscissae and the resample error at the output nodes is zero. + nff = sol === nothing ? size(y_out, 1) : mtheta + 1 + ff_x_nodes = acquire!(pool, Float64, nff) + ff_fs_nodes = acquire!(pool, Float64, nff, 4) + + if sol === nothing + @. ff_x_nodes = @view(y_out[:, 5]) / y_out[end, 5] + @. ff_fs_nodes[:, 1] = @view(y_out[:, 3])^2 + @. ff_fs_nodes[:, 2] = @view(y_out[:, 1]) / (2π) - ff_x_nodes + @. ff_fs_nodes[:, 3] = bfield.f * (@view(y_out[:, 4]) - ff_x_nodes * y_out[end, 4]) + @. ff_fs_nodes[:, 4] = @view(y_out[:, 2]) / y_out[end, 2] - ff_x_nodes + else + total_x = y_out[end, 5] + for itheta in 1:(mtheta+1) + x = theta_nodes[itheta] + eta = eta_at_sfl_angle(sol, y_out, x, total_x) + u = sol(eta) + ff_x_nodes[itheta] = x + ff_fs_nodes[itheta, 1] = u[2]^2 + ff_fs_nodes[itheta, 2] = eta / (2π) - x + ff_fs_nodes[itheta, 3] = bfield.f * (u[3] - x * y_out[end, 4]) + ff_fs_nodes[itheta, 4] = u[1] / y_out[end, 2] - x + end + end ff_fs_nodes[end, :] .= ff_fs_nodes[1, :] # enforce periodic endpoint diff --git a/src/Equilibrium/DirectEquilibriumArcLength.jl b/src/Equilibrium/DirectEquilibriumArcLength.jl index 32b93a0ac..ec4c661d3 100644 --- a/src/Equilibrium/DirectEquilibriumArcLength.jl +++ b/src/Equilibrium/DirectEquilibriumArcLength.jl @@ -75,7 +75,7 @@ outboard midplane (Z = zo, R > ro) after a minimum arc-length guard. """ @with_pool pool function arclength_fieldline_int( psifac::Float64, raw_profile::DirectRunInput, ro::Float64, zo::Float64, rs2::Float64 -)::Tuple{Matrix{Float64},DirectBField} +)::Tuple{Matrix{Float64},DirectBField,Nothing} psi0_guess = raw_profile.psio * (1.0 - psifac) r = ro + sqrt(psifac) * (rs2 - ro) @@ -142,6 +142,6 @@ outboard midplane (Z = zo, R > ro) after a minimum arc-length guard. end # bfield at the starting point carries F and P for the surface-averaged quantities - return y_out, bfield + return y_out, bfield, nothing end From aefa3cf1ff153c4c5ecc69ae2b6e81a9f25901ca Mon Sep 17 00:00:00 2001 From: logan-nc Date: Tue, 18 Aug 2026 08:44:07 -0400 Subject: [PATCH 2/4] 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 3/4] 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 4/4] 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