PerturbedEquilibrium - FEATURE - Add the psi-resolved torque response matrices T_xe and T_coil

[logan-nc]

Release note

• Audience: users
• Numerical impact: none — new output only; every previously tracked quantity is bit-identical, and the solovev_kinetic_calculated example deck gains coil forcing so three torque routes are now tracked there (harness @ 8a60480)
• 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_coil(ψ) in coil-set space, 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.

Regression report

regress --cases solovev_kinetic_calculated,diiid_n1 --refs develop,feature/pe-torque-response-matrix from the main repo root (job 126616):

Regression Report: diiid_n1
===================================================================================================================
Ref 1: develop  @ 9ff10c711 (2026-09-25)
       env: julia 1.11.7, x86_64-linux-gnu, manifest 06e27666 (pinned), 16 threads/16 BLAS
Ref 2: feature/pe-torque-response-matrix  @ 8a6048000 (2026-09-25)
       env: julia 1.11.7, x86_64-linux-gnu, manifest 06e27666 (pinned), 16 threads/16 BLAS
-------------------------------------------------------------------------------------------------------------------
Quantity                                      develop          feature/pe-torque-response-matrix  Diff       Status
-------------------------------------------------------------------------------------------------------------------
total energy Re(et[1])                        8.012318e-01     8.012318e-01                       0.0e+00    OK
total energy Im(et[1])                        4.142529e-05     4.142529e-05                       0.0e+00    OK
plasma energy Re(ep[1])                       -1.348486e+00    -1.348486e+00                      0.0e+00    OK
vacuum energy Re(ev[1])                       2.149718e+00     2.149718e+00                       0.0e+00    OK
vacuum matrix min eigenvalue                  1.873976e-01     1.873976e-01                       0.0e+00    OK
plasma energy (all)                           [35 elem]        [35 elem]                          0.0e+00    OK
vacuum energy (all)                           [35 elem]        [35 elem]                          0.0e+00    OK
total energy (all)                            [35 elem]        [35 elem]                          0.0e+00    OK
ODE steps (saved)                             2576             2576                               0.0e+00    OK
ODE steps (total)                             4572             4572                               0.0e+00    OK
q0                                            1.204212e+00     1.204212e+00                       0.0e+00    OK
q95                                           4.781723e+00     4.781723e+00                       0.0e+00    OK
beta_t                                        1.327024e-02     1.327024e-02                       0.0e+00    OK
beta_n                                        1.372511e+00     1.372511e+00                       0.0e+00    OK
internal inductance li1                       8.842392e-01     8.842392e-01                       0.0e+00    OK
internal inductance li2                       7.080847e-01     7.080847e-01                       0.0e+00    OK
internal inductance li3                       7.304433e-01     7.304433e-01                       0.0e+00    OK
poloidal beta betap1                          6.680744e-01     6.680744e-01                       0.0e+00    OK
poloidal beta betap2                          5.349834e-01     5.349834e-01                       0.0e+00    OK
poloidal beta betap3                          5.518761e-01     5.518761e-01                       0.0e+00    OK
# singular surfaces                           5                5                                  0.0e+00    OK
singular psi locations                        [5 elem]         [5 elem]                           0.0e+00    OK
singular q values                             [5 elem]         [5 elem]                           0.0e+00    OK
current beta betaj                            4.236478e-01     4.236478e-01                       0.0e+00    OK
plasma volume                                 1.829472e+01     1.829472e+01                       0.0e+00    OK
plasma current                                1.152130e+00     1.152130e+00                       0.0e+00    OK
mpert                                         35               35                                 0.0e+00    OK
npert                                         1                1                                  0.0e+00    OK
toroidal field bt0                            2.006573e+00     2.006573e+00                       0.0e+00    OK
wall field bwall                              3.880145e-01     3.880145e-01                       0.0e+00    OK
aspect ratio                                  2.845746e+00     2.845746e+00                       0.0e+00    OK
elongation kappa                              1.708350e+00     1.708350e+00                       0.0e+00    OK
q profile (checksum)                          0cd285cea88d...  0cd285cea88d...                    identical  OK
pressure profile (checksum)                   a1c48b266622...  a1c48b266622...                    identical  OK
Mercier D_I profile (checksum)                eeb06744e795...  eeb06744e795...                    identical  OK
resistive interchange D_R profile (checksum)  fa37296851f6...  fa37296851f6...                    identical  OK
ballooning Delta' profile (checksum)          6eb0ea075ccf...  6eb0ea075ccf...                    identical  OK
island half-widths                            [5 elem]         [5 elem]                           0.0e+00    OK
Chirikov parameter                            [5 elem]         [5 elem]                           0.0e+00    OK
||resonant area-weighted field||              5.189179e-04     5.189179e-04                       0.0e+00    OK
dominant-coupling singular values             [3 elem]         [3 elem]                           0.0e+00    OK
|forcing overlap with dominant mode|          1.415683e-04     1.415683e-04                       0.0e+00    OK
|delta_nominal| of coil set 1                 7.055228e-05     7.055228e-05                       0.0e+00    OK
||ddelta/d(shift)|| over coil sets            1.246233e-04     1.246233e-04                       0.0e+00    OK
||ddelta/d(tilt)|| over coil sets             4.180643e-07     4.180643e-07                       0.0e+00    OK
PE plasma energy                              3.422677e+00     3.422677e+00                       0.0e+00    OK
PE vacuum energy                              3.174510e+00     3.174510e+00                       0.0e+00    OK
PE surface energy                             5.826684e+00     5.826684e+00                       0.0e+00    OK
PE toroidal torque                            5.087465e-02     5.087465e-02                       0.0e+00    OK
NTV torque FGAR [N·m]                         5.296762e-01     5.296762e-01                       0.0e+00    OK
NTV kinetic energy dW FGAR [J]                7.924971e-02     7.924971e-02                       0.0e+00    OK
Runtime (s)                                   280.4s           272.3s                                        --
||forcing b~|| (root-area-weighted)           4.683328e-04     4.683328e-04                       0.0e+00    OK
resonant area-weighted field b^r              [5 elem]         [5 elem]                           0.0e+00    OK
===================================================================================================================
Summary: 53 unchanged

The main root's case file predates this branch, so the Solovev report there lists only the 14 pre-existing quantities (all unchanged). Re-extracted with this branch's case file (--db pointed at the same cache, --force, job 126618), showing the three new torque quantities:

Regression Report: solovev_kinetic_calculated
====================================================================================================
Ref 1: develop  @ 9ff10c711 (2026-09-25)
       env: julia 1.11.7, x86_64-linux-gnu, manifest 06e27666 (pinned), 16 threads/16 BLAS
Ref 2: feature/pe-torque-response-matrix  @ 8a6048000 (2026-09-25)
       env: julia 1.11.7, x86_64-linux-gnu, manifest 06e27666 (pinned), 16 threads/16 BLAS
----------------------------------------------------------------------------------------------------
Quantity                           develop        feature/pe-torque-response-matrix  Diff     Status
----------------------------------------------------------------------------------------------------
vacuum energy Re(ev[1])            1.043757e+01   1.043757e+01                       0.0e+00  OK
q0                                 1.900003e+00   1.900003e+00                       0.0e+00  OK
NTV fgar quadrature torque [N·m]   N/A            7.117756e+00                       N/A      N/A
total energy Re(et[1])             1.873745e+00   1.873745e+00                       0.0e+00  OK
total energy Im(et[1])             -1.434288e+00  -1.434288e+00                      0.0e+00  OK
plasma energy Re(ep[1])            -8.563829e+00  -8.563829e+00                      0.0e+00  OK
ODE steps (saved)                  599            599                                0.0e+00  OK
# singular surfaces                2              2                                  0.0e+00  OK
boundary-response torque [N·m]     N/A            4.129557e-02                       N/A      N/A
singular psi locations             [2 elem]       [2 elem]                           0.0e+00  OK
total energy (all)                 [32 elem]      [32 elem]                          0.0e+00  OK
ODE steps (total)                  734            734                                0.0e+00  OK
mpert                              32             32                                 0.0e+00  OK
q95                                3.147422e+00   3.147422e+00                       0.0e+00  OK
Runtime (s)                        213.6s         247.6s                                      --
T_xe applied-forcing torque [N·m]  N/A            4.129557e-02                       N/A      N/A
singular q values                  [2 elem]       [2 elem]                           0.0e+00  OK
npert                              1              1                                  0.0e+00  OK
====================================================================================================
Summary: 14 unchanged, 3 missing/N/A

Every pre-existing quantity is unchanged. The solovev_kinetic_calculated case deck now carries two analytic window-pane coil sets and [PerturbedEquilibrium], so it gains three tracked torques (TorqueResponse/total_torque, Energies/toroidal_torque, KineticForces/fgar/total_torque) and its runtime grows by the response and NTV-quadrature stages.

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_coil = Kc†·[U₁†U₂]·Kc·2n·i/(2μ₀) with Kc = 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.
• Coil spectrum. 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_coil)/2 == T_applied at every node. Driving the sets at scale factors s gives s†·T_coil·s/2. Left empty for file forcing and RMPFieldSum.
• 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): 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_coil vs M̃†T_xe M̃, shapes vs Info/mn_index, and the HDF5 metadata contract on the new group.

Benchmark against the NTV quadrature (examples/DIIID-like_error_field_example, C-coil only, kinetic_source="calculated", kinetic_factor=1, forward, grid_type="ldp", mpsi=64, kinetic profiles built pressure-consistent with the geqdsk; scratch deck, not committed):

domain	T_xe-contracted torque	boundary-response torque	KineticForces fgar	T_xe / fgar
psihigh=0.985 (ends between q=5.2 and 6), tol 1e-10	0.4342 N·m	0.4342 N·m	0.3843 N·m	1.130 (13 %)

Open question: the two routes differ by 13 %. They are two evaluations of the same torque from the same response and should agree better than that; I have not resolved why they do not, and this PR does not claim they agree. Best hypothesis: discretization. T_xe inherits the torque from the Euler-Lagrange kinetic matrices — the drift-kinetic operator projected onto the (m, m′) Fourier basis, evaluated on the equilibrium ψ grid (mpsi=64) and splined — whereas fgar re-integrates the full ψ–energy–pitch integrand of the response ξ with its own adaptive quadrature, whose outer tolerance defaults to rtol_psi = 1e-2 and whose ψ panels are placed independently of the solve's grid. Near the edge the kinetic profile is steep (n floored at 1e17 m⁻³ at ψ_N = 1) and the response is dominated by the q = 5 resonance at ψ_N ≈ 0.974, four grid intervals from ψ_lim, so a few-percent difference in how each route resolves that region is plausible; the run's "W inverse matrix was non-Hermitian beyond tolerance at 1497 integration steps" warning is the kinetic solve's expected non-Hermiticity and is not by itself evidence for either route. The check that would settle it is a convergence sweep in mpsi and rtol_psi on this deck; that is a follow-up. The Hermitian T_xe(ψ_lim) is positive definite there (λ ∈ [2.6, 6.0e7] N·m/T²).

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.