From c8b1f6afa7c4738df659e17d9cbb0b2bfa33a700 Mon Sep 17 00:00:00 2001 From: satelite2517 Date: Wed, 1 Jul 2026 14:33:21 +0900 Subject: [PATCH 1/5] KineticForces - BUG FIX - #269 kinetic dW inflated: don't over-scale xi^alpha by chi1 set_perturbation_data! multiplied clebsch_alpha by chi1 to form the tangential displacement fed to the NTV operator (both xs_m[3] and the JBB deweighting xms). But clebsch_alpha from FieldReconstruction is already the physical xi^alpha = xms/chi1 -- exactly what Fortran gpout_xclebsch writes (xmsout = xmsout/chi1) and what PENTRC consumes directly. The extra *chi1 (~1.65x) over-scaled the tangential displacement, inflating the kinetic dW while leaving the torque comparatively OK. Verified on the DIII-D kinetic case by feeding a matched displacement through the operator (amplitude-independent dW/T; Fortran ref -0.1823): the real-pipeline dW/T moves from -0.117 (36% off) to -0.168 (7.8% off) after the fix. The ~8% residual is grid/accuracy (mpsi, xi differences) plus the separate loglam log10 issue. Also removes the now-obsolete compensating ./chi1 in benchmark_diiid_ideal_ntv_torque.jl, which wrongly assumed gpec_xclebsch stored raw xi^alpha (the file already stores xms/chi1), so its matched test stays correct. Co-Authored-By: Claude Opus 4.8 (1M context) Claude-Session: https://claude.ai/code/session_01RLSHXUFe1LSNX3keRUMWnC --- benchmarks/benchmark_diiid_ideal_ntv_torque.jl | 7 +++---- src/KineticForces/KineticForcesStructs.jl | 9 ++++----- 2 files changed, 7 insertions(+), 9 deletions(-) diff --git a/benchmarks/benchmark_diiid_ideal_ntv_torque.jl b/benchmarks/benchmark_diiid_ideal_ntv_torque.jl index 32aa4c51..80465daa 100644 --- a/benchmarks/benchmark_diiid_ideal_ntv_torque.jl +++ b/benchmarks/benchmark_diiid_ideal_ntv_torque.jl @@ -264,15 +264,14 @@ function run_benchmark(fortran_dir::String=DEFAULT_FORTRAN_DIR) _p("\n--- Build PE state (Fortran Clebsch → dbob_m/divx_m) ---") t1 = time() - # Construct PerturbedEquilibriumState with Fortran Clebsch data - # Fortran xclebsch stores ξ^α directly; PE convention is ξ^α/χ₁ - chi1 = 2π * equil.psio + # Construct PerturbedEquilibriumState with Fortran Clebsch data. + # Fortran gpec_xclebsch already writes ξ^α = xms/χ₁, exactly what the KF operator consumes — feed directly. pe_state = PE.PerturbedEquilibriumState(; psi_grid=psi_grid_f, xi_modes=( clebsch_psi=clebsch_psi, clebsch_psi1=clebsch_psi1, - clebsch_alpha=clebsch_alpha ./ chi1 # Store as ξ^α/χ₁ per PE convention + clebsch_alpha=clebsch_alpha ) ) diff --git a/src/KineticForces/KineticForcesStructs.jl b/src/KineticForces/KineticForcesStructs.jl index 5e235812..070e602d 100644 --- a/src/KineticForces/KineticForcesStructs.jl +++ b/src/KineticForces/KineticForcesStructs.jl @@ -309,16 +309,15 @@ function set_perturbation_data!(kf_intr::KineticForcesInternal, pe_state, ffs_in psi_grid = pe_state.psi_grid npsi = length(psi_grid) mpert = ffs_intr.mpert - chi1 = kf_intr.chi1 # Build xs_m: 3 CubicSeriesInterpolants from Clebsch displacement matrices # xs_m[1] = ξ^ψ (unregularized), xs_m[2] = ∂ξ^ψ/∂ψ (regularized), xs_m[3] = ξ^α - # Note: clebsch_alpha is stored as ξ^α/χ₁, multiply by chi1 to get ξ^α itp_opts = (; extrap=ExtendExtrap()) + # clebsch_alpha is already the physical ξ^α = xms/χ₁ (as Fortran gpout_xclebsch writes it); use directly. + clebsch_alpha_mat = xi_modes.clebsch_alpha xs_m_1 = cubic_interp(psi_grid, Series(xi_modes.clebsch_psi); itp_opts...) xs_m_2 = cubic_interp(psi_grid, Series(xi_modes.clebsch_psi1); itp_opts...) - clebsch_alpha_raw = xi_modes.clebsch_alpha .* chi1 - xs_m_3 = cubic_interp(psi_grid, Series(clebsch_alpha_raw); itp_opts...) + xs_m_3 = cubic_interp(psi_grid, Series(clebsch_alpha_mat); itp_opts...) kf_intr.xs_m = [xs_m_1, xs_m_2, xs_m_3] # Build geometric matrices (S,T,X,Y,Z) for JBB deweighting @@ -355,7 +354,7 @@ function set_perturbation_data!(kf_intr::KineticForcesInternal, pe_state, ffs_in # Get Clebsch displacement vectors at this ψ xsp = view(xi_modes.clebsch_psi, ipsi, :) # ξ^ψ [mpert] xmp1 = view(xi_modes.clebsch_psi1, ipsi, :) # ∂ξ^ψ/∂ψ [mpert] - xms = view(clebsch_alpha_raw, 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) From 0858f63be9ed628cec0ac13323a93fe3cff4d0f2 Mon Sep 17 00:00:00 2001 From: satelite2517 Date: Thu, 2 Jul 2026 14:02:41 +0900 Subject: [PATCH 2/5] KineticForces - BUG FIX - #269 use axis toroidal field for bo, matching Fortran PENTRC intr.bo was equil.params.b0 = bt0, the edge-extrapolated vacuum field at rmean (equil_out.f:268). Fortran PENTRC instead uses the field at the magnetic axis: bo = |sq%f(1)|/(2pi*ro) at psi=0 (dcon_interface.f:1414). On DIII-D 147131 these differ by 2.2%, and since bo normalizes the pitch variable (lambda = mu*bo/E, v_par = 1 - lambda*B/bo) the trapped/passing boundary and every bounce-averaged quantity were shifted by that scale: wb(lambda) disagreed with Fortran by up to several percent mid-passing and ~3.5% trapped. With this fix wb/wd match Fortran to ~1e-4 across the lambda range at psi=0.3/0.5/0.7. DIIID ideal-NTV benchmark vs fresh Fortran reference (with the companion theta-quadrature fix): FGAR dW err 7.23% -> 2.94%, T err 2.66% -> 1.31%; TGAR dW 4.60% -> 4.54%, T 1.71% -> 1.65%. All within the 5% acceptance. Co-Authored-By: Claude Fable 5 --- src/KineticForces/KineticForcesStructs.jl | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/KineticForces/KineticForcesStructs.jl b/src/KineticForces/KineticForcesStructs.jl index 070e602d..15371158 100644 --- a/src/KineticForces/KineticForcesStructs.jl +++ b/src/KineticForces/KineticForcesStructs.jl @@ -254,9 +254,15 @@ the equilibrium geometry parameters needed for NTV calculations. function KineticForcesInternal(equil; verbose::Bool=false) mthsurf = length(equil.rzphi_ys) - 1 nth = mthsurf + 1 + # Fortran PENTRC (dcon_interface.f set_geom): bo = |sq%f(1)|/(2π·ro) at ψ=0 — + # the toroidal field AT THE MAGNETIC AXIS, from F extrapolated to the axis. + # NOT params.b0 (= bt0, the edge-extrapolated vacuum field at rmean, ~2% + # different on DIII-D). bo normalizes λ = μ·bo/E and v_par = 1−λB/bo, so this + # offset shifted the trapped/passing boundary and all bounce integrals (#269). + bo_axis = abs(equil.profiles.F_spline(0.0)) / (2π * equil.ro) KineticForcesInternal(; ro = equil.ro, - bo = equil.params.b0, + bo = bo_axis, chi1 = 2π * equil.psio, mthsurf, tpsi_xs = collect(range(0.0, 1.0, length=nth)), From 8c0721f154e4c00d918051e26eabdc34f8a69472 Mon Sep 17 00:00:00 2001 From: satelite2517 Date: Thu, 2 Jul 2026 14:03:05 +0900 Subject: [PATCH 3/5] KineticForces - BUG FIX - #269 port Fortran cubic-spline theta-quadrature into bounce integrals The bounce integrals (omega_b, omega_D, the pl phase cumulative, the action integral, and the matrix-path W_mu/W_E bounce averages) were trapezoid/Riemann sums. Fortran torque.F90 instead fits the samples on the uniform unit grid with spline_fit(...,"extrap") (C^2 cubic, 4-point Lagrange endpoint derivatives) and integrates the cubic exactly via spline_int. With the 1/sqrt(v_par) integrable singularities at the bounce points, the scheme difference was a leading-order systematic (dJdJ off by up to ~20% at intermediate trapped lambda). - Add _fortran_spline_derivs!/_fortran_spline_integral/_fortran_spline_cumint! porting spline_fac + spline_fit_ahg("extrap") + spline_int on a uniform grid (real and complex), verified against a dense solve to ~1e-15. - Restructure _bounce_integrate to sample arrays with Fortran's exact vpar<=0 zero-fill/hold-fill/backfill semantics (torque.F90:674-735), bug-for-bug including the wb->wd slot copy in the hold branch. - pl phase now uses the cumulative spline integral (bspl%fsi), not a running trapezoid; drop the spurious +eps in the 1/pl division (|pl|=1). - nlmda default 64 -> 128 on the torque path, matching the Fortran pentrc default and the matrix path. - fbnce lambda-interpolant BC ZeroCurvBC -> CubicFit, the analogue of Fortran cspline_fit(fbnce,'extrap'). With the companion bo axis-field fix, bounce-averaged wb/wd match Fortran to ~1e-4 and normalized dJdJ to <1% away from the trapped/ passing boundary; DIIID benchmark FGAR dW error 7.23% -> 2.94%. Co-Authored-By: Claude Fable 5 --- src/KineticForces/BounceAveraging.jl | 227 +++++++++++++++++++-------- src/KineticForces/Torque.jl | 11 +- 2 files changed, 167 insertions(+), 71 deletions(-) diff --git a/src/KineticForces/BounceAveraging.jl b/src/KineticForces/BounceAveraging.jl index 3a75b003..439c83e6 100644 --- a/src/KineticForces/BounceAveraging.jl +++ b/src/KineticForces/BounceAveraging.jl @@ -141,6 +141,94 @@ function _powspace_antideriv(x::Vector{Float64}, pow::Int) end +# ============================================================================ +# Fortran spline-quadrature helpers (uniform grid) +# ============================================================================ +# The Fortran bounce integrals are NOT Riemann/trapezoid sums: torque.F90 fits +# the θ-samples with equil/spline.f `spline_fit(...,"extrap")` (C² cubic with +# endpoint derivatives from a 4-point Lagrange fit) and integrates the fitted +# cubic exactly via `spline_int`: ∫ over interval = h/12·(6(f_i+f_{i+1}) + +# h·(d_i−d_{i+1})). The 1/√v_par integrand near bounce points makes the +# quadrature scheme a leading-order effect (trapezoid gave a systematic ~3.5% +# inflation of ∫J·B/√v_par in the trapped region → wb 3.5% low, dW off; see +# issue #269), so we reproduce the Fortran scheme exactly. + +""" + _fortran_spline_derivs!(d, f, h) + +Nodal first derivatives of the C² cubic spline through `f` on a uniform grid +of spacing `h`, with Fortran `"extrap"` endpoint conditions (derivative of the +cubic Lagrange polynomial through the first/last 4 points). Ports +`equil/spline.f` `spline_fit_ahg`/`spline_fac` specialized to a uniform grid. +Works for real and complex `f`. Requires `length(f) ≥ 4`. +""" +function _fortran_spline_derivs!(d::AbstractVector{T}, f::AbstractVector{T}, h::Float64) where {T} + n = length(f) + @assert length(d) == n && n >= 4 + # Endpoint derivatives: 4-point Lagrange (uniform grid): [-11/6, 3, -3/2, 1/3]/h + d[1] = (-11 * f[1] + 18 * f[2] - 9 * f[3] + 2 * f[4]) / (6h) + d[n] = (11 * f[n] - 18 * f[n-1] + 9 * f[n-2] - 2 * f[n-3]) / (6h) + # Interior: d[i-1] + 4·d[i] + d[i+1] = 3(f[i+1]-f[i-1])/h, i = 2..n-1, + # with known d[1], d[n] moved to the RHS. Thomas algorithm on (1,4,1). + # Sweep storage reuses d[2:n-1] for RHS. + cp = Vector{Float64}(undef, n - 2) # modified super-diagonal + @inbounds begin + rhs2 = 3 * (f[3] - f[1]) / h - d[1] + cp[1] = 1.0 / 4.0 + d[2] = rhs2 / 4.0 + for i in 3:n-1 + rhs = 3 * (f[i+1] - f[i-1]) / h + if i == n - 1 + rhs -= d[n] + end + denom = 4.0 - cp[i-2] + cp[i-1] = 1.0 / denom + d[i] = (rhs - d[i-1]) / denom + end + for i in n-2:-1:2 + d[i] -= cp[i-1] * d[i+1] + end + end + return d +end + +""" + _fortran_spline_integral(f, d, h) → total + +Exact integral of the cubic-Hermite spline with node values `f` and node +derivatives `d` on a uniform grid of spacing `h`. Ports `spline_int`: +per-interval `h/12·(6(f_i+f_{i+1}) + h·(d_i−d_{i+1}))`. +""" +function _fortran_spline_integral(f::AbstractVector{T}, d::AbstractVector{T}, h::Float64) where {T} + total = zero(T) + @inbounds for i in 1:length(f)-1 + total += (h / 12) * (6 * (f[i] + f[i+1]) + h * (d[i] - d[i+1])) + end + return total +end + +""" + _fortran_spline_cumint!(fsi, f, d, h) + +Cumulative nodal integrals (`fsi[i] = ∫₀^{x_i}` of the fitted spline), +matching Fortran `spline_int`'s `fsi` accumulation. `fsi[1] = 0`. +""" +function _fortran_spline_cumint!(fsi::AbstractVector{T}, f::AbstractVector{T}, d::AbstractVector{T}, h::Float64) where {T} + fsi[1] = zero(T) + @inbounds for i in 1:length(f)-1 + fsi[i+1] = fsi[i] + (h / 12) * (6 * (f[i] + f[i+1]) + h * (d[i] - d[i+1])) + end + return fsi +end + +"""Convenience: spline-fit ("extrap") + integrate samples `f` on a uniform grid.""" +function _fortran_spline_quad(f::AbstractVector{T}, h::Float64) where {T} + d = Vector{T}(undef, length(f)) + _fortran_spline_derivs!(d, f, h) + return _fortran_spline_integral(f, d, h) +end + + # ============================================================================ # Core bounce averaging # ============================================================================ @@ -149,7 +237,7 @@ end compute_bounce_data(psi, n, l, q, bo, bmax, bmin, ibmax, theta_bmax, tspl, mfac, chi1, ro, dbob_m_f, divx_m_f, divxfac, wdfac, mass, chrg, T_s, method; - nlmda=64, ntheta=128, + nlmda=128, ntheta=128, smat=nothing, tmat=nothing, xmat=nothing, ymat=nothing, zmat=nothing) → BounceData @@ -181,7 +269,7 @@ Ports Fortran torque.F90 lines 530-816 (GAR branch). - `method`: Method string (first char: f/t/p determines λ range) # Keyword Arguments -- `nlmda`: Number of pitch angle grid points (default 64) +- `nlmda`: Number of pitch angle grid points (default 128, matching Fortran pentrc nlmda) - `ntheta`: Number of poloidal grid points per bounce (default 128) - `smat, tmat, xmat, ymat, zmat`: Geometric matrices (mpert×mpert) for kinetic matrix path """ @@ -194,7 +282,7 @@ function compute_bounce_data( divxfac::Float64, wdfac::Float64, mass::Float64, chrg::Float64, T_s::Float64, method::String; - nlmda::Int=64, ntheta::Int=128, + nlmda::Int=128, ntheta::Int=128, smat::Union{Nothing,Matrix{ComplexF64}}=nothing, tmat::Union{Nothing,Matrix{ComplexF64}}=nothing, xmat::Union{Nothing,Matrix{ComplexF64}}=nothing, @@ -456,14 +544,12 @@ function _bounce_integrate( ntheta = length(tdt_pts) theta0 = tdt_pts[1] - # Cumulative bounce integrals - cum_wb = 0.0 - cum_wd = 0.0 - - # θ-scratch array allocated fresh each call. Pool-based reuse was tried - # (AdaptiveArrayPools) but showed no speedup and introduced severe slowdowns - # at 2+ threads; plain allocations match Fortran baseline behavior. - cum_wb_arr = zeros(Float64, ntheta) + # θ-sample arrays (Fortran bspl%fs/jvtheta; sample i ↔ Fortran node i-1). + # Allocated fresh each call. Pool-based reuse was tried (AdaptiveArrayPools) + # but showed no speedup and introduced severe slowdowns at 2+ threads; + # plain allocations match Fortran baseline behavior. + g_wb = zeros(Float64, ntheta) # bspl%fs(:,1): J·B/√v_par · dθ/dx + g_wd = zeros(Float64, ntheta) # bspl%fs(:,2): drift integrand · dθ/dx # Action integrand jvtheta = zeros(ComplexF64, ntheta) @@ -493,12 +579,23 @@ function _bounce_integrate( vpar = 1.0 - (lmda / bo) * B_val if vpar <= 0 - # Zero crossing near bounce points — handle like Fortran lines 678-697 + # Zero crossing near bounce points — Fortran torque.F90:678-697 fill + # semantics (sample i ↔ Fortran node i-1): if i < ntheta ÷ 2 - # Before midpoint: zero out previous entries + # Before midpoint: restart — zero everything up to and including + # this sample (Fortran bspl%fs(:i-1,:)=0; jvtheta(:i)=0). + fill!(view(g_wb, 1:i), 0.0) + fill!(view(g_wd, 1:i), 0.0) + fill!(view(jvtheta, 1:i), ComplexF64(0.0)) continue else - # After midpoint: clamp remaining entries + # After midpoint: hold the previous sample for the rest of the + # grid. Fortran fills BOTH quantities from fs(i-2,1) — the wb + # integrand — including the wd slot (bspl%fs(i-1:,2)=bspl%fs(i-2,1)); + # reproduced verbatim for parity. + fill!(view(g_wb, i:ntheta), g_wb[i-1]) + fill!(view(g_wd, i:ntheta), g_wb[i-1]) + fill!(view(jvtheta, i:ntheta), jvtheta[i-1]) break end end @@ -506,15 +603,9 @@ function _bounce_integrate( sqrt_vpar = sqrt(vpar) # Bounce integrands (Fortran lines 698-701) - wb_integrand = dt * jac * B_val / sqrt_vpar - wd_integrand = dt * jac * dBdpsi * (1.0 - 1.5 * lmda * B_val / bo) / sqrt_vpar + - dt * djdpsi * B_val * sqrt_vpar - - cum_wb += wb_integrand - cum_wd += wd_integrand - # Trapezoidal cumulative (matches Fortran spline_int semantics on linear fn): - # bspl%fsi(j)/Δx = g_1 + ... + g_{j-1} + g_j/2, so subtract half the current sample. - cum_wb_arr[i] = cum_wb - wb_integrand / 2 + g_wb[i] = dt * jac * B_val / sqrt_vpar + g_wd[i] = dt * jac * dBdpsi * (1.0 - 1.5 * lmda * B_val / bo) / sqrt_vpar + + dt * djdpsi * B_val * sqrt_vpar # Fourier modes at this θ (Fortran lines 702-708) — write into pre-allocated # expm buffer using the ORIGINAL expression order to preserve bit-level parity. @@ -548,15 +639,26 @@ function _bounce_integrate( wen_mt[mi, i] = wen_pre * expm[mi] / (B_val * sqrt_vpar) * phase / (2 * chi1) end end + + # Smooth backfill for points zeroed before a restart (Fortran lines 730-734, + # index ranges preserved verbatim: g fills samples 3..i-1, jv fills 2..i-1). + if i >= 3 && g_wb[i-1] == 0.0 + fill!(view(g_wb, 3:i-1), g_wb[i]) + fill!(view(g_wd, 3:i-1), g_wd[i]) + fill!(view(jvtheta, 2:i-1), jvtheta[i]) + end end - # Total bounce integrals — Fortran splines over bspl%xs = linspace(0,1,ntheta) - # and integrates via spline_int, which is ≈ (1/(ntheta-1)) × Σ f_i. The tdt(2,i) - # weights contain dθ/dx so Σ tdt·f is raw Riemann; divide by (ntheta-1) to get - # the integral over the unit linear space [0,1] that Fortran produces. - nrm = 1.0 / (ntheta - 1) - total_wb = cum_wb * nrm - total_wd = cum_wd * nrm + # Total bounce integrals — Fortran fits bspl%xs = linspace(0,1,ntheta) with + # spline_fit("extrap") and integrates the cubic exactly (spline_int). The + # tdt(2,i) weights contain dθ/dx so the samples live on the unit x-grid. + h = 1.0 / (ntheta - 1) + d_wb = Vector{Float64}(undef, ntheta) + _fortran_spline_derivs!(d_wb, g_wb, h) + fsi_wb = Vector{Float64}(undef, ntheta) + _fortran_spline_cumint!(fsi_wb, g_wb, d_wb, h) + total_wb = fsi_wb[ntheta] + total_wd = _fortran_spline_quad(g_wd, h) if total_wb ≈ 0.0 # Degenerate case — return zeros @@ -567,35 +669,22 @@ function _bounce_integrate( wbbar = ro * twopi / ((2 - sigma) * total_wb) wdbar = ro^2 * bo * wdfac * wbbar * 2 * (2 - sigma) * total_wd - # Phase factor (Fortran line 750). Ratio cum_wb_arr[i]/total_wb is dimensionless — - # (ntheta-1) cancels, no scaling. For do_matrices we keep `pl` as a Vector because - # it is referenced below; otherwise we fuse pl + bjspl → bj_integral in a single - # pass, avoiding two Vector{ComplexF64}(ntheta) allocations. - # Trapezoidal quadrature: boundary samples weighted by 0.5. jvtheta is zero at - # i=1 and i=ntheta (loop above runs 2:ntheta-1) so the boundary terms contribute - # nothing in practice, but writing the weights explicitly keeps the integration - # self-correct if the boundary handling ever changes. + # Phase factor (Fortran line 750): pl_i = exp(-2πi·lnq·h(θ_i)) with + # h = fsi_wb/((2-σ)·total_wb) the cumulative spline integral — NOT a + # running trapezoid — matching bspl%fsi exactly. pl_denom = (2 - sigma) * total_wb one_minus_sigma = 1 - sigma - bj_integral = ComplexF64(0.0) - if do_matrices - pl = Vector{ComplexF64}(undef, ntheta) - @inbounds for i in 1:ntheta - pl[i] = cis(-twopi * lnq * cum_wb_arr[i] * nrm / pl_denom) - end - @inbounds for i in 1:ntheta - w = (i == 1 || i == ntheta) ? 0.5 : 1.0 - bj_integral += w * conj(jvtheta[i]) * (pl[i] + one_minus_sigma / (pl[i] + SINGULAR_EPS)) - end - else - pl = nothing - @inbounds for i in 1:ntheta - pli = cis(-twopi * lnq * cum_wb_arr[i] * nrm / pl_denom) - w = (i == 1 || i == ntheta) ? 0.5 : 1.0 - bj_integral += w * conj(jvtheta[i]) * (pli + one_minus_sigma / (pli + SINGULAR_EPS)) - end + pl = Vector{ComplexF64}(undef, ntheta) + @inbounds for i in 1:ntheta + pl[i] = cis(-twopi * lnq * fsi_wb[i] / pl_denom) + end + + # Action bounce integral (Fortran bjspl: cspline_fit("extrap") + cspline_int). + bj_samples = Vector{ComplexF64}(undef, ntheta) + @inbounds for i in 1:ntheta + bj_samples[i] = conj(jvtheta[i]) * (pl[i] + one_minus_sigma / pl[i]) end - bj_integral *= nrm + bj_integral = _fortran_spline_quad(bj_samples, h) # |δJ|² (Fortran line 756) — division by 2 corrects quadratic form dJdJ_val = wbbar * abs(bj_integral)^2 / 2.0 / ro^2 @@ -603,21 +692,25 @@ function _bounce_integrate( # Matrix path: bounce-average W vectors and form outer products (Fortran lines 759-793) wmats_lmda = nothing if do_matrices - # Bounce-average W_μ and W_E vectors (Fortran lines 762-767). - # Trapezoidal quadrature: boundary samples weighted by 0.5 (wmu_mt and wen_mt - # are zero at i=1 and i=ntheta from the 2:ntheta-1 population loop above). + # Bounce-average W_μ and W_E vectors (Fortran lines 762-767): per mode, + # cspline_fit("extrap") + cspline_int of conj(W_m(θ))·(pl + (1-σ)/pl), + # matching the bjspl quadrature above. wmu_ba = zeros(ComplexF64, mpert) wen_ba = zeros(ComplexF64, mpert) - @inbounds for i in 1:ntheta - w = (i == 1 || i == ntheta) ? 0.5 : 1.0 - factor = w * (pl[i] + one_minus_sigma / (pl[i] + SINGULAR_EPS)) - for mi in 1:mpert - wmu_ba[mi] += conj(wmu_mt[mi, i]) * factor - wen_ba[mi] += conj(wen_mt[mi, i]) * factor + wsamp = Vector{ComplexF64}(undef, ntheta) + wder = Vector{ComplexF64}(undef, ntheta) + @inbounds for mi in 1:mpert + for i in 1:ntheta + wsamp[i] = conj(wmu_mt[mi, i]) * (pl[i] + one_minus_sigma / pl[i]) + end + _fortran_spline_derivs!(wder, wsamp, h) + wmu_ba[mi] = _fortran_spline_integral(wsamp, wder, h) + for i in 1:ntheta + wsamp[i] = conj(wen_mt[mi, i]) * (pl[i] + one_minus_sigma / pl[i]) end + _fortran_spline_derivs!(wder, wsamp, h) + wen_ba[mi] = _fortran_spline_integral(wsamp, wder, h) end - wmu_ba .*= nrm - wen_ba .*= nrm # Reshape as 1×mpert for matrix multiply (Fortran lines 771-772) wmmt = reshape(wmu_ba, 1, mpert) diff --git a/src/KineticForces/Torque.jl b/src/KineticForces/Torque.jl index 35d55a00..06c57e6c 100644 --- a/src/KineticForces/Torque.jl +++ b/src/KineticForces/Torque.jl @@ -422,7 +422,7 @@ function calculate_gar(psi, n, l, q, epsr, wdian, wdiat, welec, nuk, bo, ibmax::Int, theta_bmax::Float64, smat=nothing, tmat=nothing, xmat=nothing, ymat=nothing, zmat=nothing, - nlmda::Int=64, ntheta::Int=128, + nlmda::Int=128, ntheta::Int=128, nutype::String="harmonic", f0type::String="maxwellian", nufac::Float64=1.0, ximag::Float64=0.0, qt::Bool=false, energy_atol::Float64=1e-7, energy_rtol::Float64=1e-5, @@ -477,8 +477,10 @@ function calculate_gar(psi, n, l, q, epsr, wdian, wdiat, welec, nuk, bo, end end - # Build CubicSeriesInterpolant on normalized data (single build) - fbnce = cubic_interp(bounce.lambda, Series(fbnce_data); bc=ZeroCurvBC()) + # Build CubicSeriesInterpolant on normalized data (single build). + # CubicFit endpoint BC = 4-point polynomial fit, matching Fortran + # cspline_fit(fbnce,'extrap') endpoint-derivative extrapolation. + fbnce = cubic_interp(bounce.lambda, Series(fbnce_data); bc=CubicFit()) # 3. Set rex/imx multipliers (Fortran lines 839-847) method_suffix = length(method) >= 4 ? method[2:4] : "" @@ -808,7 +810,8 @@ function kinetic_energy_matrices_for_euler_lagrange!( end end - fbnce = cubic_interp(bounce.lambda, Series(fbnce_data); bc=ZeroCurvBC()) + # CubicFit endpoint BC matches Fortran cspline_fit(fbnce,'extrap'). + fbnce = cubic_interp(bounce.lambda, Series(fbnce_data); bc=CubicFit()) # Dual-output pitch integration: one energy sweep per (λ, E) produces both # halves. Packed return is [wmm | tmm], each length nqty. From 9eea561745839cf4b1c6b7fd6a755450f6f42d27 Mon Sep 17 00:00:00 2001 From: satelite2517 Date: Thu, 2 Jul 2026 14:03:23 +0900 Subject: [PATCH 4/5] KineticForces - MINOR - benchmark: drop stale thmax0 key, tolerate fgar-only references The generated gpec.toml carried thmax0, which ForceFreeStatesControl on this branch does not accept, so the benchmark failed at startup. Also read missing T/dW_total_* reference attributes as NaN so an fgar-only Fortran pentrc reference still benchmarks the fgar method. Co-Authored-By: Claude Fable 5 --- benchmarks/benchmark_diiid_ideal_ntv_torque.jl | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/benchmarks/benchmark_diiid_ideal_ntv_torque.jl b/benchmarks/benchmark_diiid_ideal_ntv_torque.jl index 80465daa..979621fa 100644 --- a/benchmarks/benchmark_diiid_ideal_ntv_torque.jl +++ b/benchmarks/benchmark_diiid_ideal_ntv_torque.jl @@ -77,11 +77,13 @@ Read `T_total_*` and `dW_total_*` scalar reference values from the Fortran function read_pentrc_reference(pentrc_nc::String) NCDatasets.Dataset(pentrc_nc, "r") do ds gatt = ds.attrib + # Reference runs may enable only a subset of methods (e.g. fgar only); missing attributes → NaN. + att(k) = haskey(gatt, k) ? Float64(gatt[k]) : NaN return ( - T_total_fgar=Float64(gatt["T_total_fgar"]), - T_total_tgar=Float64(gatt["T_total_tgar"]), - dW_total_fgar=Float64(gatt["dW_total_fgar"]), - dW_total_tgar=Float64(gatt["dW_total_tgar"]) + T_total_fgar=att("T_total_fgar"), + T_total_tgar=att("T_total_tgar"), + dW_total_fgar=att("dW_total_fgar"), + dW_total_tgar=att("dW_total_tgar") ) end end @@ -142,7 +144,6 @@ nn_high = 1 delta_mlow = 8 delta_mhigh = 8 mthvac = 512 -thmax0 = 1 kinetic_source = "fixed" kinetic_factor = 0.0 From 942ba2a6182d958d4b35a91ba1f217fc396f3c8f Mon Sep 17 00:00:00 2001 From: satelite2517 Date: Thu, 2 Jul 2026 23:14:29 +0900 Subject: [PATCH 5/5] KineticForces - BUG FIX - #269 match Fortran vspl: extrap v_par, knot transit start, complete descending roots Three residual bounce-averaging deviations from Fortran torque.F90, found by adversarial line-level verification after the bo/theta-quadrature fixes: 1. v_par source: Fortran builds a SEPARATE spline vspl = 1-(lambda/bo)*B fit "extrap" and takes both the 1/sqrt(v_par) factor and the bounce-point roots from it (torque.F90:599-602, :677 "more consistent w/ bnce pts than direct from tspl"). Julia computed v_par directly from the periodic tspl - exactly the path that comment avoids. The difference is purely the theta-interp end condition (extrap vs periodic): identical at nodes, divergent off-node near the theta=0/1 seam and at turning points where 1/sqrt(v_par) amplifies it. Fix: build B_extrap = cubic_interp(xs, B_vals; bc=CubicFit()) once per surface and evaluate v_par = 1-(lambda/bo)*B_extrap(theta) everywhere (valid because the extrap spline operator is affine in the nodal data, so the per-lambda vspl equals the affine image of a single lambda-independent B cubic). Numerator fields (B, J, dB/dpsi, dJ/dpsi) still come from the periodic tspl, matching Fortran. 2. Transit start: the passing and no-root branches used the refined dB/dtheta=0 extremum for t1/t2; Fortran uses the nodal-argmax knot tspl%xs(ibmax) (torque.F90:610-611, :648-649) - the refined extremum is diagnostic-only. The extrema loops keep refining the bmax/bmin VALUES (and theta_bmin, which Fortran does use), but no longer move theta_bmax. 3. Bounce-point roots: the fixed 256-point sign-change scan could miss a near-marginal root pair that Fortran's analytic per-interval spline_roots resolves, and returned roots ascending while spline_roots returns them descending (spline.f:1782-1785) - an ordering the deepest-well wrap logic is sensitive to for >=3 roots. Fix: Roots.find_zeros on the same B_extrap cubic, sorted descending; _bisect_vpar deleted. DIIID ideal-NTV benchmark vs fresh Fortran reference (with the companion chi1/bo/theta-quadrature fixes on this branch): FGAR dW err 7.23% -> 3.04%, T 2.66% -> 1.28%; TGAR dW 4.60% -> 4.54%, T 1.71% -> 1.65% - all four within the 5% acceptance (issue #269 reported 16.94%/18.70%). Kinetic-DCON eigenvalue unregressed: et1 = 1.046386-0.302321i vs Fortran 1.050683-0.301242i (0.41%/0.36%). Also tested and ruled out the issue's 083b0cc0 hypothesis: the theta-endpoint snap is a 3-ulp data correction whose measured effect on the interpolants is 3.1e-15 relative. Co-Authored-By: Claude Fable 5 --- src/KineticForces/BounceAveraging.jl | 108 ++++++++++----------------- src/KineticForces/Torque.jl | 23 ++++-- 2 files changed, 55 insertions(+), 76 deletions(-) diff --git a/src/KineticForces/BounceAveraging.jl b/src/KineticForces/BounceAveraging.jl index 439c83e6..1803c056 100644 --- a/src/KineticForces/BounceAveraging.jl +++ b/src/KineticForces/BounceAveraging.jl @@ -235,7 +235,7 @@ end """ compute_bounce_data(psi, n, l, q, bo, bmax, bmin, ibmax, theta_bmax, - tspl, mfac, chi1, ro, dbob_m_f, divx_m_f, + tspl, B_extrap, 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, @@ -277,7 +277,7 @@ function compute_bounce_data( psi::Float64, n::Int, l::Int, q::Float64, bo::Float64, bmax::Float64, bmin::Float64, ibmax::Int, theta_bmax::Float64, - tspl, mfac::Vector{Int}, chi1::Float64, ro::Float64, + tspl, B_extrap, 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, @@ -325,13 +325,13 @@ function compute_bounce_data( # Find bounce points and build θ sub-grid _, _, tdt_pts, tdt_wts = _find_bounce_points_and_grid( - lmda, bo, sigma, tspl, ibmax, theta_bmax, + lmda, bo, sigma, B_extrap, ibmax, theta_bmax, lmdatpb, lmdamax, 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, chi1, ro, mfac, dbob_m_f, divx_m_f, divxfac, wdfac, + tspl, B_extrap, chi1, ro, mfac, dbob_m_f, divx_m_f, divxfac, wdfac, do_matrices, mpert, smat, tmat, xmat, ymat, zmat) # Physical frequencies (Fortran lines 744-745) @@ -389,32 +389,42 @@ function _build_lambda_grid(method_char::Char, lmdatpb::Float64, lmdamax::Float6 end +""" + _vpar_from_extrap(B_extrap, lmda, bo, θ) → v_par + +Parallel-velocity factor `v_par = 1 − (λ/bo)·B(θ)` evaluated from the **extrap** +cubic of B (`B_extrap`), matching Fortran's separate `vspl` (torque.F90:599-600, +`vpar = vspl%f(1)` at :677 — "more consistent w/ bnce pts than direct from tspl"). +Because the "extrap" endpoint derivatives are linear in the nodal data, the extrap +cubic of `1−(λ/bo)B` equals `1−(λ/bo)·(extrap cubic of B)`, so we build `B_extrap` +once per surface (CubicFit ≡ Fortran extrap) and reuse it for the vpar factor, the +bounce-point roots, and the deepest-well test — NOT the periodic `tspl`. +""" +@inline _vpar_from_extrap(B_extrap, lmda::Float64, bo::Float64, θ::Float64) = + 1.0 - (lmda / bo) * B_extrap(mod(θ, 1.0)) + + """ Find bounce points for trapped/passing particles and build θ sub-grid. Returns (t1, t2, theta_points, theta_weights). """ function _find_bounce_points_and_grid( lmda::Float64, bo::Float64, sigma::Int, - tspl, ::Int, theta_bmax::Float64, + B_extrap, ::Int, theta_bmax::Float64, ::Float64, ::Float64, psi::Float64, ntheta::Int ) if sigma == 0 # trapped - # Build v_par(θ) = 1 - (λ/bo)*B(θ) and find roots - # Use a dense θ grid to find zero crossings - nfine = 256 - theta_fine = range(0.0, 1.0, length=nfine+1) - vpar_fine = [1.0 - (lmda / bo) * tspl(mod(θ, 1.0))[1] for θ in theta_fine] - - # Find zero crossings - bpts = Float64[] - for i in 1:nfine - if vpar_fine[i] * vpar_fine[i+1] < 0 - # Bisect for better accuracy - θ_root = _bisect_vpar(tspl, lmda, bo, theta_fine[i], theta_fine[i+1]) - push!(bpts, θ_root) - end - end + # Bounce points = roots of v_par(θ) = 1 − (λ/bo)·B_extrap(θ) in (0,1). + # Fortran fits vspl "extrap" and calls spline_roots (torque.F90:600-602), which + # solves each equilibrium interval's cubic analytically. Roots.find_zeros + # adaptively finds ALL roots of the SAME extrap cubic (B_extrap), so the root + # values match spline_roots (robust to the near-boundary interior dip a fixed + # sign-change scan could miss). Sorted DESCENDING to match spline_roots order + # (spline.f:1782-1785), which the marginally-trapped and deepest-well wrap + # logic below assume. + vpar_fn = θ -> _vpar_from_extrap(B_extrap, lmda, bo, θ) + bpts = sort!(Roots.find_zeros(vpar_fn, 0.0, 1.0); rev=true) nbpts = length(bpts) if nbpts < 1 @@ -427,7 +437,7 @@ function _find_bounce_points_and_grid( t2 = bpts[1] + 1.0 else # Find deepest potential well (Fortran lines 616-639) - t1, t2 = _find_deepest_well(bpts, tspl, lmda, bo) + t1, t2 = _find_deepest_well(bpts, B_extrap, lmda, bo) end # Power-law grid refined near bounce points @@ -443,56 +453,11 @@ function _find_bounce_points_and_grid( end -""" - _bisect_vpar(tspl, lmda, bo, θa, θb; tol=1e-12, maxiter=50) → θ - -Bisect on `θ ∈ [θa, θb]` to find a bounce point — the angle where the parallel -velocity vanishes, `v_par(θ) = 1 - (λ/bo) · B(θ) = 0`. The caller is responsible -for supplying a bracket where `v_par(θa)` and `v_par(θb)` have opposite signs; -this routine does not check, it just halves toward the sign change. - -# Arguments -- `tspl`: 1D θ-interpolant returning at least `[B(θ), …]` at first index; - `tspl(mod(θ, 1.0))[1]` extracts the local field magnitude. -- `lmda`: pitch-angle parameter λ = μ·bo / E. -- `bo`: on-axis toroidal field used to normalise λ. -- `θa, θb`: bracket endpoints (normalised θ ∈ [0,1)); typically straddling a B-peak. - -# Keyword arguments -- `tol`: termination tolerance — exits when either `|v_par(θ_mid)| < tol` or the - bracket width `θb − θa < tol`. Default `1e-12` is tight enough that the - residual `v_par` is dominated by the spline-evaluation roundoff of `tspl`. -- `maxiter`: hard iteration cap. Default `50` halves the bracket by ~10⁻¹⁵, well - past `tol` on a unit-scale bracket; the cap is a runaway guard rather than the - expected exit. On exhaustion the midpoint of the final bracket is returned. - -# Returns -- `θ::Float64`: the converged (or capped) bounce-point angle. -""" -function _bisect_vpar(tspl, lmda::Float64, bo::Float64, θa::Float64, θb::Float64; tol=1e-12, maxiter=50) - va = 1.0 - (lmda / bo) * tspl(mod(θa, 1.0))[1] - for _ in 1:maxiter - θm = 0.5 * (θa + θb) - vm = 1.0 - (lmda / bo) * tspl(mod(θm, 1.0))[1] - if abs(vm) < tol || (θb - θa) < tol - return θm - end - if va * vm < 0 - θb = θm - else - θa = θm - va = vm - end - end - return 0.5 * (θa + θb) -end - - """ Find deepest potential well among bounce point pairs. Ports Fortran lines 616-639. """ -function _find_deepest_well(bpts::Vector{Float64}, tspl, lmda::Float64, bo::Float64) +function _find_deepest_well(bpts::Vector{Float64}, B_extrap, lmda::Float64, bo::Float64) nbpts = length(bpts) best_vpar = 0.0 best_t1 = 0.0 @@ -506,7 +471,7 @@ function _find_deepest_well(bpts::Vector{Float64}, tspl, lmda::Float64, bo::Floa else θmid = 0.5 * (bpts[i] + bpts[j]) end - vpar_mid = 1.0 - (lmda / bo) * tspl(mod(θmid, 1.0))[1] + vpar_mid = _vpar_from_extrap(B_extrap, lmda, bo, θmid) if vpar_mid > best_vpar best_t1 = bpts[i] best_t2 = bpts[j] @@ -535,7 +500,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, chi1::Float64, ro::Float64, + tspl, B_extrap, 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, @@ -576,7 +541,10 @@ function _bounce_integrate( jac = tspl_f[4] djdpsi = tspl_f[5] - vpar = 1.0 - (lmda / bo) * B_val + # vpar from the extrap cubic of B (Fortran vspl, torque.F90:677), NOT the + # periodic tspl B_val — which stays as the numerator field below (Fortran + # bspl uses tspl%f(1)/(4) for the numerator, vspl%f(1) only for 1/√vpar). + vpar = 1.0 - (lmda / bo) * B_extrap(θmod) if vpar <= 0 # Zero crossing near bounce points — Fortran torque.F90:678-697 fill diff --git a/src/KineticForces/Torque.jl b/src/KineticForces/Torque.jl index 06c57e6c..671e7071 100644 --- a/src/KineticForces/Torque.jl +++ b/src/KineticForces/Torque.jl @@ -100,6 +100,10 @@ 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()) + # Extrap cubic of B for the v_par factor + bounce points (Fortran vspl, fit + # "extrap"); CubicFit ≡ Fortran extrap. Kept separate from the periodic tspl + # so v_par matches Fortran torque.F90:599-600,677 exactly (see _vpar_from_extrap). + B_extrap = cubic_interp(xs, B_vals; bc=CubicFit()) bmax = maximum(B_vals) ibmax = argmax(B_vals) @@ -144,7 +148,9 @@ function tpsi!(tpsi_var::Ref{ComplexF64}, psi::Float64, n::Int, l::Int, end if Bθ > bmax bmax = Bθ - theta_bmax = θn + # theta_bmax intentionally NOT refined: Fortran uses the nodal-argmax + # knot xs[ibmax] as the transit start t1/t2 (torque.F90:610-611,648-649); + # the refined extremum (:301) is diagnostic-only. Keep theta_bmax = xs[ibmax]. end end @@ -241,6 +247,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, ibmax=ibmax, theta_bmax=theta_bmax, + 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, @@ -419,7 +426,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, - ibmax::Int, theta_bmax::Float64, + ibmax::Int, theta_bmax::Float64, B_extrap, smat=nothing, tmat=nothing, xmat=nothing, ymat=nothing, zmat=nothing, nlmda::Int=128, ntheta::Int=128, @@ -433,7 +440,7 @@ function calculate_gar(psi, n, l, q, epsr, wdian, wdiat, welec, nuk, bo, # 1. Compute bounce-averaged quantities bounce = compute_bounce_data( psi, n, l, q, bo, bmax, bmin, ibmax, theta_bmax, - tspl, mfac, chi1, ro, dbob_m_f, divx_m_f, divxfac, wdfac, + tspl, B_extrap, mfac, chi1, ro, dbob_m_f, divx_m_f, divxfac, wdfac, mass, chrg, T_s, method; nlmda, ntheta, smat, tmat, xmat, ymat, zmat) @@ -638,6 +645,8 @@ function _setup_surface_state( end tspl = cubic_interp(xs, Series(hcat(B_vals, dBdpsi_vals, dBdtheta_vals, jac_vals, djdpsi_vals)); bc=PeriodicBC()) + # Extrap cubic of B for v_par + bounce points (Fortran vspl "extrap"); see _vpar_from_extrap. + B_extrap = cubic_interp(xs, B_vals; bc=CubicFit()) bmax = maximum(B_vals) ibmax = argmax(B_vals) @@ -671,7 +680,9 @@ function _setup_surface_state( end if Bθ > bmax bmax = Bθ - theta_bmax = θn + # theta_bmax intentionally NOT refined: Fortran uses the nodal-argmax + # knot xs[ibmax] for t1/t2 (torque.F90:610-611,648-649); refined extremum + # is diagnostic-only. Keep theta_bmax = xs[ibmax]. end end end @@ -708,7 +719,7 @@ function _setup_surface_state( return (; chrg, mass, - tspl, bmax, bmin, ibmax, theta_bmax, + tspl, B_extrap, bmax, bmin, ibmax, theta_bmax, q, n_s, T_s, welec, wdian, wdiat, wtran, wgyro, nuk, epsr, @@ -773,7 +784,7 @@ function kinetic_energy_matrices_for_euler_lagrange!( bounce = compute_bounce_data( psi, n, l, state.q, bo, state.bmax, state.bmin, state.ibmax, state.theta_bmax, - state.tspl, mfac, chi1, ro, dbob_m_f, divx_m_f, 1.0, wdfac, + state.tspl, state.B_extrap, 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)