diff --git a/docs/development/hdf5-conventions.md b/docs/development/hdf5-conventions.md index 3ec3bb633..31a05e3e0 100644 --- a/docs/development/hdf5-conventions.md +++ b/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/` | `/` (torque/energy profiles, `dTdpsi_ell` resolved by bounce harmonic on the `ell` scale, `EnergyIntegrals/`, `KineticMatrices/`); multi-ion runs add `PerSpecies///` 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_/` | diff --git a/docs/src/perturbed_equilibrium.md b/docs/src/perturbed_equilibrium.md index 5e079cfee..4bbdd6081 100644 --- a/docs/src/perturbed_equilibrium.md +++ b/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 diff --git a/examples/Solovev_kinetic_calculated_example/gpec.toml b/examples/Solovev_kinetic_calculated_example/gpec.toml index a6bce31b9..c663ab844 100644 --- a/examples/Solovev_kinetic_calculated_example/gpec.toml +++ b/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) diff --git a/regression-harness/cases/solovev_kinetic_calculated.toml b/regression-harness/cases/solovev_kinetic_calculated.toml index db6e7f8a0..05875c406 100644 --- a/regression-harness/cases/solovev_kinetic_calculated.toml +++ b/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" diff --git a/src/ForceFreeStates/CoreTypes.jl b/src/ForceFreeStates/CoreTypes.jl index ef384af7b..731f27f3b 100644 --- a/src/ForceFreeStates/CoreTypes.jl +++ b/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). diff --git a/src/KineticForces/BounceAveraging.jl b/src/KineticForces/BounceAveraging.jl index 0f974a09c..b04d0c471 100644 --- a/src/KineticForces/BounceAveraging.jl +++ b/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) diff --git a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl index e1e7703cc..2e0fb8ee6 100644 --- a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl +++ b/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") && diff --git a/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl b/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl index 18cd783f3..3353f919c 100644 --- a/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl +++ b/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 """ diff --git a/src/PerturbedEquilibrium/TorqueResponse.jl b/src/PerturbedEquilibrium/TorqueResponse.jl new file mode 100644 index 000000000..fd6cc7d70 --- /dev/null +++ b/src/PerturbedEquilibrium/TorqueResponse.jl @@ -0,0 +1,171 @@ +""" +ψ-resolved torque response matrices of a kinetic Euler-Lagrange solve. + +The complex plasma energy inside a flux surface ψ of the kinetic solve is the quadratic form +δW(ψ) = ξ†·U₂(ψ)·U₁(ψ)⁻¹·ξ/(2μ₀) of the displacement at ψ; its anti-Hermitian part is the +neoclassical toroidal viscosity torque [Logan 2013, eq. 19]. Referring the displacement back +to the applied control-surface field turns it into a matrix over the applied spectrum, the +torque response matrix, whose Hermitian part is the torque quadratic form and whose +eigenvectors are the applied spectra of extremal torque [Logan 2015, ch. 7]. +""" + +""" + torque_response_matrices(ffs, permeability, flux_conform; + forcing_b_rootarea=ComplexF64[], source_flux=zeros(ComplexF64, N, 0), source_names=String[]) -> TorqueResponse + +Build the ψ-resolved torque response of a forward kinetic solve from its stored solution. +`permeability` is the flux-space plasma response operator `P` (Φ_tot = P·Φ_x) and `flux_conform` +the b̃ → Φ operator `R = S·A`, both at the control surface; `source_flux` holds one unit-norm +Φ_x column per forcing source and `forcing_b_rootarea` the run's applied b̃ for `T_applied`. + +Per stored node `k`, with `U₁`, `U₂` the displacement and momentum blocks of `u_store` and +`U₁ₗ` the displacement block at the control surface: + + T_xe(ψ_k) = K†·[U₁(ψ_k)†·U₂(ψ_k)]·K · 2n·i/(2μ₀), K = U₁ₗ⁻¹·D·P·R, D = diag(1/(χ₁·(m − n·q_lim)·2π·i)) + +`K` maps an applied b̃ to the coefficients of the solution basis (b̃ → Φ_x → Φ_tot → boundary ξ +→ coefficients), so `T_xe` is the Fortran `gresp` chain (`dwks`, `bsurfmat`, `gind`, `gres`, +`ptof`) with the single edge inverse the stored, globally consistent basis allows. +`T_source = Ks†·[U₁†U₂]·Ks · 2n·i/(2μ₀)` with `Ks = U₁ₗ⁻¹·D·P·source_flux` is the Fortran `gcoil`. + +Errors unless the solve is a single-`n` forward integration with kinetic matrices; the ideal +Euler-Lagrange energy is Hermitian and carries no torque. Warns when the kinetic matrices are +the synthetic `"fixed"` test matrices. +""" +function torque_response_matrices( + ffs::ForceFreeStatesResult, + permeability::AbstractMatrix{ComplexF64}, + flux_conform::AbstractMatrix{ComplexF64}; + forcing_b_rootarea::AbstractVector{ComplexF64}=ComplexF64[], + source_flux::AbstractMatrix{ComplexF64}=zeros(ComplexF64, ffs.numpert_total, 0), + source_names::Vector{String}=String[] +)::TorqueResponse + ffs.integrator === :forward || + throw(ArgumentError("torque response matrices need the dense forward-integrated solution; this result came from the $(ffs.integrator) integrator")) + ffs.solution === nothing && throw(ArgumentError("torque response matrices need a ξ solution; the forward integrator produced none")) + ffs.mats.kinetic === nothing && + throw(ArgumentError("torque response matrices need a kinetic solve (kinetic_factor > 0); the ideal Euler-Lagrange energy is Hermitian and carries no torque")) + ffs.npert == 1 || throw(ArgumentError("torque response matrices are defined per toroidal harmonic; this solve stacks $(ffs.npert) (nn_low != nn_high)")) + ffs.control.kinetic_source == "fixed" && + @warn "kinetic_source = \"fixed\": the torque response is built from synthetic test matrices and has no physical meaning" + size(source_flux, 2) == length(source_names) || + throw(DimensionMismatch("source_flux has $(size(source_flux, 2)) columns but $(length(source_names)) source names were given")) + + solution = ffs.solution + N = ffs.numpert_total + npsi = solution.step + nn = ffs.nlow + μ₀ = 4π * 1e-7 + chi1 = 2π * ffs.equil.psio + m_modes = [(i - 1) % ffs.mpert + ffs.mlow for i in 1:N] + n_modes = fill(nn, N) + + # Boundary flux → boundary displacement, Φ_tot = χ₁·(m − n·q_lim)·2π·i·ξ (the same factor Λ carries). + t = [-im / (chi1 * (m_modes[i] - nn * ffs.qlim) * 2π) for i in 1:N] + U1_lim = solution.u_store[:, :, 1, npsi] + boundary_xi = Diagonal(t) * permeability + K = U1_lim \ (boundary_xi * flux_conform) + Ks = U1_lim \ (boundary_xi * source_flux) + + nsource = size(source_flux, 2) + T_xe = Array{ComplexF64}(undef, N, N, npsi) + T_source = Array{ComplexF64}(undef, nsource, nsource, npsi) + prefactor = 2 * nn * im / (2 * μ₀) + W = Matrix{ComplexF64}(undef, N, N) + for k in 1:npsi + mul!(W, adjoint(@view(solution.u_store[:, :, 1, k])), @view(solution.u_store[:, :, 2, k])) + W .*= prefactor + T_xe[:, :, k] = K' * W * K + nsource > 0 && (T_source[:, :, k] = Ks' * W * Ks) + end + + T_applied = isempty(forcing_b_rootarea) ? ComplexF64[] : torque_profile(T_xe, forcing_b_rootarea) + return TorqueResponse(solution.psi_store[1:npsi], m_modes, n_modes, T_xe, source_names, T_source, T_applied) +end + +""" + torque_profile(T_xe::AbstractArray{ComplexF64,3}, b̃) -> Vector{ComplexF64} + torque_profile(tr::TorqueResponse, b̃) -> Vector{ComplexF64} + +Contract an applied root-area-weighted spectrum `b̃` with a torque response, giving the +cumulative complex torque `T(ψ_k) = b̃†·T_xe(ψ_k)·b̃/2` at every stored node: `real` is the +toroidal torque on the plasma inside ψ_k in N·m and `imag` is `2n·δW` of that plasma in J. + +The factor 1/2 converts the stored "2×" matrix (the Fortran GPEC `T_xe` array, built on the +quadratic form of the `+n` harmonic alone) to the physical torque of the real field; the +last entry then equals the boundary-response `toroidal_torque` of the same solve. Only the +Hermitian part of `T_xe` contributes to the real part, so the torque is the same whether or +not the matrix is symmetrized first. +""" +function torque_profile(T_xe::AbstractArray{ComplexF64,3}, b̃::AbstractVector{<:Number}) + N = size(T_xe, 1) + length(b̃) == N || throw(DimensionMismatch("applied spectrum has $(length(b̃)) entries; the response matrix has $N rows")) + x = Vector{ComplexF64}(b̃) + return [dot(x, @view(T_xe[:, :, k]), x) / 2 for k in 1:size(T_xe, 3)] +end + +torque_profile(tr::TorqueResponse, b̃::AbstractVector{<:Number}) = torque_profile(tr.T_xe, b̃) + +""" + coil_flux_spectra(ffs, coil_sets, cfg) -> Matrix{ComplexF64} + +Unit-norm control-surface flux spectrum `Φ_x` of every coil set as built, one column per set +on the solve's (m, n) ordering `[numpert_total × ncoil_set]`: Biot-Savart on the control +surface `ffs.psilim` sampled with the run's own `cfg` grid, so the columns sum to the run's +coil forcing. +""" +function coil_flux_spectra(ffs::ForceFreeStatesResult, coil_sets::Vector{CoilSet}, cfg::ForcingTerms.CoilConfig) + grids = [(n, ForcingTerms.CoilForcingGrid(ffs.equil, cfg, n; psi=ffs.psilim)) for n in ffs.nlow:ffs.nhigh] + M = zeros(ComplexF64, ffs.numpert_total, length(coil_sets)) + for (j, cs) in enumerate(coil_sets) + modes = ForcingMode[] + for (n, grid) in grids + append!(modes, ForcingTerms.coil_forcing_modes(cs, grid, n, ffs.mlow, ffs.mhigh)) + end + M[:, j] = map_forcing_to_eigenmodes(modes, ffs) + end + return M +end + +""" + compute_torque_response!(state, ffs, forcing, intr, ctrl) + +Build the torque response of the solve from the plasma response already on `state` and +`intr` and store it in `state.torque_response`. Every forcing source is one column of +`T_source`: each coil set (`intr.coil_sets`, on the grid `forcing` carries) when the forcing is +coils, otherwise the run's applied spectrum as a single source named after its forcing file. +Errors when the plasma response was not computed. +""" +function compute_torque_response!( + state::PerturbedEquilibriumState, + ffs::ForceFreeStatesResult, + forcing::Union{ForcingTerms.ForcingTermsControl,ForcingTerms.RMPField}, + intr::PerturbedEquilibriumInternal, + ctrl::PerturbedEquilibriumControl +) + isempty(intr.plasma_response) && + throw(ArgumentError("compute_torque_response needs the plasma response (compute_response = true)")) + ctrl.verbose && @info "Computing torque response matrices" + flux_conform = state.rootarea_to_area_weight .* state.surface_area + + ft_ctrl = forcing isa ForcingTerms.ForcingTermsControl ? forcing : + forcing isa ForcingTerms.RMPSource ? forcing.ctrl : nothing + source_flux = zeros(ComplexF64, ffs.numpert_total, 0) + source_names = String[] + if !isempty(intr.coil_sets) && ft_ctrl !== nothing && ft_ctrl.forcing_data_format == "coil" + source_flux = coil_flux_spectra(ffs, intr.coil_sets, ForcingTerms.CoilConfig(ft_ctrl)) + source_names = [cs.name for cs in intr.coil_sets] + elseif !isempty(state.forcing_b_rootarea) + source_flux = reshape(flux_conform * state.forcing_b_rootarea, :, 1) + source_names = [ft_ctrl === nothing ? "forcing" : basename(ft_ctrl.forcing_data_file)] + end + + state.torque_response = torque_response_matrices(ffs, intr.plasma_response, flux_conform; + forcing_b_rootarea=state.forcing_b_rootarea, source_flux=source_flux, source_names=source_names) + + if ctrl.verbose + total = state.torque_response.T_applied[end] + @info "Torque response complete: applied-forcing torque $(@sprintf("%.4e", real(total))) N·m, 2n·δW $(@sprintf("%.4e", imag(total))) J" + end + return state +end diff --git a/src/PerturbedEquilibrium/Utils.jl b/src/PerturbedEquilibrium/Utils.jl index d2064aee4..e82a62d39 100644 --- a/src/PerturbedEquilibrium/Utils.jl +++ b/src/PerturbedEquilibrium/Utils.jl @@ -101,11 +101,19 @@ PerturbedEquilibrium/ │ ├── left_singular_vectors # U: resonant-field patterns over the retained surfaces [n_retained × rank] │ ├── rational_index # rows of rational_* retained in the ψ_N window [n_retained] │ └── forcing_overlap # Vᴴ·b̃_x, applied forcing's coefficient on each mode [rank] -└── Energies/ - ├── vacuum_energy - ├── surface_energy - ├── plasma_energy - └── toroidal_torque +├── Energies/ +│ ├── vacuum_energy +│ ├── surface_energy +│ ├── plasma_energy +│ └── toroidal_torque +└── TorqueResponse/ # kinetic forward solves with compute_torque_response; b̃ basis, cumulative in ψ + ├── psi # ψ_N of each stored node [npsi] + ├── m, n # (m, n) of each row/column of T_xe [numpert_total] + ├── T_xe # torque response on b̃, N·m/T² (ComplexF64 [numpert_total, numpert_total, npsi]) + ├── source_name # forcing-source names: the coil sets, or the forcing file [nsource] + ├── T_source # torque response in forcing-source space (ComplexF64 [nsource, nsource, npsi]) + ├── T_applied # b̃ₓ†·T_xe(ψ)·b̃ₓ/2: Re = torque inside ψ, Im = 2n·δW inside ψ (ComplexF64 [npsi]) + └── total_torque # Re T_applied at the control surface [N·m] ``` """ function write_outputs_to_HDF5( @@ -247,6 +255,24 @@ function write_outputs_to_HDF5( energy_group["plasma_energy"] = state.plasma_energy energy_group["toroidal_torque"] = state.toroidal_torque + # Torque response matrices (kinetic forward solves only) + tr = state.torque_response + if tr !== nothing + tr_group = haskey(pe_group, "TorqueResponse") ? pe_group["TorqueResponse"] : create_group(pe_group, "TorqueResponse") + tr_group["psi"] = tr.psi + tr_group["m"] = tr.m_modes + tr_group["n"] = tr.n_modes + tr_group["T_xe"] = tr.T_xe + if !isempty(tr.source_names) + tr_group["source_name"] = tr.source_names + tr_group["T_source"] = tr.T_source + end + if !isempty(tr.T_applied) + tr_group["T_applied"] = tr.T_applied + tr_group["total_torque"] = real(tr.T_applied[end]) + end + end + annotate_pe!(pe_group) end end @@ -387,6 +413,26 @@ const PE_H5_ANNOTATIONS = [ "Energies/toroidal_torque" => (; long_name="boundary-response toroidal torque −2n·Im⟨Φ_tot,Λ⁻¹Φ_tot⟩/4: equals the volume-integrated Euler-Lagrange kinetic torque for converged self-consistent solutions; distinct construction from the KineticForces NTV torque", units="N*m" + ), + "TorqueResponse/psi" => (; long_name="normalized poloidal flux ψ_N of each stored node of the solve; the torque response matrices are cumulative inside it", scale="psi"), + "TorqueResponse/m" => (; long_name="poloidal mode number of each row and column of T_xe", scale="m"), + "TorqueResponse/n" => (; long_name="toroidal mode number of each row and column of T_xe", dims=("mode",), attach=(1 => "TorqueResponse/m",)), + "TorqueResponse/T_xe" => + (; + long_name="torque response matrix on the applied root-area-weighted field b̃, cumulative inside ψ_N: b̃†·T_xe·b̃/2 is the complex torque (Re = toroidal torque, Im = 2n·δW ideal+kinetic); the Hermitian part is the torque quadratic form", + units="N*m/T^2", dims=("mode_row", "mode_col", "psi"), attach=(1 => "TorqueResponse/m", 2 => "TorqueResponse/m", 3 => "TorqueResponse/psi")), + "TorqueResponse/source_name" => (; long_name="name of each forcing source of T_source (coil set or forcing file)"), + "TorqueResponse/T_source" => + (; + long_name="torque response matrix in forcing-source space, cumulative inside ψ_N: s†·T_source·s/2 is the complex torque of the sources driven at s times their strength as built", + units="N*m", dims=("source_row", "source_col", "psi"), attach=(3 => "TorqueResponse/psi",)), + "TorqueResponse/T_applied" => + (; + long_name="run's own forcing contracted with the torque response, b̃ₓ†·T_xe(ψ)·b̃ₓ/2: Re = toroidal torque on the plasma inside ψ_N, Im = 2n·δW (ideal + kinetic) of that plasma", + units="N*m", dims=("psi",), attach=(1 => "TorqueResponse/psi",)), + "TorqueResponse/total_torque" => (; + long_name="toroidal torque of the run's own forcing from the torque response at the control surface, Re T_applied(ψ_lim); equals Energies/toroidal_torque", + units="N*m" ) ] diff --git a/test/runtests.jl b/test/runtests.jl index 005d65ff1..c9393875b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -64,6 +64,7 @@ else include("./runtests_kinetic_validity.jl") include("./runtests_multiion.jl") include("./runtests_fullruns.jl") + include("./runtests_torque_response.jl") include("./runtests_coils.jl") include("./runtests_imas.jl") include("./runtests_rerun_from_h5.jl") diff --git a/test/runtests_torque_response.jl b/test/runtests_torque_response.jl new file mode 100644 index 000000000..2e827efbf --- /dev/null +++ b/test/runtests_torque_response.jl @@ -0,0 +1,166 @@ +using Test +using HDF5 +using LinearAlgebra + +# ψ-resolved torque response matrices of a kinetic forward solve: the Solovev calculated-kinetic +# fixture driven by two analytic window-pane coil sets, so T_xe, T_source and the run's own +# forcing are all exercised on one (cheap, mpsi=16) solve. + +const _GPE_TR = GeneralizedPerturbedEquilibrium +const _PE_TR = _GPE_TR.PerturbedEquilibrium + +include("h5_metadata_check.jl") + +const _TORQUE_RESPONSE_SECTIONS = """ + +[ForcingTerms] +forcing_data_format = "coil" # Format: "ascii", "hdf5", or "coil" (Biot-Savart from 3D wires) +mtheta_coil = 128 # Poloidal grid for boundary B·n̂ evaluation +nzeta_coil = 32 # Toroidal grid (0 = auto = 32·n) + +[[ForcingTerms.coil_set]] +name = "pane_upper" # Analytic window-pane array above the midplane +source = "window_pane" # Rectangular picture-frame coils placed by explicit corners +n_coils = 6 # Number of frames around the torus +rz_corners = [[1.5, 0.05], [1.5, 0.45]] # Frame corners (R, Z) [m] +gap_fraction = 0.15 # Toroidal gap between neighbouring frames as a fraction of the pitch +currents = [1000.0, 500.0, -500.0, -1000.0, -500.0, 500.0] # n=1 cosine phasing, 1 kA peak [A] + +[[ForcingTerms.coil_set]] +name = "pane_lower" # Analytic window-pane array below the midplane +source = "window_pane" # Rectangular picture-frame coils placed by explicit corners +n_coils = 6 # Number of frames around the torus +rz_corners = [[1.5, -0.45], [1.5, -0.05]] # Frame corners (R, Z) [m] +gap_fraction = 0.15 # Toroidal gap between neighbouring frames as a fraction of the pitch +currents = [1000.0, 0.0, -1000.0, -1000.0, 0.0, 1000.0] # n=1 phasing shifted by 60 degrees, 1 kA peak [A] + +[PerturbedEquilibrium] +compute_response = true # Compute plasma response to forcing +compute_singular_coupling = false # 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 +""" + +@testset "Torque response matrices" begin + template = joinpath(@__DIR__, "test_data", "regression_solovev_kinetic_calculated") + + mktempdir() do dir + for name in readdir(template) + cp(joinpath(template, name), joinpath(dir, name)) + end + toml = joinpath(dir, "gpec.toml") + write(toml, read(toml, String) * _TORQUE_RESPONSE_SECTIONS) + + res = _GPE_TR.main([dir]) + ffs = res.ffs + pe = res.pe + tr = pe.torque_response + sol = ffs.solution + N = ffs.numpert_total + npsi = sol.step + nn = ffs.nlow + μ₀ = 4π * 1e-7 + + @testset "shapes and mode ordering" begin + @test tr isa _PE_TR.TorqueResponse + @test size(tr.T_xe) == (N, N, npsi) + @test size(tr.T_source) == (2, 2, npsi) + @test tr.source_names == ["pane_upper", "pane_lower"] + @test length(tr.psi) == npsi && tr.psi == sol.psi_store[1:npsi] + @test tr.psi[end] ≈ ffs.psilim + h5open(joinpath(dir, "gpec.h5"), "r") do h5 + mn = read(h5["Info/mn_index"]) + @test tr.m_modes == mn[:, 1] + @test tr.n_modes == mn[:, 2] + g = h5["PerturbedEquilibrium/TorqueResponse"] + @test read(g["m"]) == tr.m_modes + @test read(g["n"]) == tr.n_modes + @test read(g["T_xe"]) == tr.T_xe + @test read(g["T_source"]) == tr.T_source + @test read(g["source_name"]) == tr.source_names + @test read(g["T_applied"]) == tr.T_applied + @test read(g["total_torque"]) == real(tr.T_applied[end]) + @test read(g["psi"]) == tr.psi + bad = _collect_metadata_violations(h5) + isempty(bad) || @error "metadata violations" bad + @test isempty(bad) + end + end + + # The applied forcing, referred back to the solution basis without the function under test: + # Φ_x = R·b̃, boundary ξ = D·P·Φ_x with the flux-space permeability P recovered from the stored b̃ form. + R = pe.rootarea_to_area_weight .* pe.surface_area + P = R * pe.permeability / R + chi1 = 2π * ffs.equil.psio + D = Diagonal([-im / (chi1 * (tr.m_modes[i] - nn * ffs.qlim) * 2π) for i in 1:N]) + b̃ = pe.forcing_b_rootarea + Φx = R * b̃ + ξ = D * P * Φx + U1l = sol.u_store[:, :, 1, npsi] + U2l = sol.u_store[:, :, 2, npsi] + T_expected = dot(ξ, (U2l / U1l) * ξ) * (2 * nn * im / (2μ₀)) / 2 + + # The run's coil sets, rebuilt the way the response stage builds them. + ft_ctrl = _GPE_TR.forcing_terms_control(_GPE_TR.TOML.parsefile(toml)) + cfg = _GPE_TR.ForcingTerms.CoilConfig(ft_ctrl) + coil_sets = _GPE_TR.ForcingTerms.load_coil_sets(cfg, nn; equil=ffs.equil) + + @testset "self-consistency with the kinetic δW of the same solve" begin + T_total = tr.T_applied[end] + @test real(T_total) ≈ real(T_expected) rtol = 1e-8 + @test imag(T_total) ≈ imag(T_expected) rtol = 1e-8 + @test real(T_total) != 0 + @test real(T_total) ≈ pe.toroidal_torque rtol = 1e-6 + @test _PE_TR.torque_profile(tr, b̃) == tr.T_applied + # Cumulative in ψ: every node's matrix is the plasma inside it, so the same + # contraction at an interior node is the torque inside that node. + k = npsi ÷ 2 + ξk = sol.u_store[:, :, 1, k] * (U1l \ ξ) + Tk = dot(ξk, (sol.u_store[:, :, 2, k] / sol.u_store[:, :, 1, k]) * ξk) * (2 * nn * im / (2μ₀)) / 2 + @test tr.T_applied[k] ≈ Tk rtol = 1e-8 + end + + @testset "Hermitian part carries the torque" begin + T_lim = tr.T_xe[:, :, end] + T_h = (T_lim + T_lim') / 2 + T_a = (T_lim - T_lim') / 2 + @test ishermitian(T_h) + @test real(dot(b̃, T_h, b̃)) / 2 ≈ real(tr.T_applied[end]) rtol = 1e-10 + @test abs(real(dot(b̃, T_a, b̃))) < 1e-10 * abs(real(tr.T_applied[end])) + @test imag(dot(b̃, T_a, b̃)) / 2 ≈ imag(tr.T_applied[end]) rtol = 1e-10 + @test all(isreal, eigvals(T_h)) + end + + @testset "forcing-source matrix: coil sets" begin + M = _PE_TR.coil_flux_spectra(ffs, coil_sets, cfg) + @test size(M) == (N, 2) + @test vec(sum(M; dims=2)) ≈ Φx rtol = 1e-10 + M̃ = R \ M + for k in (1, npsi ÷ 2, npsi) + @test tr.T_source[:, :, k] ≈ M̃' * tr.T_xe[:, :, k] * M̃ rtol = 1e-10 + end + @test sum(tr.T_source[:, :, end]) / 2 ≈ tr.T_applied[end] rtol = 1e-10 + end + + @testset "forcing-source matrix: a forcing file is one source" begin + file_state = deepcopy(pe) + file_intr = _PE_TR.PerturbedEquilibriumInternal(; plasma_response=P) + file_ctrl = _GPE_TR.ForcingTerms.ForcingTermsControl(; forcing_data_file="forcing.dat") + _PE_TR.compute_torque_response!(file_state, ffs, file_ctrl, file_intr, _PE_TR.PerturbedEquilibriumControl()) + ts = file_state.torque_response + @test ts.source_names == ["forcing.dat"] + @test size(ts.T_source) == (1, 1, npsi) + @test ts.T_source[1, 1, :] ./ 2 ≈ tr.T_applied rtol = 1e-8 + end + + @testset "argument checks" begin + flux_conform = pe.rootarea_to_area_weight .* pe.surface_area + @test_throws DimensionMismatch _PE_TR.torque_response_matrices(ffs, P, flux_conform; source_flux=zeros(ComplexF64, N, 1)) + @test_throws DimensionMismatch _PE_TR.torque_profile(tr, ComplexF64[1.0]) + empty_state = _PE_TR.PerturbedEquilibriumState() + empty_intr = _PE_TR.PerturbedEquilibriumInternal() + @test_throws ArgumentError _PE_TR.compute_torque_response!(empty_state, ffs, ft_ctrl, empty_intr, _PE_TR.PerturbedEquilibriumControl()) + end + end +end