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

docs/development/hdf5-conventions.md

@@ -43,7 +43,7 @@ Top level (11 groups):
 | `ForceFreeStates/` | `Solutions/ForwardIntegration/` (u-solutions), `Solutions/GalerkinIntegration/` (closed ξ profiles in the shared layout, the gal surface list, debug-gated `Basis/`), `EulerLagrangeMatrices/{Ideal,Kinetic}`, `FreeBoundaryStability/`, `EdgeScan/` |
 | `LocalStability/` | Mercier `D_I`, resistive interchange `D_R`, `ballooning_Delta_prime` on `psi`; the ballooning α boundary on `ballooning_psi` |
 | `SingularSurfaces/` | Per-rational-surface data: `rational_psi`/`rational_q`/`rational_m`/`rational_n`, GGJ coefficients, `Delta_prime_matrix`/`Delta_prime_raw`/`Delta_coil`/`pest3_A`/`pest3_B`/`pest3_Gamma` (Riccati or Galerkin alike), the inner-layer `Match/` results, `Kinetic/` |
-| `PerturbedEquilibrium/` | `ForcingModes/`, `Response/`, `ResponseMatrices/`, `SingularCoupling/`, `Energies/`, control-surface spectra |
+| `PerturbedEquilibrium/` | `ForcingModes/`, `Response/`, `ResponseMatrices/`, `SingularCoupling/`, `Energies/`, `TorqueResponse/` (kinetic forward solves: ψ-resolved `T_xe`, `T_source`), control-surface spectra |
 | `KineticForces/` | `<method>/` (torque/energy profiles, `dTdpsi_ell` resolved by bounce harmonic on the `ell` scale, `EnergyIntegrals/`, `KineticMatrices/`); multi-ion runs add `PerSpecies/<species>/<method>/` with the same per-method layout, summing to the top-level total |
 | `ErrorFields/` | `CoilSensitivities/` (per-coil-set control-surface spectra and their rigid shift/tilt derivatives, `DominantMode/` full-window projection); `MonteCarlo/` (intrinsic and corrected `\|δ\|` histograms over the sampled tolerances, per batch and averaged); `Risk/` (threshold density, `P(lock\|δ)`, locking probabilities, `ToleranceScan/`); `NTV/` (correction-coil overlap and NTV torque couplings per kAt) |
 | `Tearing/` | `PerSurface/` (+ `DpMatrix/`), `Roots/`, `LayerWidths/`, `Diagnostics/{ValidRoots,Poles,FilteredRoots}`, `Scan/Surface_<k>/` |

docs/src/perturbed_equilibrium.md

@@ -81,6 +81,55 @@ GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.check_mode_basis
 GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.compute_dominant_coupling!
 ```
 
+## Torque response matrices
+
+A kinetic forward solve (`integrator = "forward"`, `kinetic_factor > 0`) carries the
+drift-kinetic response inside its Euler-Lagrange energy: the complex plasma energy inside a
+flux surface, `δW(ψ) = ξ†·U₂(ψ)·U₁(ψ)⁻¹·ξ/(2μ₀)`, has an anti-Hermitian part that is the
+neoclassical toroidal viscosity torque (Logan & Park 2013, eq. 19). With
+`compute_torque_response = true` the response stage refers that quadratic form back to the
+applied root-area-weighted control-surface field b̃ — through the permeability, the
+boundary flux-to-displacement factor `χ₁(m − n·q_lim)·2πi` and the b̃ → Φ conform `R` — and
+stores the result at every stored node as the **torque response matrix** `T_xe(ψ)`, cumulative
+in ψ, in `PerturbedEquilibrium/TorqueResponse/`. It is the Fortran GPEC `T_xe`
+(`gpout_dw_matrix`), in the same basis as the stored response matrices, with units N·m/T².
+
+Two conventions are worth stating explicitly:
+
+  - **Factor of 2.** `T_xe` matches the Fortran array, whose quadratic form is built on the `+n`
+    harmonic alone and is therefore twice the physical value. The physical cumulative complex
+    torque of an applied spectrum is `T(ψ) = b̃†·T_xe(ψ)·b̃/2` — `torque_profile` applies the
+    1/2 — with `real(T)` the toroidal torque on the plasma inside ψ and `imag(T)` its `2n·δW`
+    (ideal plus kinetic). At the control surface `real(T)` equals `Energies/toroidal_torque`.
+  - **Hermitian part.** Only the Hermitian part `(T_xe + T_xe†)/2` contributes to the torque;
+    its eigenvectors are the applied spectra of extremal torque and its eigenvalues their torque
+    per unit `‖b̃‖²/2`. The anti-Hermitian part carries the energy.
+
+The forcing-source form `T_source(ψ) = M†·T·M` is written alongside, with one column of `M` per
+forcing source of the run's `[ForcingTerms]`, each as built: every coil set at its deck currents,
+or the forcing file's spectrum as a single source. Driving the sources at scale factors `s` gives
+`s†·T_source·s/2`.
+
+```julia
+using GeneralizedPerturbedEquilibrium.PerturbedEquilibrium
+using LinearAlgebra
+
+tr = pe_state.torque_response                       # a TorqueResponse
+T = torque_profile(tr, pe_state.forcing_b_rootarea) # cumulative complex torque of the run's forcing
+real(T[end])                                        # total toroidal torque [N·m]
+
+T_h = Hermitian((tr.T_xe[:, :, end] + tr.T_xe[:, :, end]') / 2)
+λ, V = eigen(T_h)                                   # V[:, end]: applied b̃ of maximum torque
+```
+
+```@docs
+GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.TorqueResponse
+GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.torque_response_matrices
+GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.torque_profile
+GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.coil_flux_spectra
+GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.compute_torque_response!
+```
+
 ## Plotting per-surface results against ψ or q
 
 `SingularCoupling/` quantities are indexed by rational-surface **index**, not by q: with

examples/Solovev_kinetic_calculated_example/gpec.toml

@@ -4,7 +4,8 @@
 # the multi-ion species loop (pure plasma, n_i = n_e, so no impurity species arises). The
 # equilibrium is generated analytically from the embedded [SOL_INPUT] section, with n·T tied
 # to the Solovev P₀≈4.27e4 Pa via the local kinetic.dat.
-# Also the only deck that runs [PerturbedEquilibrium] on a kinetic solve.
+# Also the only deck that runs [PerturbedEquilibrium] on a kinetic solve, so it builds the
+# ψ-resolved torque response matrices (T_xe on the applied field, T_source for the forcing file).
 
 [Equilibrium]
 eq_type = "sol"                         # Type of the input 2D equilibrium file
@@ -37,6 +38,7 @@ fixed_boundary = false                  # Use fixed boundary conditions
 output_eigenmodes = true                # Output eigenmode fields as b-fields
 compute_response = true                 # Compute plasma response to forcing
 compute_singular_coupling = true        # Compute singular layer coupling metrics
+compute_torque_response = true          # Build the ψ-resolved torque response matrices T_xe and T_source
 verbose = false                         # Enable verbose logging
 write_outputs_to_HDF5 = true            # Write perturbed equilibrium outputs to HDF5
 regularization_width = 0.05             # Regularization width for singular surfaces (0 = disabled)

regression-harness/cases/solovev_kinetic_calculated.toml

@@ -113,6 +113,32 @@ extract = "value"
 label = "npert"
 noise_threshold = 0
 
+# Torque of the run's own forcing, three routes on one kinetic solve: the ψ-resolved torque
+# response matrix T_xe contracted with the applied field (its Hermitian part carries the torque),
+# the boundary-response torque from the plasma inductance (same number, different algebra), and
+# the KineticForces NTV quadrature of the same displacement (an independent evaluation that
+# agrees only to the level the kinetic-matrix and quadrature discretizations do).
+[quantities.torque_response_total]
+h5path = "PerturbedEquilibrium/TorqueResponse/total_torque"
+type = "real_scalar"
+extract = "value"
+label = "T_xe applied-forcing torque [N·m]"
+noise_threshold = 1e-10
+
+[quantities.boundary_torque]
+h5path = "PerturbedEquilibrium/Energies/toroidal_torque"
+type = "real_scalar"
+extract = "value"
+label = "boundary-response torque [N·m]"
+noise_threshold = 1e-10
+
+[quantities.ntv_fgar_torque]
+h5path = "KineticForces/fgar/total_torque"
+type = "real_scalar"
+extract = "value"
+label = "NTV fgar quadrature torque [N·m]"
+noise_threshold = 1e-10
+
 # Perturbed equilibrium — Clebsch displacements and regularized b^θ from the kinetic solve
 [quantities.pe_clebsch_alpha_norm]
 h5path = "PerturbedEquilibrium/Response/xi_clebsch_alpha"

src/ForceFreeStates/CoreTypes.jl

@@ -138,7 +138,7 @@ gpec.toml.
   - `rational_offset_fraction::Float64` - Distance beyond last rational surface (normalised ∈ [0,1) in units of 1/n). Only used when `truncate_at_rational_offset` is true. Fortran STRIDE convention is 0.2 (truncate 20 % of one rational-surface spacing above the last surface), retained here.
   - `frobenius_order::Int` - Order of singular layer (Frobenius) expansion at rational surfaces. Default 6 (Fortran STRIDE convention for Δ' calculations; lower values trade accuracy for speed).
   - `qhigh::Float64` - Integration terminated at q limit determined by minimum of qhigh and qa from equil
-  - `kinetic_source::String` - Kinetic matrix source: "fixed" (X-shaped test matrices scaled by kinetic_factor relative to ideal matrix Frobenius norms; Ak, Dk, Hk Hermitian, Bk, Ck, Ek non-Hermitian), "calculated" (PENTRC — not yet implemented)
+  - `kinetic_source::String` - Kinetic matrix source: "fixed" (X-shaped test matrices scaled by kinetic_factor relative to ideal matrix Frobenius norms; Ak, Dk, Hk Hermitian, Bk, Ck, Ek non-Hermitian), "calculated" (the drift-kinetic matrices of the KineticForces model, evaluated on the `[KineticForces]` profiles; needs `integrator = "forward"`)
   - `kinetic_factor::Float64` - Dimensionless scaling factor for kinetic matrices. Zero (the default) disables the kinetic path; any positive value enables it and scales the kinetic matrices: when kinetic_source="fixed", scales X-shaped test matrices relative to ideal matrix norms; when kinetic_source="calculated", applied as uniform post-hoc multiplier to W and T components.
   - `qlow::Float64` - Integration terminated at q limit determined by minimum of qlow and q0 from equil
   - `dW_edge_scan_start::Float64` - If less than psilim, records a dW(ψ) diagnostic scan over [dW_edge_scan_start, psilim] on odet.edge_scan. The integration domain (psilim) is always controlled by qhigh / psihigh and is not modified by this scan (unless `truncate_at_dW_peak=true`, see caveats below).

src/KineticForces/BounceAveraging.jl

@@ -1024,11 +1024,11 @@ function _bounce_integrate(
             for i in 1:ntheta
                 wsamp[i] = conj(wmu_mt[mi, i]) * (pl[i] + one_minus_sigma / pl[i])
             end
-            wmu_ba[mi] = dot(scr.int_w, wsamp)
+            wmu_ba[mi] = conj(dot(scr.int_w, wsamp))
             for i in 1:ntheta
                 wsamp[i] = conj(wen_mt[mi, i]) * (pl[i] + one_minus_sigma / pl[i])
             end
-            wen_ba[mi] = dot(scr.int_w, wsamp)
+            wen_ba[mi] = conj(dot(scr.int_w, wsamp))
         end
 
         # Reshape as 1×mpert for matrix multiply (Fortran lines 771-772)

src/PerturbedEquilibrium/PerturbedEquilibrium.jl

@@ -22,6 +22,7 @@ import DelimitedFiles: readdlm
 # Include module files
 include("PerturbedEquilibriumStructs.jl")
 include("ResponseMatrices.jl")
+include("TorqueResponse.jl")
 include("FieldReconstruction.jl")
 include("Response.jl")
 include("SingularCoupling.jl")
@@ -38,6 +39,7 @@ export ForcingMode
 export compute_perturbed_equilibrium
 export write_outputs_to_HDF5
 export ResonantCoupling, DominantCoupling, dominant_coupling, rootarea_field, coupling_overlap, check_mode_basis, CORE_PSI_LOW, CORE_PSI_HIGH
+export TorqueResponse, torque_response_matrices, torque_profile
 
 """
     compute_perturbed_equilibrium(ffs, forcing, ctrl, intr; runtimes=nothing)::PerturbedEquilibriumState
@@ -102,6 +104,14 @@ function compute_perturbed_equilibrium(
         compute_plasma_response!(state, equil, solution, ffs.free_boundary.wt0, mtheta_vacuum, ffs, intr, ctrl, ffs.metric, mats)
     end
 
+    # Torque response matrices of a kinetic forward solve, on the plasma response just computed
+    if ctrl.compute_torque_response
+        ctrl.compute_response || throw(ArgumentError("compute_torque_response needs compute_response = true"))
+        ForceFreeStates.require_solution(ffs, "torque response matrices") ||
+            throw(ArgumentError("compute_torque_response needs the dense forward-integrated ξ solution"))
+        compute_torque_response!(state, ffs, forcing, intr, ctrl)
+    end
+
     # Step 3: Compute singular coupling metrics
     if ctrl.compute_singular_coupling &&
        ForceFreeStates.require(ffs, :free_boundary, "singular coupling calculation") &&

src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl

@@ -1,3 +1,38 @@
+"""
+    TorqueResponse
+
+ψ-resolved torque response of a forward kinetic Euler-Lagrange solve, expressed on the
+root-area-weighted control-surface field b̃ (the basis of `forcing_b_rootarea` and of the
+stored response matrices). The matrices are cumulative in ψ: `T_xe[:, :, k]` is the response
+of the plasma inside `psi[k]`, and the last slice is the total.
+
+The quadratic form `b̃†·T_xe(ψ)·b̃` is the complex "2×" torque, matching the Fortran GPEC
+`T_xe` array: the physical torque and energy of the applied spectrum `b̃` are recovered with a
+factor 1/2 (`torque_profile`), `T(ψ) = b̃†·T_xe(ψ)·b̃/2`, whose real part is the toroidal torque
+on the plasma inside ψ and whose imaginary part is `2n·δW` of that plasma (ideal plus kinetic).
+
+## Fields
+
+  - `psi::Vector{Float64}` - ψ_N of each stored node of the solve `[npsi]`
+  - `m_modes`, `n_modes::Vector{Int}` - (m, n) of each row and column of `T_xe` `[numpert_total]`
+  - `T_xe::Array{ComplexF64,3}` - torque response on b̃, N·m/T² `[numpert_total × numpert_total × npsi]`
+  - `source_names::Vector{String}` - name of each forcing source of `T_source`: the coil sets, or the forcing file `[nsource]`
+  - `T_source::Array{ComplexF64,3}` - torque response in forcing-source space `[nsource × nsource × npsi]`,
+    N·m per squared scale factor of each source as built: driving the sources at `s` times their deck
+    currents or file amplitudes gives the complex 2× torque `s†·T_source·s`. Empty without forcing.
+  - `T_applied::Vector{ComplexF64}` - the run's own forcing contracted, `b̃ₓ†·T_xe(ψ)·b̃ₓ/2` `[npsi]`;
+    empty when no forcing spectrum was supplied
+"""
+struct TorqueResponse
+    psi::Vector{Float64}
+    m_modes::Vector{Int}
+    n_modes::Vector{Int}
+    T_xe::Array{ComplexF64,3}
+    source_names::Vector{String}
+    T_source::Array{ComplexF64,3}
+    T_applied::Vector{ComplexF64}
+end
+
 """
     PerturbedEquilibriumControl
 
@@ -13,6 +48,7 @@ High Priority (MWE):
   - `output_eigenmodes::Bool` - Output mode fields as b-fields (default: true)
   - `compute_response::Bool` - Compute plasma response (default: true)
   - `compute_singular_coupling::Bool` - Compute singular coupling metrics (default: true)
+  - `compute_torque_response::Bool` - Build the ψ-resolved torque response matrices `T_xe` and `T_source` (default: false). Needs a forward-integrated kinetic solve (`integrator = "forward"`, `kinetic_factor > 0`) and `compute_response`; errors otherwise
   - `verbose::Bool` - Enable verbose logging (default: true)
 
 Output Settings:
@@ -37,6 +73,7 @@ Regularization:
     output_eigenmodes::Bool = true
     compute_response::Bool = true
     compute_singular_coupling::Bool = true
+    compute_torque_response::Bool = false
     verbose::Bool = true
 
     # Output settings
@@ -163,6 +200,11 @@ well-conditioned flux-space inductances L, Λ:
     the boundary-response torque, zero for ideal
     (Hermitian) runs. Equals the volume-integrated Euler-Lagrange kinetic torque only for converged
     self-consistent solutions, and is a distinct construction from the KineticForces NTV torque.
+
+Torque response (kinetic forward solves with `compute_torque_response`):
+
+  - `torque_response::Union{Nothing,TorqueResponse}` - the ψ-resolved torque response matrices on the
+    applied b̃ and on the coil sets, with the run's own forcing contracted; see [`TorqueResponse`](@ref)
 """
 @kwdef mutable struct PerturbedEquilibriumState
     # Radial grid (FFS ODE integration ψ_n values) [npsi]
@@ -229,6 +271,9 @@ well-conditioned flux-space inductances L, Λ:
     surface_energy::Float64 = 0.0
     plasma_energy::Float64 = 0.0
     toroidal_torque::Float64 = 0.0
+
+    # ψ-resolved torque response matrices of a kinetic forward solve
+    torque_response::Union{Nothing,TorqueResponse} = nothing
 end
 
 """