[logan-nc]

Release note

• Audience: users
• Numerical impact: none (harness @ 9aeee4d)
• Migration: none

KineticForces/<method>/ now writes dTdpsi_ell, the torque density resolved by bounce harmonic ℓ,
on a new ell dimension scale. The values were already computed one harmonic at a time and then
summed; only the sum was kept. Multi-code NTV benchmarking compares the ℓ breakdown, and the Fortran
parent already emits it.

Regression report

Regression Report: solovev_kinetic_ntv
========================================================================================================
Ref 1: develop  @ 84b7e4800 (2026-10-20)
       env: julia 1.11.7, x86_64-linux-gnu, manifest 7027ae20 (pinned), 16 threads/16 BLAS, jit build
Ref 2: feature/kf-per-harmonic-torque  @ 9aeee4d67 (2026-10-20)
       env: julia 1.11.7, x86_64-linux-gnu, manifest 7027ae20 (pinned), 16 threads/16 BLAS, jit build
--------------------------------------------------------------------------------------------------------
Quantity                                   develop       feature/kf-per-harmonic-torque  Diff     Status
--------------------------------------------------------------------------------------------------------
NTV torque fgar [Re, Im]                   [1 elem]      [1 elem]                        0.0e+00  OK    
NTV ψ quadrature evaluations               45            45                              0.0e+00  OK    
root-area-weighted total energy Re(et[1])  1.090570e-02  1.090570e-02                    0.0e+00  OK    
# singular surfaces                        2             2                               0.0e+00  OK    
singular psi locations                     [2 elem]      [2 elem]                        0.0e+00  OK    
q0                                         1.900006e+00  1.900006e+00                    0.0e+00  OK    
Runtime (s)                                229.3s        228.8s                                   --    
========================================================================================================
Summary: 6 unchanged
================================================================
Case: solovev_kinetic_multiion — Solovev n=1, ideal + PE + multi-ion NTV (D + T + electron summed)
================================================================

Regression Report: solovev_kinetic_multiion
===========================================================================================
Ref 1: develop  @ 84b7e4800 (2026-10-20)
       env: julia 1.11.7, x86_64-linux-gnu, manifest 7027ae20 (pinned), 16 threads/16 BLAS, jit build
Ref 2: feature/kf-per-harmonic-torque  @ 9aeee4d67 (2026-10-20)
       env: julia 1.11.7, x86_64-linux-gnu, manifest 7027ae20 (pinned), 16 threads/16 BLAS, jit build
-------------------------------------------------------------------------------------------
Quantity                     develop        feature/kf-per-harmonic-torque  Diff     Status
-------------------------------------------------------------------------------------------
NTV total (D+T+imp+e) [N·m]  2.978404e-01   2.978404e-01                    0.0e+00  OK    
NTV Deuterium [N·m]          1.874428e-01   1.874428e-01                    0.0e+00  OK    
NTV Tritium [N·m]            1.943137e-01   1.943137e-01                    0.0e+00  OK    
NTV electron [N·m]           -8.391612e-02  -8.391612e-02                   0.0e+00  OK    
total energy Re(et[1])       1.090570e-02   1.090570e-02                    0.0e+00  OK    
# singular surfaces          2              2                               0.0e+00  OK    
Runtime (s)                  234.0s         236.5s                                   --    
===========================================================================================
Summary: 6 unchanged

Run after merging origin/develop into the branch (it was 195 commits behind, which would have attributed every intervening develop change to this PR). The merge was clean and the report confirms it: nothing moves against current develop.

Full harness on ITER SDCC (Julia 1.12.6, all 19 cases, same refs 84b7e4800 vs 9aeee4d67): 17 unchanged; gal_resistive_pe N/A on both refs; one changed row, diiid_slayer_n1 SLAYER γ_Hz at 0.02%, which is that case's known run-to-run scatter (this PR does not touch SLAYER).

Review package

https://claude.ai/artifact/4iEvyhb9RQ3S7Xv1gxmVLB

Four figures, so what this adds can be judged without reading the ℓ loop: the decomposition itself
per species and phasing; the sum rule with its 1e-16 residual, which is what makes it a
decomposition rather than a second estimate of the torque; the share of each harmonic's integral
carried by its three largest ψ nodes, which is the honest limit on quoting the high harmonics; and
the per-ℓ comparison across the 195-commit develop merge. All drawn from the feature's own output.
The link is internal and may need access granted — ask @logan-nc, or ask and the page can be sent as
a file instead.

Notes for reviewers

Why

The torque density is already computed harmonic by harmonic — the ℓ loop in integrate_psi_quadgk
fills harm_vals and then reduces it — but only the sum was ever kept.

Multi-code NTV benchmarking wants the breakdown. The ITPA TC-24 activity is comparing the
contributions from ℓ = 0, ±1 … ±5 across GPEC, NTVTOK and MARS-K, and the Fortran parent already
writes the ℓ-resolved profile (dTdpsi_<method> in PENTRC's NetCDF output), so the Julia code could
not answer that request at all.

What

KineticForces/<method>/ gains dTdpsi_ell on a new ell dimension scale. Shape is (npsi, nell)
in Julia, following the repo's (psi, m) axis convention, so a C-ordered reader (h5py, xarray) sees
(nell, npsi). This is not the Fortran orientation: PENTRC defines dTdpsi_<method> as
(i, ell, psi) in Fortran order (pentrc/torque.F90), which those readers see as (psi, ell, i)
with real and imaginary parts split on i. A Python comparison against Fortran output therefore
needs a transpose; the attached psi and ell dimension scales let xarray align the axes by name.
Summing over ℓ reproduces dTdpsi. Per-species groups get it with no extra code.

Three things that needed care

• The per-ℓ log is permuted by the same sortperm as psi. QuadGK visits nodes in refinement
  order, not in ψ order, so an unsorted log would misalign with the grid silently rather than
  erroring.
• harm_vals is reused at every ψ node, so the log stores a copy.
• combine_species_states interpolates onto a union grid, and multi-ion runs need per-ℓ on the
  summed total as well as per species. The interpolation is now a shared helper used by both the
  total and each ℓ column, and the combined dtdpsi is taken as the row sum of the per-ℓ matrix —
  which makes the sum rule exact by construction rather than true to interpolation roundoff. A
  species without per-ℓ data leaves the combined matrix empty instead of erroring.

Caveats

Quadrature adapts on the summed integrand, so a sharply resonant individual harmonic can be less
well resolved than the total. Panel boundaries sit at every ℓ's kinetic-resonance surface, which
limits but does not remove this. Profiles still come from the first toroidal mode only, so for
npert > 1 the sum over ℓ matches dtdpsi but neither integrates to total_torque.

Testing and verification — complete; ready to review and merge

In-tree tests. New assertions in runtests_multiion.jl and runtests_kinetic.jl cover the sum
rule through combine_species_states, the on-disk shape and orientation, the ℓ-sum on disk, the
dimension-scale attachment, and the mixed case where one species lacks per-ℓ data. The on-disk
orientation was checked from both sides: HDF5.jl reads (npsi, nell) and h5py reads (nell, npsi),
with values indexed as intended.

(An earlier version of this description quoted pass counts of 64/64 and 1845/1845. Those were
measured at the original base, before this branch was brought up to develop, so they are no longer
the right numbers to cite — CI below is the live evidence.)

CI. All checks green at the head of this branch, including runtests on both Julia 1.11 and
1.x. This is the load-bearing result for the develop merge described next.

Branch is current with develop. It was 195 commits behind, which would have made a
--refs develop,<branch> harness comparison attribute every intervening develop change to this
PR. Merged (cleanly) before the report above was run, so that report isolates this PR's changes.

The per-ℓ output was verified to survive the merge, on a production case. The sum rule alone
could not establish this: combine_species_states now derives the total as the per-ℓ row sum, so
Σ_ℓ ≡ total holds by construction and cannot fail. Nor could the regression harness — its kinetic
cases do not set bounce_harmonic_max, so they run the default of 1, covering ℓ ∈ {−1,0,+1} and
tracking only summed and per-species totals. So an ITER-scale n=3 case from the ITPA TC-24 activity
was rerun at this branch's head and its ℓ-resolved integrals compared against the same case on the
pre-merge build:

	D	T	e	C
ℓ ordering preserved	yes	yes	yes	yes
all harmonic signs preserved	yes	yes	yes	yes
species total, rel. diff	5.7e-4	7.9e-4	5.1e-5	1.8e-3
ℓ=0, rel. diff	6.3e-4	5.8e-4	5.2e-5	1.9e-4
worst single harmonic	8.0e-3	8.3e-3	6.2e-3	4.3e-3

Ordering and signs are the invariants worth gating here, because they are grid-independent and are
what an ℓ reversal, an index shift or a one-sided ±ℓ symmetrization would break. They hold exactly.

The sub-percent movement comes from develop, not from this PR. The ξ grid that b_n is built on
went from 1642 surfaces to 886, so the fgar quadrature placed its 345 nodes almost entirely
differently, while the field itself agrees to 0.1% at matched ψ_n. That grid is the saved
Euler-Lagrange step record (psi_store), and the drop is the expected effect of #480 (per-column ODE
absolute tolerance, save_interval 3 → 1). Its harness report shows the same ratio on diiid_n1
(2675 → 1440 saved steps, 0.54 against 0.54 here), and its release note documents NTV moving at the
percent level.

That rerun and its comparison script live in the TC-24 benchmark repository, since they need that
activity's ITER equilibrium — so they are not reproducible in-tree, and the numbers are offered as
evidence rather than as a check you can repeat here.

Nothing outstanding from my side. The only thing gating merge is review approval.

---

🤖 Generated with Claude Code

[matt-pharr]

Looks good @logan-nc, I pushed a fix for one small bug that claude found with the envelopes changing between calls of the envelope function, and edited the PR description to accurately reflect the structure of the output from the fortran code.

[matt-pharr]

One note for anyone using dTdpsi_ell for the TC-24 comparison, from checking the per-ℓ logic against PENTRC (pentrc/torque.F90, tintgrnd/tintgrl_lsode):

The per-harmonic quantity, its ℓ ordering (−nl first), and the total (a plain sum over ℓ) match Fortran exactly. Both write tpsi for each ℓ with no per-ℓ factors. The one difference is error control. PENTRC integrates each ℓ as its own LSODE component (neq = 2(2nl+1), per-component rtol/atol), so every harmonic is error-controlled. Julia's QuadGK adapts only on the summed integrand. Panel boundaries at each ℓ's kinetic-resonance surface limit the effect, but a weak or sharply resonant individual harmonic can be less well resolved in Julia than in Fortran even when the totals agree. So per-ℓ disagreements between the codes on high harmonics should be checked against this before being read as physics. A tighter rtol_psi is the knob if a specific harmonic needs it.

Not blocking. It's already in the caveats; this just records where it comes from in the Fortran.