From 2afc6c3ebe732b2f99083253c12a48bb9aef6c3f Mon Sep 17 00:00:00 2001 From: logan-nc Date: Sat, 15 Aug 2026 13:58:02 -0400 Subject: [PATCH 01/14] KineticForces - BUG FIX - fit the parallel-velocity spline periodically The bounce-point spline samples 1-(lmda/bo)*B(theta), which is periodic on the closed poloidal interval, but was fitted with non-periodic endpoint conditions while the adjacent equilibrium B spline (tspl) was already periodic. The resulting fit is only C0 at the theta=0/1 seam, so it can manufacture false near-seam extrema and bounce-root pairs. On a synthetic shaped field the seam derivative mismatch drops ~5x (1.4e-4 -> 2.6e-5) when the fit is made periodic. Fixed at both construction sites (the tpsi! quadrature path and the _setup_surface_state kinetic-matrix path). Renamed B_extrap -> B_vpar and _vpar_from_extrap -> _vpar_from_spline, since after this change the old names describe a fit that no longer exists. Julia analog of Fortran GPEC PR #284, which fixed the same spline_fit(vspl, "extrap") -> "periodic" defect in pentrc/torque.F90 after it produced non-finite omega_b/omega_D and an LSODE failure on an ITER case. Co-Authored-By: Claude Fable 5 --- src/KineticForces/BounceAveraging.jl | 38 ++++++++++++++-------------- src/KineticForces/Torque.jl | 25 ++++++++++-------- 2 files changed, 33 insertions(+), 30 deletions(-) diff --git a/src/KineticForces/BounceAveraging.jl b/src/KineticForces/BounceAveraging.jl index 63a59aedc..9db0033fe 100644 --- a/src/KineticForces/BounceAveraging.jl +++ b/src/KineticForces/BounceAveraging.jl @@ -147,7 +147,7 @@ end """ compute_bounce_data(psi, n, l, q, bo, bmax, bmin, theta_bmax, - tspl, B_extrap, mfac, chi1, ro, dbob_m_f, divx_m_f, + tspl, B_vpar, mfac, chi1, ro, dbob_m_f, divx_m_f, divxfac, wdfac, mass, chrg, T_s, method; nlmda=128, ntheta=128, smat=nothing, tmat=nothing, xmat=nothing, @@ -168,8 +168,8 @@ Ports Fortran torque.F90 lines 530-816 (GAR branch). - `bmax, bmin`: Max/min of B(θ) at this ψ - `theta_bmax`: θ location of Bmax (nodal knot; the passing-transit start) - `tspl`: Periodic poloidal interpolant: tspl(θ) → [B, dB/dψ, dB/dθ, J, dJ/dψ] -- `B_extrap`: Endpoint-fit (non-periodic) cubic of B(θ) used for v_par and the - bounce-point roots (the Fortran `vspl` equivalent) +- `B_vpar`: Periodic cubic of B(θ) used for v_par and the bounce-point roots + (the Fortran `vspl` equivalent) - `mfac`: Poloidal mode numbers [mlow:mhigh] - `chi1`: 2π·ψ₀ flux normalization - `ro`: Major radius [m] @@ -190,7 +190,7 @@ function compute_bounce_data( psi::Float64, n::Int, l::Int, q::Float64, bo::Float64, bmax::Float64, bmin::Float64, theta_bmax::Float64, - tspl, B_extrap, mfac::Vector{Int}, chi1::Float64, ro::Float64, + tspl, B_vpar, mfac::Vector{Int}, chi1::Float64, ro::Float64, dbob_m_f::Vector{ComplexF64}, divx_m_f::Vector{ComplexF64}, divxfac::Float64, wdfac::Float64, mass::Float64, chrg::Float64, @@ -241,12 +241,12 @@ function compute_bounce_data( # Find bounce points and build θ sub-grid _, _, tdt_pts, tdt_wts = _find_bounce_points_and_grid( - lmda, bo, sigma, B_extrap, theta_bmax, psi, ntheta) + lmda, bo, sigma, B_vpar, theta_bmax, psi, ntheta) # Bounce integrals over θ (Fortran lines 674-735) wbbar, wdbar, dJdJ_val, wmats_lmda = _bounce_integrate( tdt_pts, tdt_wts, lmda, lnq, sigma, n, q, bo, - tspl, B_extrap, chi1, ro, mfac, dbob_m_f, divx_m_f, divxfac, wdfac, + tspl, B_vpar, chi1, ro, mfac, dbob_m_f, divx_m_f, divxfac, wdfac, do_matrices, mpert, smat, tmat, xmat, ymat, zmat, scr) # Physical frequencies (Fortran lines 744-745) @@ -409,12 +409,12 @@ end """ -Parallel-velocity factor `v_par = 1 − (λ/bo)·B(θ)` from the endpoint-fit cubic of B -(`B_extrap`, built where the surface interpolants are constructed), keeping v_par +Parallel-velocity factor `v_par = 1 − (λ/bo)·B(θ)` from the periodic cubic of B +(`B_vpar`, built where the surface interpolants are constructed), keeping v_par consistent with the bounce-point roots as in Fortran's `vspl`. """ -@inline _vpar_from_extrap(B_extrap, lmda::Float64, bo::Float64, θ::Float64) = - 1.0 - (lmda / bo) * B_extrap(mod(θ, 1.0)) +@inline _vpar_from_spline(B_vpar, lmda::Float64, bo::Float64, θ::Float64) = + 1.0 - (lmda / bo) * B_vpar(mod(θ, 1.0)) """ @@ -423,14 +423,14 @@ Returns (t1, t2, theta_points, theta_weights). """ function _find_bounce_points_and_grid( lmda::Float64, bo::Float64, sigma::Int, - B_extrap, theta_bmax::Float64, psi::Float64, + B_vpar, theta_bmax::Float64, psi::Float64, ntheta::Int ) if sigma == 0 # trapped - # Bounce points: all roots of v_par(θ) = 1 − (λ/bo)·B_extrap(θ) in (0,1), + # Bounce points: all roots of v_par(θ) = 1 − (λ/bo)·B_vpar(θ) in (0,1), # sorted descending — the same order as Fortran spline_roots, which the # marginally-trapped and deepest-well wrap logic below assume. - vpar_fn = θ -> _vpar_from_extrap(B_extrap, lmda, bo, θ) + vpar_fn = θ -> _vpar_from_spline(B_vpar, lmda, bo, θ) bpts = sort!(Roots.find_zeros(vpar_fn, 0.0, 1.0); rev=true) nbpts = length(bpts) @@ -443,7 +443,7 @@ function _find_bounce_points_and_grid( t1 = bpts[1] t2 = bpts[1] + 1.0 else - t1, t2 = _find_deepest_well(bpts, B_extrap, lmda, bo) + t1, t2 = _find_deepest_well(bpts, B_vpar, lmda, bo) end # Power-law grid refined near bounce points @@ -463,7 +463,7 @@ end Find the deepest potential well (largest midpoint v_par) among bounce-point pairs, handling pairs that wrap through θ = 0/1. """ -function _find_deepest_well(bpts::Vector{Float64}, B_extrap, lmda::Float64, bo::Float64) +function _find_deepest_well(bpts::Vector{Float64}, B_vpar, lmda::Float64, bo::Float64) nbpts = length(bpts) best_vpar = 0.0 best_t1 = 0.0 @@ -477,7 +477,7 @@ function _find_deepest_well(bpts::Vector{Float64}, B_extrap, lmda::Float64, bo:: else θmid = 0.5 * (bpts[i] + bpts[j]) end - vpar_mid = _vpar_from_extrap(B_extrap, lmda, bo, θmid) + vpar_mid = _vpar_from_spline(B_vpar, lmda, bo, θmid) if vpar_mid > best_vpar best_t1 = bpts[i] best_t2 = bpts[j] @@ -506,7 +506,7 @@ Ports Fortran torque.F90 lines 674-793. function _bounce_integrate( tdt_pts::Vector{Float64}, tdt_wts::Vector{Float64}, lmda::Float64, lnq::Float64, sigma::Int, n::Int, q::Float64, bo::Float64, - tspl, B_extrap, chi1::Float64, ro::Float64, + tspl, B_vpar, chi1::Float64, ro::Float64, mfac::Vector{Int}, dbob_m_f::Vector{ComplexF64}, divx_m_f::Vector{ComplexF64}, divxfac::Float64, wdfac::Float64, do_matrices::Bool, mpert::Int, @@ -550,9 +550,9 @@ function _bounce_integrate( jac = tspl_f[4] djdpsi = tspl_f[5] - # v_par from the endpoint-fit cubic (consistent with the bounce points); + # v_par from the periodic cubic (consistent with the bounce points); # the periodic tspl B_val remains the numerator field in the integrands. - vpar = 1.0 - (lmda / bo) * B_extrap(θmod) + vpar = 1.0 - (lmda / bo) * B_vpar(θmod) if vpar <= 0 # Negative v_par near a bounce point: same fill rules as the Fortran bounce loop. diff --git a/src/KineticForces/Torque.jl b/src/KineticForces/Torque.jl index f9af86485..e4bbd3c35 100644 --- a/src/KineticForces/Torque.jl +++ b/src/KineticForces/Torque.jl @@ -94,10 +94,12 @@ function tpsi!(tpsi_var::Ref{ComplexF64}, psi::Float64, n::Int, l::Int, # Create periodic interpolant for poloidal quantities tspl = cubic_interp(xs, Series(hcat(B_vals, dBdpsi_vals, dBdtheta_vals, jac_vals, djdpsi_vals)); bc=PeriodicBC()) - # v_par and the bounce points use a separate endpoint-fit (non-periodic) cubic - # of B, like Fortran's vspl. The endpoint fit of 1−(λ/bo)B equals 1−(λ/bo) - # times the fit of B, so one B_extrap per surface serves every λ. - B_extrap = cubic_interp(xs, B_vals; bc=CubicFit()) + # v_par and the bounce points use a separate scalar cubic of B (Fortran's vspl). + # 1−(λ/bo)B is periodic on the closed θ interval, so the fit must be periodic too: + # a non-periodic endpoint fit is only C⁰ at the θ=0/1 seam and manufactures false + # near-seam extrema and root pairs there. The fit of 1−(λ/bo)B equals 1−(λ/bo) + # times the fit of B, so one B_vpar per surface serves every λ. + B_vpar = cubic_interp(xs, B_vals; bc=PeriodicBC()) bmax = maximum(B_vals) ibmax = argmax(B_vals) @@ -226,7 +228,7 @@ function tpsi!(tpsi_var::Ref{ComplexF64}, psi::Float64, n::Int, l::Int, method, op_wmats; chi1=intr.chi1, ro=intr.ro, mfac=intr.mfac, mpert=intr.mpert, theta_bmax=theta_bmax, - B_extrap=B_extrap, + B_vpar=B_vpar, smat=smat_f, tmat=tmat_f, xmat=xmat_f, ymat=ymat_f, zmat=zmat_f, energy_atol=atol_xlmda, energy_rtol=rtol_xlmda, @@ -401,7 +403,7 @@ function calculate_gar(psi, n, l, q, epsr, wdian, wdiat, welec, nuk, bo, bmax, bmin, n_s::Float64, T_s::Float64, mass, chrg, tspl, dbob_m_f, divx_m_f, divxfac, wdfac, method, op_wmats; chi1::Float64, ro::Float64, mfac::Vector{Int}, mpert::Int, - theta_bmax::Float64, B_extrap, + theta_bmax::Float64, B_vpar, smat=nothing, tmat=nothing, xmat=nothing, ymat=nothing, zmat=nothing, nlmda::Int=128, ntheta::Int=128, @@ -415,7 +417,7 @@ function calculate_gar(psi, n, l, q, epsr, wdian, wdiat, welec, nuk, bo, # Bounce-averaged quantities per pitch angle bounce = compute_bounce_data( psi, n, l, q, bo, bmax, bmin, theta_bmax, - tspl, B_extrap, mfac, chi1, ro, dbob_m_f, divx_m_f, divxfac, wdfac, + tspl, B_vpar, mfac, chi1, ro, dbob_m_f, divx_m_f, divxfac, wdfac, mass, chrg, T_s, method; nlmda, ntheta, smat, tmat, xmat, ymat, zmat) @@ -615,8 +617,9 @@ function _setup_surface_state( end tspl = cubic_interp(xs, Series(hcat(B_vals, dBdpsi_vals, dBdtheta_vals, jac_vals, djdpsi_vals)); bc=PeriodicBC()) - # Endpoint-fit (non-periodic) cubic of B for v_par and bounce points (Fortran vspl equivalent). - B_extrap = cubic_interp(xs, B_vals; bc=CubicFit()) + # Periodic cubic of B for v_par and bounce points (Fortran vspl equivalent); 1−(λ/bo)B + # is periodic on the closed θ interval, so a non-periodic fit would break at the seam. + B_vpar = cubic_interp(xs, B_vals; bc=PeriodicBC()) bmax = maximum(B_vals) ibmax = argmax(B_vals) @@ -685,7 +688,7 @@ function _setup_surface_state( return (; chrg, mass, - tspl, B_extrap, bmax, bmin, theta_bmax, + tspl, B_vpar, bmax, bmin, theta_bmax, q, n_s, T_s, welec, wdian, wdiat, wtran, wgyro, nuk, epsr, @@ -750,7 +753,7 @@ function kinetic_energy_matrices_for_euler_lagrange!( bounce = compute_bounce_data( psi, n, l, state.q, bo, state.bmax, state.bmin, state.theta_bmax, - state.tspl, state.B_extrap, mfac, chi1, ro, dbob_m_f, divx_m_f, 1.0, wdfac, + state.tspl, state.B_vpar, mfac, chi1, ro, dbob_m_f, divx_m_f, 1.0, wdfac, state.mass, state.chrg, state.T_s, "fwmm"; nlmda, ntheta, smat=smat_f, tmat=tmat_f, xmat=xmat_f, ymat=ymat_f, zmat=zmat_f) From 53ae808102d7312187c65f0dc7c930b225611639 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Sat, 15 Aug 2026 14:01:44 -0400 Subject: [PATCH 02/14] KineticForces - BUG FIX - drop the spurious major radius from the omega_D prefactor The general-aspect-ratio precession prefactor carried ro^2 where it should carry ro. wbbar = ro*2pi/((2-sigma)*I1) already contains one factor of ro that its own normalization bhat = sqrt(2T/m)/ro cancels; reusing wbbar inside wdbar imports that ro a third time, while dhat = (T/q)/(bo*ro^2) removes only the two written explicitly. Dimensionally, both bounce integrals carry the J*b*dtheta = dl length, so I1 is a length and I2/I1 is 1/Wb. With T/q in volts and V/Wb = 1/s, the prefactor 4*pi*wdfac*(I2/I1)*(T/q) is already a frequency, and the surviving ro left omega_D in m/s. Verified by holding the physics fixed and varying only the machine size: the old form scales as ro (ratio 2.0 when ro is doubled), the corrected form is ro-invariant (ratio 1.0). omega_b, built from the same I1, is untouched and keeps its correct v_th/L scaling. Julia analog of Fortran GPEC PR #281, which measured the same correction against an independently validated guiding-centre precession operator: the least-squares slope of omega_D against the reference moves from -6.49 (= -ro on that ITER equilibrium) to -1.02, with omega_b agreeing to within [0.9964, 1.0011] both before and after as the control. This changes omega_D by a factor of ro on every machine, so it moves the resonance denominator and all NTV torque. Only the magnitude is addressed here; the overall sign convention is a separate question and is not touched. Co-Authored-By: Claude Fable 5 --- src/KineticForces/BounceAveraging.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/KineticForces/BounceAveraging.jl b/src/KineticForces/BounceAveraging.jl index 63a59aedc..372bcdaa8 100644 --- a/src/KineticForces/BounceAveraging.jl +++ b/src/KineticForces/BounceAveraging.jl @@ -632,9 +632,13 @@ function _bounce_integrate( return 0.0, 0.0, 0.0, nothing end - # Bounce-averaged frequencies + # Bounce-averaged frequencies. wbbar already carries one factor of ro that its own + # normalization bhat = sqrt(2T/m)/ro cancels; reusing it inside wdbar imports that ro + # a third time while dhat = (T/q)/(bo·ro²) removes only the two written explicitly, so + # the drift prefactor takes ro, not ro². (Otherwise ω_D = wdbar·dhat carries a surplus + # length: 4π·(I₂/I₁)·(T/q) is already V/Wb = 1/s, so the extra ro leaves m/s.) wbbar = ro * twopi / ((2 - sigma) * total_wb) - wdbar = ro^2 * bo * wdfac * wbbar * 2 * (2 - sigma) * total_wd + wdbar = ro * bo * wdfac * wbbar * 2 * (2 - sigma) * total_wd # Phase factor pl_i = exp(-2πi·lnq·fsi_wb(θ_i)/((2-σ)·total_wb)), using the # cumulative spline integral of the bounce action. From b69e40ce5e17a023285372fb19572e3e3c907c52 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Sat, 15 Aug 2026 14:07:57 -0400 Subject: [PATCH 03/14] KineticForces - IMPROVEMENT - give the energy integration its own tolerances The energy (x) integration was handed the same atol_xlmda/rtol_xlmda pair as the pitch (lambda) integration that encloses it. Since the pitch integrand IS the energy integral, the outer integrator was asked to resolve its integrand to the same tolerance to which that integrand was itself computed, so it chases the inner integrator's quadrature noise instead of converging. Adds atol_x/rtol_x, defaulting (NaN sentinel) to nested_tolerance_margin = 1e-2 times the pitch tolerances, which extends the nesting rule the struct already documents one level up for rtol_psi vs rtol_xlmda. Wired through both the psi-quadrature path (tpsi!) and the kinetic-matrix path. This deliberately changes default numerical behaviour: shipped decks now integrate the energy variable to atol 1e-10 / rtol 1e-7 rather than 1e-8 / 1e-5, and pay for it in runtime. A deck can set atol_x/rtol_x explicitly to override the derived values, or widen nested_tolerance_margin to recover the old cost. Julia analog of the second commit of Fortran GPEC PR #280. Co-Authored-By: Claude Fable 5 --- .../CalculatedKineticMatrices.jl | 4 +++- src/KineticForces/Compute.jl | 4 +++- src/KineticForces/KineticForcesStructs.jl | 19 ++++++++++++----- src/KineticForces/Torque.jl | 21 +++++++++++++++---- 4 files changed, 37 insertions(+), 11 deletions(-) diff --git a/src/KineticForces/CalculatedKineticMatrices.jl b/src/KineticForces/CalculatedKineticMatrices.jl index d36fead99..66e969287 100644 --- a/src/KineticForces/CalculatedKineticMatrices.jl +++ b/src/KineticForces/CalculatedKineticMatrices.jl @@ -121,7 +121,9 @@ function compute_calculated_kinetic_matrices( kf_ctrl.zi, kf_ctrl.mi, kf_ctrl.wdfac, kf_ctrl.divxfac, kf_ctrl.electron, equil, intr_t, kinetic_profiles; nutype=kf_ctrl.nutype, f0type=kf_ctrl.f0type, nufac=kf_ctrl.nufac, - atol_xlmda=kf_ctrl.atol_xlmda, rtol_xlmda=kf_ctrl.rtol_xlmda + atol_xlmda=kf_ctrl.atol_xlmda, rtol_xlmda=kf_ctrl.rtol_xlmda, + atol_x=kf_ctrl.atol_x, rtol_x=kf_ctrl.rtol_x, + nested_tolerance_margin=kf_ctrl.nested_tolerance_margin ) full_w .+= block_w full_t .+= block_t diff --git a/src/KineticForces/Compute.jl b/src/KineticForces/Compute.jl index 4175a62b4..7a316a5cd 100644 --- a/src/KineticForces/Compute.jl +++ b/src/KineticForces/Compute.jl @@ -135,7 +135,9 @@ function integrate_psi_quadgk( 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) + atol_xlmda=ctrl.atol_xlmda, rtol_xlmda=ctrl.rtol_xlmda, + atol_x=ctrl.atol_x, rtol_x=ctrl.rtol_x, + nested_tolerance_margin=ctrl.nested_tolerance_margin) harm_vals[ell_idx] = thread_tpsi[tid][] is_matrix_method && (harm_elems[ell_idx] .= w) end diff --git a/src/KineticForces/KineticForcesStructs.jl b/src/KineticForces/KineticForcesStructs.jl index 6df40e6a8..b153e39ae 100644 --- a/src/KineticForces/KineticForcesStructs.jl +++ b/src/KineticForces/KineticForcesStructs.jl @@ -95,11 +95,20 @@ ctrl = KineticForcesControl(; (Symbol(k) => v for (k, v) in inputs["KineticForce nn::Int = 1 # Toroidal mode number nl::Int = 1 # Bounce harmonic number - # Tolerances. - # *_xlmda: shared tolerances for inner λ (pitch) and x (energy) integrations - # *_psi: tolerances for outer ψ quadrature - atol_xlmda::Float64 = 1e-8 # Absolute tolerance for inner pitch + energy integrations - rtol_xlmda::Float64 = 1e-5 # Relative tolerance for inner pitch + energy integrations + # 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 + 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. + 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). diff --git a/src/KineticForces/Torque.jl b/src/KineticForces/Torque.jl index f9af86485..5fd0c9475 100644 --- a/src/KineticForces/Torque.jl +++ b/src/KineticForces/Torque.jl @@ -35,7 +35,14 @@ function tpsi!(tpsi_var::Ref{ComplexF64}, psi::Float64, n::Int, l::Int, op_wmats::Union{Nothing,Array{ComplexF64,3}}=nothing, rex_override::Union{Nothing,Float64}=nothing, imx_override::Union{Nothing,Float64}=nothing, - atol_xlmda::Float64=1e-9, rtol_xlmda::Float64=1e-6) + atol_xlmda::Float64=1e-9, rtol_xlmda::Float64=1e-6, + atol_x::Float64=NaN, rtol_x::Float64=NaN, + nested_tolerance_margin::Float64=1e-2) + + # The pitch integrand is itself the energy integral, so the energy level is resolved + # tighter than the pitch level (explicit atol_x/rtol_x override the derived values). + atol_energy = isnan(atol_x) ? nested_tolerance_margin * atol_xlmda : atol_x + rtol_energy = isnan(rtol_x) ? nested_tolerance_margin * rtol_xlmda : rtol_x # Enforce bounds if psi > 1 @@ -229,7 +236,7 @@ function tpsi!(tpsi_var::Ref{ComplexF64}, psi::Float64, n::Int, l::Int, B_extrap=B_extrap, smat=smat_f, tmat=tmat_f, xmat=xmat_f, ymat=ymat_f, zmat=zmat_f, - energy_atol=atol_xlmda, energy_rtol=rtol_xlmda, + energy_atol=atol_energy, energy_rtol=rtol_energy, pitch_atol=atol_xlmda, pitch_rtol=rtol_xlmda, rex_override=rex_override, imx_override=imx_override) end @@ -892,7 +899,9 @@ function compute_kinetic_matrices_at_psi!( electron::Bool, equil, intr::KineticForcesInternal, kinetic_profiles::Equilibrium.KineticProfileSplines; nutype::String="harmonic", f0type::String="maxwellian", nufac::Float64=1.0, - atol_xlmda::Float64=1e-9, rtol_xlmda::Float64=1e-6) + atol_xlmda::Float64=1e-9, rtol_xlmda::Float64=1e-6, + atol_x::Float64=NaN, rtol_x::Float64=NaN, + nested_tolerance_margin::Float64=1e-2) # Bypass ψ > 1 (no kinetic contribution outside plasma) if psi > 1 @@ -901,13 +910,17 @@ function compute_kinetic_matrices_at_psi!( return nothing end + # See tpsi!: the energy integral is the pitch integrand, so it is resolved tighter. + atol_energy = isnan(atol_x) ? nested_tolerance_margin * atol_xlmda : atol_x + rtol_energy = isnan(rtol_x) ? nested_tolerance_margin * rtol_xlmda : rtol_x + state = _setup_surface_state(psi, zi, mi, electron, equil, intr, kinetic_profiles) kinetic_energy_matrices_for_euler_lagrange!( kwmat, ktmat, state, psi, n, l, wdfac, intr; nutype, f0type, nufac, - energy_atol=atol_xlmda, energy_rtol=rtol_xlmda, + energy_atol=atol_energy, energy_rtol=rtol_energy, pitch_atol=atol_xlmda, pitch_rtol=rtol_xlmda) return nothing From 09cb856d3eb07464e01f24b505ef4297bc85738b Mon Sep 17 00:00:00 2001 From: logan-nc Date: Thu, 20 Aug 2026 11:47:01 -0400 Subject: [PATCH 04/14] 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 05/14] 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 06/14] 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 07/14] 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 08/14] 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 09/14] 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 10/14] 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 11/14] 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 12/14] 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 13/14] 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 14/14] 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