[logan-nc]

Release note

• Audience: users
• Numerical impact: self-consistent kinetic runs (kinetic_source = "calculated") on up-down-asymmetric equilibria move: the kinetic Euler-Lagrange matrices change, and with them the kinetic eigenfunctions, plasma response, δW and torques. Up-down-symmetric equilibria are unaffected up to round-off; the two tracked Solovev kinetic cases move by ~1e-7 relative in the energies and by 0.02–0.16 % in fgar and the ξ/b^θ norms after one adaptive ODE step flips. The solovev_kinetic_calculated deck also turns on compute_torque_response, so three torque routes are now tracked there (harness @ 7333488)
• Migration: none

A forward kinetic solve (integrator = "forward", kinetic_factor > 0) can now write the ψ-resolved torque response matrices of Fortran GPEC's gpout_dw_matrix: T_xe(ψ, m, m′) on the applied root-area-weighted control-surface field and T_source(ψ) in forcing-source space (one column per coil set, or the forcing file as a single source), under PerturbedEquilibrium/TorqueResponse/, gated by the new [PerturbedEquilibrium] compute_torque_response control. The run's own forcing is contracted into a cumulative torque profile, and torque_profile does the same for any applied spectrum.

This PR also fixes a conjugation error in the bounce-averaged W vectors that build the kinetic Euler-Lagrange matrices. Because of it, self-consistent kinetic runs on up-down-asymmetric equilibria responded to the up-down-mirrored perturbation. With the fix, the torque from T_xe and the KineticForces fgar quadrature agree to 0.01–0.05 % on converged grids. Fortran GPEC has the same error, fixed in PrincetonUniversity/GPEC#357.

Regression report

Full harness, all 19 cases, develop @ 1cb9bf0 vs this branch @ 7333488 (develop merged in; SDCC regress_fanout.sh, harness and case files from this branch, Julia 1.12.6, CI manifest pin):

Result	Cases
Unchanged	16 cases, 229 quantities. Only nonzero diffs: diiid_slayer_n1 SLAYER roots (Q ≤ 3.8e-07, ω/γ ≤ 0.014 Hz), below its noise thresholds, OK
Changed	solovev_kinetic_calculated (12 of 19 quantities) and solovev_kinetic_nuzero (5 of 14), tables below
New quantity	solovev_kinetic_calculated: T_xe applied-forcing torque (N/A on develop)
N/A on both refs	gal_resistive_pe, all 8 quantities; same on develop, not from this PR

Both changed cases are Solovev, which is up-down symmetric. There the conjugation fix is an identity in exact arithmetic, and the energies move only at ~1e-7 relative. In solovev_kinetic_calculated that is enough to flip one adaptive ODE step (606 → 607 saved), which accounts for the larger moves in the ψ-resolved ξ and b^θ norms (0.11–0.16 %) and in fgar (0.02 %). The new TorqueResponse/total_torque (2.285706e-03 N·m) equals the tracked boundary-response torque Energies/toroidal_torque to every printed digit. On this deck it is 0.8 % above fgar, as expected at mpsi = 16 (see the convergence notes below).

Regression Report: solovev_kinetic_calculated
=====================================================================================================================
Ref 1: develop  @ 1cb9bf099 (2026-10-25)
       env: julia 1.12.6, x86_64-linux-gnu, manifest 4c9d54fa (pinned), 16 threads/8 BLAS, jit build
Ref 2: feature/pe-torque-response-matrix  @ 733348854 (2026-10-25)
       env: julia 1.12.6, x86_64-linux-gnu, manifest 4c9d54fa (pinned), 16 threads/8 BLAS, jit build
---------------------------------------------------------------------------------------------------------------------
Quantity                           develop        feature/pe-torque-response-matrix  Diff               Status
---------------------------------------------------------------------------------------------------------------------
vacuum energy Re(ev[1])            1.054773e+01   1.054774e+01                       7.595e-07 (0.00%)  ** CHANGED **
q0                                 1.900003e+00   1.900003e+00                       0.0e+00            OK
NTV fgar quadrature torque [N·m]   2.267374e-03   2.266826e-03                       5.487e-07 (0.02%)  ** CHANGED **
||xi_clebsch_alpha||               6.040373e-01   6.033791e-01                       6.582e-04 (0.11%)  ** CHANGED **
total energy Re(et[1])             1.828343e+00   1.828343e+00                       5.759e-07 (0.00%)  ** CHANGED **
total energy Im(et[1])             -1.326721e+00  -1.326722e+00                      3.332e-07 (0.00%)  ** CHANGED **
plasma energy Re(ep[1])            -8.719392e+00  -8.719393e+00                      1.335e-06 (0.00%)  ** CHANGED **
ODE steps (saved)                  606            607                                1.000e+00 (0.17%)  ** CHANGED **
# singular surfaces                2              2                                  0.0e+00            OK
boundary-response torque [N·m]     2.285706e-03   2.285706e-03                       5.832e-10 (0.00%)  ** CHANGED **
singular psi locations             [2 elem]       [2 elem]                           0.0e+00            OK
total energy (all)                 [32 elem]      [32 elem]                          8.922e-06 (0.00%)  ** CHANGED **
ODE steps (total)                  605            606                                1.000e+00 (0.17%)  ** CHANGED **
mpert                              32             32                                 0.0e+00            OK
q95                                3.147422e+00   3.147422e+00                       0.0e+00            OK
||Jb_theta_reg||                   3.939787e-03   3.946156e-03                       6.369e-06 (0.16%)  ** CHANGED **
Runtime (s)                        199.9s         199.0s                                                --
T_xe applied-forcing torque [N·m]  N/A            2.285706e-03                       N/A                N/A
||dxi_clebsch_psi/dpsi||           3.790879e+00   3.786745e+00                       4.134e-03 (0.11%)  ** CHANGED **
singular q values                  [2 elem]       [2 elem]                           0.0e+00            OK
npert                              1              1                                  0.0e+00            OK
=====================================================================================================================
Summary: 12 changed, 7 unchanged, 1 missing/N/A

Regression Report: solovev_kinetic_nuzero
===========================================================================================================
Ref 1: develop  @ 1cb9bf099 (2026-10-25)
       env: julia 1.12.6, x86_64-linux-gnu, manifest 4c9d54fa (pinned), 16 threads/8 BLAS, jit build
Ref 2: feature/pe-torque-response-matrix  @ 733348854 (2026-10-25)
       env: julia 1.12.6, x86_64-linux-gnu, manifest 4c9d54fa (pinned), 16 threads/8 BLAS, jit build
-----------------------------------------------------------------------------------------------------------
Quantity                 develop        feature/pe-torque-response-matrix  Diff               Status
-----------------------------------------------------------------------------------------------------------
vacuum energy Re(ev[1])  1.058724e+01   1.058724e+01                       1.280e-07 (0.00%)  ** CHANGED **
q0                       1.900003e+00   1.900003e+00                       0.0e+00            OK
total energy Re(et[1])   2.290787e+00   2.290786e+00                       4.112e-07 (0.00%)  ** CHANGED **
total energy Im(et[1])   -1.564787e+00  -1.564787e+00                      5.225e-08 (0.00%)  ** CHANGED **
plasma energy Re(ep[1])  -8.296450e+00  -8.296451e+00                      5.392e-07 (0.00%)  ** CHANGED **
ODE steps (saved)        608            608                                0.0e+00            OK
# singular surfaces      2              2                                  0.0e+00            OK
singular psi locations   [2 elem]       [2 elem]                           0.0e+00            OK
total energy (all)       [32 elem]      [32 elem]                          3.863e-06 (0.00%)  ** CHANGED **
ODE steps (total)        607            607                                0.0e+00            OK
mpert                    32             32                                 0.0e+00            OK
q95                      3.147422e+00   3.147422e+00                       0.0e+00            OK
Runtime (s)              201.1s         198.9s                                                --
singular q values        [2 elem]       [2 elem]                           0.0e+00            OK
npert                    1              1                                  0.0e+00            OK
===========================================================================================================
Summary: 5 changed, 9 unchanged

No tracked case is up-down asymmetric with calculated kinetic matrices, so the harness cannot show the fix's real effect. The evidence for that is in the notes below.

Notes for reviewers

Construction. Per stored node k, T_xe(ψ_k) = K†·[U₁(ψ_k)†U₂(ψ_k)]·K · 2n·i/(2μ₀) with K = U₁ₗ⁻¹·D·P·R (D = diag(1/(χ₁(m − n·q_lim)·2πi)), P the flux-space permeability, R = S·A the b̃ → Φ conform). This is the Fortran dwks → bsurfmat → gind → gres → ptof chain collapsed to one edge inverse, which the globally consistent u_store basis allows; a scratch check against the literal Fortran route (U₂U₁⁻¹, U₁/U₁ₗ, per-step) agreed to 2e-16 at a mid-ψ node. T_source = Ks†·[U₁†U₂]·Ks·2n·i/(2μ₀) with Ks = U₁ₗ⁻¹·D·P·M.

Decisions the Fortran left open, and what was chosen:

• Factor of 2. T_xe matches the Fortran array (no /2; "general (2×) inductance"). The physical torque is b̃†·T_xe·b̃/2, applied in torque_profile and in the stored T_applied / total_torque. With this convention real(T_applied[end]) equals Energies/toroidal_torque exactly (both reduce to −n·Im(ξ†U₂U₁⁻¹ξ)/(2μ₀); the vacuum part of wt0 is Hermitian and drops out). Note the OMFIT consumer ntv_limits_of_efc.py omits the /2 and therefore reports twice the torque; pypec/post.py optimize_torque keeps it.
• Basis. T_xe acts on b̃ (forcing_b_rootarea, the basis of the stored response matrices), not on unit-norm Φ_x — R†·gres·R is exactly Fortran's ptof†·gres·ptof since ptof = Σ√A = R.
• Forcing sources. Every forcing source is one column of M. With coil forcing, M[:, c] is each coil set's unit-norm Φ_x as built, at its deck currents, on the run's own boundary grid (Fortran coilmn from field_bs_psi carries coil%cur too), so Σ_c M[:, c] reproduces the run's forcing to 1e-10 and sum(T_source)/2 == T_applied at every node. Any other forcing is a single source: the run's applied Φ_x, named after its forcing file (forcing for an RMPFieldSum), so its 1×1 T_source/2 == T_applied. Driving the sources at scale factors s gives s†·T_source·s/2. Reading several forcing files as separate sources is left for later.
• Imaginary part. imag(b̃†T_xe b̃/2) = 2n·δW of the whole Euler-Lagrange energy (ideal + kinetic), not the KineticForces total_energy (kinetic only) — named and annotated accordingly.
• Single n only. Errors on nn_low != nn_high; the 2n prefactor is per harmonic.

Self-consistency test (test/runtests_torque_response.jl, Solovev calculated-kinetic fixture with two window-pane coil sets; numbers from a3301e5, before the conjugation fix): the run's forcing referred to the solution basis directly from u_store gives T = 1.1060988185e-3 + 2.6342971530e-2 i N·m; T_applied[end] = 1.1060988185e-3 + 2.6342971530e-2 i (agrees to 1e-15), Energies/toroidal_torque = 1.1060987849e-3 (3e-8 relative, the Λ-route algebra). Also checked: the same contraction at an interior node, Hermitian/anti-Hermitian split, T_source vs M̃†T_xe M̃, the single-source file case, shapes vs Info/mn_index, and the HDF5 metadata contract on the new group.

Benchmark against the NTV quadrature: the 13 % gap is resolved. This body used to report T_xe/fgar = 1.130 on examples/DIIID-like_error_field_example at mpsi = 64 and left the cause open. It had two causes, and neither is in this PR's construction.

1. Unconverged grid. At mpsi = 64 neither route is converged. fgar alone moves 21 % between mpsi 64 and 128, and the gap changes sign. On every deck tried, T_xe/fgar scatters by about ±4 % below mpsi ≈ 160. It is meaningful from mpsi ≈ 192 and converged by 576–1024. Compare total torques only at converged mpsi, or compare the ψ-resolved densities.
2. A conjugation error in the kinetic matrices, fixed here in 92e2c95. On converged grids a ~2 % total gap remained. It sat in a 13 % local excess of T_xe over fgar across ψ ∈ [0.905, 0.935], around the kinetic resonance at ψ ≈ 0.92, and no resolution or tolerance setting moved it. The bounce-averaged W vectors that build the six kinetic Euler-Lagrange matrices were formed as w_m = ∫ conj(f_m)·pl, which pairs the D_m amplitudes with conjugated phases. The matrices therefore encoded the drift-kinetic action of the up-down-mirrored perturbation, while fgar's scalar |∫ jvtheta·conj(pl)| was already physical. Conjugating the bounce averages (wmu_ba, wen_ba in src/KineticForces/BounceAveraging.jl) makes the matrix route's complex δJ equal fgar's for each pitch angle and bounce harmonic: δJ / conj(bj) has a median of 0.999987 + 0.000001i over 4992 samples, where the old pairing scattered across the complex plane. Up-down-symmetric equilibria such as Solovev are immune, which is why no regression case caught it.

DIII-D-like C-coil deck, calculated kinetic matrices, forward integrator (scratch decks on SDCC, not committed). It differs from the deck above in two ways: it uses the example's own kinetic profiles, and the q = 5.2 truncation puts ψ_lim at 0.978.

grid	build	T_xe [N·m]	fgar [N·m]	T_xe / fgar	same ratio on ψ ∈ [0.905, 0.935]
mpsi = 256	before fix	1.7828e-4	1.7512e-4	1.0180	1.128
mpsi = 256	with fix	1.7148e-4	1.7155e-4	0.9996	0.995
grid_type = "auto"	before fix	1.7856e-4	1.7562e-4	1.0168	1.128
grid_type = "auto"	with fix	1.7160e-4	1.7185e-4	0.9985	0.9985

Both torques move with the fix, because the self-consistent ξ changes too.

Two fixed-boundary TokaMaker families test whether the gap tracks up-down asymmetry (R₀ = 1.7 m, a = 0.6 m, κ = 1.7, mean δ = 0.35, β_N = 1.5). One varies δ_upper − δ_lower from 0 to 0.5 and the other varies κ_upper − κ_lower, with 11 shapes each. Every shape ran at 26 grids, mpsi 128–1024, with and without the fix: 1144 runs, none failed. At converged grids (mean over mpsi 576–1024) the gap before the fix grows with the asymmetry, to between −0.57 % and +1.09 % across the Δδ family and to −0.28 % across the Δκ family. With the fix it is 0.01–0.05 % for every shape. The symmetric shapes give the same result with and without the fix. The Δδ plots are in the comment thread above.

Fortran shows the same thing (PrincetonUniversity/GPEC#357). On its DIII-D-like run_kinetic_example at mpsi = 256, fgar ÷ the self-consistent torque goes from 0.995 to 1.0016 with the fix. On its Solovev example the output is identical with and without the fix, to every printed digit.

Caveat observed on the way: with the example's default psihigh=0.995 (domain ending at q=6.2, just past the q=6 rational; cond(U₁) ~ 1e7–1e9 over the last nodes, tol 1e-8) the identity broke — T_xe route −110 N·m vs boundary route 1.95 N·m vs fgar 1739 N·m, and U₂U₁⁻¹ at the last stored node differed from the free-boundary W_plasma by 300 % and from the previous node by 116 %. The construction is only as good as the edge solution it is built from; when Energies/toroidal_torque and TorqueResponse/total_torque disagree, the solve's edge is ill-conditioned. This is a property of the forward solve near a rational, not of this PR, but worth knowing before trusting either torque on such a deck.

Also in passing: the stale kinetic_source docstring in ForceFreeStatesControl ("calculated (PENTRC — not yet implemented)") now describes the implemented path.

---

🤖 Generated with Claude Code

https://claude.ai/code/session_01BzKDmtKYDotJWrcxWf1oJu

[logan-nc]

Not urgent — @matt-pharr, this is a self-contained addition (new output only, gated off by default) and can wait behind your current work. The one thing worth a skeptical eye when you get to it is the 13 % gap between the T_xe-contracted torque and the fgar quadrature, stated as an open question in the body.

[matt-pharr]

@logan-nc I have been in the process of investigating the T_xe vs. fgar difference you mentioned for the past ~week, and I think it's a bug in the kinetic-matrix bounce averaging (pentrc/torque.F90, the wmmt/wemt bounce-averaged W vectors), not something from this PR. Claude thinks "the W vectors are conjugated inconsistently with the scalar |δJ|² route, so the self-consistent kinetic matrices effectively see the up-down-mirrored perturbation." The missing conjugation also exists in fortran GPEC and only affects non-up-down-symmetric equilibria, so Solovev doesn't show it. I am doing a bunch of fortran and julia runs both at the moment. It doesn't show up in the d3d example because you get lucky with a spurious cancellation that happens in 147131 specifically. You can hold 1.6.0 on fortran on this if you want but otherwise we can make it the first hotfix. I think it will still take another day or two to finish hunting down just because the runs to check it at higher and higher resolution take a while and I don't want to bog down feynman (especially since I graduated :^) )

[logan-nc]

Thanks @matt-pharr ! This would be a great bug hunt-down if it works out.

Don't you have big new fancy ITER computers you can bog down instead of our little old feynman thats doing all it can before APS?

[matt-pharr]

I am currently running a convergence scan on 20 nodes on ITER's SDCC cluster, which has pretty low usage on a saturday night :) I am going to be running all of my harness/GPEC stuff there from now on I think, so long as everything is nice and effectively parallelized!

[matt-pharr]

@logan-nc

Based on this scan using different upper and lower deltas, can confirm that the torque gap between FGAR and T_xe is only present for up-down-asymmetric equilibria, and is fixed by the aforementioned complex conjugate. The difference in torques between the two paths grows the greater a difference in upper and lower triangularity the equilibrium has. This bug is present in fortran as well, and I will open a very small PR there to fix it before 1.6. In the second figure you can see that with the fix we converge cleanly to the two disagreeing by ~0.01%, which could be due to some inherent difference in the two treatments still but regardless is much smaller than the present error, and is the same for all levels of up-down asymmetry. We can do more convergence testing for 2.0.0 while we are writing the paper IMO and that could expose more errors, that is the sort of thing has exposed all these little numerical improvements for Delta' calculations we have merged in recently, and I think a series of convergence plots for key quantities would make some nice appendix plots for the paper.

[logan-nc]

Great sleuthing @matt-pharr!
Will you implement the fix here or in a separate PR?

[matt-pharr]

@logan-nc I can do so here since I have already posted the plots here. I opened a separate PR in fortran to correct this.

[matt-pharr]

Went through and checked everything for correctness, some docs are a bit egregious in terms of length and verbosity but I think it is not worth holding this on that since we want to do a docstring pass once the dust has settled from 2.0.0 API changes. Merging this now