From de0db4f551b4774df959109e146dfffd6a5532b3 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Sun, 16 Aug 2026 18:28:59 -0400 Subject: [PATCH 01/23] 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 02/23] 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 03/23] 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 04/23] 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 From 09cb856d3eb07464e01f24b505ef4297bc85738b Mon Sep 17 00:00:00 2001 From: logan-nc Date: Thu, 20 Aug 2026 11:47:01 -0400 Subject: [PATCH 05/23] KF - NEW FEATURE - Suppress kinetic terms where the zero-orbit-width ordering fails near the axis Physics ruling (issue #376 DIII-D kinetic pathology): the drift-kinetic model loses validity where thermal ion orbit widths reach the local minor radius. psi_c = outermost crossing of by max(potato width (q^2 rho^2 R0)^(1/3), banana width q rho/sqrt(eps), poloidal gyroradius q rho/eps), computed from the equilibrium and kinetic profiles at runtime -- no user tuning parameters (the Fortran ktanh_flag precedent needed four). A C2 quintic envelope zeroes the calculated kinetic increments below psi_c (kernel evaluation skipped) and rises to 1 at 2 psi_c; the same boundary and envelope apply to the NTV torque psi quadrature (one source of truth). One Bool (axis_validity_suppression, default true) to disable for debugging. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01LzbLFQKyuRE5DYZmLokKmk --- examples/DIIID-like_ideal_example/gpec.toml | 1 + .../Solovev_kinetic_NTV_example/gpec.toml | 1 + examples/a10_kinetic_example/gpec.toml | 1 + .../CalculatedKineticMatrices.jl | 68 +++++++++----- src/KineticForces/Compute.jl | 61 ++++++++----- src/KineticForces/KineticForcesStructs.jl | 88 ++++++++++--------- src/KineticForces/Utils.jl | 56 +++++++++++- 7 files changed, 188 insertions(+), 88 deletions(-) diff --git a/examples/DIIID-like_ideal_example/gpec.toml b/examples/DIIID-like_ideal_example/gpec.toml index b84e12872..1edc25231 100644 --- a/examples/DIIID-like_ideal_example/gpec.toml +++ b/examples/DIIID-like_ideal_example/gpec.toml @@ -87,5 +87,6 @@ f0type = "maxwellian" # Equilibrium distribution moment = "pressure" # Pressure-moment NTV torque atol_xlmda = 1e-9 # Absolute tolerance for inner pitch + energy integrations rtol_xlmda = 1e-5 # Relative tolerance for inner pitch + energy integrations +axis_validity_suppression = true # Suppress kinetic terms where the zero-orbit-width ordering fails near the axis (profile-derived boundary, no tuning parameters) write_outputs_to_HDF5 = true # Write outputs to the HDF5 file verbose = true # Enable verbose logging diff --git a/examples/Solovev_kinetic_NTV_example/gpec.toml b/examples/Solovev_kinetic_NTV_example/gpec.toml index e52b0d537..7a54a0397 100644 --- a/examples/Solovev_kinetic_NTV_example/gpec.toml +++ b/examples/Solovev_kinetic_NTV_example/gpec.toml @@ -87,3 +87,4 @@ f0fac = 1 # Scale toroidal field at constant pressure (β, q change; Φ, p, [KineticForces] kinetic_file = "kinetic.dat" # Kinetic profile file: psi_n, n_i, n_e, T_i, T_e, omega_E columns +axis_validity_suppression = true # Suppress kinetic terms where the zero-orbit-width ordering fails near the axis (profile-derived boundary, no tuning parameters) diff --git a/examples/a10_kinetic_example/gpec.toml b/examples/a10_kinetic_example/gpec.toml index fd8881122..ae602799f 100644 --- a/examples/a10_kinetic_example/gpec.toml +++ b/examples/a10_kinetic_example/gpec.toml @@ -59,3 +59,4 @@ nutype = "harmonic" # Collision operator (zero, small, krook, harmoni f0type = "maxwellian" # Distribution function (maxwellian, jkp, cgl) atol_xlmda = 1e-9 # Absolute tolerance for inner pitch + energy integrations rtol_xlmda = 1e-5 # Relative tolerance for inner pitch + energy integrations +axis_validity_suppression = true # Suppress kinetic terms where the zero-orbit-width ordering fails near the axis (profile-derived boundary, no tuning parameters) diff --git a/src/KineticForces/CalculatedKineticMatrices.jl b/src/KineticForces/CalculatedKineticMatrices.jl index d36fead99..9d8a19def 100644 --- a/src/KineticForces/CalculatedKineticMatrices.jl +++ b/src/KineticForces/CalculatedKineticMatrices.jl @@ -31,22 +31,25 @@ tracked as follow-up work blocked on PR #196 — see the plan's "Out of scope" section. # Arguments -- `ffs_ctrl`: ForceFreeStatesControl (carries `kinetic_factor`, `kinetic_source`) -- `equil`: PlasmaEquilibrium with 2D interpolants and named profile/geometry splines -- `ffs_intr`: ForceFreeStatesInternal (mode indexing) -- `metric`: MetricData (provides ψ grid via `metric.xs`) -- `ffit`: FourFitVars (used only for `numpert_total` cross-check) + + - `ffs_ctrl`: ForceFreeStatesControl (carries `kinetic_factor`, `kinetic_source`) + - `equil`: PlasmaEquilibrium with 2D interpolants and named profile/geometry splines + - `ffs_intr`: ForceFreeStatesInternal (mode indexing) + - `metric`: MetricData (provides ψ grid via `metric.xs`) + - `ffit`: FourFitVars (used only for `numpert_total` cross-check) # Keyword arguments -- `kf_ctrl`: KineticForcesControl, defaults to `KineticForcesControl()`. Used to - carry NTV-specific knobs (nl, zi, mi, wdfac, divxfac, electron) that the - KineticForces kernel needs but ForceFreeStatesControl does not expose. -- `kinetic_profiles::Equilibrium.KineticProfileSplines`: Required. Named kinetic- - profile splines loaded via `Equilibrium.load_kinetic_profiles`. + + - `kf_ctrl`: KineticForcesControl, defaults to `KineticForcesControl()`. Used to + carry NTV-specific knobs (nl, zi, mi, wdfac, divxfac, electron) that the + KineticForces kernel needs but ForceFreeStatesControl does not expose. + - `kinetic_profiles::Equilibrium.KineticProfileSplines`: Required. Named kinetic- + profile splines loaded via `Equilibrium.load_kinetic_profiles`. # Returns -- `kw_flat::Array{ComplexF64,3}`: Energy matrices, shape `(mpsi, np^2, 6)` -- `kt_flat::Array{ComplexF64,3}`: Torque matrices, shape `(mpsi, np^2, 6)` + + - `kw_flat::Array{ComplexF64,3}`: Energy matrices, shape `(mpsi, np^2, 6)` + - `kt_flat::Array{ComplexF64,3}`: Torque matrices, shape `(mpsi, np^2, 6)` """ function compute_calculated_kinetic_matrices( _ffs_ctrl, @@ -54,10 +57,13 @@ function compute_calculated_kinetic_matrices( ffs_intr, metric, ffit; - kf_ctrl::KineticForcesControl = KineticForcesControl(), + kf_ctrl::KineticForcesControl=KineticForcesControl(), kinetic_profiles::Equilibrium.KineticProfileSplines, + psis::Union{Nothing,Vector{Float64}}=nothing ) - xs = metric.xs + # The kernel is a pure function of psi (it evaluates equilibrium splines), so it can be + # driven over any knot list; default is the full equilibrium grid. + xs = psis === nothing ? metric.xs : psis mpsi = length(xs) mpert = ffs_intr.mpert npert = ffs_intr.npert @@ -93,22 +99,36 @@ function compute_calculated_kinetic_matrices( # ipsi row of kw_flat/kt_flat. Per-thread copies of kf_intr provide isolated # tpsi_* θ-grid buffers and interpolant hint refs; geometric/profile splines # are read-only and safely shared through deepcopy semantics. + # Near-axis validity envelope: suppress the drift-kinetic increments where the + # zero-orbit-width ordering fails; kernel evaluation is skipped where it is identically 0. + env = ones(Float64, mpsi) + if kf_ctrl.axis_validity_suppression + psi_c = kinetic_axis_validity_psi(kinetic_profiles, equil; + zi=kf_ctrl.zi, mi=kf_ctrl.mi, electron=kf_ctrl.electron) + if psi_c > 0 + env .= kinetic_axis_validity_envelope.(xs, psi_c) + @info "Kinetic axis-validity suppression: psi_c=$(round(psi_c; sigdigits=3)), envelope reaches 1 at " * + "psi=$(round(2 * psi_c; sigdigits=3)); $(count(iszero, env)) of $mpsi surfaces skipped" + end + end + nl = kf_ctrl.nl nthreads = Threads.maxthreadid() - thread_intrs = [deepcopy(kf_intr) for _ in 1:nthreads] - thread_full_w = [zeros(ComplexF64, mpert, mpert, 6) for _ in 1:nthreads] - thread_full_t = [zeros(ComplexF64, mpert, mpert, 6) for _ in 1:nthreads] + thread_intrs = [deepcopy(kf_intr) for _ in 1:nthreads] + thread_full_w = [zeros(ComplexF64, mpert, mpert, 6) for _ in 1:nthreads] + thread_full_t = [zeros(ComplexF64, mpert, mpert, 6) for _ in 1:nthreads] thread_block_w = [zeros(ComplexF64, mpert, mpert, 6) for _ in 1:nthreads] thread_block_t = [zeros(ComplexF64, mpert, mpert, 6) for _ in 1:nthreads] Threads.@threads for ipsi in 1:mpsi - tid = Threads.threadid() - intr_t = thread_intrs[tid] - full_w = thread_full_w[tid] - full_t = thread_full_t[tid] + tid = Threads.threadid() + intr_t = thread_intrs[tid] + full_w = thread_full_w[tid] + full_t = thread_full_t[tid] block_w = thread_block_w[tid] block_t = thread_block_t[tid] - psi = xs[ipsi] + psi = xs[ipsi] + env[ipsi] == 0.0 && continue for in_idx in 1:npert n = ffs_intr.nlow + in_idx - 1 fill!(full_w, 0) @@ -131,8 +151,8 @@ function compute_calculated_kinetic_matrices( row_offset = (in_idx - 1) * mpert for k in 1:6, j in 1:mpert, i in 1:mpert idx = (row_offset + j - 1) * np + (row_offset + i) - kw_flat[ipsi, idx, k] = full_w[i, j, k] - kt_flat[ipsi, idx, k] = full_t[i, j, k] + kw_flat[ipsi, idx, k] = env[ipsi] * full_w[i, j, k] + kt_flat[ipsi, idx, k] = env[ipsi] * full_t[i, j, k] end end end diff --git a/src/KineticForces/Compute.jl b/src/KineticForces/Compute.jl index 4175a62b4..3d9c3d8cf 100644 --- a/src/KineticForces/Compute.jl +++ b/src/KineticForces/Compute.jl @@ -41,7 +41,7 @@ Warn when the ψ torque quadrature terminated without satisfying the requested t silent-garbage scenario for weak applied fields, since NTV scales as δB². """ function check_psi_quadrature_convergence(total::ComplexF64, quad_err::Float64, - ctrl::KineticForcesControl, method::String) + ctrl::KineticForcesControl, method::String) if quad_err > max(ctrl.atol_psi, ctrl.rtol_psi * abs(total)) @warn "ψ torque quadrature ($method) did not converge within maxevals_psi=$(ctrl.maxevals_psi): " * "error estimate $quad_err vs |T|=$(abs(total)) N·m. Raise maxevals_psi or loosen rtol_psi." @@ -63,14 +63,16 @@ diagnostic T(ψ) profile at no extra cost (the values are computed anyway — we just keep them). # Returns + NamedTuple with: -- `total::ComplexF64`: Total integrated torque -- `torque_profile`: NamedTuple of (psi, dtdpsi, t_cumulative) from evaluation points -- `matrix_integrated`: Trapezoidal-integrated mpert×mpert×6 matrix (if matrix method) -- `psi_nsteps::Int`: Number of integrand evaluations -- `psi_quad_error::Float64`: Quadrature error estimate for the total torque -- `panel_psis::Vector{Float64}`: Quadrature panel boundaries actually used (bounds + interior resonant surfaces) -- `resonance_psis::Vector{Float64}`: Located kinetic-resonance ψ surfaces (Ω_ℓ(x=1)=0), for diagnostics/plotting + + - `total::ComplexF64`: Total integrated torque + - `torque_profile`: NamedTuple of (psi, dtdpsi, t_cumulative) from evaluation points + - `matrix_integrated`: Trapezoidal-integrated mpert×mpert×6 matrix (if matrix method) + - `psi_nsteps::Int`: Number of integrand evaluations + - `psi_quad_error::Float64`: Quadrature error estimate for the total torque + - `panel_psis::Vector{Float64}`: Quadrature panel boundaries actually used (bounds + interior resonant surfaces) + - `resonance_psis::Vector{Float64}`: Located kinetic-resonance ψ surfaces (Ω_ℓ(x=1)=0), for diagnostics/plotting The integral is paneled at the rational-surface ψ locations (`intr.sing_psis`) and capped at `ctrl.maxevals_psi` evaluations; a warning is emitted if the quadrature fails to reach @@ -95,9 +97,21 @@ function integrate_psi_quadgk( if x0 >= xout return (total=ComplexF64(0.0), torque_profile=nothing, matrix_integrated=nothing, psi_nsteps=0, psi_quad_error=0.0, - panel_psis=Float64[], resonance_psis=Float64[]) + panel_psis=Float64[], resonance_psis=Float64[]) end + # Near-axis validity suppression: same boundary and envelope as the EL kinetic + # matrices (one source of truth), applied to the torque density; the quadrature + # domain starts at the boundary since the integrand is identically zero below it. + psi_c = ctrl.axis_validity_suppression ? + kinetic_axis_validity_psi(kinetic_profiles, equil; zi=zi, mi=mi, electron=electron) : 0.0 + x0 = max(x0, psi_c) + if x0 >= xout + return (total=ComplexF64(0.0), torque_profile=nothing, matrix_integrated=nothing, psi_nsteps=0, psi_quad_error=0.0, + panel_psis=Float64[], resonance_psis=Float64[]) + end + psi_c > 0 && @info "Kinetic axis-validity suppression in ψ torque quadrature: domain starts at psi_c=$(round(psi_c; sigdigits=3))" + # The outer ψ-integral (the QuadGK batch / ψ-node loop) stays serial: QuadGK's refine # loop invokes the callback with small batches (~15 nodes), so threading it is # fork-join-bound. Instead thread the inner bounce-harmonic loop (2·nl+1 harmonics), @@ -133,9 +147,9 @@ function integrate_psi_quadgk( w = is_matrix_method ? thread_wtw[tid] : nothing is_matrix_method && fill!(w, 0) tpsi!(thread_tpsi[tid], psi, n, l, zi, mi, wdfac, divxfac, - electron, method, equil, thread_intrs[tid], kinetic_profiles; - op_wmats=w, - atol_xlmda=ctrl.atol_xlmda, rtol_xlmda=ctrl.rtol_xlmda) + electron, method, equil, thread_intrs[tid], kinetic_profiles; + op_wmats=w, + atol_xlmda=ctrl.atol_xlmda, rtol_xlmda=ctrl.rtol_xlmda) harm_vals[ell_idx] = thread_tpsi[tid][] is_matrix_method && (harm_elems[ell_idx] .= w) end @@ -145,6 +159,7 @@ function integrate_psi_quadgk( for ell_idx in 1:nharm total += harm_vals[ell_idx] end + total *= kinetic_axis_validity_envelope(psi, psi_c) y[k] = total push!(logged_psi, psi) @@ -154,6 +169,7 @@ function integrate_psi_quadgk( for ell_idx in 1:nharm elems_accum .+= harm_elems[ell_idx] end + elems_accum .*= kinetic_axis_validity_envelope(psi, psi_c) push!(logged_elems, elems_accum) end end @@ -162,7 +178,7 @@ function integrate_psi_quadgk( # Panels at the rational surfaces the run resolved plus the kinetic-resonance # surfaces (thermal-energy Ω_ℓ = 0 for ℓ ∈ -nl:nl) — both are torque-density peaks. resonance_psis = kinetic_resonance_psi_nodes(kinetic_profiles, equil; n, nl, zi, mi, electron, wdfac) - pts = psi_panel_points(vcat(intr.sing_psis, resonance_psis), x0, xout) + pts = psi_panel_points(vcat(intr.sing_psis, resonance_psis, [2 * psi_c]), x0, xout) bi = QuadGK.BatchIntegrand(psi_batch!, ComplexF64[], Float64[]) total, quad_err = quadgk(bi, pts...; atol=ctrl.atol_psi, rtol=ctrl.rtol_psi, maxevals=ctrl.maxevals_psi) @@ -205,7 +221,7 @@ function integrate_psi_quadgk( "$(length(pts) - 1) panels ($n_rational rational + $n_resonance kinetic resonance surfaces)" return (total=total, torque_profile=torque_profile, matrix_integrated=matrix_integrated, psi_nsteps=npts, psi_quad_error=quad_err, - panel_psis=pts, resonance_psis=sort(resonance_psis)) + panel_psis=pts, resonance_psis=sort(resonance_psis)) end @@ -224,15 +240,16 @@ For multi-n calculations, loops over toroidal mode numbers and assembles block-diagonal kinetic matrices. # Arguments -- `state::KineticForcesState`: Accumulates results for all methods -- `intr::KineticForcesInternal`: Internal state with equilibrium data -- `ctrl::KineticForcesControl`: Control parameters specifying which methods to run -- `equil`: PlasmaEquilibrium with 2D interpolants -- `kinetic_profiles::Equilibrium.KineticProfileSplines`: Named kinetic-profile splines + + - `state::KineticForcesState`: Accumulates results for all methods + - `intr::KineticForcesInternal`: Internal state with equilibrium data + - `ctrl::KineticForcesControl`: Control parameters specifying which methods to run + - `equil`: PlasmaEquilibrium with 2D interpolants + - `kinetic_profiles::Equilibrium.KineticProfileSplines`: Named kinetic-profile splines """ function compute_torque_all_methods!(state::KineticForcesState, intr::KineticForcesInternal, - ctrl::KineticForcesControl, equil, - kinetic_profiles::Equilibrium.KineticProfileSplines) + ctrl::KineticForcesControl, equil, + kinetic_profiles::Equilibrium.KineticProfileSplines) for entry in METHOD_REGISTRY getfield(ctrl, entry.flag) || continue @@ -311,7 +328,7 @@ function compute_torque_all_methods!(state::KineticForcesState, intr::KineticFor t_cumulative=t_cum_out, psi_nsteps=psi_nsteps_total, panel_psis=panel_psis_out, - resonance_psis=resonance_psis_out, + resonance_psis=resonance_psis_out ) state.method_results[method] = result_entry diff --git a/src/KineticForces/KineticForcesStructs.jl b/src/KineticForces/KineticForcesStructs.jl index 0d1b46464..f3e1fbf9b 100644 --- a/src/KineticForces/KineticForcesStructs.jl +++ b/src/KineticForces/KineticForcesStructs.jl @@ -3,15 +3,16 @@ Single source of truth for the NTV calculation methods. Each entry is a NamedTuple `(name, flag, kind, doc)`: -- `name` — short method identifier used as the HDF5 group key and in `intr.method` -- `flag` — the `KineticForcesControl` field symbol that enables the method -- `kind` — dispatch routing tag consumed by `method_kind` / `Torque.jl` - (`:gar` for the GAR/matrix family, `:fcgl`/`:rlar`/`:clar` for the - three special-cased methods) -- `doc` — one-line description printed in verbose output + + - `name` — short method identifier used as the HDF5 group key and in `intr.method` + - `flag` — the `KineticForcesControl` field symbol that enables the method + - `kind` — dispatch routing tag consumed by `method_kind` / `Torque.jl` + (`:gar` for the GAR/matrix family, `:fcgl`/`:rlar`/`:clar` for the + three special-cased methods) + - `doc` — one-line description printed in verbose output The method names/docs and the `Compute.jl` enable list are all derived from this -tuple, and `Torque.jl` routes on `kind`, so the methods are enumerated in one place. +tuple, and `Torque.jl` routes on `kind`, so the methods are enumerated in one place. To add a method: append an entry here and add the matching `*_flag` field to `KineticForcesControl`. """ @@ -56,6 +57,7 @@ User-facing control parameters from the TOML `[KineticForces]` section. Configures which NTV methods to run, species parameters, tolerances, and output options. Constructed via keyword arguments or from a TOML dict: + ```julia ctrl = KineticForcesControl(; (Symbol(k) => v for (k, v) in inputs["KineticForces"])...) ``` @@ -123,6 +125,7 @@ ctrl = KineticForcesControl(; (Symbol(k) => v for (k, v) in inputs["KineticForce nufac::Float64 = 1.0 # Collisionality scaling divxfac::Float64 = 1.0 # div(xi_perp) scaling + axis_validity_suppression::Bool = true # Suppress kinetic terms where the zero-orbit-width ordering fails near the axis (see kinetic_axis_validity_psi); envelope and boundary are profile-derived, no tuning parameters # Energy integration parameters nutype::String = "harmonic" # Collision operator: "zero", "small", "krook", "harmonic" @@ -176,12 +179,13 @@ Internal working state for KineticForces calculations. Holds equilibrium-derived quantities, profile interpolants, and integration results. Fields replacing former module-level globals: -- `ro`, `bo`, `chi1`: Equilibrium geometry parameters -- `mthsurf`, `mfac`: Poloidal grid info -- `dbob_m`, `divx_m`: Perturbation mode interpolants -- `sing_psis`: Rational-surface ψ locations (sorted, from the stability analysis), used as - panel boundaries for the outer ψ torque quadrature so the resonant peaks fall on - Gauss-Kronrod interval endpoints instead of driving deep adaptive bisection + + - `ro`, `bo`, `chi1`: Equilibrium geometry parameters + - `mthsurf`, `mfac`: Poloidal grid info + - `dbob_m`, `divx_m`: Perturbation mode interpolants + - `sing_psis`: Rational-surface ψ locations (sorted, from the stability analysis), used as + panel boundaries for the outer ψ torque quadrature so the resonant peaks fall on + Gauss-Kronrod interval endpoints instead of driving deep adaptive bisection Equilibrium and kinetic profile data are read directly from the `PlasmaEquilibrium` (`equil.profiles`, `equil.geometry`) and the @@ -271,17 +275,17 @@ function KineticForcesInternal(equil; verbose::Bool=false) # Axis toroidal field F(0)/ro that normalizes λ = μ·bo/E; F_spline stores 2πF. bo_axis = abs(equil.profiles.F_spline(0.0)) / (2π * equil.ro) KineticForcesInternal(; - ro = equil.ro, - bo = bo_axis, - chi1 = 2π * equil.psio, + ro=equil.ro, + bo=bo_axis, + chi1=2π * equil.psio, mthsurf, - tpsi_xs = collect(range(0.0, 1.0, length=nth)), - tpsi_B = Vector{Float64}(undef, nth), - tpsi_dBdpsi = Vector{Float64}(undef, nth), - tpsi_dBdtheta = Vector{Float64}(undef, nth), - tpsi_jac = Vector{Float64}(undef, nth), - tpsi_djdpsi = Vector{Float64}(undef, nth), - verbose, + tpsi_xs=collect(range(0.0, 1.0; length=nth)), + tpsi_B=Vector{Float64}(undef, nth), + tpsi_dBdpsi=Vector{Float64}(undef, nth), + tpsi_dBdtheta=Vector{Float64}(undef, nth), + tpsi_jac=Vector{Float64}(undef, nth), + tpsi_djdpsi=Vector{Float64}(undef, nth), + verbose ) end @@ -291,19 +295,21 @@ end Populate perturbation data from PerturbedEquilibriumState into KineticForcesInternal. Builds three interpolant sets from PE Clebsch displacements: -1. `xs_m` — [ξ^ψ, ∂ξ^ψ/∂ψ, ξ^α] CubicSeriesInterpolants over ψ -2. `dbob_m` — δB/B Fourier modes via JBB deweighting (Fortran set_peq) -3. `divx_m` — ∇·ξ⊥ Fourier modes via JBB deweighting + + 1. `xs_m` — [ξ^ψ, ∂ξ^ψ/∂ψ, ξ^α] CubicSeriesInterpolants over ψ + 2. `dbob_m` — δB/B Fourier modes via JBB deweighting (Fortran set_peq) + 3. `divx_m` — ∇·ξ⊥ Fourier modes via JBB deweighting The JBB deweighting algorithm (Fortran pentrc/inputs.f90:828-868): -1. Apply geometric matrices S,T,X,Y,Z in m-space -2. Inverse DFT to θ-space -3. Divide by J·B² at each θ -4. Forward DFT back to m-space + + 1. Apply geometric matrices S,T,X,Y,Z in m-space + 2. Inverse DFT to θ-space + 3. Divide by J·B² at each θ + 4. Forward DFT back to m-space """ function set_perturbation_data!(kf_intr::KineticForcesInternal, pe_state, ffs_intr, - equil::Equilibrium.PlasmaEquilibrium, - metric::ForceFreeStates.MetricData) + equil::Equilibrium.PlasmaEquilibrium, + metric::ForceFreeStates.MetricData) # Copy mode numbers from FFS kf_intr.mlow = ffs_intr.mlow kf_intr.mhigh = ffs_intr.mhigh @@ -372,9 +378,9 @@ function set_perturbation_data!(kf_intr::KineticForcesInternal, pe_state, ffs_in psi = psi_grid[ipsi] # Get Clebsch displacement vectors at this ψ - xsp = view(xi_modes.clebsch_psi, ipsi, :) # ξ^ψ [mpert] + xsp = view(xi_modes.clebsch_psi, ipsi, :) # ξ^ψ [mpert] xmp1 = view(xi_modes.clebsch_psi1, ipsi, :) # ∂ξ^ψ/∂ψ [mpert] - xms = view(clebsch_alpha_mat, ipsi, :) # ξ^α [mpert] + xms = view(clebsch_alpha_mat, ipsi, :) # ξ^α [mpert] # Evaluate geometric matrices at ψ → mpert² flat vectors, reshape to mpert×mpert geom_mats.smats(smat_flat, psi; hint=hint_s) @@ -394,8 +400,8 @@ function set_perturbation_data!(kf_intr::KineticForcesInternal, pe_state, ffs_in mul!(jbb_kapx, smat, xsp) mul!(jbb_kapx, tmat, xms, 1.0 + 0.0im, 1.0 + 0.0im) # += tmat * xms mul!(jbb_divx, xmat, xmp1) - mul!(jbb_divx, ymat, xsp, 1.0 + 0.0im, 1.0 + 0.0im) # += ymat * xsp - mul!(jbb_divx, zmat, xms, 1.0 + 0.0im, 1.0 + 0.0im) # += zmat * xms + mul!(jbb_divx, ymat, xsp, 1.0 + 0.0im, 1.0 + 0.0im) # += ymat * xsp + mul!(jbb_divx, zmat, xms, 1.0 + 0.0im, 1.0 + 0.0im) # += zmat * xms @. jbb_dbob = -(jbb_divx + jbb_kapx) # Inverse DFT to θ-space, divide by J·B², forward DFT back @@ -422,9 +428,9 @@ Matches Fortran set_peq lines 859-868: transforms JBB-weighted m-space data to θ-space, removes the J·B² weighting at each poloidal angle, and transforms back. """ function _jbb_deweight!(out::AbstractVector{ComplexF64}, jbb_modes::Vector{ComplexF64}, - ft::Utilities.FourierTransforms.FourierTransform, - psi::Float64, equil::Equilibrium.PlasmaEquilibrium, - mthsurf::Int, theta_buf::Vector{ComplexF64}) + ft::Utilities.FourierTransforms.FourierTransform, + psi::Float64, equil::Equilibrium.PlasmaEquilibrium, + mthsurf::Int, theta_buf::Vector{ComplexF64}) # Inverse DFT: m-space → θ-space theta_buf .= Utilities.FourierTransforms.inverse(ft, jbb_modes) @@ -494,8 +500,8 @@ Accumulated results from all KineticForces computations. Written to gpec.h5 under the "KineticForces" group. """ @kwdef mutable struct KineticForcesState - method_results::Dict{String, MethodResult} = Dict{String, MethodResult}() + method_results::Dict{String,MethodResult} = Dict{String,MethodResult}() # Block-diagonal kinetic matrices: key=method, value=(numpert_total, numpert_total, 6) - kinetic_matrices::Dict{String, Array{ComplexF64,3}} = Dict{String, Array{ComplexF64,3}}() + kinetic_matrices::Dict{String,Array{ComplexF64,3}} = Dict{String,Array{ComplexF64,3}}() completed::Bool = false end diff --git a/src/KineticForces/Utils.jl b/src/KineticForces/Utils.jl index 51c0cd3a8..97cfa4665 100644 --- a/src/KineticForces/Utils.jl +++ b/src/KineticForces/Utils.jl @@ -71,7 +71,7 @@ the estimate degenerates. Panel placement only needs ~peak-width accuracy, so th estimates (single spline evaluations) are sufficient and no bounce averaging is performed. """ function kinetic_resonance_psi_nodes(kinetic_profiles::Equilibrium.KineticProfileSplines, equil; - n::Int, nl::Int, zi::Int=1, mi::Int=2, electron::Bool=false, wdfac::Float64=1.0, xeval::Float64=2.5) + n::Int, nl::Int, zi::Int=1, mi::Int=2, electron::Bool=false, wdfac::Float64=1.0, xeval::Float64=2.5) chrg = electron ? -e : zi * e mass = electron ? me : mi * mp T_spline = electron ? kinetic_profiles.Te_spline : kinetic_profiles.Ti_spline @@ -89,3 +89,57 @@ function kinetic_resonance_psi_nodes(kinetic_profiles::Equilibrium.KineticProfil grid = filter(x -> x > 0, kinetic_profiles.xs) return _resonance_nodes_from_frequencies(wbhat_f, kinetic_profiles.omegaE_spline, wdhat_f, grid; n, nl, xeval) end + +""" + kinetic_axis_validity_psi(kinetic_profiles, equil; zi=1, mi=2, electron=false) → Float64 + +ψ_N below which the zero-orbit-width drift-kinetic ordering fails: the outermost ψ on the +kinetic-profile grid where any of the three thermal orbit-width scales reaches the local minor +radius ⟨r⟩ — potato width `(q²ρ²R₀)^(1/3)`, banana width `q·ρ/√ε`, and poloidal gyroradius +`q·ρ/ε` (ε = ⟨r⟩/⟨R⟩, clamped as in `kinetic_resonance_psi_nodes`; thermal gyroradius +ρ = √(2·m·T)·/(Z·e·B₀) of the computed species). Inside this boundary trapped bananas become +potato orbits with width comparable to r itself, so the bounce-averaged kinetic response is +evaluated outside its validity domain (measured consequence: diverging kinetic increments and a +pathological EL step count). Returns 0.0 when no criterion is met anywhere. The Fortran precedent +(`ktanh_flag`, dcon/fourfit.F) suppressed the same region with four hand-tuned knobs; here the +boundary is derived from the profiles with no user parameters. +""" +function kinetic_axis_validity_psi(kinetic_profiles::Equilibrium.KineticProfileSplines, equil; + zi::Int=1, mi::Int=2, electron::Bool=false) + chrg = electron ? e : zi * e + mass = electron ? me : mi * mp + T_spline = electron ? kinetic_profiles.Te_spline : kinetic_profiles.Ti_spline + q_spline = equil.profiles.q_spline + avg_r = equil.geometry.avg_r_spline + avg_R = equil.geometry.avg_R_spline + ro = abs(equil.ro) + bo = abs(equil.params.b0) + psi_c = 0.0 + for psi in kinetic_profiles.xs + psi <= 0 && continue + r = avg_r(psi) + r <= 0 && continue + eps = max(r / avg_R(psi), 1e-6) + rho = mass * sqrt(2 * T_spline(psi) / mass) / (abs(chrg) * bo) + q = abs(q_spline(psi)) + w_orbit = max(cbrt(q^2 * rho^2 * ro), q * rho / sqrt(eps), q * rho / eps) + w_orbit >= r && (psi_c = max(psi_c, psi)) + end + return psi_c +end + +""" + kinetic_axis_validity_envelope(psi, psi_c) → Float64 + +C² quintic smoothstep for the near-axis kinetic suppression: 0 for ψ ≤ ψ_c (drift-kinetic model +invalid; kernel evaluation may be skipped), rising over `[ψ_c, 2ψ_c]`, 1 above. The transition +width is tied to ψ_c itself, so there is no independent width parameter. `psi_c ≤ 0` returns 1 +(no suppression). +""" +function kinetic_axis_validity_envelope(psi::Float64, psi_c::Float64) + psi_c <= 0 && return 1.0 + t = (psi - psi_c) / psi_c + t <= 0 && return 0.0 + t >= 1 && return 1.0 + return t^3 * (10 + t * (6 * t - 15)) +end From f95230d38c1e5e0e7b8d61712f49fcbc3ae5b047 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Thu, 20 Aug 2026 11:56:29 -0400 Subject: [PATCH 06/23] KF - NEW FEATURE - Output kinetic-model validity profiles under KineticForces/Validity Whenever kinetic profiles are used (self-consistent matrices or NTV post-processing), write the thermal orbit-width scales (rho_i, rho_banana, rho_theta, w_potato), the local geometry (r_minor, d_separatrix), the profile gradient lengths (L_p, L_q), the near-axis boundary psi_c with its applied envelope, and an is_valid array (orbit width < r, rho_banana < L_p and L_q, orbit width < distance to separatrix). Validity outside the near-axis envelope is flagged, never suppressed -- the far edge can dominate the physical NTV. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01LzbLFQKyuRE5DYZmLokKmk --- src/GeneralizedPerturbedEquilibrium.jl | 39 +++++++++-- src/HDF5Schema.jl | 38 +++++++++-- .../CalculatedKineticMatrices.jl | 2 +- src/KineticForces/Utils.jl | 67 +++++++++++++++++++ 4 files changed, 132 insertions(+), 14 deletions(-) diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 7b17683cf..d1d85d22c 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -76,7 +76,7 @@ using .ForceFreeStates: eulerlagrange_integration, free_run, normalize_eigenfunc using .ForceFreeStates: galerkin_solve, write_galerkin!, GalerkinResult, gal_matched_odestate const _DEPRECATED_FFS_KEYS = ("mer_flag", "force_wv_symmetry", "ode_flag", "cyl_flag", "mat_flag", "reform_eq_with_psilim", - "use_riccati", "use_parallel", "parallel_threads", "populate_dense_xi") + "use_riccati", "use_parallel", "parallel_threads", "populate_dense_xi") const _DEPRECATED_EQUIL_KEYS = ("power_bp", "power_b", "power_r", "power_rc") # Drop deprecated keys from a parsed gpec.toml section so legacy files keep parsing @@ -419,10 +419,11 @@ function main_from_inputs( # Inject the KineticForces callback so the "calculated" source can # invoke compute_calculated_kinetic_matrices without ForceFreeStates # importing KineticForces (which would invert the load order). - calculated_cb = (c, e, i, m, f) -> - KineticForces.compute_calculated_kinetic_matrices( - c, e, i, m, f; - kf_ctrl=kf_ctrl, kinetic_profiles=kinetic_profiles) + calculated_cb = + (c, e, i, m, f; psis=nothing) -> + KineticForces.compute_calculated_kinetic_matrices( + c, e, i, m, f; + kf_ctrl=kf_ctrl, kinetic_profiles=kinetic_profiles, psis=psis) make_kinetic_matrix(ctrl, equil, ffit, intr, metric; calculated_source=calculated_cb) @@ -495,6 +496,8 @@ function main_from_inputs( inputs, forcing_modes_snapshot, gal_data; + kinetic_profiles=kinetic_profiles, + kf_ctrl=kf_ctrl, locstab=locstab, ballooning_boundary=ballooning_boundary ) @@ -705,7 +708,9 @@ function write_outputs_to_HDF5( forcing_modes::Union{Nothing,Vector{ForcingTerms.ForcingMode}}=nothing, gal_data::Union{GalerkinResult,Nothing}=nothing; locstab::Union{FastInterpolations.CubicSeriesInterpolant,Nothing}=nothing, - ballooning_boundary=(psi=Float64[], alpha=Float64[], alpha_critical=Float64[]) + ballooning_boundary=(psi=Float64[], alpha=Float64[], alpha_critical=Float64[]), + kinetic_profiles=nothing, + kf_ctrl=nothing ) # Idempotent: already done if a PerturbedEquilibrium stage ran. Leaves the stores empty @@ -821,6 +826,28 @@ function write_outputs_to_HDF5( out_h5["$fwd/xi_s"] = odet.xi_s_store out_h5["$fwd/crit"] = odet.crit_store + # Kinetic-model validity diagnostics: orbit-width scales vs local geometry and profile + # gradient lengths, whenever kinetic profiles were used (self-consistent matrices or NTV + # post-processing). Diagnostic only — nothing outside the near-axis envelope is suppressed. + if kinetic_profiles !== nothing + kfc = kf_ctrl === nothing ? KineticForces.KineticForcesControl() : kf_ctrl + vp = KineticForces.kinetic_validity_profiles(kinetic_profiles, equil; + zi=kfc.zi, mi=kfc.mi, electron=kfc.electron) + vg = "KineticForces/Validity" + out_h5["$vg/psi"] = vp.psi + out_h5["$vg/rho_i"] = vp.rho_i + out_h5["$vg/rho_banana"] = vp.rho_banana + out_h5["$vg/rho_theta"] = vp.rho_theta + out_h5["$vg/w_potato"] = vp.w_potato + out_h5["$vg/r_minor"] = vp.r_minor + out_h5["$vg/L_p"] = vp.L_p + out_h5["$vg/L_q"] = vp.L_q + out_h5["$vg/d_separatrix"] = vp.d_separatrix + out_h5["$vg/psi_c"] = vp.psi_c + out_h5["$vg/envelope"] = kfc.axis_validity_suppression ? vp.envelope : ones(length(vp.psi)) + out_h5["$vg/is_valid"] = Int8.(vp.is_valid) + end + # Write edge stability scan data (only present when psiedge < psilim). # Generalized (W, N) pencil energies — power-normalized, Jacobian-invariant; these are # the values findmax_dW_edge! uses to choose the truncation point. diff --git a/src/HDF5Schema.jl b/src/HDF5Schema.jl index a1a47ef44..06487b7ea 100644 --- a/src/HDF5Schema.jl +++ b/src/HDF5Schema.jl @@ -141,6 +141,25 @@ const MAIN_H5_ANNOTATIONS = [ "LocalStability/alpha_critical" => (; long_name="critical normalized pressure gradient α for first ballooning stability", dims=("psi_ballooning",), attach=(1 => "LocalStability/ballooning_psi",)), + # --- KineticForces/Validity/ --- + "KineticForces/Validity/psi" => (; long_name="normalized poloidal flux ψ_N of the kinetic validity profiles", scale="psi"), + "KineticForces/Validity/rho_i" => (; long_name="thermal ion gyroradius √(2mT)/(Z·e·B₀)", units="m", attach=(1 => "KineticForces/Validity/psi",)), + "KineticForces/Validity/rho_banana" => (; long_name="thermal banana orbit width q·ρ_i/√ε", units="m", attach=(1 => "KineticForces/Validity/psi",)), + "KineticForces/Validity/rho_theta" => (; long_name="thermal poloidal gyroradius q·ρ_i/ε", units="m", attach=(1 => "KineticForces/Validity/psi",)), + "KineticForces/Validity/w_potato" => (; long_name="potato orbit width (q²ρ_i²R₀)^(1/3)", units="m", attach=(1 => "KineticForces/Validity/psi",)), + "KineticForces/Validity/r_minor" => (; long_name="surface-average minor radius ⟨r⟩", units="m", attach=(1 => "KineticForces/Validity/psi",)), + "KineticForces/Validity/L_p" => (; long_name="pressure gradient scale length |p|/|dp/dr|", units="m", attach=(1 => "KineticForces/Validity/psi",)), + "KineticForces/Validity/L_q" => (; long_name="safety-factor gradient scale length |q|/|dq/dr|", units="m", attach=(1 => "KineticForces/Validity/psi",)), + "KineticForces/Validity/d_separatrix" => (; long_name="distance to the separatrix ⟨r⟩(1) − ⟨r⟩(ψ)", units="m", attach=(1 => "KineticForces/Validity/psi",)), + "KineticForces/Validity/psi_c" => (; long_name="near-axis kinetic validity boundary: outermost ψ_N where a thermal orbit width reaches ⟨r⟩"), + "KineticForces/Validity/envelope" => + (; long_name="near-axis suppression envelope applied to the calculated kinetic terms (1 = unsuppressed)", units="1", attach=(1 => "KineticForces/Validity/psi",)), + "KineticForces/Validity/is_valid" => (; + long_name="1 where every zero-orbit-width ordering holds: max orbit width < ⟨r⟩, ρ_banana < L_p and L_q, max orbit width < d_separatrix", + units="1", + attach=(1 => "KineticForces/Validity/psi",) + ), + # --- ForceFreeStates/Solutions/ForwardIntegration/ --- "ForceFreeStates/Solutions/ForwardIntegration/nstep" => (; long_name="number of saved solution snapshots"), "ForceFreeStates/Solutions/ForwardIntegration/nstep_total" => (; long_name="total ODE solver steps taken"), @@ -239,7 +258,7 @@ const MAIN_H5_ANNOTATIONS = [ "SurfaceGeometries/Plasma/z" => (; long_name="Cartesian z of plasma-surface point cloud", units="m"), "SurfaceGeometries/Wall/x" => (; long_name="Cartesian x of wall point cloud", units="m"), "SurfaceGeometries/Wall/y" => (; long_name="Cartesian y of wall point cloud", units="m"), - "SurfaceGeometries/Wall/z" => (; long_name="Cartesian z of wall point cloud", units="m"), + "SurfaceGeometries/Wall/z" => (; long_name="Cartesian z of wall point cloud", units="m") ] # Euler-Lagrange operator matrices: same wording per letter, Ideal/ and Kinetic/ variants. @@ -252,7 +271,7 @@ const _ELM_IDEAL_LETTERS = [ ("H", "Euler-Lagrange primitive coefficient matrix H"), ("F", "Euler-Lagrange derived coefficient matrix F"), ("K", "Euler-Lagrange derived coefficient matrix K"), - ("G", "Euler-Lagrange derived coefficient matrix G"), + ("G", "Euler-Lagrange derived coefficient matrix G") ] # The kinetic branch overwrites only A, B, C, K, G and adds f0; D, E, H, F are shared # unchanged from the ideal set and are not re-emitted. @@ -262,14 +281,19 @@ const _ELM_KINETIC_LETTERS = [ ("C", "Euler-Lagrange primitive coefficient matrix C"), ("K", "Euler-Lagrange derived coefficient matrix K"), ("G", "Euler-Lagrange derived coefficient matrix G"), - ("f0", "raw kinetic component matrix f0"), + ("f0", "raw kinetic component matrix f0") ] const ELM_H5_ANNOTATIONS = vcat( ["ForceFreeStates/EulerLagrangeMatrices/psi" => (; long_name="normalized poloidal flux ψ_N grid of the operator matrices", scale="psi")], - ["ForceFreeStates/EulerLagrangeMatrices/Ideal/$l" => - (; long_name="ideal " * d, dims=("psi", "mode_row", "mode_col"), attach=(1 => "ForceFreeStates/EulerLagrangeMatrices/psi",)) for (l, d) in _ELM_IDEAL_LETTERS], - ["ForceFreeStates/EulerLagrangeMatrices/Kinetic/$l" => - (; long_name="kinetic-modified " * d, dims=("psi", "mode_row", "mode_col"), attach=(1 => "ForceFreeStates/EulerLagrangeMatrices/psi",)) for (l, d) in _ELM_KINETIC_LETTERS] + [ + "ForceFreeStates/EulerLagrangeMatrices/Ideal/$l" => + (; long_name="ideal " * d, dims=("psi", "mode_row", "mode_col"), attach=(1 => "ForceFreeStates/EulerLagrangeMatrices/psi",)) for (l, d) in _ELM_IDEAL_LETTERS + ], + [ + "ForceFreeStates/EulerLagrangeMatrices/Kinetic/$l" => + (; long_name="kinetic-modified " * d, dims=("psi", "mode_row", "mode_col"), attach=(1 => "ForceFreeStates/EulerLagrangeMatrices/psi",)) for + (l, d) in _ELM_KINETIC_LETTERS + ] ) """ diff --git a/src/KineticForces/CalculatedKineticMatrices.jl b/src/KineticForces/CalculatedKineticMatrices.jl index 9d8a19def..492a62403 100644 --- a/src/KineticForces/CalculatedKineticMatrices.jl +++ b/src/KineticForces/CalculatedKineticMatrices.jl @@ -108,7 +108,7 @@ function compute_calculated_kinetic_matrices( if psi_c > 0 env .= kinetic_axis_validity_envelope.(xs, psi_c) @info "Kinetic axis-validity suppression: psi_c=$(round(psi_c; sigdigits=3)), envelope reaches 1 at " * - "psi=$(round(2 * psi_c; sigdigits=3)); $(count(iszero, env)) of $mpsi surfaces skipped" + "psi=$(round(2 * psi_c; sigdigits=3)) (kernel evaluation skipped below psi_c)" maxlog = 1 end end diff --git a/src/KineticForces/Utils.jl b/src/KineticForces/Utils.jl index 97cfa4665..e2ab60c43 100644 --- a/src/KineticForces/Utils.jl +++ b/src/KineticForces/Utils.jl @@ -143,3 +143,70 @@ function kinetic_axis_validity_envelope(psi::Float64, psi_c::Float64) t >= 1 && return 1.0 return t^3 * (10 + t * (6 * t - 15)) end + +""" + kinetic_validity_profiles(kinetic_profiles, equil; zi=1, mi=2, electron=false) → NamedTuple + +Radial profiles of the drift-kinetic validity diagnostics, on the kinetic-profile ψ grid: +`psi`; the thermal orbit-width scales `rho_i` (gyroradius √(2mT)/(Z·e·B₀)), `rho_banana` +(q·ρ/√ε), `rho_theta` (poloidal gyroradius q·ρ/ε), `w_potato` ((q²ρ²R₀)^(1/3)); the local +geometry `r_minor` (⟨r⟩) and `d_separatrix` (⟨r⟩(1) − ⟨r⟩(ψ)); the profile gradient lengths +`L_p` and `L_q` (|X|·|dr/dψ|/|dX/dψ|, Inf where the profile is flat); the near-axis suppression +boundary `psi_c` (`kinetic_axis_validity_psi`) with its `envelope`; and `is_valid` — true where +every zero-orbit-width ordering holds at coefficient 1: max orbit width < r, ρ_banana < L_p and +L_q, and max(ρ_i, ρ_banana, ρ_θ) < d_separatrix. Validity is diagnostic only — nothing outside +the near-axis envelope is suppressed (the far edge and steep-gradient regions are flagged, not +zeroed, since they can dominate the physical NTV). +""" +function kinetic_validity_profiles(kinetic_profiles::Equilibrium.KineticProfileSplines, equil; + zi::Int=1, mi::Int=2, electron::Bool=false) + chrg = electron ? e : zi * e + mass = electron ? me : mi * mp + T_spline = electron ? kinetic_profiles.Te_spline : kinetic_profiles.Ti_spline + q_spline, q_deriv = equil.profiles.q_spline, equil.profiles.q_deriv + P_spline, P_deriv = equil.profiles.P_spline, equil.profiles.P_deriv + avg_r = equil.geometry.avg_r_spline + avg_R = equil.geometry.avg_R_spline + r_deriv = deriv1(avg_r) + ro = abs(equil.ro) + bo = abs(equil.params.b0) + r_sep = avg_r(1.0) + psi_c = kinetic_axis_validity_psi(kinetic_profiles, equil; zi=zi, mi=mi, electron=electron) + + psi = [x for x in kinetic_profiles.xs if 0 < x <= 1] + n = length(psi) + rho_i = zeros(n) + rho_banana = zeros(n) + rho_theta = zeros(n) + w_potato = zeros(n) + r_minor = zeros(n) + L_p = zeros(n) + L_q = zeros(n) + d_separatrix = zeros(n) + envelope = zeros(n) + is_valid = falses(n) + for (i, x) in pairs(psi) + r = avg_r(x) + eps = max(r / avg_R(x), 1e-6) + rho = mass * sqrt(2 * T_spline(x) / mass) / (abs(chrg) * bo) + q = abs(q_spline(x)) + drdpsi = abs(r_deriv(x)) + rho_i[i] = rho + rho_banana[i] = q * rho / sqrt(eps) + rho_theta[i] = q * rho / eps + w_potato[i] = cbrt(q^2 * rho^2 * ro) + r_minor[i] = r + d_separatrix[i] = max(r_sep - r, 0.0) + dP = abs(P_deriv(x)) + L_p[i] = dP > 0 ? abs(P_spline(x)) * drdpsi / dP : Inf + dq = abs(q_deriv(x)) + L_q[i] = dq > 0 ? q * drdpsi / dq : Inf + envelope[i] = kinetic_axis_validity_envelope(x, psi_c) + w_orbit = max(w_potato[i], rho_banana[i], rho_theta[i]) + is_valid[i] = w_orbit < r && rho_banana[i] < L_p[i] && rho_banana[i] < L_q[i] && + max(rho_i[i], rho_banana[i], rho_theta[i]) < d_separatrix[i] + end + return (psi=psi, rho_i=rho_i, rho_banana=rho_banana, rho_theta=rho_theta, w_potato=w_potato, + r_minor=r_minor, L_p=L_p, L_q=L_q, d_separatrix=d_separatrix, + psi_c=psi_c, envelope=envelope, is_valid=is_valid) +end From 3f1dbfef948d2b07511f9195e412bf0710767d95 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Thu, 20 Aug 2026 12:45:05 -0400 Subject: [PATCH 07/23] KF - BUGFIX - Resolve the validity-envelope band on coarse grids; open-or-create the KineticForces group The envelope has structure on the psi_c scale; coarse kinetic decks (m16) cannot represent env*(increment) and the spline overshoot can land on a rational surface inside the transition band, corrupting the eigenvalues. Augment the kernel grid with knots across [psi_c, 2 psi_c] (band ends pinned -- the smoothstep is only C2 there) on the full-grid path and seed them on the certified path. Also open-or-create KineticForces in the NTV writer, which collided with the Validity group. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01LzbLFQKyuRE5DYZmLokKmk --- src/ForceFreeStates/Kinetic.jl | 27 ++++++++++++++++++++------ src/GeneralizedPerturbedEquilibrium.jl | 10 +++++++++- src/KineticForces/Output.jl | 25 +++++++++++++----------- 3 files changed, 44 insertions(+), 18 deletions(-) diff --git a/src/ForceFreeStates/Kinetic.jl b/src/ForceFreeStates/Kinetic.jl index da1202bae..fa9f954bb 100644 --- a/src/ForceFreeStates/Kinetic.jl +++ b/src/ForceFreeStates/Kinetic.jl @@ -26,11 +26,19 @@ function make_kinetic_matrix( ffit::FourFitVars, intr::ForceFreeStatesInternal, metric::MetricData; - calculated_source::Union{Nothing,Function}=nothing + calculated_source::Union{Nothing,Function}=nothing, + axis_validity_psi_c::Float64=0.0 ) xs = metric.xs mpsi = length(xs) + # The near-axis validity envelope (KineticForces) has structure on the scale of the + # suppression boundary; coarse equilibrium grids cannot represent env·(increment), and the + # spline overshoot can land on a rational surface. Pin the band ends (the smoothstep is + # only C² there) and resolve the transition with a fixed set of knots. + band_knots(lo, hi) = axis_validity_psi_c > 0 ? + [x for x in range(axis_validity_psi_c, 2 * axis_validity_psi_c; length=9) if lo < x < hi] : Float64[] + # Get raw kinetic matrices (scaling is baked into each source) if ctrl.kinetic_source == "fixed" kw_flat, kt_flat = fixed_kinetic_matrices(intr.mpert, mpsi, ctrl.kinetic_factor, intr.mlow, ffit, xs) @@ -41,7 +49,14 @@ function make_kinetic_matrix( "calling make_kinetic_matrix directly, or pass " * "`calculated_source=KineticForces.compute_calculated_kinetic_matrices` explicitly." ) - kw_flat, kt_flat = calculated_source(ctrl, equil, intr, metric, ffit) + band = band_knots(xs[1], xs[end]) + if isempty(band) + kw_flat, kt_flat = calculated_source(ctrl, equil, intr, metric, ffit) + else + xs = sort!(unique!(vcat(collect(xs), band))) + mpsi = length(xs) + kw_flat, kt_flat = calculated_source(ctrl, equil, intr, metric, ffit; psis=xs) + end kw_flat .*= ctrl.kinetic_factor kt_flat .*= ctrl.kinetic_factor else @@ -55,13 +70,13 @@ function make_kinetic_matrix( end # Pre-compute FKG derived matrices (corresponds to Fortran method=0) - _compute_fkg_matrices!(ffit, equil, intr, metric, kw_flat, kt_flat) + _compute_fkg_matrices!(ffit, equil, intr, metric, kw_flat, kt_flat; xs=xs) return nothing end """ - _compute_fkg_matrices!(ffit, equil, intr, metric, kw_flat, kt_flat) + _compute_fkg_matrices!(ffit, equil, intr, metric, kw_flat, kt_flat; xs=xs) Pre-compute the derived F, K, G kinetic matrices at each ψ grid point and store as splines. This corresponds to `fourfit_kinetic_matrix` method=0 in the Fortran code (Fortran `fourfit.F` lines 1170-1260). @@ -78,9 +93,9 @@ function _compute_fkg_matrices!( intr::ForceFreeStatesInternal, metric::MetricData, kw_flat::Array{ComplexF64,3}, - kt_flat::Array{ComplexF64,3} + kt_flat::Array{ComplexF64,3}; + xs::Vector{Float64}=metric.xs ) - xs = metric.xs mpsi = length(xs) np = intr.numpert_total mpert = intr.mpert diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index d1d85d22c..52d255a07 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -424,8 +424,16 @@ function main_from_inputs( KineticForces.compute_calculated_kinetic_matrices( c, e, i, m, f; kf_ctrl=kf_ctrl, kinetic_profiles=kinetic_profiles, psis=psis) + # Near-axis validity boundary for the calculated kinetic matrices: the envelope band + # must be resolved by the kernel grid, so make_kinetic_matrix needs its location. + axis_psi_c = 0.0 + if ctrl.kinetic_source == "calculated" && kf_ctrl.axis_validity_suppression && kinetic_profiles !== nothing + axis_psi_c = KineticForces.kinetic_axis_validity_psi( + kinetic_profiles, equil; + zi=kf_ctrl.zi, mi=kf_ctrl.mi, electron=kf_ctrl.electron) + end make_kinetic_matrix(ctrl, equil, ffit, intr, metric; - calculated_source=calculated_cb) + calculated_source=calculated_cb, axis_validity_psi_c=axis_psi_c) # 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; diff --git a/src/KineticForces/Output.jl b/src/KineticForces/Output.jl index dc613c04f..975ed6f03 100644 --- a/src/KineticForces/Output.jl +++ b/src/KineticForces/Output.jl @@ -12,14 +12,15 @@ then write to gpec.h5 in a single pass. Write all KineticForces results to the "KineticForces" group in gpec.h5. # Arguments -- `h5file::HDF5.File`: Open HDF5 file handle -- `state::KineticForcesState`: Accumulated computation results -- `dVdpsi_spline`: Optional dV/dψ_N profile interpolant; when given, dV/dψ_N is - written at the quadrature points so the torque density dT/dV = (dT/dψ)/(dV/dψ) - is directly available + + - `h5file::HDF5.File`: Open HDF5 file handle + - `state::KineticForcesState`: Accumulated computation results + - `dVdpsi_spline`: Optional dV/dψ_N profile interpolant; when given, dV/dψ_N is + written at the quadrature points so the torque density dT/dV = (dT/dψ)/(dV/dψ) + is directly available """ function write_to_hdf5!(h5file::HDF5.File, state::KineticForcesState; dVdpsi_spline=nothing) - g = create_group(h5file, "KineticForces") + g = haskey(h5file, "KineticForces") ? h5file["KineticForces"] : create_group(h5file, "KineticForces") for (method_name, result) in state.method_results mg = create_group(g, method_name) @@ -113,8 +114,9 @@ Write variable-length integration trajectory records using offset-indexed concat This is the standard HDF5 ragged array pattern for storing variable-length data. # Arguments -- `mg::HDF5.Group`: HDF5 group for this method -- `records::Vector{EnergyIntegrationResult}`: Integration records to write + + - `mg::HDF5.Group`: HDF5 group for this method + - `records::Vector{EnergyIntegrationResult}`: Integration records to write """ function write_integration_records!(mg::HDF5.Group, records::Vector{EnergyIntegrationResult}) rg = create_group(mg, "EnergyIntegrals") @@ -146,13 +148,14 @@ end Print a summary of KineticForces results to stdout. # Arguments -- `state::KineticForcesState`: Accumulated computation results -- `verbose::Bool`: Print detailed per-surface results + + - `state::KineticForcesState`: Accumulated computation results + - `verbose::Bool`: Print detailed per-surface results """ function print_summary(state::KineticForcesState; verbose::Bool=false) for (method_name, result) in state.method_results @printf("%-8s T_phi = %11.3e 2n*dW_k = %11.3e\n", - method_name, real(result.total_torque), imag(result.total_torque)) + method_name, real(result.total_torque), imag(result.total_torque)) end if verbose for (method_name, _) in state.kinetic_matrices From 12fca467cc7c9e0aeb4022c42a7caedc90fda655 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Thu, 20 Aug 2026 13:41:47 -0400 Subject: [PATCH 08/23] KF - MINOR - Add units to the Validity psi_c annotation (schema metadata contract) Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01LzbLFQKyuRE5DYZmLokKmk --- src/HDF5Schema.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/HDF5Schema.jl b/src/HDF5Schema.jl index 06487b7ea..e1a73bb4e 100644 --- a/src/HDF5Schema.jl +++ b/src/HDF5Schema.jl @@ -151,7 +151,7 @@ const MAIN_H5_ANNOTATIONS = [ "KineticForces/Validity/L_p" => (; long_name="pressure gradient scale length |p|/|dp/dr|", units="m", attach=(1 => "KineticForces/Validity/psi",)), "KineticForces/Validity/L_q" => (; long_name="safety-factor gradient scale length |q|/|dq/dr|", units="m", attach=(1 => "KineticForces/Validity/psi",)), "KineticForces/Validity/d_separatrix" => (; long_name="distance to the separatrix ⟨r⟩(1) − ⟨r⟩(ψ)", units="m", attach=(1 => "KineticForces/Validity/psi",)), - "KineticForces/Validity/psi_c" => (; long_name="near-axis kinetic validity boundary: outermost ψ_N where a thermal orbit width reaches ⟨r⟩"), + "KineticForces/Validity/psi_c" => (; long_name="near-axis kinetic validity boundary: outermost ψ_N where a thermal orbit width reaches ⟨r⟩", units="1"), "KineticForces/Validity/envelope" => (; long_name="near-axis suppression envelope applied to the calculated kinetic terms (1 = unsuppressed)", units="1", attach=(1 => "KineticForces/Validity/psi",)), "KineticForces/Validity/is_valid" => (; From 0def3533b66d90310fcb2eede0dca669cd3a045d Mon Sep 17 00:00:00 2001 From: logan-nc Date: Fri, 21 Aug 2026 17:52:15 -0400 Subject: [PATCH 09/23] EQUIL - NEW FEATURE - Pin located kinetic-resonance surfaces into the auto psi grid The two-pass auto grid's criterion is ideal-driven and knows nothing about kinetic resonance locations. When a run builds calculated kinetic matrices, locate the Omega_l = 0 surfaces (same locator as the NTV quadrature paneling) and insert them as plain knots via merge_mandatory_nodes -- knot-at-node, no cleared zone, inserted before rational bracketing so the Delta-prime clean-interval treatment wins locally. Nodes inside the near-axis validity region are suppressed anyway and not pinned. DIII-D: +2 net knots, et[1] unchanged to 2e-6, EL steps drop 14%. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01LzbLFQKyuRE5DYZmLokKmk --- src/Equilibrium/GridRefinement.jl | 8 +++++++- src/GeneralizedPerturbedEquilibrium.jl | 26 +++++++++++++++++++++++++- 2 files changed, 32 insertions(+), 2 deletions(-) diff --git a/src/Equilibrium/GridRefinement.jl b/src/Equilibrium/GridRefinement.jl index e58eb81dc..d9441bbe8 100644 --- a/src/Equilibrium/GridRefinement.jl +++ b/src/Equilibrium/GridRefinement.jl @@ -418,7 +418,7 @@ Build the refined pass-2 ψ grid from a formed pass-1 equilibrium: measured-curv density (`_knot_density`), equidistribution, a global minimum-spacing floor (`enforce_min_spacing`), and rational-surface bracketing (`bracket_mandatory_nodes`). `tau` is the target interpolation accuracy (`psi_accuracy`); `kin` optionally supplies kinetic profiles -whose pedestal gradients attract knots; `mandatory` lists rational-surface ψ values to bracket; +whose pedestal gradients attract knots; `mandatory` lists rational-surface ψ values to bracket; `pinned` lists ψ values inserted as plain knots without a cleared zone (kinetic-resonance surfaces); `singfac_min` and `n_min` (smallest |n| in the run) set each surface's matching half-stencil `dpsi = singfac_min/(n_min·|q′|)`, and the bracket half-width is `bracket_coef·dpsi` (floored at `min_spacing`). Rational surfaces are bracketed, not pinned: a knot on the surface would make the @@ -428,6 +428,7 @@ function refined_psi_grid(equil::PlasmaEquilibrium; tau::Float64, kin::Union{Nothing,KineticProfileSplines}=nothing, mandatory::Vector{Float64}=Float64[], + pinned::Vector{Float64}=Float64[], singfac_min::Float64=1e-4, n_min::Int=1, bracket_coef::Float64=BRACKET_COEF, @@ -449,6 +450,11 @@ function refined_psi_grid(equil::PlasmaEquilibrium; N == N_cap && M_total > N_cap && @warn "refined_psi_grid: knot count capped at $N_cap (density integral wants $(ceil(Int, M_total))); psi_accuracy=$tau may not be attainable" grid = enforce_min_spacing(_equidistribute(xs, rho, N), min_spacing) + # Pinned knots (e.g. kinetic-resonance surfaces): knot-at-node semantics via + # merge_mandatory_nodes — no cleared zone. Inserted before rational bracketing, so a + # pinned node inside a rational's bracket zone is cleared by it (the Δ′ clean-interval + # requirement wins locally; the rational's own dense floor resolves that neighbourhood). + grid = isempty(pinned) ? grid : merge_mandatory_nodes(grid, pinned) isempty(mandatory) && return grid min_half_widths = [max(bracket_coef * singfac_min / (n_min * abs(equil.profiles.q_deriv(m))), min_spacing) for m in mandatory] return bracket_mandatory_nodes(grid, mandatory, min_half_widths, min_spacing) diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 52d255a07..b95d2f263 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -249,8 +249,32 @@ function main_from_inputs( # Smallest |n| in the run sets the widest matching half-stencil dpsi = singfac_min/(n_min·|q′|), # so the rational-surface brackets clear a zone large enough for every mode. n_min = minimum(abs(n) for n in intr.nlow:intr.nhigh if n != 0) + # Pin the located kinetic-resonance surfaces (Ω_ℓ = 0) as plain knots when the run + # builds calculated kinetic matrices — the ideal-driven grid criterion knows nothing + # about kinetic resonance locations. Nodes inside the near-axis validity region are + # suppressed anyway and not pinned. + pinned = Float64[] + if ctrl.kinetic_factor > 0 && ctrl.kinetic_source == "calculated" && kinetic_profiles !== nothing + for n_res in intr.nlow:intr.nhigh + n_res == 0 && continue + append!( + pinned, + 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 + if kf_ctrl.axis_validity_suppression + psi_c_grid = KineticForces.kinetic_axis_validity_psi(kinetic_profiles, equil; + zi=kf_ctrl.zi, mi=kf_ctrl.mi, electron=kf_ctrl.electron) + filter!(p -> p > psi_c_grid, pinned) + end + isempty(pinned) || + @info "Pinning $(length(pinned)) kinetic-resonance surfaces into the ψ grid: $(round.(sort(pinned); digits=3))" + end psi_nodes = Equilibrium.refined_psi_grid(equil; - tau=eq_config.psi_accuracy, kin=kinetic_profiles, mandatory=mandatory, + tau=eq_config.psi_accuracy, kin=kinetic_profiles, mandatory=mandatory, pinned=pinned, singfac_min=ctrl.singfac_min, n_min=n_min) rerun_input = if additional_input !== nothing # Analytic *Config, IMAS dd, or prebuilt RunInput — all re-formable. The IMAS From 168857b8bd7f102bb8056534e553476a6ba30889 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Sat, 22 Aug 2026 13:45:59 -0400 Subject: [PATCH 10/23] PE/KF - BUGFIX! - Disable ideal-singularity regularization in self-consistent kinetic runs reg_spot smooths the ideal 1/(m-nq) divergence of the displacements before they drive the NTV integrand. The self-consistent kinetic Euler-Lagrange operator has no such divergence -- det(F-bar) is complex and nonzero at the rationals (Park & Logan, Phys. Plasmas 24, 032505 (2017) Eq. 70) -- so regularizing there suppresses a finite physical response, and inconsistently: xi^psi is never regularized, so damping the other two breaks their near-resonance cancellation in dB/B. Measured (DIII-D-like, n=1) against the EL solution's own dissipation: NTV torque 0.1655 vs 0.1322 N*m with reg_spot=0.05, and 0.1324 vs 0.1322 (0.15%) with it off. Kinetic runs now force reg_spot=0 and log the override; ideal runs are unchanged, where without it the displacement and torque diverge by four orders of magnitude. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01LzbLFQKyuRE5DYZmLokKmk --- docs/src/kinetic_forces.md | 36 +++++++++++++++++++ src/GeneralizedPerturbedEquilibrium.jl | 13 ++++++- .../FieldReconstruction.jl | 13 ++++--- .../PerturbedEquilibriumStructs.jl | 22 ++++++------ 4 files changed, 68 insertions(+), 16 deletions(-) diff --git a/docs/src/kinetic_forces.md b/docs/src/kinetic_forces.md index bfdf825ca..e81d59aff 100644 --- a/docs/src/kinetic_forces.md +++ b/docs/src/kinetic_forces.md @@ -89,6 +89,42 @@ terms respectively, and do not modify the stored kinetic profile splines. when computing the `toroidal_rotation_factor` back-solve. Julia uses a clean reimplementation with consistent pre-scaling derivatives throughout. +## Regularization: why kinetic runs set `reg_spot = 0` + +`[PerturbedEquilibrium] reg_spot` smooths the displacements before they drive the NTV +integrand, multiplying ``\xi^{\psi\prime}`` and ``\xi^\alpha`` by +``Q^2/(Q^2 + \mathrm{reg\_spot}^2)`` with ``Q = m - nq``. It exists because **ideal** MHD is +singular at the rationals: ``\bar F_\mathrm{ideal} = Q F Q`` has ``\det \bar F = 0`` there, so +those two components diverge as ``1/Q`` and the torque integral does not converge. + +The **self-consistent kinetic** workflow has no such singularity. Park & Logan +([Phys. Plasmas 24, 032505 (2017)](https://doi.org/10.1063/1.4978562), §III D) decompose the +kinetic composite matrix as ``F_k = Q \bar F_k Q - P_l^\dagger Q - Q P_u + R_1`` where +``R_1 \neq 0`` at ``Q = 0``; with finite torque ``\det \bar F`` is complex, its zeros leave the +real ``\psi`` axis, and the singularity is removed from both the solution and the torque +integral. (Torque-free kinetic energy principles instead *shift and split* the zeros to +``\psi_r \mp r_{L,R}``, where the singularity is logarithmic and integrable — still not a case +for smoothing.) + +GPEC therefore **forces `reg_spot = 0` whenever `kinetic_factor > 0`**, logging the override. +Leaving it on suppresses a finite physical response and does so inconsistently — ``\xi^\psi`` +is never regularized, so damping the other two breaks their near-resonance cancellation in +``\delta B/B`` and leaves a spurious residue driving the NTV integrand. + +Measured on the DIII-D-like H-mode case (n = 1, C-coil drive), comparing the NTV torque against +the Euler–Lagrange solution's own dissipation ``-2n\,\mathrm{Im}\langle \xi, u_2\rangle/4\mu_0`` +— two independent calculations of the same quantity: + +| configuration | max ``|\xi^\alpha|`` | NTV torque [N·m] | EL dissipation [N·m] | +|---|---|---|---| +| ideal, `reg_spot = 0` | 643.8 | 6074.3 | — | +| ideal, `reg_spot = 0.05` | 0.058 | 0.554 | — | +| kinetic, `reg_spot = 0.05` | 0.055 | 0.1655 | 0.1322 | +| kinetic, `reg_spot = 0` | 0.059 | **0.1324** | **0.1322** | + +The ideal rows show why the knob exists; the kinetic rows show why it must be off there — the +two independent torques agree to 0.15 % with no regularization, and to 20 % with it. + ## HDF5 outputs: complex torque convention and the EnergyIntegrals layout The method level of `KineticForces//` reports the two physical scalars a user diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index b95d2f263..9dbd50752 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -603,8 +603,19 @@ function main_from_inputs( ft_ctrl = ForcingTerms.ForcingTermsControl() # Use defaults end + pe_raw_in = inputs["PerturbedEquilibrium"] + # reg_spot smooths the ideal 1/(m−nq) divergence of ξ^ψ′ and ξ^α before they reach the NTV + # integrand. The self-consistent kinetic Euler-Lagrange operator has no such divergence — + # det(F̄) is complex and nonzero at the rationals — so regularizing there suppresses a finite + # physical response, and does so inconsistently (ξ^ψ is left unregularized, breaking the + # near-resonance cancellation in δB/B). Park & Logan, Phys. Plasmas 24, 032505 (2017) §III D. + if ctrl.kinetic_factor > 0 && get(pe_raw_in, "reg_spot", 0.0) != 0 + @info "Self-consistent kinetic run: overriding reg_spot=$(pe_raw_in["reg_spot"]) with 0 " * + "(the kinetic terms remove the ideal resonant singularity; see docs/src/kinetic_forces.md)" + pe_raw_in = merge(pe_raw_in, Dict("reg_spot" => 0.0)) + end pe_ctrl = PerturbedEquilibrium.PerturbedEquilibriumControl(; - (Symbol(k) => v for (k, v) in inputs["PerturbedEquilibrium"])... + (Symbol(k) => v for (k, v) in pe_raw_in)... ) pe_intr = PerturbedEquilibrium.PerturbedEquilibriumInternal(; dir_path=intr.dir_path) diff --git a/src/PerturbedEquilibrium/FieldReconstruction.jl b/src/PerturbedEquilibrium/FieldReconstruction.jl index e491477c4..1b0e5f8e3 100644 --- a/src/PerturbedEquilibrium/FieldReconstruction.jl +++ b/src/PerturbedEquilibrium/FieldReconstruction.jl @@ -19,6 +19,7 @@ where χ₁ = 2π·Ψ₀ [Park Phys. Plasmas 14, 052110 (2007) eq. 8-10]. Clebsch displacement components for PENTRC (matches Fortran gpout_xclebsch): ξ^ψ = xsp_mn (unregularized) ∂ξ^ψ/∂ψ = xmp1_mn (regularized: xsp1 * singfac²/(singfac² + reg_spot²)) + reg_spot = 0 in kinetic runs — no ideal singularity to smooth ξ^α = xms_mn (regularized: -A⁻¹(B·xmp1 + C·xsp), divided by χ₁ in output) Contravariant displacement from Jacobian convolution (matches Fortran gpeq_contra): @@ -404,7 +405,7 @@ function compute_clebsch_displacements( # xms = -(A\B)*xmp1 - (A\C)*xsp xsp_vec = view(xi_psi_modes, ipsi, :) mul!(xms_vec, bmat, xmp1_vec) # xms = B*xmp1 - mul!(xms_vec, cmat_buf, xsp_vec, 1.0+0.0im, 1.0+0.0im) # xms += C*xsp + mul!(xms_vec, cmat_buf, xsp_vec, 1.0 + 0.0im, 1.0 + 0.0im) # xms += C*xsp # amat is positive-definite by construction (Newcomb kinetic-energy form), so cholesky is # safe. cholesky! factorizes in place (amat is a per-thread scratch buffer, refilled by # ffit.amats each surface), avoiding a fresh factorization allocation per surface. @@ -980,10 +981,12 @@ function _apply_rzphi_transform( # Per-thread scratch (the immutable `ft` functor and `geom` are shared read-only): θ-space # transform inputs/outputs (length mtheta) and mode-space forward-DFT outputs (length mpert), # so the DFTs run in place with no per-surface allocation. - bufs = [(R=zeros(ComplexF64, mtheta), Z=zeros(ComplexF64, mtheta), P=zeros(ComplexF64, mtheta), - psi=zeros(ComplexF64, mtheta), th=zeros(ComplexF64, mtheta), ze=zeros(ComplexF64, mtheta), - Ro=zeros(ComplexF64, mpert), Zo=zeros(ComplexF64, mpert), Po=zeros(ComplexF64, mpert)) - for _ in 1:Threads.maxthreadid()] + bufs = [ + (R=zeros(ComplexF64, mtheta), Z=zeros(ComplexF64, mtheta), P=zeros(ComplexF64, mtheta), + psi=zeros(ComplexF64, mtheta), th=zeros(ComplexF64, mtheta), ze=zeros(ComplexF64, mtheta), + Ro=zeros(ComplexF64, mpert), Zo=zeros(ComplexF64, mpert), Po=zeros(ComplexF64, mpert)) + for _ in 1:Threads.maxthreadid() + ] Threads.@threads :static for ipsi in 1:npsi buf = bufs[Threads.threadid()] diff --git a/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl b/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl index 1ce0895d6..e44edb522 100644 --- a/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl +++ b/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl @@ -26,8 +26,9 @@ Medium Priority (defer for MWE): - `singular_point_method::String` - Method for singular point treatment (default: "standard") Regularization: - # High Priority (MWE) - - `reg_spot::Float64` - Regularization width for singular surface smoothing (default: 0.05). Set to 0 to disable. Must be ≥ 0. +# High Priority (MWE) + + - `reg_spot::Float64` - Regularization width for singular surface smoothing (default: 0.05). Set to 0 to disable. Must be ≥ 0. Forced to 0 in self-consistent kinetic runs, whose Euler-Lagrange operator has no resonant singularity to smooth. """ @kwdef struct PerturbedEquilibriumControl # High Priority (MWE) @@ -121,6 +122,7 @@ Metadata [n_rational] — identifies each (surface, n) row: Control-surface forcing/response spectra [numpert_total], in the three Pharr (2026) field representations (all tesla; no flux/weber is stored): + - `forcing_b`/`response_b` - bare normal field b (Σ⁻¹·b̃) - `forcing_b_rootarea`/`response_b_rootarea` - root-area-weighted field b̃ (coordinate-invariant) - `forcing_b_area`/`response_b_area` - area-weighted field b̄ (= S·b̃; flux is Φ = A·b̄) @@ -185,18 +187,18 @@ well-conditioned flux-space inductances L, Λ: rational_surface_idx::Vector{Int} = Int[] # Control-surface forcing/response spectra in the three weightings of field representations [numpert_total], tesla - forcing_b::Vector{ComplexF64} = ComplexF64[] # bare normal field b (forcing Φ_x) - forcing_b_rootarea::Vector{ComplexF64} = ComplexF64[] # root-area-weighted field b̃ (coordinate-invariant) - forcing_b_area::Vector{ComplexF64} = ComplexF64[] # area-weighted field b̄ - response_b::Vector{ComplexF64} = ComplexF64[] # bare normal field b (response Φ_tot = P·Φ_x) + forcing_b::Vector{ComplexF64} = ComplexF64[] # bare normal field b (forcing Φ_x) + forcing_b_rootarea::Vector{ComplexF64} = ComplexF64[] # root-area-weighted field b̃ (coordinate-invariant) + forcing_b_area::Vector{ComplexF64} = ComplexF64[] # area-weighted field b̄ + response_b::Vector{ComplexF64} = ComplexF64[] # bare normal field b (response Φ_tot = P·Φ_x) response_b_rootarea::Vector{ComplexF64} = ComplexF64[] # root-area-weighted field b̃ - response_b_area::Vector{ComplexF64} = ComplexF64[] # area-weighted field b̄ + response_b_area::Vector{ComplexF64} = ComplexF64[] # area-weighted field b̄ # Control surface matrices [numpert_total × numpert_total], root-area-weighted field (b̃) space - plasma_inductance::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) # Λ̃ (field space) + plasma_inductance::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) # Λ̃ (field space) surface_inductance::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) # L̃ (field space) - permeability::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) # P̃ = R⁻¹·Λ·L⁻¹·R - reluctance::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) # ϱ̃ = R†·L⁻¹·(Λ−L)·L⁻¹·R + permeability::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) # P̃ = R⁻¹·Λ·L⁻¹·R + reluctance::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) # ϱ̃ = R†·L⁻¹·(Λ−L)·L⁻¹·R rootarea_to_area_weight::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) # S = Σ/√A at psilim: b̃→b̄ recovery operator surface_area::Float64 = 0.0 # scalar control-surface area A = ∫J|∇ψ|dθ (flux: Φ = A·b̄; conform R = S·A) From 5aa11223c7123820f07aa73fe800d12987eb7967 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Sat, 22 Aug 2026 16:35:33 -0400 Subject: [PATCH 11/23] FFS - IMPROVEMENT - Report sub-threshold near-singular kinetic F-bar structure The cond(F-bar) scan that locates kinetic singular surfaces already sweeps 2000 points and is written to SingularSurfaces/Kinetic/scan_cond, but only peaks above the 1e8 singular threshold were surfaced. On a DIII-D-like case the strongest peaks sit at 3e5-5e6 -- real shifted/split resonance structure (Park & Logan Eq. 70) that stayed silent, and which tracks the NTV torque-density peaks at low collisionality/rotation. Report the strongest few. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01LzbLFQKyuRE5DYZmLokKmk --- src/ForceFreeStates/Sing.jl | 27 +++++++++++++++++++++------ 1 file changed, 21 insertions(+), 6 deletions(-) diff --git a/src/ForceFreeStates/Sing.jl b/src/ForceFreeStates/Sing.jl index f41e6dd2f..363c094d1 100644 --- a/src/ForceFreeStates/Sing.jl +++ b/src/ForceFreeStates/Sing.jl @@ -123,9 +123,11 @@ function sing_lim!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl, # strategy. Multi-n runs are not supported — the "outermost rational + dmlim/n" cutoff depends # on which n is used — and fall back to qhigh / psihigh truncation with a warning. if ctrl.set_psilim_via_dmlim && intr.nlow <= 0 - error("sing_lim!: set_psilim_via_dmlim = true requires a resolved toroidal range, but got intr.nlow=$(intr.nlow). " * - "Assign intr.nlow / intr.nhigh (from ctrl.nn_low / ctrl.nn_high) before calling sing_lim!, " * - "or set set_psilim_via_dmlim = false to truncate via qhigh / psihigh instead.") + error( + "sing_lim!: set_psilim_via_dmlim = true requires a resolved toroidal range, but got intr.nlow=$(intr.nlow). " * + "Assign intr.nlow / intr.nhigh (from ctrl.nn_low / ctrl.nn_high) before calling sing_lim!, " * + "or set set_psilim_via_dmlim = false to truncate via qhigh / psihigh instead." + ) elseif ctrl.set_psilim_via_dmlim && intr.nlow != intr.nhigh @warn "set_psilim_via_dmlim = true is ignored for multi-n runs (nn_low=$(intr.nlow), nn_high=$(intr.nhigh)); falling back to qhigh / psihigh truncation." elseif ctrl.set_psilim_via_dmlim @@ -266,7 +268,7 @@ function compute_sing_asymptotics( # This is the parameter α but for all modes - α = 0 for non-resonant modes power[ipert_res] .= -alpha - power[ipert_res .+ intr.numpert_total] .= alpha + power[ipert_res.+intr.numpert_total] .= alpha # Zeroth-order non-resonant solutions for ipert in 1:intr.numpert_total @@ -303,7 +305,7 @@ function compute_sing_asymptotics( msg *= @sprintf(" m0mat(1,2)= %+.12e %+.12ei\n", real(m0mat[1, 2]), imag(m0mat[1, 2])) msg *= @sprintf(" m0mat(2,1)= %+.12e %+.12ei\n", real(m0mat[2, 1]), imag(m0mat[2, 1])) msg *= @sprintf(" m0mat(2,2)= %+.12e %+.12ei\n", real(m0mat[2, 2]), imag(m0mat[2, 2])) - di = m0mat[1, 1]*m0mat[2, 2] - m0mat[2, 1]*m0mat[1, 2] + di = m0mat[1, 1] * m0mat[2, 2] - m0mat[2, 1] * m0mat[1, 2] msg *= @sprintf(" di= %+.12e, alpha= %+.12e %+.12ei\n", real(di), real(alpha[1]), imag(alpha[1])) msg *= @sprintf(" psifac= %+.12e, r1=%d, ipert0=%d\n", singp.psifac, r1[1], ipert0) msg *= @sprintf(" vmat(ip,ip,2,0)= %+.8e %+.8ei\n", real(vmat[ipert0, ipert0, 2, 1]), imag(vmat[ipert0, ipert0, 2, 1])) @@ -771,7 +773,7 @@ function sing_get_ua(sing_asymp::SingAsymptotics, dpsi::Float64) # Restore powers (unshear v→u) — matches Fortran STRIDE sing_get_ua for i in eachindex(r1) - pfac = pfac_base ^ sing_asymp.alpha[i] # dpsi^α + pfac = pfac_base^sing_asymp.alpha[i] # dpsi^α ua[:, r2[2*i-1], :] ./= pfac # big solution column: /dpsi^α ua[:, r2[2*i], :] .*= pfac # small solution column: *dpsi^α ua[r1[i], :, 1] ./= sqrtfac # resonant row ξ: /√dpsi @@ -1432,6 +1434,19 @@ function find_kinetic_singular_surfaces!(ffit::FourFitVars, equil::Equilibrium.P end end + # Peaks below the threshold are not singular surfaces, but they mark where the kinetic F̄ comes + # closest to singular — the shifted/split resonances of Park & Logan Eq. (70). Report the + # strongest few so sharp kinetic structure is visible rather than silent (on a DIII-D-like case + # these track the NTV torque-density peaks at low collisionality/rotation). + subthreshold = [i for i in 2:(ngrid-1) if cond_vals[i] > cond_vals[i-1] && cond_vals[i] > cond_vals[i+1] && + cond_threshold / 100 < cond_vals[i] <= cond_threshold] + if !isempty(subthreshold) + top = sort(subthreshold; by=i -> -cond_vals[i])[1:min(3, length(subthreshold))] + @info "Kinetic F̄ near-singular structure below the singular threshold at " * + join(["ψ=$(round(psi_grid[i]; digits=4)) (cond=$(round(cond_vals[i]; sigdigits=3)))" for i in top], ", ") * + " — full scan in SingularSurfaces/Kinetic/scan_cond; check the ψ grid resolves these if results look grid-sensitive" + end + # Refine each peak to find the precise ψ location kinsing_surfaces = SingType[] for idx in peak_indices From b3d7ab9a83760861b9db2f649f89b95a6aef247a Mon Sep 17 00:00:00 2001 From: logan-nc Date: Sat, 22 Aug 2026 17:30:39 -0400 Subject: [PATCH 12/23] FFS - NEW FEATURE - Add kinetic knots across unresolved near-singular F-bar structure The cond(F-bar) scan already locates shifted/split kinetic resonances (Park & Logan Eq. 70) but only reported them. Measure each sub-threshold peak's FWHM and, where the grid puts fewer than three knots inside it, evaluate the kernel at a few targeted psi and splice them in -- existing values are reused, so the cost is one kernel call per added knot, not a re-formation. Respects MIN_KNOT_SPACING, the near-axis validity band, and a 24-knot cap; a resolved grid inserts nothing. Measured on DIII-D: no insertions on the nominal or low-collisionality cases; on the collisionless slow-rotation case two peaks (psi=0.53, 0.51, FWHM 1.5e-3 and 2.5e-3) had zero knots inside them, and one of the two coincides with a top NTV torque-density peak. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01LzbLFQKyuRE5DYZmLokKmk --- src/ForceFreeStates/Kinetic.jl | 91 ++++++++++++++++++++++++++++++++++ 1 file changed, 91 insertions(+) diff --git a/src/ForceFreeStates/Kinetic.jl b/src/ForceFreeStates/Kinetic.jl index fa9f954bb..a1c070baa 100644 --- a/src/ForceFreeStates/Kinetic.jl +++ b/src/ForceFreeStates/Kinetic.jl @@ -1,3 +1,69 @@ + +""" + refine_grid_at_fbar_peaks(xs, kw, kt, evaluate, ffit, equil, intr, psi_c; + ngrid=1000, relaxed_frac=0.01, target=3, max_add=24) → (xs, kw, kt) + +Insert kinetic evaluation knots across near-singular structure of F̄ that the grid does not +resolve. Scans cond(F̄) (the same operator `find_kinetic_singular_surfaces!` searches — Park & +Logan Eq. 70, so shifted and split resonances are included), takes peaks between +`relaxed_frac`·threshold and the singular threshold, measures each peak's FWHM, and adds knots +only where fewer than `target` knots lie inside it. New points respect `MIN_KNOT_SPACING`, stay +above the near-axis validity band, and are capped at `max_add`; each costs one kernel evaluation +and the existing values are reused. A well-resolved grid inserts nothing. +""" +function refine_grid_at_fbar_peaks(xs::Vector{Float64}, kw::Array{ComplexF64,3}, kt::Array{ComplexF64,3}, + evaluate::Function, ffit::FourFitVars, equil::Equilibrium.PlasmaEquilibrium, + intr::ForceFreeStatesInternal, psi_c::Float64; + ngrid::Int=1000, relaxed_frac::Float64=0.01, target::Int=3, max_add::Int=24, + cond_threshold::Float64=1e8) + + lo, hi = xs[1], xs[end] + scan = collect(range(lo, hi; length=ngrid)) + hint = Ref(1) + cond_vals = [ + try + evaluate_fbar_condition(x, ffit, equil, intr; hint=hint) + catch + Inf + end for x in scan + ] + + add = Float64[] + for i in 2:(ngrid-1) + c = cond_vals[i] + (c > cond_vals[i-1] && c > cond_vals[i+1] && relaxed_frac * cond_threshold < c <= cond_threshold) || continue + l = i + while l > 1 && cond_vals[l] > c / 2 + l -= 1 + end + r = i + while r < ngrid && cond_vals[r] > c / 2 + r += 1 + end + inside = count(x -> scan[l] <= x <= scan[r], xs) + inside >= target && continue + w = (scan[r] - scan[l]) / 3 + for x in (scan[i], scan[i] - w, scan[i] + w) + (lo < x < hi && x > 2 * psi_c) || continue + minimum(abs.(xs .- x)) < Equilibrium.MIN_KNOT_SPACING && continue + isempty(add) || minimum(abs.(add .- x)) >= Equilibrium.MIN_KNOT_SPACING || continue + push!(add, x) + end + end + isempty(add) && return xs, kw, kt + length(add) > max_add && (add = sort(add)[1:max_add]) + + sort!(add) + @info "Kinetic grid: $(length(add)) knot(s) added across unresolved near-singular F̄ structure at " * + "ψ=$(round.(add; digits=4)) (cond peaks below the singular threshold)" + kw_new, kt_new = evaluate(add) + allxs = vcat(xs, add) + perm = sortperm(allxs) + kw_all = cat(kw, kw_new; dims=1)[perm, :, :] + kt_all = cat(kt, kt_new; dims=1)[perm, :, :] + return allxs[perm], kw_all, kt_all +end + """ make_kinetic_matrix(ctrl, equil, ffit, intr, metric; calculated_source=nothing) @@ -72,6 +138,31 @@ function make_kinetic_matrix( # Pre-compute FKG derived matrices (corresponds to Fortran method=0) _compute_fkg_matrices!(ffit, equil, intr, metric, kw_flat, kt_flat; xs=xs) + # The FKG splines now exist, so F̄ can be scanned: add knots only where near-singular + # structure (shifted/split kinetic resonances) falls in an interval that does not resolve it. + if ctrl.kinetic_source == "calculated" && calculated_source !== nothing + xs2, kw_flat, kt_flat = refine_grid_at_fbar_peaks( + collect(xs), kw_flat, kt_flat, + psis -> begin + kwn, ktn = calculated_source(ctrl, equil, intr, metric, ffit; psis=psis) + (kwn .* ctrl.kinetic_factor, ktn .* ctrl.kinetic_factor) + end, + ffit, equil, intr, axis_validity_psi_c) + if length(xs2) != length(xs) + xs = xs2 + # _compute_fkg_matrices! folds the kinetic increments into amats/bmats/cmats (saving the + # ideal copies), so restore those before recomputing or the increments are added twice. + ffit.amats = ffit.amats_ideal + ffit.bmats = ffit.bmats_ideal + ffit.cmats = ffit.cmats_ideal + for ic in 1:6 + ffit.kwmats[ic] = cubic_interp(xs, Series(@view(kw_flat[:, :, ic])); ffit.itp_opts...) + ffit.ktmats[ic] = cubic_interp(xs, Series(@view(kt_flat[:, :, ic])); ffit.itp_opts...) + end + _compute_fkg_matrices!(ffit, equil, intr, metric, kw_flat, kt_flat; xs=xs) + end + end + return nothing end From f587e9becea1da59e5785866dd5f986c59b5c830 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Sat, 22 Aug 2026 18:31:50 -0400 Subject: [PATCH 13/23] PE - BUGFIX - Force reg_spot=0 in kinetic runs even when the deck omits the key The override tested get(inputs,"reg_spot",0.0)!=0, so it only fired when a deck set reg_spot explicitly; a deck relying on the struct default (0.05) silently kept regularization on in a self-consistent kinetic run -- exactly the case the change exists to prevent. Compare against the struct default and always force 0, logging whenever the prior effective value was nonzero. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01LzbLFQKyuRE5DYZmLokKmk --- src/GeneralizedPerturbedEquilibrium.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 9dbd50752..debbed227 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -609,8 +609,9 @@ function main_from_inputs( # det(F̄) is complex and nonzero at the rationals — so regularizing there suppresses a finite # physical response, and does so inconsistently (ξ^ψ is left unregularized, breaking the # near-resonance cancellation in δB/B). Park & Logan, Phys. Plasmas 24, 032505 (2017) §III D. - if ctrl.kinetic_factor > 0 && get(pe_raw_in, "reg_spot", 0.0) != 0 - @info "Self-consistent kinetic run: overriding reg_spot=$(pe_raw_in["reg_spot"]) with 0 " * + if ctrl.kinetic_factor > 0 + prev = get(pe_raw_in, "reg_spot", PerturbedEquilibrium.PerturbedEquilibriumControl().reg_spot) + prev == 0 || @info "Self-consistent kinetic run: overriding reg_spot=$prev with 0 " * "(the kinetic terms remove the ideal resonant singularity; see docs/src/kinetic_forces.md)" pe_raw_in = merge(pe_raw_in, Dict("reg_spot" => 0.0)) end From 37c4eebb4eefd362e140773323bd192bc108864a Mon Sep 17 00:00:00 2001 From: logan-nc Date: Fri, 4 Sep 2026 15:03:54 -0400 Subject: [PATCH 14/23] KF - BUGFIX! - Keep the validity envelope clear of rational surfaces A rational surface inside the envelope's transition band gets its near-singular kinetic increments multiplied by a rapidly varying, near-zero envelope, which the matrix splines cannot represent; the overshoot propagates NaNs into the stability solve. Caught by the Solovev D-T deck, where the widest-orbit species (tritium) puts psi_c at 0.12 and the rational sits at 0.1221. Push the boundary past any rational whose window the band would cut through -- where orbit widths already reach the resonance is not trustworthy anyway -- and compute that boundary once per run, threading it to the kernel, the band knots and the Validity output so all three agree. Also from the clean-code review: single orbit_widths helper feeding both the boundary and the diagnostic profiles, and Validity's psi_c/envelope/is_valid now describe one species (the widest-orbit one) instead of contradicting each other. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01LzbLFQKyuRE5DYZmLokKmk --- src/GeneralizedPerturbedEquilibrium.jl | 16 +++- .../CalculatedKineticMatrices.jl | 5 +- src/KineticForces/Output.jl | 12 ++- src/KineticForces/Utils.jl | 93 +++++++++++++++---- 4 files changed, 97 insertions(+), 29 deletions(-) diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 021ae81d0..574fa4d9f 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -583,14 +583,18 @@ function prepare_force_free_states!( # Inject the KineticForces callback so the "calculated" source can # invoke compute_calculated_kinetic_matrices without ForceFreeStates # importing KineticForces (which would invert the load order). + # One boundary for the whole run: the widest-orbit species, moved clear of any rational + # surface whose window the transition band would otherwise cut through. + axis_psi_c = KineticForces.axis_validity_boundary(kf_ctrl, species, kinetic_profiles, equil, + Float64[sng.psifac for sng in intr.sing]) calculated_cb = (c, e, i, m, f; psis::Vector{Float64}=Float64[]) -> KineticForces.compute_calculated_kinetic_matrices( c, e, i, m, f; - kf_ctrl=kf_ctrl, kinetic_profiles=kinetic_profiles, species=species, psis=psis) + kf_ctrl=kf_ctrl, kinetic_profiles=kinetic_profiles, species=species, psis=psis, + axis_psi_c=axis_psi_c) mats = build_kinetic_matrix_splines(ctrl, equil, mats, intr, metric; - calculated_source=calculated_cb, - axis_validity_psi_c=KineticForces.axis_validity_boundary(kf_ctrl, species, kinetic_profiles, equil)) + calculated_source=calculated_cb, axis_validity_psi_c=axis_psi_c) # 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; @@ -974,7 +978,8 @@ function run_kinetic_forces( h5open(joinpath(result.dir_path, kf_ctrl.HDF5_filename), "cw") do h5file KineticForces.write_to_hdf5!(h5file, kf_state; dVdpsi_spline=result.equil.profiles.dVdpsi_spline) - KineticForces.write_validity!(h5file, kf_ctrl, species, kinetic_profiles, result.equil) + KineticForces.write_validity!(h5file, kf_ctrl, species, kinetic_profiles, result.equil, + Float64[sng.psifac for sng in result.surfaces]) end end else @@ -1000,7 +1005,8 @@ function run_kinetic_forces( end KineticForces.write_to_hdf5!(h5file, kf_state; dVdpsi_spline=result.equil.profiles.dVdpsi_spline) - KineticForces.write_validity!(h5file, kf_ctrl, species, kinetic_profiles, result.equil) + KineticForces.write_validity!(h5file, kf_ctrl, species, kinetic_profiles, result.equil, + Float64[sng.psifac for sng in result.surfaces]) end end end diff --git a/src/KineticForces/CalculatedKineticMatrices.jl b/src/KineticForces/CalculatedKineticMatrices.jl index 5c1a1bc90..a698c848d 100644 --- a/src/KineticForces/CalculatedKineticMatrices.jl +++ b/src/KineticForces/CalculatedKineticMatrices.jl @@ -61,7 +61,8 @@ function compute_calculated_kinetic_matrices( kf_ctrl::KineticForcesControl=KineticForcesControl(), kinetic_profiles::Equilibrium.KineticProfileSplines, species::Union{Nothing,AbstractVector{<:Equilibrium.ResolvedNTVSpecies}}=nothing, - psis::Vector{Float64}=Float64[] + psis::Vector{Float64}=Float64[], + axis_psi_c::Float64=0.0 ) # The kernel is a pure function of psi (it evaluates equilibrium splines), so it can be # driven over any knot list; default is the full equilibrium grid. @@ -119,7 +120,7 @@ function compute_calculated_kinetic_matrices( # ordering fails, taking the widest-orbit species. Kernel evaluation is skipped where it is 0. env = ones(Float64, mpsi) if kf_ctrl.axis_validity_suppression - psi_c = kinetic_axis_validity_psi(splist, equil) + psi_c = axis_psi_c if psi_c > 0 env .= kinetic_axis_validity_envelope.(xs, psi_c) @info "Kinetic axis-validity suppression: psi_c=$(round(psi_c; sigdigits=3)), envelope reaches 1 at " * diff --git a/src/KineticForces/Output.jl b/src/KineticForces/Output.jl index 91f16959b..771873251 100644 --- a/src/KineticForces/Output.jl +++ b/src/KineticForces/Output.jl @@ -200,11 +200,15 @@ geometry, profile gradient lengths, the near-axis boundary and its applied envel far edge and steep-gradient regions are flagged rather than zeroed. """ function write_validity!(h5file::HDF5.File, kf_ctrl::KineticForcesControl, species, - kinetic_profiles, equil) + kinetic_profiles, equil, rationals::Vector{Float64}=Float64[]) kinetic_profiles === nothing && return nothing - vp = kinetic_validity_profiles(kinetic_profiles, equil; - zi=kf_ctrl.zi, mi=kf_ctrl.mi, electron=kf_ctrl.electron) - psi_c = axis_validity_boundary(kf_ctrl, species, kinetic_profiles, equil) + # One species sets all three of psi_c, envelope and is_valid: the widest-orbit species when a + # resolved set is available, else the control's single species. + vp = + species === nothing ? + kinetic_validity_profiles(kinetic_profiles, equil; zi=kf_ctrl.zi, mi=kf_ctrl.mi, electron=kf_ctrl.electron) : + kinetic_validity_profiles(species, equil) + psi_c = axis_validity_boundary(kf_ctrl, species, kinetic_profiles, equil, rationals) root = haskey(h5file, "KineticForces") ? h5file["KineticForces"] : create_group(h5file, "KineticForces") haskey(root, "Validity") && return nothing g = create_group(root, "Validity") diff --git a/src/KineticForces/Utils.jl b/src/KineticForces/Utils.jl index 2c922e1f2..80a04f1bf 100644 --- a/src/KineticForces/Utils.jl +++ b/src/KineticForces/Utils.jl @@ -90,6 +90,26 @@ function kinetic_resonance_psi_nodes(kinetic_profiles::Equilibrium.KineticProfil return _resonance_nodes_from_frequencies(wbhat_f, kinetic_profiles.omegaE_spline, wdhat_f, grid; n, nl, xeval) end +""" + orbit_widths(psi, T_spline, q_spline, avg_r, avg_R, mass, chrg, ro, bo) + +Thermal orbit-width scales of one species at ψ: the minor radius ⟨r⟩ and inverse aspect ratio ε, +the gyroradius ρ = √(2mT)/(Z·e·B₀), and the three widths that bound the zero-orbit-width ordering — +potato (q²ρ²R₀)^(1/3), banana qρ/√ε, and poloidal gyroradius qρ/ε. Single source of truth for both +the validity boundary and the diagnostic profiles. +""" +@inline function orbit_widths(psi::Float64, T_spline, q_spline, avg_r, avg_R, + mass::Float64, chrg::Float64, ro::Float64, bo::Float64) + r = avg_r(psi) + eps = max(r / avg_R(psi), 1e-6) + rho = sqrt(2 * mass * T_spline(psi)) / (abs(chrg) * bo) + q = abs(q_spline(psi)) + return (; r, eps, rho, + w_potato=cbrt(q^2 * rho^2 * ro), + rho_banana=q * rho / sqrt(eps), + rho_theta=q * rho / eps) +end + """ kinetic_axis_validity_psi(kinetic_profiles, equil; zi=1, mi=2, electron=false) → Float64 @@ -117,13 +137,9 @@ function kinetic_axis_validity_psi(kinetic_profiles::Equilibrium.KineticProfileS psi_c = 0.0 for psi in kinetic_profiles.xs psi <= 0 && continue - r = avg_r(psi) - r <= 0 && continue - eps = max(r / avg_R(psi), 1e-6) - rho = mass * sqrt(2 * T_spline(psi) / mass) / (abs(chrg) * bo) - q = abs(q_spline(psi)) - w_orbit = max(cbrt(q^2 * rho^2 * ro), q * rho / sqrt(eps), q * rho / eps) - w_orbit >= r && (psi_c = max(psi_c, psi)) + w = orbit_widths(psi, T_spline, q_spline, avg_r, avg_R, mass, chrg, ro, bo) + w.r <= 0 && continue + max(w.w_potato, w.rho_banana, w.rho_theta) >= w.r && (psi_c = max(psi_c, psi)) end return psi_c end @@ -203,14 +219,13 @@ function kinetic_validity_profiles(kinetic_profiles::Equilibrium.KineticProfileS is_valid = falses(n) for (i, x) in pairs(psi) r = avg_r(x) - eps = max(r / avg_R(x), 1e-6) - rho = mass * sqrt(2 * T_spline(x) / mass) / (abs(chrg) * bo) + w = orbit_widths(x, T_spline, q_spline, avg_r, avg_R, mass, chrg, ro, bo) + rho_i[i] = w.rho + rho_banana[i] = w.rho_banana + rho_theta[i] = w.rho_theta + w_potato[i] = w.w_potato q = abs(q_spline(x)) drdpsi = abs(r_deriv(x)) - rho_i[i] = rho - rho_banana[i] = q * rho / sqrt(eps) - rho_theta[i] = q * rho / eps - w_potato[i] = cbrt(q^2 * rho^2 * ro) r_minor[i] = r d_separatrix[i] = max(r_sep - r, 0.0) dP = abs(P_deriv(x)) @@ -227,6 +242,26 @@ function kinetic_validity_profiles(kinetic_profiles::Equilibrium.KineticProfileS psi_c=psi_c, envelope=envelope, is_valid=is_valid) end +""" + clear_rational_windows(psi_c, rationals) → Float64 + +Move the near-axis validity boundary outward until no rational surface lies inside the envelope's +transition band `[ψ_c, 2ψ_c]`. A rational sitting in the band gets its (near-singular) kinetic +increments multiplied by a rapidly varying, near-zero envelope, which the matrix splines cannot +represent — the resulting overshoot propagates NaNs into the stability solve. Where the orbit width +already reaches ⟨r⟩ the resonance is not trustworthy anyway, so the boundary is pushed past the +surface (suppressing it wholly) rather than cutting through it. +""" +function clear_rational_windows(psi_c::Float64, rationals::Vector{Float64})::Float64 + psi_c <= 0 && return psi_c + for _ in 1:length(rationals) + inside = filter(r -> psi_c - Equilibrium.RATIONAL_RES_RADIUS <= r <= 2 * psi_c, rationals) + isempty(inside) && break + psi_c = maximum(inside) + Equilibrium.RATIONAL_RES_RADIUS + end + return psi_c +end + """ axis_validity_boundary(kf_ctrl, species, kinetic_profiles, equil) → Float64 @@ -234,11 +269,14 @@ end disabled or no kinetic profiles were loaded. Falls back to `kf_ctrl.zi`/`mi` when the resolved species set is unavailable. """ -function axis_validity_boundary(kf_ctrl::KineticForcesControl, species, kinetic_profiles, equil)::Float64 +function axis_validity_boundary(kf_ctrl::KineticForcesControl, species, kinetic_profiles, equil, + rationals::Vector{Float64}=Float64[])::Float64 (kf_ctrl.axis_validity_suppression && kinetic_profiles !== nothing) || return 0.0 - species === nothing && return kinetic_axis_validity_psi(kinetic_profiles, equil; - zi=kf_ctrl.zi, mi=kf_ctrl.mi, electron=kf_ctrl.electron) - return kinetic_axis_validity_psi(species, equil) + psi_c = + species === nothing ? + kinetic_axis_validity_psi(kinetic_profiles, equil; zi=kf_ctrl.zi, mi=kf_ctrl.mi, electron=kf_ctrl.electron) : + kinetic_axis_validity_psi(species, equil) + return clear_rational_windows(psi_c, rationals) end """ @@ -260,7 +298,26 @@ function resonance_grid_nodes(ctrl, kf_ctrl::KineticForcesControl, kinetic_profi n=n_res, nl=kf_ctrl.nl, zi=kf_ctrl.zi, mi=kf_ctrl.mi, electron=kf_ctrl.electron, wdfac=kf_ctrl.wdfac)) end - psi_c = axis_validity_boundary(kf_ctrl, species, kinetic_profiles, equil) + psi_c = axis_validity_boundary(kf_ctrl, species, kinetic_profiles, equil, + Float64[sng.psifac for sng in intr.sing]) psi_c > 0 && filter!(p -> p > psi_c, nodes) return nodes end + +""" + kinetic_validity_profiles(species, equil) → NamedTuple + +Validity diagnostics for a resolved multi-species set, computed for the species that sets the +boundary (largest ψ_c). Reporting one species' profiles alongside another's ψ_c would make +`is_valid` contradict `envelope` in the same output group. +""" +function kinetic_validity_profiles(species::AbstractVector{<:Equilibrium.ResolvedNTVSpecies}, equil) + isempty(species) && error("kinetic_validity_profiles: empty species set") + widest = species[1] + psi_c = -1.0 + for sp in species + c = kinetic_axis_validity_psi(sp.profiles, equil; zi=sp.z, mi=sp.m, electron=sp.electron) + c > psi_c && (psi_c = c; widest = sp) + end + return kinetic_validity_profiles(widest.profiles, equil; zi=widest.z, mi=widest.m, electron=widest.electron) +end From 4dd32b280a86659d66d2ef871cd0131bb849b79c Mon Sep 17 00:00:00 2001 From: logan-nc Date: Fri, 4 Sep 2026 15:25:23 -0400 Subject: [PATCH 15/23] FFS/KF/PE - MINOR - Clean-code review follow-ups Single named constants for the kinetic singular threshold and its relaxed reporting fraction (they were duplicated literals in two files that could drift), one local-maxima pass instead of two, the regularization policy moved into PerturbedEquilibrium which owns reg_spot, the band-knot count named and justified, the struct field documented in the docstring per convention, and two allocating minimum(abs.(...)) checks replaced with any(). Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01LzbLFQKyuRE5DYZmLokKmk --- src/ForceFreeStates/Fourfit.jl | 1 + src/ForceFreeStates/Kinetic.jl | 15 +++++++---- src/ForceFreeStates/Surfaces/Finding.jl | 26 ++++++++++++------- src/GeneralizedPerturbedEquilibrium.jl | 22 +--------------- src/KineticForces/KineticForcesStructs.jl | 15 ++++++++--- .../PerturbedEquilibriumStructs.jl | 21 +++++++++++++++ 6 files changed, 60 insertions(+), 40 deletions(-) diff --git a/src/ForceFreeStates/Fourfit.jl b/src/ForceFreeStates/Fourfit.jl index 0f8095e3c..57cc28128 100644 --- a/src/ForceFreeStates/Fourfit.jl +++ b/src/ForceFreeStates/Fourfit.jl @@ -443,6 +443,7 @@ rational's resolution window, preserving the Δ′-stencil structure the equilib (`Equilibrium.RATIONAL_RES_RADIUS`). """ function core_capped_knots(xs::Vector{Float64}, rationals::Vector{Float64})::Vector{Int} + length(xs) < 3 && return collect(eachindex(xs)) cap_edge = 0.1 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) diff --git a/src/ForceFreeStates/Kinetic.jl b/src/ForceFreeStates/Kinetic.jl index 9cbb8001f..e057276e1 100644 --- a/src/ForceFreeStates/Kinetic.jl +++ b/src/ForceFreeStates/Kinetic.jl @@ -1,4 +1,9 @@ +# Knots laid across the envelope's [ψ_c, 2ψ_c] transition. The envelope is a quintic smoothstep, +# so a cubic spline needs several interior knots to follow it without overshoot; nine keeps the +# residual below the kernel's own tolerance even on coarse decks. +const BAND_KNOTS = 9 + """ refine_grid_at_fbar_peaks(xs, kw, kt, evaluate, kin, equil, intr, psi_c; ngrid=1000, relaxed_frac=0.01, target=3, max_add=24) → (xs, kw, kt) @@ -14,8 +19,8 @@ and the existing values are reused. A well-resolved grid inserts nothing. function refine_grid_at_fbar_peaks(xs::Vector{Float64}, kw::Array{ComplexF64,3}, kt::Array{ComplexF64,3}, evaluate::Function, kin::KineticMatrices, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal, psi_c::Float64; - ngrid::Int=1000, relaxed_frac::Float64=0.01, target::Int=3, max_add::Int=24, - cond_threshold::Float64=1e8) + ngrid::Int=1000, relaxed_frac::Float64=KINETIC_RELAXED_FRAC, target::Int=3, max_add::Int=24, + cond_threshold::Float64=KINETIC_SINGULAR_COND) lo, hi = xs[1], xs[end] scan = collect(range(lo, hi; length=ngrid)) @@ -45,8 +50,8 @@ function refine_grid_at_fbar_peaks(xs::Vector{Float64}, kw::Array{ComplexF64,3}, w = (scan[r] - scan[l]) / 3 for x in (scan[i], scan[i] - w, scan[i] + w) (lo < x < hi && x > 2 * psi_c) || continue - minimum(abs.(xs .- x)) < Equilibrium.MIN_KNOT_SPACING && continue - isempty(add) || minimum(abs.(add .- x)) >= Equilibrium.MIN_KNOT_SPACING || continue + any(y -> abs(y - x) < Equilibrium.MIN_KNOT_SPACING, xs) && continue + any(y -> abs(y - x) < Equilibrium.MIN_KNOT_SPACING, add) && continue push!(add, x) end end @@ -104,7 +109,7 @@ function build_kinetic_matrix_splines( # spline overshoot can land on a rational surface. Pin the band ends (the smoothstep is # only C² there) and resolve the transition with a fixed set of knots. band_knots(lo, hi) = axis_validity_psi_c > 0 ? - [x for x in range(axis_validity_psi_c, 2 * axis_validity_psi_c; length=9) if lo < x < hi] : Float64[] + [x for x in range(axis_validity_psi_c, 2 * axis_validity_psi_c; length=BAND_KNOTS) if lo < x < hi] : Float64[] # Get raw kinetic matrices (scaling is baked into each source) if ctrl.kinetic_source == "fixed" diff --git a/src/ForceFreeStates/Surfaces/Finding.jl b/src/ForceFreeStates/Surfaces/Finding.jl index 7c17bd5dc..9cc7b6112 100644 --- a/src/ForceFreeStates/Surfaces/Finding.jl +++ b/src/ForceFreeStates/Surfaces/Finding.jl @@ -241,6 +241,11 @@ function evaluate_fbar_condition(psi::Float64, kin::KineticMatrices, equil::Equi return cond(fbar) end +# Kinetic F̄ counts as singular above this condition number; structure within +# KINETIC_RELAXED_FRAC of it is reported as near-singular (the shifted/split resonances). +const KINETIC_SINGULAR_COND = 1.0e8 +const KINETIC_RELAXED_FRAC = 0.01 + """ find_kinetic_singular_surfaces!(mats, equil, intr; ngrid=2000, cond_threshold=1e8) @@ -259,7 +264,13 @@ Algorithm: 3. Refine each peak with golden-section minimization of -cond 4. Filter by threshold and resonance condition """ -function find_kinetic_singular_surfaces!(mats::MatrixSplines, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal; ngrid::Int=2000, cond_threshold::Float64=1e8) +function find_kinetic_singular_surfaces!( + mats::MatrixSplines, + equil::Equilibrium.PlasmaEquilibrium, + intr::ForceFreeStatesInternal; + ngrid::Int=2000, + cond_threshold::Float64=KINETIC_SINGULAR_COND +) kin = mats.kinetic kin === nothing && error("find_kinetic_singular_surfaces! requires a kinetic fit; call build_kinetic_matrix_splines first") psilow = equil.profiles.xs[1] @@ -284,21 +295,16 @@ function find_kinetic_singular_surfaces!(mats::MatrixSplines, equil::Equilibrium intr.kinsing_scan_threshold = cond_threshold # Find local maxima of cond(F̄): points where cond increases then decreases - peak_indices = Int[] - for i in 2:(ngrid-1) - if cond_vals[i] > cond_vals[i-1] && cond_vals[i] > cond_vals[i+1] && cond_vals[i] > cond_threshold - push!(peak_indices, i) - end - end + local_maxima = [i for i in 2:(ngrid-1) if cond_vals[i] > cond_vals[i-1] && cond_vals[i] > cond_vals[i+1]] + peak_indices = filter(i -> cond_vals[i] > cond_threshold, local_maxima) # Peaks below the threshold are not singular surfaces, but they mark where the kinetic F̄ comes # closest to singular — the shifted/split resonances of Park & Logan Eq. (70). Report the # strongest few so sharp kinetic structure is visible rather than silent (on a DIII-D-like case # these track the NTV torque-density peaks at low collisionality/rotation). - subthreshold = [i for i in 2:(ngrid-1) if cond_vals[i] > cond_vals[i-1] && cond_vals[i] > cond_vals[i+1] && - cond_threshold / 100 < cond_vals[i] <= cond_threshold] + subthreshold = filter(i -> KINETIC_RELAXED_FRAC * cond_threshold < cond_vals[i] <= cond_threshold, local_maxima) if !isempty(subthreshold) - top = sort(subthreshold; by=i -> -cond_vals[i])[1:min(3, length(subthreshold))] + top = sort(subthreshold; by=i -> cond_vals[i], rev=true)[1:min(3, length(subthreshold))] @info "Kinetic F̄ near-singular structure below the singular threshold at " * join(["ψ=$(round(psi_grid[i]; digits=4)) (cond=$(round(cond_vals[i]; sigdigits=3)))" for i in top], ", ") * " — full scan in SingularSurfaces/Kinetic/scan_cond; check the ψ grid resolves these if results look grid-sensitive" diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 574fa4d9f..aea43f817 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -855,26 +855,6 @@ function run_perturbed_equilibrium( return pe_state end -""" - kinetic_regularization_kwargs(ffs, kwargs) -> NamedTuple - -`reg_spot` smooths the **ideal** 1/(m−nq) divergence of ξ^ψ′ and ξ^α before they drive the NTV -integrand. A self-consistent kinetic solve has no such divergence — det(F̄) is complex and nonzero -at the rationals (Park & Logan, Phys. Plasmas 24, 032505 (2017) §III D) — so regularizing there -suppresses a finite physical response, and does so inconsistently, since ξ^ψ is never regularized. -Force `reg_spot = 0` for kinetic solves and log the override; ideal solves keep their setting. -""" -function kinetic_regularization_kwargs(ffs::ForceFreeStatesResult, kwargs) - base = values(kwargs) - default_reg = PerturbedEquilibrium.PerturbedEquilibriumControl().reg_spot - prev = get(base, :reg_spot, default_reg) - kinetic = ffs.mats.kinetic !== nothing - kinetic && prev != 0 && - @info "Self-consistent kinetic run: overriding reg_spot=$prev with 0 " * - "(the kinetic terms remove the ideal resonant singularity; see docs/src/kinetic_forces.md)" - return merge(base, (; reg_spot=kinetic ? 0.0 : prev)) -end - """ perturbed_equilibrium(ffs, rmp; forcing_modes=nothing, coil_sets=nothing, kwargs...) -> PerturbedEquilibriumState @@ -901,7 +881,7 @@ function perturbed_equilibrium( kwargs... ) ctrl = ffs.control - pe_ctrl = PerturbedEquilibrium.PerturbedEquilibriumControl(; kinetic_regularization_kwargs(ffs, kwargs)...) + pe_ctrl = PerturbedEquilibrium.PerturbedEquilibriumControl(; PerturbedEquilibrium.kinetic_regularization_kwargs(ffs, kwargs)...) pe_intr = PerturbedEquilibrium.PerturbedEquilibriumInternal(; dir_path=ffs.dir_path) # Inner-layer penetrated resonant field; zeros under ideal closure. diff --git a/src/KineticForces/KineticForcesStructs.jl b/src/KineticForces/KineticForcesStructs.jl index 3966d7410..a1a7f0a90 100644 --- a/src/KineticForces/KineticForcesStructs.jl +++ b/src/KineticForces/KineticForcesStructs.jl @@ -86,6 +86,13 @@ ctrl = KineticForcesControl(; (Symbol(k) => v for (k, v) in inputs["KineticForce Immutable: vary a field by building a new control rather than assigning to one (the multi-species loop does this per species, and `check_psi_quadrature_convergence`'s test builds a second control for its differing tolerance). + + - `axis_validity_suppression::Bool` - Suppress the calculated kinetic terms where the + zero-orbit-width ordering fails near the axis: below ψ_c — the outermost ψ at which a thermal + orbit width of the widest-orbit species reaches ⟨r⟩, moved clear of any rational surface the + transition would otherwise cut through — the increments are zeroed and the kernel is skipped, + rising to full strength at 2ψ_c through a C² envelope. Boundary and envelope come from the + equilibrium and kinetic profiles; there are no tuning parameters. Default `true`. """ @kwdef struct KineticForcesControl # Moment type @@ -151,7 +158,7 @@ builds a second control for its differing tolerance). nufac::Float64 = 1.0 # Collisionality scaling divxfac::Float64 = 1.0 # div(xi_perp) scaling - axis_validity_suppression::Bool = true # Suppress kinetic terms where the zero-orbit-width ordering fails near the axis (see kinetic_axis_validity_psi); envelope and boundary are profile-derived, no tuning parameters + axis_validity_suppression::Bool = true # documented in the `## Fields` docstring # Energy integration parameters nutype::String = "harmonic" # Collision operator: "zero", "small", "krook", "harmonic" @@ -334,9 +341,9 @@ The JBB deweighting algorithm (Fortran pentrc/inputs.f90:828-868): 4. Forward DFT back to m-space """ function set_perturbation_data!(kf_intr::KineticForcesInternal, pe_state::PerturbedEquilibrium.PerturbedEquilibriumState, - ffs::ForceFreeStates.ForceFreeStatesResult, - equil::Equilibrium.PlasmaEquilibrium, - metric::ForceFreeStates.MetricData) + ffs::ForceFreeStates.ForceFreeStatesResult, + equil::Equilibrium.PlasmaEquilibrium, + metric::ForceFreeStates.MetricData) # Copy mode numbers from FFS kf_intr.mlow = ffs.mlow kf_intr.mhigh = ffs.mhigh diff --git a/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl b/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl index ee0d23a85..52a97be33 100644 --- a/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl +++ b/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl @@ -26,6 +26,7 @@ Medium Priority (defer for MWE): - `singular_point_method::String` - Method for singular point treatment (default: "standard") Regularization: + # High Priority (MWE) - `reg_spot::Float64` - Regularization width for singular surface smoothing (default: 0.05). Set to 0 to disable. Must be ≥ 0. Forced to 0 in self-consistent kinetic runs, whose Euler-Lagrange operator has no resonant singularity to smooth. @@ -210,3 +211,23 @@ well-conditioned flux-space inductances L, Λ: plasma_energy::Float64 = 0.0 toroidal_torque::Float64 = 0.0 end + +""" + kinetic_regularization_kwargs(ffs, kwargs) -> NamedTuple + +`reg_spot` smooths the **ideal** 1/(m−nq) divergence of ξ^ψ′ and ξ^α before they drive the NTV +integrand. A self-consistent kinetic solve has no such divergence — det(F̄) is complex and nonzero +at the rationals (Park & Logan, Phys. Plasmas 24, 032505 (2017) §III D) — so regularizing there +suppresses a finite physical response, and does so inconsistently, since ξ^ψ is never regularized. +Force `reg_spot = 0` for kinetic solves and log the override; ideal solves keep their setting. +""" +function kinetic_regularization_kwargs(ffs::ForceFreeStatesResult, kwargs) + base = values(kwargs) + default_reg = PerturbedEquilibrium.PerturbedEquilibriumControl().reg_spot + prev = get(base, :reg_spot, default_reg) + kinetic = ffs.mats.kinetic !== nothing + kinetic && prev != 0 && + @info "Self-consistent kinetic run: overriding reg_spot=$prev with 0 " * + "(the kinetic terms remove the ideal resonant singularity; see docs/src/kinetic_forces.md)" + return merge(base, (; reg_spot=kinetic ? 0.0 : prev)) +end From 2dcbbf17079547eea6ed1cc1ebb955dd00d894f1 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Sat, 19 Sep 2026 09:11:29 -0400 Subject: [PATCH 16/23] KineticForces - BUGFIX! - Use one validity boundary and bound its rational clearing Independent review of #414 found the near-axis suppression boundary computed three different ways, and two ways it can run away. One boundary per run. The NTV psi torque quadrature recomputed psi_c from the single species it was integrating and never called clear_rational_windows, while the EL kinetic matrices and the Validity output used the widest-orbit species' cleared boundary. Validity/psi_c in gpec.h5 therefore described a different domain from the one the torque was integrated over, and the standalone torque kept near-axis contributions the EL build had already zeroed - which is exactly the cross-check the PE-vs-NTV agreement rests on. axis_validity_boundary is now the only entry point: run_kinetic_forces computes it once and threads it through compute_torque_all_methods! into integrate_psi_quadgk. Contiguity. kinetic_axis_validity_psi took the outermost trigger anywhere on the profile grid. q is not monotonic, so an edge q rise or a pedestal gradient could re-trigger the criterion far from the axis and drag the suppressed band across a large fraction of the plasma. The scan now stops at the first surface where the ordering holds and warns if the criterion fires again further out. Cascade. clear_rational_windows rescans [psi_c, 2 psi_c] after each move, so the window widens as psi_c grows; with rationals at 0.15/0.25/0.45 a psi_c of 0.10 walked out to 0.455. The move is capped at VALIDITY_CLEAR_MAX_FACTOR times the orbit-width boundary, past which the uncleared boundary is kept with a warning - a band cutting a rational is resolved by the band knots, whereas suppressing half the plasma is not recoverable. Grid pinning. resonance_grid_nodes cleared its boundary against intr.sing, which sing_find! does not populate until long after the grid is built, so the clearing was silently a no-op. The nodes are now located inside maybe_reform_equilibrium against that pass's own rational list. Validity coverage. write_validity! was reachable only from run_kinetic_forces, whose body is gated on a perturbed-equilibrium state, so a calculated-kinetic FFS-only run wrote no Validity group at all. It is written from the ForceFreeStates output stage instead, wherever kinetic profiles were loaded. kinetic_regularization_kwargs is duck-typed on ffs so the override is testable, and its message no longer says "self-consistent" - it covers both kinetic sources. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01VdRsGGg4YRG6E9fgLGbpdF --- src/GeneralizedPerturbedEquilibrium.jl | 35 ++++++++---- .../CalculatedKineticMatrices.jl | 5 +- src/KineticForces/Compute.jl | 18 ++++--- src/KineticForces/Utils.jl | 54 ++++++++++++++----- .../PerturbedEquilibriumStructs.jl | 20 +++---- 5 files changed, 90 insertions(+), 42 deletions(-) diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index aea43f817..df8a0678e 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -211,7 +211,7 @@ function main_from_inputs( kf_ctrl, kinetic_profiles, kf_species = load_kinetic_context(inputs, intr, ctrl, equil) equil = maybe_reform_equilibrium(equil, eq_config, additional_input, intr, ctrl, kinetic_profiles; - pinned=KineticForces.resonance_grid_nodes(ctrl, kf_ctrl, kinetic_profiles, kf_species, equil, intr)) + kf_ctrl=kf_ctrl, kf_species=kf_species) @info "Equilibrium construction completed in $(@sprintf("%.3f", time() - equil_start)) s" @@ -257,6 +257,15 @@ function main_from_inputs( @info "Results written to $(ctrl.HDF5_filename)" end + # Validity describes the kinetic model, not the torque contraction, so it is written wherever + # kinetic profiles were loaded — including FFS-only runs that never reach run_kinetic_forces. + if ctrl.write_outputs_to_HDF5 && kinetic_profiles !== nothing + h5open(joinpath(ffs_result.dir_path, ctrl.HDF5_filename), "cw") do h5file + KineticForces.write_validity!(h5file, kf_ctrl, kf_species, kinetic_profiles, equil, + Float64[sng.psifac for sng in intr.sing]) + end + end + @info "Force-Free States completed in $(@sprintf("%.3f", time() - ffs_start)) s" # Early exit if user only requested force-free states (SLAYER still runs). @@ -380,7 +389,8 @@ function load_kinetic_context( end """ - maybe_reform_equilibrium(equil, eq_config, additional_input, intr, ctrl, kinetic_profiles) -> equil + maybe_reform_equilibrium(equil, eq_config, additional_input, intr, ctrl, kinetic_profiles; + kf_ctrl=nothing, kf_species=nothing) -> equil Two-pass auto grid: measure the pass-1 equilibrium's curvature (profiles, geometry, kinetic profiles), pin knots on rational surfaces, and re-form on the refined grid from the in-memory @@ -393,11 +403,16 @@ function maybe_reform_equilibrium( intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl, kinetic_profiles; - pinned::Vector{Float64}=Float64[] + kf_ctrl::Union{Nothing,KineticForces.KineticForcesControl}=nothing, + kf_species=nothing ) Equilibrium.wants_two_pass(eq_config) || return equil mandatory = ForceFreeStates.rational_psi_nodes(equil; nlow=intr.nlow, nhigh=intr.nhigh) + # Kinetic-resonance nodes are located here rather than at the call site so the validity + # boundary that filters them is cleared against this pass's rationals (intr.sing is empty yet). + pinned = kf_ctrl === nothing ? Float64[] : + KineticForces.resonance_grid_nodes(ctrl, kf_ctrl, kinetic_profiles, kf_species, equil, intr; rationals=mandatory) # Smallest |n| in the run sets the widest matching half-stencil dpsi = singfac_min/(n_min·|q′|), # so the rational-surface brackets clear a zone large enough for every mode. n_min = minimum(abs(n) for n in intr.nlow:intr.nhigh if n != 0) @@ -942,6 +957,10 @@ function run_kinetic_forces( @info "\n KineticForces\n$_SECTION" kf_start = time() + # One boundary for the whole run, shared with the EL kinetic matrices and the Validity output. + axis_psi_c = KineticForces.axis_validity_boundary(kf_ctrl, species, kinetic_profiles, result.equil, + Float64[sng.psifac for sng in result.surfaces]) + # Standalone NTV torque diagnostics need a PE state (they contract kinetic operators # against ξ). The self-consistent kinetic_source="calculated" path produces none — skip. if pe_state === nothing @@ -953,13 +972,12 @@ function run_kinetic_forces( if species === nothing kf_state = KineticForces.KineticForcesState() - KineticForces.compute_torque_all_methods!(kf_state, kf_intr, kf_ctrl, result.equil, kinetic_profiles) + KineticForces.compute_torque_all_methods!(kf_state, kf_intr, kf_ctrl, result.equil, kinetic_profiles; + axis_psi_c=axis_psi_c) if kf_ctrl.write_outputs_to_HDF5 h5open(joinpath(result.dir_path, kf_ctrl.HDF5_filename), "cw") do h5file KineticForces.write_to_hdf5!(h5file, kf_state; dVdpsi_spline=result.equil.profiles.dVdpsi_spline) - KineticForces.write_validity!(h5file, kf_ctrl, species, kinetic_profiles, result.equil, - Float64[sng.psifac for sng in result.surfaces]) end end else @@ -972,7 +990,8 @@ function run_kinetic_forces( (f => getfield(kf_ctrl, f) for f in fieldnames(KineticForces.KineticForcesControl))..., zi=sp.z, mi=sp.m, electron=sp.electron, ion_species=KineticForces.IonSpecies[]) st = KineticForces.KineticForcesState() - KineticForces.compute_torque_all_methods!(st, kf_intr, sctrl, result.equil, sp.profiles) + KineticForces.compute_torque_all_methods!(st, kf_intr, sctrl, result.equil, sp.profiles; + axis_psi_c=axis_psi_c) push!(states, st) push!(labels, sp.label) end @@ -985,8 +1004,6 @@ function run_kinetic_forces( end KineticForces.write_to_hdf5!(h5file, kf_state; dVdpsi_spline=result.equil.profiles.dVdpsi_spline) - KineticForces.write_validity!(h5file, kf_ctrl, species, kinetic_profiles, result.equil, - Float64[sng.psifac for sng in result.surfaces]) end end end diff --git a/src/KineticForces/CalculatedKineticMatrices.jl b/src/KineticForces/CalculatedKineticMatrices.jl index de916b0ac..d6436b710 100644 --- a/src/KineticForces/CalculatedKineticMatrices.jl +++ b/src/KineticForces/CalculatedKineticMatrices.jl @@ -116,8 +116,9 @@ function compute_calculated_kinetic_matrices( [Equilibrium.ResolvedNTVSpecies(kf_ctrl.zi, kf_ctrl.mi, kf_ctrl.electron, "single", kinetic_profiles)] : species - # Near-axis validity envelope: suppress the drift-kinetic increments where the zero-orbit-width - # ordering fails, taking the widest-orbit species. Kernel evaluation is skipped where it is 0. + # Near-axis validity envelope: suppresses the drift-kinetic increments where the zero-orbit-width + # ordering fails (kernel evaluation skipped where it is 0). The same real scalar must scale kw and + # kt at the same ψ — that is what keeps the blocks Hermitian and the energy/torque relation intact. env = ones(Float64, mpsi) if kf_ctrl.axis_validity_suppression psi_c = axis_psi_c diff --git a/src/KineticForces/Compute.jl b/src/KineticForces/Compute.jl index 1ad347561..3ed4fdd9d 100644 --- a/src/KineticForces/Compute.jl +++ b/src/KineticForces/Compute.jl @@ -83,7 +83,7 @@ function integrate_psi_quadgk( wdfac::Float64, divxfac::Float64, electron::Bool, method::String, equil, intr::KineticForcesInternal, ctrl::KineticForcesControl, kinetic_profiles::Equilibrium.KineticProfileSplines; - psi_min::Float64=0.0, psi_max::Float64=1.0 + psi_min::Float64=0.0, psi_max::Float64=1.0, axis_psi_c::Float64=0.0 ) is_matrix_method = occursin("mm", method) mpert = intr.mpert @@ -100,11 +100,10 @@ function integrate_psi_quadgk( panel_psis=Float64[], resonance_psis=Float64[]) end - # Near-axis validity suppression: same boundary and envelope as the EL kinetic - # matrices (one source of truth), applied to the torque density; the quadrature - # domain starts at the boundary since the integrand is identically zero below it. - psi_c = ctrl.axis_validity_suppression ? - kinetic_axis_validity_psi(kinetic_profiles, equil; zi=zi, mi=mi, electron=electron) : 0.0 + # Near-axis validity suppression on the torque density. `axis_psi_c` is the run's single + # boundary (widest-orbit species, cleared of rationals), so the quadrature domain matches the + # EL kinetic matrices and the Validity output; below it the integrand is identically zero. + psi_c = ctrl.axis_validity_suppression ? axis_psi_c : 0.0 x0 = max(x0, psi_c) if x0 >= xout return (total=ComplexF64(0.0), torque_profile=nothing, matrix_integrated=nothing, psi_nsteps=0, psi_quad_error=0.0, @@ -305,10 +304,13 @@ block-diagonal kinetic matrices. - `ctrl::KineticForcesControl`: Control parameters specifying which methods to run - `equil`: PlasmaEquilibrium with 2D interpolants - `kinetic_profiles::Equilibrium.KineticProfileSplines`: Named kinetic-profile splines + - `axis_psi_c`: the run's near-axis validity boundary, suppressing the torque density below it + (`KineticForces.axis_validity_boundary`); 0 disables the suppression """ function compute_torque_all_methods!(state::KineticForcesState, intr::KineticForcesInternal, ctrl::KineticForcesControl, equil, - kinetic_profiles::Equilibrium.KineticProfileSplines) + kinetic_profiles::Equilibrium.KineticProfileSplines; + axis_psi_c::Float64=0.0) for entry in METHOD_REGISTRY getfield(ctrl, entry.flag) || continue @@ -351,7 +353,7 @@ function compute_torque_all_methods!(state::KineticForcesState, intr::KineticFor n, ctrl.nl, ctrl.zi, ctrl.mi, ctrl.wdfac, ctrl.divxfac, ctrl.electron, method, equil, intr, ctrl, kinetic_profiles; - psi_min=ctrl.psilims[1], psi_max=ctrl.psilims[2]) + psi_min=ctrl.psilims[1], psi_max=ctrl.psilims[2], axis_psi_c=axis_psi_c) total_torque += result.total psi_nsteps_total += result.psi_nsteps diff --git a/src/KineticForces/Utils.jl b/src/KineticForces/Utils.jl index 80a04f1bf..d2281de21 100644 --- a/src/KineticForces/Utils.jl +++ b/src/KineticForces/Utils.jl @@ -117,10 +117,9 @@ end kinetic-profile grid where any of the three thermal orbit-width scales reaches the local minor radius ⟨r⟩ — potato width `(q²ρ²R₀)^(1/3)`, banana width `q·ρ/√ε`, and poloidal gyroradius `q·ρ/ε` (ε = ⟨r⟩/⟨R⟩, clamped as in `kinetic_resonance_psi_nodes`; thermal gyroradius -ρ = √(2·m·T)·/(Z·e·B₀) of the computed species). Inside this boundary trapped bananas become -potato orbits with width comparable to r itself, so the bounce-averaged kinetic response is -evaluated outside its validity domain (measured consequence: diverging kinetic increments and a -pathological EL step count). Returns 0.0 when no criterion is met anywhere. The Fortran precedent +ρ = √(2·m·T)·/(Z·e·B₀) of the computed species). Inside this boundary trapped bananas become potato +orbits of width comparable to r itself, so the bounce-averaged kinetic response is evaluated outside +its validity domain. Returns 0.0 when no criterion is met anywhere. The Fortran precedent (`ktanh_flag`, dcon/fourfit.F) suppressed the same region with four hand-tuned knobs; here the boundary is derived from the profiles with no user parameters. """ @@ -134,13 +133,23 @@ function kinetic_axis_validity_psi(kinetic_profiles::Equilibrium.KineticProfileS avg_R = equil.geometry.avg_R_spline ro = abs(equil.ro) bo = abs(equil.params.b0) + # Only the region contiguous with the axis counts: q is not monotonic, so an edge q rise or a + # pedestal gradient can re-trigger the criterion far out, where the model has not lost validity. psi_c = 0.0 + psi_outer = 0.0 + contiguous = true for psi in kinetic_profiles.xs psi <= 0 && continue w = orbit_widths(psi, T_spline, q_spline, avg_r, avg_R, mass, chrg, ro, bo) w.r <= 0 && continue - max(w.w_potato, w.rho_banana, w.rho_theta) >= w.r && (psi_c = max(psi_c, psi)) + if max(w.w_potato, w.rho_banana, w.rho_theta) >= w.r + contiguous ? (psi_c = psi) : (psi_outer = max(psi_outer, psi)) + else + contiguous = false + end end + psi_outer > 0 && @warn "kinetic_axis_validity_psi: orbit widths also reach ⟨r⟩ out to ψ=$(round(psi_outer; sigdigits=3)), " * + "outside the axis-contiguous region ending at ψ_c=$(round(psi_c; sigdigits=3)); only the latter is suppressed" return psi_c end @@ -242,22 +251,39 @@ function kinetic_validity_profiles(kinetic_profiles::Equilibrium.KineticProfileS psi_c=psi_c, envelope=envelope, is_valid=is_valid) end +# Bound on how far clearing rationals out of the envelope band may push ψ_c, as a multiple of the +# orbit-width boundary; beyond it the uncleared boundary is kept (see `clear_rational_windows`). +const VALIDITY_CLEAR_MAX_FACTOR = 1.5 + """ clear_rational_windows(psi_c, rationals) → Float64 Move the near-axis validity boundary outward until no rational surface lies inside the envelope's -transition band `[ψ_c, 2ψ_c]`. A rational sitting in the band gets its (near-singular) kinetic -increments multiplied by a rapidly varying, near-zero envelope, which the matrix splines cannot -represent — the resulting overshoot propagates NaNs into the stability solve. Where the orbit width -already reaches ⟨r⟩ the resonance is not trustworthy anyway, so the boundary is pushed past the -surface (suppressing it wholly) rather than cutting through it. +transition band, tested over `[ψ_c - RATIONAL_RES_RADIUS, 2ψ_c]` so a surface just inside the band +edge is cleared too. A rational sitting in the band gets its (near-singular) kinetic increments +multiplied by a rapidly varying, near-zero envelope, which the matrix splines cannot represent — +the resulting overshoot propagates NaNs into the stability solve. Where the orbit width already +reaches ⟨r⟩ the resonance is not trustworthy anyway, so the boundary is pushed past the surface +(suppressing it wholly) rather than cutting through it. + +The band widens as ψ_c moves, so clearing can cascade through a rational-dense core. The move is +capped at `VALIDITY_CLEAR_MAX_FACTOR`× the orbit-width boundary; past that the uncleared boundary +is kept with a warning, since suppressing a large fraction of the plasma is worse than a band that +cuts a rational (which the band knots resolve). """ function clear_rational_windows(psi_c::Float64, rationals::Vector{Float64})::Float64 psi_c <= 0 && return psi_c + psi_c0 = psi_c + limit = VALIDITY_CLEAR_MAX_FACTOR * psi_c0 for _ in 1:length(rationals) inside = filter(r -> psi_c - Equilibrium.RATIONAL_RES_RADIUS <= r <= 2 * psi_c, rationals) isempty(inside) && break psi_c = maximum(inside) + Equilibrium.RATIONAL_RES_RADIUS + if psi_c > limit + @warn "clear_rational_windows: clearing the envelope band of rationals would push ψ_c from " * + "$(round(psi_c0; sigdigits=3)) past $(round(limit; sigdigits=3)); keeping the orbit-width boundary" + return psi_c0 + end end return psi_c end @@ -286,10 +312,11 @@ end located Ω_ℓ = 0 resonance surfaces for every toroidal mode in the run, outside the near-axis validity region (nodes there are suppressed anyway). Empty for ideal runs — the ideal grid criterion knows nothing about kinetic resonances, so this is the only thing that puts knots on -them. +them. `rationals` is the pass-1 rational-surface list the boundary is cleared against; the caller +supplies it because `intr.sing` is not populated until after the grid is built. """ function resonance_grid_nodes(ctrl, kf_ctrl::KineticForcesControl, kinetic_profiles, species, - equil, intr)::Vector{Float64} + equil, intr; rationals::Vector{Float64}=Float64[])::Vector{Float64} nodes = Float64[] (ctrl.kinetic_factor > 0 && ctrl.kinetic_source == "calculated" && kinetic_profiles !== nothing) || return nodes for n_res in intr.nlow:intr.nhigh @@ -298,8 +325,7 @@ function resonance_grid_nodes(ctrl, kf_ctrl::KineticForcesControl, kinetic_profi n=n_res, nl=kf_ctrl.nl, zi=kf_ctrl.zi, mi=kf_ctrl.mi, electron=kf_ctrl.electron, wdfac=kf_ctrl.wdfac)) end - psi_c = axis_validity_boundary(kf_ctrl, species, kinetic_profiles, equil, - Float64[sng.psifac for sng in intr.sing]) + psi_c = axis_validity_boundary(kf_ctrl, species, kinetic_profiles, equil, rationals) psi_c > 0 && filter!(p -> p > psi_c, nodes) return nodes end diff --git a/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl b/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl index 52a97be33..2a53285e3 100644 --- a/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl +++ b/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl @@ -27,9 +27,9 @@ Medium Priority (defer for MWE): Regularization: -# High Priority (MWE) - - - `reg_spot::Float64` - Regularization width for singular surface smoothing (default: 0.05). Set to 0 to disable. Must be ≥ 0. Forced to 0 in self-consistent kinetic runs, whose Euler-Lagrange operator has no resonant singularity to smooth. + - `reg_spot::Float64` - Regularization width for singular surface smoothing (default: 0.05). Set to 0 + to disable. Must be ≥ 0. Forced to 0 in kinetic runs, whose Euler-Lagrange operator has no + resonant singularity to smooth. """ @kwdef struct PerturbedEquilibriumControl # High Priority (MWE) @@ -216,18 +216,20 @@ end kinetic_regularization_kwargs(ffs, kwargs) -> NamedTuple `reg_spot` smooths the **ideal** 1/(m−nq) divergence of ξ^ψ′ and ξ^α before they drive the NTV -integrand. A self-consistent kinetic solve has no such divergence — det(F̄) is complex and nonzero -at the rationals (Park & Logan, Phys. Plasmas 24, 032505 (2017) §III D) — so regularizing there -suppresses a finite physical response, and does so inconsistently, since ξ^ψ is never regularized. -Force `reg_spot = 0` for kinetic solves and log the override; ideal solves keep their setting. +integrand. A kinetic solve has no such divergence — det(F̄) is complex and nonzero at the rationals +(Park & Logan, Phys. Plasmas 24, 032505 (2017) §III D) — so regularizing there suppresses a finite +physical response, and does so inconsistently, since ξ^ψ is never regularized. Force +`reg_spot = 0` whenever a kinetic fit is present (either `kinetic_source`, both of which remove the +resonant singularity) and log the override; ideal solves keep their setting. `ffs` is duck-typed: +only `ffs.mats.kinetic` is read, so the override is unit-testable without a full run. """ -function kinetic_regularization_kwargs(ffs::ForceFreeStatesResult, kwargs) +function kinetic_regularization_kwargs(ffs, kwargs) base = values(kwargs) default_reg = PerturbedEquilibrium.PerturbedEquilibriumControl().reg_spot prev = get(base, :reg_spot, default_reg) kinetic = ffs.mats.kinetic !== nothing kinetic && prev != 0 && - @info "Self-consistent kinetic run: overriding reg_spot=$prev with 0 " * + @info "Kinetic run: overriding reg_spot=$prev with 0 " * "(the kinetic terms remove the ideal resonant singularity; see docs/src/kinetic_forces.md)" return merge(base, (; reg_spot=kinetic ? 0.0 : prev)) end From 292c50b24bb50133701278ca305f8cdb5a90ceeb Mon Sep 17 00:00:00 2001 From: logan-nc Date: Sat, 19 Sep 2026 09:11:37 -0400 Subject: [PATCH 17/23] KineticForces - TEST - Cover the near-axis validity boundary and regularization override The validity machinery shipped with no unit coverage. Adds runtests_kinetic_validity.jl (registered in runtests.jl) over the pieces the review found defects in: - kinetic_axis_validity_envelope: endpoint values, the midpoint, monotonicity, and a C2 check that the second difference straddling each band end vanishes linearly in h - clear_rational_windows: empty and far-rational no-ops, a single in-band rational, and the cascade construction that must hit the cap and keep the orbit-width boundary - kinetic_axis_validity_psi: a q spike far from the axis must not move the boundary - kinetic_regularization_kwargs: ideal pass-through, kinetic forcing 0, and the case where the deck omits reg_spot entirely and must pick up the struct default The equilibrium is a duck-typed stand-in - the boundary scan reads only q, , , ro and b0 - so the tests run in seconds without building a PlasmaEquilibrium. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01VdRsGGg4YRG6E9fgLGbpdF --- test/runtests.jl | 1 + test/runtests_kinetic_validity.jl | 114 ++++++++++++++++++++++++++++++ 2 files changed, 115 insertions(+) create mode 100644 test/runtests_kinetic_validity.jl diff --git a/test/runtests.jl b/test/runtests.jl index 7845eea60..5fa49ba86 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -48,6 +48,7 @@ else include("./runtests_dispersion_polish.jl") include("./runtests_slayer_runner.jl") include("./runtests_kinetic.jl") + include("./runtests_kinetic_validity.jl") include("./runtests_multiion.jl") include("./runtests_fullruns.jl") include("./runtests_coils.jl") diff --git a/test/runtests_kinetic_validity.jl b/test/runtests_kinetic_validity.jl new file mode 100644 index 000000000..c84b1d217 --- /dev/null +++ b/test/runtests_kinetic_validity.jl @@ -0,0 +1,114 @@ +@testset "Kinetic near-axis validity boundary" begin + + GPE = GeneralizedPerturbedEquilibrium + KF = GPE.KineticForces + EQ = GPE.Equilibrium + PE = GPE.PerturbedEquilibrium + + # Duck-typed stand-in for PlasmaEquilibrium: the boundary scan reads only these fields. + mock_equil(q, r, R; ro=1.7, b0=2.0) = + (; profiles=(; q_spline=q), geometry=(; avg_r_spline=r, avg_R_spline=R), ro=ro, params=(; b0=b0)) + + # Kinetic profiles carrying a prescribed Ti(ψ) in Joules; only xs and Ti_spline are read. + function mock_profiles(xs, Ti_eV) + ones_ = ones(length(xs)) + Ti = Ti_eV .* 1.602176634e-19 + return EQ.KineticProfileSplines(xs, 1e19 .* ones_, 1e19 .* ones_, Ti, Ti, + zero(ones_), 17.0 .* ones_, ones_, ones_, ones_) + end + + @testset "kinetic_axis_validity_envelope" begin + psi_c = 0.04 + @test KF.kinetic_axis_validity_envelope(0.0, psi_c) == 0.0 + @test KF.kinetic_axis_validity_envelope(psi_c, psi_c) == 0.0 + @test KF.kinetic_axis_validity_envelope(2 * psi_c, psi_c) == 1.0 + @test KF.kinetic_axis_validity_envelope(1.0, psi_c) == 1.0 + @test KF.kinetic_axis_validity_envelope(1.5 * psi_c, psi_c) ≈ 0.5 atol = 1e-12 + + # Disabled boundary leaves everything unsuppressed. + @test KF.kinetic_axis_validity_envelope(0.0, 0.0) == 1.0 + @test KF.kinetic_axis_validity_envelope(0.5, -1.0) == 1.0 + + # Monotone across the transition. + band = range(psi_c, 2 * psi_c; length=101) + vals = [KF.kinetic_axis_validity_envelope(x, psi_c) for x in band] + @test all(diff(vals) .>= -1e-15) + + # C² at the band ends: the second difference straddling each end vanishes linearly in h. + # A C¹-but-not-C² join would leave it constant, and a slope kink would blow it up as 1/h. + d2(x, h) = (KF.kinetic_axis_validity_envelope(x + h, psi_c) - 2 * KF.kinetic_axis_validity_envelope(x, psi_c) + + KF.kinetic_axis_validity_envelope(x - h, psi_c)) / h^2 + for edge in (psi_c, 2 * psi_c) + h = psi_c * 1e-2 + @test abs(d2(edge, h) / d2(edge, h / 2)) ≈ 2 atol = 0.2 + end + end + + @testset "clear_rational_windows" begin + R = EQ.RATIONAL_RES_RADIUS + + @test KF.clear_rational_windows(0.04, Float64[]) == 0.04 + @test KF.clear_rational_windows(0.0, [0.05]) == 0.0 + + # A rational well outside the band is left alone. + @test KF.clear_rational_windows(0.04, [0.5]) == 0.04 + + # A rational inside [ψ_c - R, 2ψ_c] pushes the boundary just past it. + @test KF.clear_rational_windows(0.10, [0.1221]) ≈ 0.1221 + R + + # Cascade guard: rationals spaced so each move widens the band into the next one. The move + # is capped at VALIDITY_CLEAR_MAX_FACTOR × ψ_c, and the orbit-width boundary is kept. + cascade = [0.15, 0.25, 0.45] + local psi_capped + @test_logs (:warn,) match_mode = :any begin + psi_capped = KF.clear_rational_windows(0.10, cascade) + end + @test psi_capped == 0.10 + @test KF.VALIDITY_CLEAR_MAX_FACTOR > 1 + end + + @testset "kinetic_axis_validity_psi is contiguous with the axis" begin + xs = collect(0.005:0.005:1.0) + # r = a√ψ, R constant: ε → 0 on axis, so ρ_θ = qρ/ε diverges there and the criterion fires. + a, R0 = 0.6, 1.7 + r_of(psi) = a * sqrt(psi) + R_of(_) = R0 + Ti = fill(2000.0, length(xs)) # 2 keV, flat + + # Smooth q: the invalid region is a single band touching the axis. + q_smooth(psi) = 1.0 + 2.0 * psi^2 + psi_c = KF.kinetic_axis_validity_psi(mock_profiles(xs, Ti), mock_equil(q_smooth, r_of, R_of)) + @test psi_c > 0 + @test psi_c < 0.5 + + # Same profile with a narrow q spike far from the axis. The spike re-triggers the criterion, + # but the model has not lost validity there, so the boundary must not jump out to it. + q_spiked(psi) = q_smooth(psi) + 400.0 * exp(-((psi - 0.9) / 0.01)^2) + local psi_c_spiked + @test_logs (:warn,) match_mode = :any begin + psi_c_spiked = KF.kinetic_axis_validity_psi(mock_profiles(xs, Ti), mock_equil(q_spiked, r_of, R_of)) + end + @test psi_c_spiked ≈ psi_c + @test psi_c_spiked < 0.5 + end + + @testset "kinetic_regularization_kwargs" begin + default_reg = PE.PerturbedEquilibriumControl().reg_spot + ideal = (; mats=(; kinetic=nothing)) + kinetic = (; mats=(; kinetic=(;))) + + # Ideal runs keep whatever the deck asked for, explicit or defaulted. + @test PE.kinetic_regularization_kwargs(ideal, pairs((; reg_spot=0.05))).reg_spot == 0.05 + @test PE.kinetic_regularization_kwargs(ideal, pairs((;))).reg_spot == default_reg + + # Kinetic runs force 0 — including when the deck omits the key entirely, which must pick up + # the struct default rather than leaving reg_spot unset. + @test PE.kinetic_regularization_kwargs(kinetic, pairs((; reg_spot=0.05))).reg_spot == 0.0 + @test PE.kinetic_regularization_kwargs(kinetic, pairs((;))).reg_spot == 0.0 + @test default_reg != 0.0 # otherwise the override above proves nothing + + # Unrelated kwargs pass through untouched. + out = PE.kinetic_regularization_kwargs(kinetic, pairs((; reg_spot=0.05, filter_modes=true))) + @test out.filter_modes == true + end +end From fd3411a55203efaa77c493bda4693218bc033deb Mon Sep 17 00:00:00 2001 From: logan-nc Date: Sat, 19 Sep 2026 09:11:37 -0400 Subject: [PATCH 18/23] ForceFreeStates/KineticForces - MINOR - Trim the review-flagged comments Review feedback on #414: the block comments are too long. Collapses the pinned-knot block in refined_psi_grid, the BAND_KNOTS and validity-envelope blocks in Kinetic.jl, and the sub-threshold reporting block in Finding.jl to one or two lines each, and rewraps the refined_psi_grid docstring line that had grown past the 180-char margin. Also folds find_kinetic_singular_surfaces!'s numbered "Algorithm:" list into prose, per the repo's no-step-numbering rule. docs/src/kinetic_forces.md keeps the regularization argument and its citations - it is user-facing reference material, not a comment - but drops the single-case DIII-D validation table, which belongs in the review package rather than in permanent docs where it goes stale. Adds the Logan (2015) Ch. 7 Eq. 7.46 reference alongside Park & Logan (2017), since it states the same result in GPEC's own notation. The three example decks get a shorter inline description for axis_validity_suppression, matching the surrounding comment column. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01VdRsGGg4YRG6E9fgLGbpdF --- docs/src/kinetic_forces.md | 21 +++++++------------ examples/DIIID-like_ideal_example/gpec.toml | 2 +- .../Solovev_kinetic_NTV_example/gpec.toml | 2 +- examples/a10_kinetic_example/gpec.toml | 2 +- src/Equilibrium/GridRefinement.jl | 9 ++++---- src/ForceFreeStates/Kinetic.jl | 11 ++++------ src/ForceFreeStates/Surfaces/Finding.jl | 14 ++++--------- 7 files changed, 22 insertions(+), 39 deletions(-) diff --git a/docs/src/kinetic_forces.md b/docs/src/kinetic_forces.md index 6cbde2a58..95cf10315 100644 --- a/docs/src/kinetic_forces.md +++ b/docs/src/kinetic_forces.md @@ -153,26 +153,19 @@ kinetic composite matrix as ``F_k = Q \bar F_k Q - P_l^\dagger Q - Q P_u + R_1`` real ``\psi`` axis, and the singularity is removed from both the solution and the torque integral. (Torque-free kinetic energy principles instead *shift and split* the zeros to ``\psi_r \mp r_{L,R}``, where the singularity is logarithmic and integrable — still not a case -for smoothing.) +for smoothing.) Logan, *Electromagnetic Torque in Tokamaks with Toroidal Asymmetries* (PhD thesis, +Princeton, 2015), Ch. 7 Eq. 7.46 states the same resonance shifting/splitting result in GPEC's own +notation. GPEC therefore **forces `reg_spot = 0` whenever `kinetic_factor > 0`**, logging the override. Leaving it on suppresses a finite physical response and does so inconsistently — ``\xi^\psi`` is never regularized, so damping the other two breaks their near-resonance cancellation in ``\delta B/B`` and leaves a spurious residue driving the NTV integrand. -Measured on the DIII-D-like H-mode case (n = 1, C-coil drive), comparing the NTV torque against -the Euler–Lagrange solution's own dissipation ``-2n\,\mathrm{Im}\langle \xi, u_2\rangle/4\mu_0`` -— two independent calculations of the same quantity: - -| configuration | max ``|\xi^\alpha|`` | NTV torque [N·m] | EL dissipation [N·m] | -|---|---|---|---| -| ideal, `reg_spot = 0` | 643.8 | 6074.3 | — | -| ideal, `reg_spot = 0.05` | 0.058 | 0.554 | — | -| kinetic, `reg_spot = 0.05` | 0.055 | 0.1655 | 0.1322 | -| kinetic, `reg_spot = 0` | 0.059 | **0.1324** | **0.1322** | - -The ideal rows show why the knob exists; the kinetic rows show why it must be off there — the -two independent torques agree to 0.15 % with no regularization, and to 20 % with it. +The check that this is right is that the NTV torque and the Euler–Lagrange solution's own +dissipation ``-2n\,\mathrm{Im}\langle \xi, u_2\rangle/4\mu_0`` — two independent calculations +of the same quantity — agree to a fraction of a percent with the regularization off, and disagree +by tens of percent with it on. ## HDF5 outputs: complex torque convention and the EnergyIntegrals layout diff --git a/examples/DIIID-like_ideal_example/gpec.toml b/examples/DIIID-like_ideal_example/gpec.toml index 148d8ebf8..9896f4c89 100644 --- a/examples/DIIID-like_ideal_example/gpec.toml +++ b/examples/DIIID-like_ideal_example/gpec.toml @@ -87,6 +87,6 @@ f0type = "maxwellian" # Equilibrium distribution moment = "pressure" # Pressure-moment NTV torque atol_xlmda = 1e-9 # Absolute tolerance for inner pitch + energy integrations rtol_xlmda = 1e-5 # Relative tolerance for inner pitch + energy integrations -axis_validity_suppression = true # Suppress kinetic terms where the zero-orbit-width ordering fails near the axis (profile-derived boundary, no tuning parameters) +axis_validity_suppression = true # Suppress kinetic terms inside the near-axis orbit-width validity limit write_outputs_to_HDF5 = true # Write outputs to the HDF5 file verbose = true # Enable verbose logging diff --git a/examples/Solovev_kinetic_NTV_example/gpec.toml b/examples/Solovev_kinetic_NTV_example/gpec.toml index e174c5ab9..0e45a65fd 100644 --- a/examples/Solovev_kinetic_NTV_example/gpec.toml +++ b/examples/Solovev_kinetic_NTV_example/gpec.toml @@ -87,4 +87,4 @@ f0fac = 1 # Scale toroidal field at constant pressure (β, q change; Φ, p, [KineticForces] kinetic_file = "kinetic.dat" # Kinetic profile file: psi_n, n_i, n_e, T_i, T_e, omega_E columns -axis_validity_suppression = true # Suppress kinetic terms where the zero-orbit-width ordering fails near the axis (profile-derived boundary, no tuning parameters) +axis_validity_suppression = true # Suppress kinetic terms inside the near-axis orbit-width validity limit diff --git a/examples/a10_kinetic_example/gpec.toml b/examples/a10_kinetic_example/gpec.toml index 7854db0e1..7eabffae0 100644 --- a/examples/a10_kinetic_example/gpec.toml +++ b/examples/a10_kinetic_example/gpec.toml @@ -59,4 +59,4 @@ nutype = "harmonic" # Collision operator (zero, small, krook, harmoni f0type = "maxwellian" # Distribution function (maxwellian, jkp, cgl) atol_xlmda = 1e-9 # Absolute tolerance for inner pitch + energy integrations rtol_xlmda = 1e-5 # Relative tolerance for inner pitch + energy integrations -axis_validity_suppression = true # Suppress kinetic terms where the zero-orbit-width ordering fails near the axis (profile-derived boundary, no tuning parameters) +axis_validity_suppression = true # Suppress kinetic terms inside the near-axis orbit-width validity limit diff --git a/src/Equilibrium/GridRefinement.jl b/src/Equilibrium/GridRefinement.jl index d9441bbe8..0b64337a7 100644 --- a/src/Equilibrium/GridRefinement.jl +++ b/src/Equilibrium/GridRefinement.jl @@ -418,7 +418,8 @@ Build the refined pass-2 ψ grid from a formed pass-1 equilibrium: measured-curv density (`_knot_density`), equidistribution, a global minimum-spacing floor (`enforce_min_spacing`), and rational-surface bracketing (`bracket_mandatory_nodes`). `tau` is the target interpolation accuracy (`psi_accuracy`); `kin` optionally supplies kinetic profiles -whose pedestal gradients attract knots; `mandatory` lists rational-surface ψ values to bracket; `pinned` lists ψ values inserted as plain knots without a cleared zone (kinetic-resonance surfaces); +whose pedestal gradients attract knots; `mandatory` lists rational-surface ψ values to bracket; +`pinned` lists ψ values inserted as plain knots without a cleared zone (kinetic-resonance surfaces); `singfac_min` and `n_min` (smallest |n| in the run) set each surface's matching half-stencil `dpsi = singfac_min/(n_min·|q′|)`, and the bracket half-width is `bracket_coef·dpsi` (floored at `min_spacing`). Rational surfaces are bracketed, not pinned: a knot on the surface would make the @@ -450,10 +451,8 @@ function refined_psi_grid(equil::PlasmaEquilibrium; N == N_cap && M_total > N_cap && @warn "refined_psi_grid: knot count capped at $N_cap (density integral wants $(ceil(Int, M_total))); psi_accuracy=$tau may not be attainable" grid = enforce_min_spacing(_equidistribute(xs, rho, N), min_spacing) - # Pinned knots (e.g. kinetic-resonance surfaces): knot-at-node semantics via - # merge_mandatory_nodes — no cleared zone. Inserted before rational bracketing, so a - # pinned node inside a rational's bracket zone is cleared by it (the Δ′ clean-interval - # requirement wins locally; the rational's own dense floor resolves that neighbourhood). + # Pinned knots use knot-at-node semantics (no cleared zone); inserted before rational bracketing + # so a rational's bracket still wins locally where the two overlap. grid = isempty(pinned) ? grid : merge_mandatory_nodes(grid, pinned) isempty(mandatory) && return grid min_half_widths = [max(bracket_coef * singfac_min / (n_min * abs(equil.profiles.q_deriv(m))), min_spacing) for m in mandatory] diff --git a/src/ForceFreeStates/Kinetic.jl b/src/ForceFreeStates/Kinetic.jl index e057276e1..d148b12e9 100644 --- a/src/ForceFreeStates/Kinetic.jl +++ b/src/ForceFreeStates/Kinetic.jl @@ -1,7 +1,6 @@ -# Knots laid across the envelope's [ψ_c, 2ψ_c] transition. The envelope is a quintic smoothstep, -# so a cubic spline needs several interior knots to follow it without overshoot; nine keeps the -# residual below the kernel's own tolerance even on coarse decks. +# Knots across the envelope's [ψ_c, 2ψ_c] quintic-smoothstep transition, enough for a cubic spline +# to follow it without overshoot even on coarse decks. const BAND_KNOTS = 9 """ @@ -104,10 +103,8 @@ function build_kinetic_matrix_splines( xs = metric.xs mpsi = length(xs) - # The near-axis validity envelope (KineticForces) has structure on the scale of the - # suppression boundary; coarse equilibrium grids cannot represent env·(increment), and the - # spline overshoot can land on a rational surface. Pin the band ends (the smoothstep is - # only C² there) and resolve the transition with a fixed set of knots. + # Coarse equilibrium grids cannot represent env·(increment) across the validity band, and the + # spline overshoot can land on a rational; pin the band ends (C² there) and resolve the middle. band_knots(lo, hi) = axis_validity_psi_c > 0 ? [x for x in range(axis_validity_psi_c, 2 * axis_validity_psi_c; length=BAND_KNOTS) if lo < x < hi] : Float64[] diff --git a/src/ForceFreeStates/Surfaces/Finding.jl b/src/ForceFreeStates/Surfaces/Finding.jl index 9cc7b6112..9edfe6e0d 100644 --- a/src/ForceFreeStates/Surfaces/Finding.jl +++ b/src/ForceFreeStates/Surfaces/Finding.jl @@ -257,12 +257,8 @@ Mirrors the intent of Fortran `ksing_find` (`sing.f:1486-1616`) which finds zero det(F̄) via adaptive bisection. Here we use condition number peaks instead of determinant zeros for better numerical robustness and scale invariance. -Algorithm: - - 1. Evaluate cond(F̄) on a dense ψ grid - 2. Find local maxima (peaks where gradient changes from + to -) - 3. Refine each peak with golden-section minimization of -cond - 4. Filter by threshold and resonance condition +Evaluates cond(F̄) on a dense ψ grid, finds its local maxima, refines each with golden-section +minimization of -cond, then filters by threshold and resonance condition. """ function find_kinetic_singular_surfaces!( mats::MatrixSplines, @@ -298,10 +294,8 @@ function find_kinetic_singular_surfaces!( local_maxima = [i for i in 2:(ngrid-1) if cond_vals[i] > cond_vals[i-1] && cond_vals[i] > cond_vals[i+1]] peak_indices = filter(i -> cond_vals[i] > cond_threshold, local_maxima) - # Peaks below the threshold are not singular surfaces, but they mark where the kinetic F̄ comes - # closest to singular — the shifted/split resonances of Park & Logan Eq. (70). Report the - # strongest few so sharp kinetic structure is visible rather than silent (on a DIII-D-like case - # these track the NTV torque-density peaks at low collisionality/rotation). + # Peaks below threshold mark near-singular structure (the shifted/split resonances of Park & + # Logan Eq. 70); report the strongest few so sharp kinetic structure is visible rather than silent. subthreshold = filter(i -> KINETIC_RELAXED_FRAC * cond_threshold < cond_vals[i] <= cond_threshold, local_maxima) if !isempty(subthreshold) top = sort(subthreshold; by=i -> cond_vals[i], rev=true)[1:min(3, length(subthreshold))] From 0f8dd766fb0a359d5cb18721e12df6b8b5c7491d Mon Sep 17 00:00:00 2001 From: logan-nc Date: Sat, 19 Sep 2026 09:11:37 -0400 Subject: [PATCH 19/23] KineticForces - MINOR - Compress the nested-tolerance commentary Review feedback on #414 flagged this block as too long. It is not actually this branch's code - it arrived on develop with the nested energy/pitch tolerances and shows in the pull request diff only because the merge base makes develop's commits look like additions - but the comment rule applies either way, so it is trimmed here rather than deferred: 27 comment lines down to 14, one line per field. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01VdRsGGg4YRG6E9fgLGbpdF --- src/KineticForces/KineticForcesStructs.jl | 25 ++++++----------------- 1 file changed, 6 insertions(+), 19 deletions(-) diff --git a/src/KineticForces/KineticForcesStructs.jl b/src/KineticForces/KineticForcesStructs.jl index 8834a689f..16d895ada 100644 --- a/src/KineticForces/KineticForcesStructs.jl +++ b/src/KineticForces/KineticForcesStructs.jl @@ -130,32 +130,19 @@ builds a second control for its differing tolerance). nn::Int = 1 # Toroidal mode number nl::Int = 1 # Bounce harmonic number - # Tolerances, outermost to innermost: ψ quadrature ⊃ λ (pitch) ⊃ x (energy). - # Each level must be resolved more tightly than the one enclosing it, or the outer - # integrator chases its integrand's own quadrature noise instead of converging. - # *_xlmda: tolerances for the λ (pitch) integration - # *_x: tolerances for the x (energy) integration nested inside it; NaN ⇒ derive as - # nested_tolerance_margin × the pitch tolerances - # *_psi: tolerances for the outer ψ quadrature + # Tolerances, outermost to innermost: ψ quadrature ⊃ λ (pitch) ⊃ x (energy); each level must + # resolve tighter than the one enclosing it, or the outer integrator chases quadrature noise. atol_xlmda::Float64 = 1e-8 # Absolute tolerance for the inner pitch integration rtol_xlmda::Float64 = 1e-5 # Relative tolerance for the inner pitch integration atol_x::Float64 = NaN # Absolute tolerance for the energy integration (NaN ⇒ derived) rtol_x::Float64 = NaN # Relative tolerance for the energy integration (NaN ⇒ derived) - # The pitch integrand IS the energy integral, so the energy level is resolved this much - # tighter than the pitch level by default. + # Energy tolerances default to this fraction of the pitch ones (the pitch integrand is itself an energy integral). nested_tolerance_margin::Float64 = 1e-2 # Factor relating derived energy tolerances to the pitch ones - # rtol_psi is the primary convergence knob: ~2 significant figures matches the validity - # of the NTV model approximations. Do not set it tighter than the noise floor of the - # inner integrals (keep rtol_psi ≳ 10 × rtol_xlmda). + # Primary convergence knob (~2 sig figs matches NTV model validity); keep rtol_psi ≳ 10 × rtol_xlmda or it chases inner-integral noise. rtol_psi::Float64 = 1e-2 # Relative tolerance for outer ψ quadrature - # atol_psi is in N·m and therefore amplitude-sensitive: NTV scales as δB², so a 10× - # weaker applied field gives a 100× smaller torque and any fixed absolute tolerance can - # silently dominate termination with O(1) relative error. Default 0 (rtol-only); a - # nonzero value is an expert opt-out for near-net-zero-torque cases and triggers a - # warning when it dominates. + # Amplitude-sensitive (NTV ∝ δB²); default 0 (rtol-only). Nonzero is an expert opt-out for near-zero-torque cases and warns when it dominates. atol_psi::Float64 = 0.0 # Absolute tolerance for outer ψ quadrature [N·m] - # Runaway guard for the outer quadrature (e.g. sign-cancelling torque density with tiny - # net torque under rtol-only control); a convergence warning fires when hit. + # Runaway guard for rtol-only control (e.g. sign-cancelling torque density); warns when hit. maxevals_psi::Int = 2000 # Max integrand evaluations for outer ψ quadrature # Scaling factors From 2abea51d292928328fcb55190c3f3bcfc4b29c95 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Sat, 19 Sep 2026 09:16:10 -0400 Subject: [PATCH 20/23] KineticForces - BUGFIX - Annotate the Validity group where it is written KineticForces/Validity carried no long_name/units at all. Its metadata table sat in HDF5Schema.MAIN_H5_ANNOTATIONS, which apply_main_h5_metadata! applies at the end of write_outputs_to_HDF5 - but write_validity! runs afterwards, in a separate h5open, so the pass never saw those datasets and the entries were dead. HDF5Schema's own header states the rule: sub-writers keep their tables next to their own writers. The table moves to Output.jl as KF_VALIDITY_H5_ANNOTATIONS and write_validity! annotates the group it just created, matching how write_to_hdf5! handles the per-method groups. The entries also gain the `dims` attribute the contract requires of rank >= 1 profiles. Corrects the rationale given in 2dcbbf170 for the cascade cap. That commit argued a transition band cutting a rational is harmless because the band knots resolve it; the history says otherwise - the band knots landed in 3f1dbfef9 and the NaN that motivated clear_rational_windows was found after them, in 37c4eebb4. The capped branch does reinstate the hazard, so the warning now says so and names the two ways out (axis_validity_suppression = false, or a finer grid there) instead of implying the knots cover it. The behaviour is unchanged: swallowing half the plasma is the worse failure. The contiguity warning gains maxlog=1 - it is reachable from four call sites, once per species - and its unit test now asserts the boundary rather than the log line, since maxlog makes the log order-dependent across the suite. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01VdRsGGg4YRG6E9fgLGbpdF --- src/HDF5Schema.jl | 19 ------------------- src/KineticForces/Output.jl | 24 ++++++++++++++++++++++++ src/KineticForces/Utils.jl | 14 +++++++++----- test/runtests_kinetic_validity.jl | 7 +++---- 4 files changed, 36 insertions(+), 28 deletions(-) diff --git a/src/HDF5Schema.jl b/src/HDF5Schema.jl index 512264af2..2f2f6f175 100644 --- a/src/HDF5Schema.jl +++ b/src/HDF5Schema.jl @@ -143,25 +143,6 @@ const MAIN_H5_ANNOTATIONS = [ "LocalStability/alpha_critical" => (; long_name="critical normalized pressure gradient α for first ballooning stability", dims=("psi_ballooning",), attach=(1 => "LocalStability/ballooning_psi",)), - # --- KineticForces/Validity/ --- - "KineticForces/Validity/psi" => (; long_name="normalized poloidal flux ψ_N of the kinetic validity profiles", scale="psi"), - "KineticForces/Validity/rho_i" => (; long_name="thermal ion gyroradius √(2mT)/(Z·e·B₀)", units="m", attach=(1 => "KineticForces/Validity/psi",)), - "KineticForces/Validity/rho_banana" => (; long_name="thermal banana orbit width q·ρ_i/√ε", units="m", attach=(1 => "KineticForces/Validity/psi",)), - "KineticForces/Validity/rho_theta" => (; long_name="thermal poloidal gyroradius q·ρ_i/ε", units="m", attach=(1 => "KineticForces/Validity/psi",)), - "KineticForces/Validity/w_potato" => (; long_name="potato orbit width (q²ρ_i²R₀)^(1/3)", units="m", attach=(1 => "KineticForces/Validity/psi",)), - "KineticForces/Validity/r_minor" => (; long_name="surface-average minor radius ⟨r⟩", units="m", attach=(1 => "KineticForces/Validity/psi",)), - "KineticForces/Validity/L_p" => (; long_name="pressure gradient scale length |p|/|dp/dr|", units="m", attach=(1 => "KineticForces/Validity/psi",)), - "KineticForces/Validity/L_q" => (; long_name="safety-factor gradient scale length |q|/|dq/dr|", units="m", attach=(1 => "KineticForces/Validity/psi",)), - "KineticForces/Validity/d_separatrix" => (; long_name="distance to the separatrix ⟨r⟩(1) − ⟨r⟩(ψ)", units="m", attach=(1 => "KineticForces/Validity/psi",)), - "KineticForces/Validity/psi_c" => (; long_name="near-axis kinetic validity boundary: outermost ψ_N where a thermal orbit width reaches ⟨r⟩", units="1"), - "KineticForces/Validity/envelope" => - (; long_name="near-axis suppression envelope applied to the calculated kinetic terms (1 = unsuppressed)", units="1", attach=(1 => "KineticForces/Validity/psi",)), - "KineticForces/Validity/is_valid" => (; - long_name="1 where every zero-orbit-width ordering holds: max orbit width < ⟨r⟩, ρ_banana < L_p and L_q, max orbit width < d_separatrix", - units="1", - attach=(1 => "KineticForces/Validity/psi",) - ), - # --- ForceFreeStates/Solutions/ForwardIntegration/ --- "ForceFreeStates/Solutions/ForwardIntegration/nstep" => (; long_name="number of saved solution snapshots"), "ForceFreeStates/Solutions/ForwardIntegration/nstep_total" => (; long_name="total ODE solver steps taken"), diff --git a/src/KineticForces/Output.jl b/src/KineticForces/Output.jl index 771873251..867fafff4 100644 --- a/src/KineticForces/Output.jl +++ b/src/KineticForces/Output.jl @@ -224,5 +224,29 @@ function write_validity!(h5file::HDF5.File, kf_ctrl::KineticForcesControl, speci g["psi_c"] = psi_c g["envelope"] = psi_c > 0 ? kinetic_axis_validity_envelope.(vp.psi, psi_c) : ones(length(vp.psi)) g["is_valid"] = Int8.(vp.is_valid) + Utilities.HDF5Annotations.annotate!(g, KF_VALIDITY_H5_ANNOTATIONS) return nothing end + +# Metadata table for KineticForces/Validity/ (paths relative to the group). It lives here, not in +# HDF5Schema, because this group is written after the main writer's annotation pass has run. +const KF_VALIDITY_H5_ANNOTATIONS = [ + "psi" => (; long_name="normalized poloidal flux ψ_N of the kinetic validity profiles", scale="psi"), + "rho_i" => (; long_name="thermal ion gyroradius √(2mT)/(Z·e·B₀)", units="m", dims=("psi",), attach=(1 => "psi",)), + "rho_banana" => (; long_name="thermal banana orbit width q·ρ_i/√ε", units="m", dims=("psi",), attach=(1 => "psi",)), + "rho_theta" => (; long_name="thermal poloidal gyroradius q·ρ_i/ε", units="m", dims=("psi",), attach=(1 => "psi",)), + "w_potato" => (; long_name="potato orbit width (q²ρ_i²R₀)^(1/3)", units="m", dims=("psi",), attach=(1 => "psi",)), + "r_minor" => (; long_name="surface-average minor radius ⟨r⟩", units="m", dims=("psi",), attach=(1 => "psi",)), + "L_p" => (; long_name="pressure gradient scale length |p|/|dp/dr|", units="m", dims=("psi",), attach=(1 => "psi",)), + "L_q" => (; long_name="safety-factor gradient scale length |q|/|dq/dr|", units="m", dims=("psi",), attach=(1 => "psi",)), + "d_separatrix" => (; long_name="distance to the separatrix ⟨r⟩(1) − ⟨r⟩(ψ)", units="m", dims=("psi",), attach=(1 => "psi",)), + "psi_c" => (; long_name="near-axis kinetic validity boundary: outermost ψ_N where a thermal orbit width reaches ⟨r⟩", units="1"), + "envelope" => + (; long_name="near-axis suppression envelope applied to the calculated kinetic terms (1 = unsuppressed)", units="1", dims=("psi",), attach=(1 => "psi",)), + "is_valid" => (; + long_name="1 where every zero-orbit-width ordering holds: max orbit width < ⟨r⟩, ρ_banana < L_p and L_q, max orbit width < d_separatrix", + units="1", + dims=("psi",), + attach=(1 => "psi",) + ) +] diff --git a/src/KineticForces/Utils.jl b/src/KineticForces/Utils.jl index d2281de21..2cf25c2f5 100644 --- a/src/KineticForces/Utils.jl +++ b/src/KineticForces/Utils.jl @@ -149,7 +149,7 @@ function kinetic_axis_validity_psi(kinetic_profiles::Equilibrium.KineticProfileS end end psi_outer > 0 && @warn "kinetic_axis_validity_psi: orbit widths also reach ⟨r⟩ out to ψ=$(round(psi_outer; sigdigits=3)), " * - "outside the axis-contiguous region ending at ψ_c=$(round(psi_c; sigdigits=3)); only the latter is suppressed" + "outside the axis-contiguous region ending at ψ_c=$(round(psi_c; sigdigits=3)); only the latter is suppressed" maxlog = 1 return psi_c end @@ -267,9 +267,11 @@ reaches ⟨r⟩ the resonance is not trustworthy anyway, so the boundary is push (suppressing it wholly) rather than cutting through it. The band widens as ψ_c moves, so clearing can cascade through a rational-dense core. The move is -capped at `VALIDITY_CLEAR_MAX_FACTOR`× the orbit-width boundary; past that the uncleared boundary -is kept with a warning, since suppressing a large fraction of the plasma is worse than a band that -cuts a rational (which the band knots resolve). +capped at `VALIDITY_CLEAR_MAX_FACTOR`× the orbit-width boundary; past that the uncleared boundary is +kept with a warning. That reinstates the hazard this function exists to remove — the band knots +alone did not prevent it — so the warning says so: suppressing a large fraction of the plasma is the +worse failure, but the capped run needs `axis_validity_suppression = false` or a finer ψ grid if the +solve then produces NaNs. """ function clear_rational_windows(psi_c::Float64, rationals::Vector{Float64})::Float64 psi_c <= 0 && return psi_c @@ -281,7 +283,9 @@ function clear_rational_windows(psi_c::Float64, rationals::Vector{Float64})::Flo psi_c = maximum(inside) + Equilibrium.RATIONAL_RES_RADIUS if psi_c > limit @warn "clear_rational_windows: clearing the envelope band of rationals would push ψ_c from " * - "$(round(psi_c0; sigdigits=3)) past $(round(limit; sigdigits=3)); keeping the orbit-width boundary" + "$(round(psi_c0; sigdigits=3)) past $(round(limit; sigdigits=3)); keeping the orbit-width boundary. " * + "A rational then sits in the transition band — if the stability solve produces NaNs, set " * + "axis_validity_suppression = false or refine the ψ grid there." return psi_c0 end end diff --git a/test/runtests_kinetic_validity.jl b/test/runtests_kinetic_validity.jl index c84b1d217..5a5eac8c2 100644 --- a/test/runtests_kinetic_validity.jl +++ b/test/runtests_kinetic_validity.jl @@ -83,11 +83,10 @@ # Same profile with a narrow q spike far from the axis. The spike re-triggers the criterion, # but the model has not lost validity there, so the boundary must not jump out to it. + # The accompanying @warn carries maxlog=1, so whether it fires here depends on what ran + # earlier in the suite; the boundary itself is the contract worth asserting. q_spiked(psi) = q_smooth(psi) + 400.0 * exp(-((psi - 0.9) / 0.01)^2) - local psi_c_spiked - @test_logs (:warn,) match_mode = :any begin - psi_c_spiked = KF.kinetic_axis_validity_psi(mock_profiles(xs, Ti), mock_equil(q_spiked, r_of, R_of)) - end + psi_c_spiked = KF.kinetic_axis_validity_psi(mock_profiles(xs, Ti), mock_equil(q_spiked, r_of, R_of)) @test psi_c_spiked ≈ psi_c @test psi_c_spiked < 0.5 end From a9bc438f777d4b99cf0244544a5ef28d25c67521 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Sat, 3 Oct 2026 18:21:22 -0400 Subject: [PATCH 21/23] Equilibrium - BUGFIX! - Sample every traced surface at common straight-fieldline angles through one tracer contract The common-abscissae fix applied only to direct_fieldline_int. efit_arclength still resampled each surface from its own adaptive solver steps, the same white-in-psi noise, and equilibrium_solver chose between the two behaviours on `sol === nothing`, so any tracer that returned no dense solution silently kept the noisy path. Sampling is now part of the tracer contract. Every tracer takes theta_nodes and returns its surface already sampled there, through one shared helper (sample_trace_at_sfl_angles, which root-solves the monotone SFL integral on the dense solution in the tracer's own integration variable). equilibrium_solver keeps a single path with no duplicated formulas. The arclength tracer now solves with dense output, brackets in arc length rather than eta, and unwraps eta at each sample against its bracketing step. The dense solution no longer leaves the tracer, so the concrete return annotation is restored and the type instability introduced earlier on this branch is gone. The root-solve brackets from the saved steps instead of re-evaluating the dense solution, and the unreachable linear-interpolation fallback is now an error. Direct-tracer results are unchanged; efit_arclength results move. Measured on the EQDSK_COCOS_02 fixture, efit and efit_arclength rzphi now agree to 2-3e-6, down from 2.4-3.9e-4. Also backports the length(xs) < 3 guard in core_capped_knots, which otherwise returns a duplicate knot on a one-knot grid. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01VdRsGGg4YRG6E9fgLGbpdF --- src/Equilibrium/DirectEquilibrium.jl | 119 +++++++++--------- src/Equilibrium/DirectEquilibriumArcLength.jl | 22 ++-- src/ForceFreeStates/Fourfit.jl | 1 + 3 files changed, 75 insertions(+), 67 deletions(-) diff --git a/src/Equilibrium/DirectEquilibrium.jl b/src/Equilibrium/DirectEquilibrium.jl index 586d81001..d22756413 100644 --- a/src/Equilibrium/DirectEquilibrium.jl +++ b/src/Equilibrium/DirectEquilibrium.jl @@ -228,33 +228,51 @@ function direct_position!(raw_profile::DirectRunInput) end """ - eta_at_sfl_angle(sol, y_out, x, total_x) -> Float64 + sfl_sample_time(sol, ts, sfl, sfl_index, x) -> (t, lo) -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). +Solver time `t` at which the normalised straight-fieldline angle reaches `x` ∈ [0, 1], and the index +`lo` of the saved step that starts its bracket. `sfl` is the SFL integral ∫jac·dl/Bp at the saved +times `ts`; it is monotone, so it brackets the root to one step, and Brent finishes on component +`sfl_index` of the dense solution. """ -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)) +function sfl_sample_time(sol, ts::Vector{Float64}, sfl::AbstractVector{Float64}, sfl_index::Int, x::Float64) + n = length(ts) + x <= 0 && return ts[1], 1 + x >= 1 && return ts[n], n - 1 + target = x * sfl[n] + hi = clamp(searchsortedfirst(sfl, target), 2, n) 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()) + flo, fhi = sfl[lo] - target, sfl[hi] - target + flo == 0 && return ts[lo], lo + fhi == 0 && return ts[hi], lo + flo < 0 < fhi || error("sfl_sample_time: SFL integral is not monotone across a solver step (x = $x)") + return find_zero(t -> sol(t)[sfl_index] - target, (ts[lo], ts[hi]), Roots.Brent()), lo end """ - direct_fieldline_int(psifac, raw_profile, ro, zo, rs2) + sample_trace_at_sfl_angles(sol, y_out, sfl_index, theta_nodes, state_row) -> Matrix{Float64} + +Sample a traced flux surface at the straight-fieldline angles `theta_nodes` (0 to 1). Row `i` is +`state_row(sol(t), t, lo)` at the time `t` where the normalised SFL angle reaches `theta_nodes[i]` +(see `sfl_sample_time`), in the 5-column `y_out` layout. The last row is the saved closed endpoint, +so the surface totals `y_out[end, :]` are exact. Every surface sharing the same abscissae keeps the +resample error from varying between neighbouring surfaces. Every field-line tracer must return its +surface through this function. +""" +function sample_trace_at_sfl_angles(sol, y_out::Matrix{Float64}, sfl_index::Int, theta_nodes, state_row::F) where {F} + ts = sol.t::Vector{Float64} + sfl = view(y_out, :, 5) + samples = Matrix{Float64}(undef, length(theta_nodes), 5) + for i in 1:(length(theta_nodes)-1) + t, lo = sfl_sample_time(sol, ts, sfl, sfl_index, Float64(theta_nodes[i])) + samples[i, :] .= state_row(sol(t), t, lo) + end + samples[end, :] .= @view y_out[end, :] + return samples +end + +""" + direct_fieldline_int(psifac, raw_profile, ro, zo, rs2, theta_nodes) -> (y_out, bfield) Performs the field-line integration for a single flux surface. This is a Julia adaptation of the Fortran `direct_fl_int` subroutine. Note that the array `y_out` is now indexed @@ -266,22 +284,22 @@ from 1:5 rather than 0:4 as in Fortran. - `raw_profile`: `DirectRunInput` object containing splines and parameters. - `ro`, `zo`: Coordinates of the magnetic axis [m]. - `rs2`: R-coordinate of the outboard separatrix [m]. + - `theta_nodes`: straight-fieldline angles (0 to 1) at which to sample the surface. ## Returns: - - `y_out`: A matrix containing the integrated quantities vs. the geometric angle `η`. + - `y_out`: A matrix containing the integrated quantities at each of `theta_nodes` + (see `sample_trace_at_sfl_angles`). - `y_out[:, 1]`: η (geometric poloidal angle) - `y_out[:, 2]`: ∫(dl/Bp) - `y_out[:, 3]`: rfac (radial distance from magnetic axis) - `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) +function direct_fieldline_int(psifac::Float64, raw_profile::DirectRunInput, ro::Float64, zo::Float64, rs2::Float64, + theta_nodes)::Tuple{Matrix{Float64},DirectBField} # Find the starting point on the flux surface (outboard midplane) psi0_guess = raw_profile.psio * (1.0 - psifac) @@ -320,12 +338,13 @@ 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) - # 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). + # Dense output so the surface can be sampled at prescribed SFL angles between solver steps. 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, sol + y_out = hcat(sol.t::Vector{Float64}, sol_matrix) + # State is [∫dl/Bp, rfac, ∫dl/(R²Bp), ∫jac·dl/Bp] and the integration variable is η itself. + return sample_trace_at_sfl_angles(sol, y_out, 4, theta_nodes, (u, eta, _) -> (eta, u[1], u[2], u[3], u[4])), bfield end """ @@ -495,6 +514,9 @@ robustness. - `raw_profile`: A `DirectRunInput` object containing the initial splines (`psi_in`, `sq_in`) and run parameters (`equil_input`). + - `fieldline_int`: the flux-surface tracer, called as `fieldline_int(psifac, raw_profile, ro, zo, + rs2, theta_nodes) -> (y_out, bfield)` with `y_out` already sampled at `theta_nodes` through + `sample_trace_at_sfl_angles`. ## Returns: @@ -520,45 +542,22 @@ robustness. mpsi = length(psi_nodes) - 1 theta_nodes = range(0.0, 1.0; length=mtheta + 1) + ff_x_nodes = collect(theta_nodes) sq_nodes = zeros!(pool, Float64, mpsi + 1, 4) rzphi_nodes = zeros!(pool, Float64, mpsi + 1, mtheta + 1, 4) ff_val = zeros!(pool, Float64, 4) ff_deriv_val = zeros!(pool, Float64, 4) for ipsi in (mpsi+1):-1:1 # outermost to innermost - y_out, bfield, sol = fieldline_int(psi_nodes[ipsi], raw_profile, ro, zo, rs2) + # Every tracer returns the surface sampled at theta_nodes, so all surfaces share SFL abscissae. + y_out, bfield = fieldline_int(psi_nodes[ipsi], raw_profile, ro, zo, rs2, theta_nodes) checkpoint!(pool, Float64) - # 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 = acquire!(pool, Float64, mtheta + 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 ff_fs_nodes[end, :] .= ff_fs_nodes[1, :] # enforce periodic endpoint diff --git a/src/Equilibrium/DirectEquilibriumArcLength.jl b/src/Equilibrium/DirectEquilibriumArcLength.jl index ec4c661d3..e88559450 100644 --- a/src/Equilibrium/DirectEquilibriumArcLength.jl +++ b/src/Equilibrium/DirectEquilibriumArcLength.jl @@ -59,10 +59,11 @@ so the solver never restricts step size for them. end """ - arclength_fieldline_int(psifac, raw_profile, ro, zo, rs2) + arclength_fieldline_int(psifac, raw_profile, ro, zo, rs2, theta_nodes) -> (y_out, bfield) Arc-length-parameterized flux surface integration. Drop-in replacement for -`direct_fieldline_int` with identical return format: +`direct_fieldline_int` with identical return format, sampled at the straight-fieldline angles +`theta_nodes` (see `sample_trace_at_sfl_angles`): - `y_out[:, 1]`: geometric angle η ∈ 0 to 2π (CCW from outboard midplane) - `y_out[:, 2]`: accumulated ∫dl/Bp @@ -74,8 +75,8 @@ The ODE is terminated by a `ContinuousCallback` that detects the return to the 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,Nothing} + psifac::Float64, raw_profile::DirectRunInput, ro::Float64, zo::Float64, rs2::Float64, theta_nodes +)::Tuple{Matrix{Float64},DirectBField} psi0_guess = raw_profile.psio * (1.0 - psifac) r = ro + sqrt(psifac) * (rs2 - ro) @@ -114,13 +115,14 @@ outboard midplane (Z = zo, R > ro) after a minimum arc-length guard. reltol_vec = [equil_config.etol, equil_config.etol, 1e20, 1e20, 1e20] abstol_vec = [1e-8, 1e-8, 1e20, 1e20, 1e20] sol = solve(prob, BS5(); callback=callback, reltol=reltol_vec, abstol=abstol_vec, - dt=2π / 200, adaptive=true, dense=false) + dt=2π / 200, adaptive=true, dense=true) n = length(sol.u) y_out = Matrix{Float64}(undef, n, 5) # Compute geometric angles and unwrap to monotone [0, 2π] - η_raw = [atan(sol.u[i][2] - zo, sol.u[i][1] - ro) for i in 1:n] + η_atan = [atan(sol.u[i][2] - zo, sol.u[i][1] - ro) for i in 1:n] + η_raw = copy(η_atan) for i in 2:n Δ = η_raw[i] - η_raw[i-1] if Δ > π @@ -141,7 +143,13 @@ outboard midplane (Z = zo, R > ro) after a minimum arc-length guard. y_out[i, 5] = sol.u[i][5] end + # The integration variable is arc length, so η at a sample is unwrapped against its bracketing step. + function state_row(u, _s, lo) + η = y_out[lo, 1] + rem2pi(atan(u[2] - zo, u[1] - ro) - η_atan[lo], RoundNearest) + return (η, u[3], sqrt((u[1] - ro)^2 + (u[2] - zo)^2), u[4], u[5]) + end + # bfield at the starting point carries F and P for the surface-averaged quantities - return y_out, bfield, nothing + return sample_trace_at_sfl_angles(sol, y_out, 5, theta_nodes, state_row), bfield end diff --git a/src/ForceFreeStates/Fourfit.jl b/src/ForceFreeStates/Fourfit.jl index 0f8095e3c..57cc28128 100644 --- a/src/ForceFreeStates/Fourfit.jl +++ b/src/ForceFreeStates/Fourfit.jl @@ -443,6 +443,7 @@ rational's resolution window, preserving the Δ′-stencil structure the equilib (`Equilibrium.RATIONAL_RES_RADIUS`). """ function core_capped_knots(xs::Vector{Float64}, rationals::Vector{Float64})::Vector{Int} + length(xs) < 3 && return collect(eachindex(xs)) cap_edge = 0.1 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) From bfeaef972180dad405bf45c778096af5ee5f138b Mon Sep 17 00:00:00 2001 From: logan-nc Date: Sat, 3 Oct 2026 18:21:22 -0400 Subject: [PATCH 22/23] Equilibrium - TEST - Check traced straight-fieldline coordinates do not depend on solver steps The common-SFL-abscissae fix had no test coverage. The existing equilibrium tests compare at 1e-3 to 10 % and are blind to the ~1e-6 white-in-psi noise it removes. Two checks of the property itself rather than the mechanism: - Solovev rzphi at off-grid points is independent of psi resolution (64 vs 128). - On the EQDSK_COCOS_02 fixture, the efit and efit_arclength tracers agree. They take different step sequences along the same surfaces, so any disagreement is step dependence. Each threshold sits about 10x or more from both sides of the measured behaviour, and was mutation-checked: all four checks fail on develop (old resampling in both tracers), and the tracer check fails on this branch before the arclength tracer was converted. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01VdRsGGg4YRG6E9fgLGbpdF --- test/runtests.jl | 1 + test/runtests_sfl_sampling.jl | 41 +++++++++++++++++++++++++++++++++++ 2 files changed, 42 insertions(+) create mode 100644 test/runtests_sfl_sampling.jl diff --git a/test/runtests.jl b/test/runtests.jl index 7845eea60..e47e01017 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -25,6 +25,7 @@ else include("./runtests_fouriertransforms.jl") include("./runtests_vacuum.jl") include("./runtests_equil.jl") + include("./runtests_sfl_sampling.jl") include("./runtests_grid_refinement.jl") include("./runtests_coordinate_invariant.jl") include("./runtests_eulerlagrange.jl") diff --git a/test/runtests_sfl_sampling.jl b/test/runtests_sfl_sampling.jl new file mode 100644 index 000000000..f08646791 --- /dev/null +++ b/test/runtests_sfl_sampling.jl @@ -0,0 +1,41 @@ +using Test +using GeneralizedPerturbedEquilibrium + +# The straight-fieldline coordinate a traced equilibrium produces must not depend on where the ODE +# solver happened to place its steps. Sampling each surface at its own solver steps and resampling +# leaves an error that varies from surface to surface (white noise in ψ). These checks fail by +# 10–1000× when that happens, against either tracer. +@testset "Straight-fieldline sampling is independent of solver steps" begin + Eq = GeneralizedPerturbedEquilibrium.Equilibrium + psis = range(0.21, 0.89; length=23) + thetas = range(0.013, 0.97; length=17) + # Max over off-grid points of |a − b|, relative to the largest |a|. + function max_rel_diff(a, b) + va = [a((p, t)) for p in psis, t in thetas] + vb = [b((p, t)) for p in psis, t in thetas] + return maximum(abs.(va .- vb)) / maximum(abs.(va)) + end + + @testset "Solovev: rzphi independent of ψ resolution" begin + sol_eq(mpsi) = Eq.equilibrium_solver(Eq.sol_run( + Eq.EquilibriumConfig(; eq_type="sol", eq_filename="unused", jac_type="pest", grid_type="ldp", + psilow=1e-4, psihigh=0.99999, mpsi=mpsi, mtheta=128), + Eq.SolovevConfig(64, 64, 64, 1.6, 0.33, 1.0, 1.9, 1.0, 1.0, 1.0))) + coarse, fine = sol_eq(64), sol_eq(128) + # ν is identically zero for this equilibrium, so it carries no signal. + for field in (:rzphi_offset, :rzphi_jac, :rzphi_rsquared) + # Measured: ≤ 2.2e-7 with common SFL abscissae, ≥ 2.1e-4 with per-surface resampling. + @test max_rel_diff(getfield(coarse, field), getfield(fine, field)) < 1e-5 + end + end + + @testset "EFIT: efit and efit_arclength tracers agree" begin + efit_file = joinpath(@__DIR__, "test_data", "CHEASE_test_data", "EQDSK_COCOS_02") + efit_eq(kind) = Eq.setup_equilibrium(Eq.EquilibriumConfig(; eq_filename=efit_file, eq_type=kind, + jac_type="boozer", grid_type="ldp", psilow=0.01, psihigh=0.994, mpsi=64, mtheta=128)) + direct, arclength = efit_eq("efit"), efit_eq("efit_arclength") + # Two integrators with different step sequences along the same surfaces. Measured on the SFL + # Jacobian: 2.3e-6 when both sample at common abscissae, ≥ 2.4e-4 if either resamples its steps. + @test max_rel_diff(direct.rzphi_jac, arclength.rzphi_jac) < 2e-5 + end +end From 2311583c6b65d7e7bce461e54a9819e855dd0a52 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Sat, 3 Oct 2026 19:18:41 -0400 Subject: [PATCH 23/23] Equilibrium - BUGFIX! - Evaluate F and P at the traced surface rather than at the solver's last evaluation direct_fieldline_int returned its ODE scratch DirectBField, which the right-hand side, the surface-refinement callback and the dense interpolant overwrite on every evaluation. F, P and everything built on them (q, nu, the F and P profiles, beta) were therefore read wherever the solver last happened to evaluate, not at the surface start point the docstring promises. The arclength tracer already kept a separate buffer for the ODE (bfield_ode) and returned the start point. The direct tracer now does the same. The change is about 1e-10 in q and nu on a fixed grid, and the sampled geometry (r^2, angle offset, Jacobian) is bit-identical. It also makes F and P independent of evaluation order: the old value moved whenever sampling order changed, which is how this was found. Applying this one fix to both the old and the refactored sampling makes their equilibria bitwise identical, so the sampling restructure in the previous commit is otherwise a pure restructure. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01VdRsGGg4YRG6E9fgLGbpdF --- src/Equilibrium/DirectEquilibrium.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/Equilibrium/DirectEquilibrium.jl b/src/Equilibrium/DirectEquilibrium.jl index d22756413..0fe1f9430 100644 --- a/src/Equilibrium/DirectEquilibrium.jl +++ b/src/Equilibrium/DirectEquilibrium.jl @@ -324,10 +324,11 @@ function direct_fieldline_int(psifac::Float64, raw_profile::DirectRunInput, ro:: u0 = zeros(Float64, 4) u0[2] = sqrt((r - ro)^2 + (z - zo)^2) - bfield = DirectBField() + # The ODE right-hand side overwrites its own field buffer, so keep `bfield` at the start point (F, P at ψ₀). + bfield_ode = DirectBField() equil_config = raw_profile.config params = FieldLineDerivParams(ro, zo, raw_profile.psi_in, raw_profile.sq_in, sq_in_deriv, raw_profile.psio, - equil_config.power_bp, equil_config.power_b, equil_config.power_r, equil_config.power_rc, bfield) + equil_config.power_bp, equil_config.power_b, equil_config.power_r, equil_config.power_rc, bfield_ode) # Use a callback to refine the solution at each step to stay on the flux surface function refine_affect!(integrator)