kinetic NTV δW inflated vs Fortran: χ₁ over-scaling + three bounce-averaging bugs#310
kinetic NTV δW inflated vs Fortran: χ₁ over-scaling + three bounce-averaging bugs#310jaesun57 wants to merge 5 commits into
Conversation
…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) <[email protected]> Claude-Session: https://claude.ai/code/session_01RLSHXUFe1LSNX3keRUMWnC
…ng 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 <[email protected]>
…ture 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 <[email protected]>
…ar-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 <[email protected]>
… 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 083b0cc 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 <[email protected]>
There was a problem hiding this comment.
Pull request overview
This PR addresses the remaining kinetic NTV δW discrepancy against the Fortran GPEC/PENTRC reference by correcting multiple issues in the KineticForces bounce-averaging chain (and related normalization choices), bringing the DIII-D benchmark metrics back within the stated acceptance targets.
Changes:
- Aligns key normalizations with Fortran (notably
boon-axis field and ξᵅ scaling expectations) to remove systematic offsets in λ-space and bounce quantities. - Reworks trapped bounce-point detection and θ-integral evaluation to match the Fortran “extrap” spline + exact-spline-integration scheme rather than trapezoid/scan-based approaches.
- Updates the DIII-D ideal NTV torque benchmark harness to tolerate missing Fortran reference attributes and to pass ξᵅ in the expected units.
Reviewed changes
Copilot reviewed 4 out of 4 changed files in this pull request and generated no comments.
| File | Description |
|---|---|
src/KineticForces/Torque.jl |
Threads B_extrap into GAR/bounce paths, increases default nlmda, and aligns λ-interpolant boundary conditions with Fortran-style “extrap”. |
src/KineticForces/KineticForcesStructs.jl |
Uses on-axis bo derived from F_spline(0) and fixes ξᵅ handling so it is not over-scaled. |
src/KineticForces/BounceAveraging.jl |
Implements Fortran-style “extrap” spline derivative/integration helpers; switches bounce roots/v∥ evaluation to the extrap spline and updates quadrature accordingly. |
benchmarks/benchmark_diiid_ideal_ntv_torque.jl |
Makes reference attribute reads robust to missing keys and adjusts ξᵅ handling to match the corrected operator expectations. |
jhalpern30
left a comment
There was a problem hiding this comment.
Thanks for the fix! The physics changes look ok to me, but I'll let @logan-nc give the signoff on that since the kinetics are out of my area of expertise.
My main comments are ironically relating to comments - LLMs have a habit of absolutely going crazy with multi-line comments that pollute the codebase, and this PR introduces a lot of those that need cleaning up. I've started marking up a few of them, but most of the comments in this PR should be looked at with this lens and it can be cut down a lot. In general try to "humanize" the comments a bit more:
- Avoid these long, multi-line comments around a single line that focus on what was changed versus what it does. LLMs have a habit of writing the logic behind what was changed, versus a simple statement describing the actual code that currently exists. The "what was changed" comments will quickly go out of date. Similarly, avoid adding a multi-line comment in place of a removed line of code
- Avoid referencing Fortran and specific lines of Fortran explicitly, which will again go out of date. Focus on describing the physics of what was implemented, and if helpful just add a "performs the same calculation as XXX in the Fortran" if that would be helpful
Feel free to ping me again once you've done a pass through, happy to take a look at the whole PR then.
| # 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. | ||
|
|
||
| """ |
There was a problem hiding this comment.
Are these really needed? We migrated from the Fortran spline package. If possible, use FastInterpolations API whenever available.
| # 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 |
There was a problem hiding this comment.
Avoid this many references to Fortran - this is a new code. Ok to provide enough for an interested party to reference the original fortran code, but summarizing like this pollute this codebase. This comment could be a line or two if we really need the spline wrappers below
| # 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 |
There was a problem hiding this comment.
This whole part of the comment can be cut off. Specific error values will quickly go out of date as well as references to github issues
| """ | ||
| 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, |
There was a problem hiding this comment.
B_extrap needs to be described in the docstring
| """ | ||
| _vpar_from_extrap(B_extrap, lmda, bo, θ) → v_par | ||
|
|
||
| Parallel-velocity factor `v_par = 1 − (λ/bo)·B(θ)` evaluated from the **extrap** |
There was a problem hiding this comment.
Way too long of a docstring for a short inline function - I think this is more suitable for where B_extrap is created/in the main docstring of tpsi!/something like that. If its this essential of a detail, shouldn't be hidden in the docstring of an inline
| push!(bpts, θ_root) | ||
| end | ||
| end | ||
| # Bounce points = roots of v_par(θ) = 1 − (λ/bo)·B_extrap(θ) in (0,1). |
There was a problem hiding this comment.
Explanation not needed. Something like "# Build v_par(θ) = 1 - (λ/bo)*B(θ) and find roots using extrapolated B field" would be sufficient
| vpar = 1.0 - (lmda / bo) * B_extrap(θmod) | ||
|
|
||
| if vpar <= 0 | ||
| # Zero crossing near bounce points — handle like Fortran lines 678-697 |
There was a problem hiding this comment.
The comments for this entire section could probably have stayed the same
Summary
Closes #269. The kinetic δW drift against Fortran GPEC (dW_fgar 16.94 %, dW_tgar 18.70 % in the
issue) is resolved: all four DIIID PENTRC benchmark metrics are within the 5 % acceptance bar
at production tolerances (converged dW_tgar sits at 5.49 %, the known t/p-boundary residual —
see the tolerance table below), and the kinetic-DCON eigenvalue path is verified unregressed at
0.4 %. The issue's leading
hypothesis (develop
083b0cc0, periodic-θ endpoint snap) is tested and ruled out; the realroot cause was a set of independent bugs in the KineticForces bounce-averaging chain.
Four bugs fixed on this branch:
set_perturbation_data!multiplied ξ^α by χ₁ althoughclebsch_alphais already the physical ξ^α (= Fortrangpec_xclebsch, which storesxmsout/chi1, gpout.f:5643)c8b1f6afbowasequil.params.b0= bt0, the edge-extrapolated vacuum field at rmean (equil_out.f:268). Fortran PENTRC uses the axis field|F(0)|/ro(dcon_interface.f:1414). Since λ = μ·bo/E, the trapped/passing boundary and every bounce average were shifted0858f63bspline_fit(...,"extrap")(C² cubic, 4-pt Lagrange endpoints) and integrates exactly viaspline_int. With the 1/√v_par integrable endpoint singularities this is a leading-order scheme difference8c0721f1vsplfit "extrap" and takes bothvpar=vspl%f(1)andspline_roots(vspl)from it (torque.F90:599-602, 677 — "more consistent w/ bnce pts than direct from tspl"). Also: transit start t1/t2 must be the nodal knotxs[ibmax](torque.F90:610-611, 648-649), not the refined dB/dθ=0 extremum; and roots must be complete + descending likespline_rootsFix #4 exploits the identity
extrap-cubic(1−(λ/bo)B) = 1−(λ/bo)·extrap-cubic(B)(extrap endpointderivatives are linear in the nodal data), so a single λ-independent
B_extrapper surfacereplaces Fortran's per-λ
vsplexactly. Root-finding isRoots.find_zeroson the same cubic(replaces a fixed 256-point sign-change scan that could miss near-marginal root pairs), sorted
descending to match
spline_roots(spline.f:1782-1785).Validation — three-way (Fortran / Julia-before / Julia-after)
Fresh Fortran reference generated with the current binaries (v1.5.5-512) on
docs/examples/DIIID_ideal_example(dcon → gpec → pentrc) andDIIID_kinetic_example(kin_flag=t dcon). "Before" = branch baseline
c8b1f6af(χ₁ fix in, bounce code original);"after" = this PR. Both Julia runs consume the same Fortran
gpec_xclebsch_n1.outξ, so thecomparison isolates the KF operator.
NTV torque and δW (DIIID, n=1, nl=4, fgar+tgar)
The after-fix dW matches PR #224's pre-merge quality (dW_fgar 2.68 %, dW_tgar 4.81 %). Underlying
bounce quantities: wb/wd match Fortran to ~1e-4 across λ at ψ = 0.3/0.5/0.7; normalized dJdJ < 1 %
away from the trapped/passing boundary. Remaining tgar ~4.5 % is the t/p-boundary dip node +
λ-integration scheme (LSODE vs QuadGK), where both codes are marginal.
Kinetic-DCON eigenvalue (kin_flag=t path, kinetic_source="calculated")
The bounce fixes do not regress the self-consistent eigenvalue path (before ≈ after to ~3e-5):
et is a ψ-integrated quantity and is insensitive to the local quadrature details that dominate the
NTV δW. Both clear the 5 % (Re) / 20 % (Im) acceptance with large margin.
ξ^α (χ₁ fix)
Before the fix the tangential displacement fed to the kinetic operator was
ξ^α · χ₁(χ₁ = 2π·ψ₀ = 1.648 on this equilibrium); after the fix it is identical to the Fortran
gpec_xclebschcolumn. The old NTV benchmark masked this with a compensating./chi1(removed in
c8b1f6af) — which is also why the bug never showed in CI (coverage gap, see below).Figures (attached)
Three-way overlays (Fortran black dashed / before orange / after crimson):
ntv_fgar_profiles.png,ntv_tgar_profiles.png— dT/dψ, cumulative T, dδW/dψ, cumulative δWkinetic_dcon_et.png— Re/Im(et₁)xi_alpha_chi1.png— Σ_m |ξ^α|(ψ): before = ×1.648 over-scale, after ≡ FortranReproducer scripts (three-way dump drivers + plotter) are kept out of this PR to keep the diff
src-only; available on the local validation setup (or on request).
The 083b0cc hypothesis — tested and ruled out
The issue's leading hypothesis blamed develop
083b0cc0("Snap periodic theta endpoint at 2Dinterpolant construction sites") for the dW drift. Three independent findings rule it out:
the strict FastInterpolations
PeriodicBCvalidator rejects the unsnapped array withy[1]=7.188034825477611vsy[end]=7.188034825477608: a relative drift of 4e-16.snapped each way (the full ambiguity the commit resolves) changes the interpolant by
3.1e-15 relative (dense-grid max). A 12–13 % dW shift would require ~4×10¹³ amplification.
083b0cc0remains in place, untouched.The issue's intuition that Im is more boundary-sensitive than Re was correct — but the sensitive
boundary is the bounce-integrand quadrature (1/√v_par endpoints), not the equilibrium
interpolant θ-endpoint.
Relation to PR #224's open items
"Matrix Path Gap" (5–50 % rel_frob on A_k…H_k vs Fortran, smooth ψ-dependent scaling
‖F‖/‖J‖ ≈ 1.01–1.35, no m-structure, not quadrature-order): the bo scale error (2.2 % on λ),
the trapezoid θ-quadrature systematic, and nlmda 64→128 are exactly this signature. The 0.4 %
et₁ agreement after the fixes suggests the plateau is largely closed; re-running
compare_kinetic_matrices.jl(validation-proofs branch) is a cheap follow-up confirmation.Tolerance sensitivity: KineticForces — Complete PE→KF pipeline and validate DIIID-PENTRC benchmark #224 documented dW error exploding as the outer-ψ tolerance loosens
past 1e-5 (2.7 → 12.8 → 53.8 → 73.7 %). Current
KineticForcesControldefaults areatol_psi=rtol_psi=1e-2. Convergence check of the after-fix numbers (ψ-steps 435 → 3165):No tolerance explosion occurs with the fixed code — dW moves < 1 pp between 1e-2 and 1e-5,
so the KineticForces — Complete PE→KF pipeline and validate DIIID-PENTRC benchmark #224 blow-up mechanism is gone (it was fed by the bounce-side errors). The converged
dW_tgar sits marginally above the 5 % bar (5.49 %), which is the known t/p-boundary residual
tracked in follow-up 3, not an outer-quadrature artifact.
Follow-ups (not blocking)
kinetic_forces/fgar/{total_torque,total_energy}through thepost-PE path — the coverage gap that let KineticForces — dW_kinetic (imaginary NTV energy) ~17–19% off Fortran on DIIID benchmark #269 ship.
atol_psi/rtol_psidefaults vs KineticForces — Complete PE→KF pipeline and validate DIIID-PENTRC benchmark #224's 1e-5 floor.🤖 Generated with Claude Code