Perturbed Equilibrium

The PerturbedEquilibrium module computes the plasma response to external magnetic perturbations.

Types

GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.PerturbedEquilibriumControl — Type

PerturbedEquilibriumControl

User-facing control parameters from TOML [PerturbedEquilibrium] section.

Fields

Note: Forcing data file settings are now in [ForcingTerms] section.

High Priority (MWE):

fixed_boundary::Bool - Fixed boundary flag (default: false)
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)
verbose::Bool - Enable verbose logging (default: true)

Output Settings:

output_filename::String - Combined output file with ForceFreeStates results (default: uses ForceFreeStates HDF5_filename)
write_outputs_to_HDF5::Bool - Write perturbed equilibrium outputs to HDF5 (default: true)

Medium Priority (defer for MWE):

filter_modes::Bool - Enable mode filtering (default: false)
singular_point_method::String - Method for singular point treatment (default: "standard")

Regularization:

High Priority (MWE)

reg_spot::Float64 - Regularization width for singular surface smoothing (default: 0.05). Set to 0 to disable. Must be ≥ 0.

source

GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.PerturbedEquilibriumInternal — Type

PerturbedEquilibriumInternal

Internal state variables for perturbed equilibrium calculations.

Fields

dir_path::String - Working directory path
forcing_modes::Vector{ForcingMode} - Loaded forcing mode data
coil_sets::Vector{CoilSet} - Coil geometry used (when forcing_data_format = "coil"); captured for the gpec.h5 rerun snapshot
plasma_response::Matrix{ComplexF64} - Plasma response matrix
singular_coupling_metrics::Dict{String,Float64} - Coupling metrics at singular surfaces
m_modes::Vector{Int} - Poloidal mode numbers for each index i in 1:numpert_total
n_modes::Vector{Int} - Toroidal mode numbers for each index i in 1:numpert_total

source

GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.PerturbedEquilibriumState — Type

PerturbedEquilibriumState

Results from perturbed equilibrium calculations.

Fields

Response fields (mode space):

xi_modes::Union{Nothing, NamedTuple} - Displacement (psi, theta, zeta) [npsi, numpert_total]
b_modes::Union{Nothing, NamedTuple} - Magnetic field; psi=b^ψ, bpsiareaweighted=b^ψ/⟨J·|∇ψ|⟩θ, theta/zeta=unregularized, thetareg/zetareg=regularized [npsi, numpert_total]
b_n_modes::Union{Nothing, Matrix{ComplexF64}} - Physical normal field bn [npsi, numperttotal]
xi_n_modes::Union{Nothing, Matrix{ComplexF64}} - Physical normal displacement xin [npsi, numperttotal]

Coupling matrices [nrational × numperttotal] — one row per resonant (surface, n) pair. Each row maps the full applied field to the resonant response at that surface. Matches Fortran C_f_x_out, C_i_x_out, etc. (shape [modeC, mout]).

C_resonant_area_weighted_field - Φ_r/A^r coupling (resonant area-weighted field b^r in tesla; singcoup row 1) [Pharr 2026]
C_resonant_current - Resonant current coupling (singcoup row 2)
C_island_width_sq - (w/2)² coupling (singcoup row 3)
C_penetrated_area_weighted_field - Penetrated area-weighted field coupling (singcoup row 4)
C_delta_prime - Δ' coupling (singcoup row 5)

Applied resonant vectors [nrational] = C · forcingamplitudes. Matches Fortran Phi_res, w_isl, K_isl, Delta.

resonant_area_weighted_field, resonant_current, island_width_sq, penetrated_area_weighted_field, delta_prime

Diagnostics [n_rational]:

island_half_width::Vector{Float64} - w/2 = sqrt(|islandwidthsq|) from applied forcing
chirikov_parameter::Vector{Float64} - Island overlap metric

Metadata [n_rational] — identifies each (surface, n) row:

rational_psi, rational_q, rational_m_res, rational_n, rational_surface_idx

Dominant resonant-coupling modes — the run summary SVD U·diag(σ)·Vᴴ of C_resonant_area_weighted_field over every resonant row (window post hoc with dominant_coupling on a ResonantCoupling):

dominant_singular_values - σ, descending [rank]
dominant_right_singular_vectors - V, applied-b̃ spectra ranked by resonant drive [numpert_total × rank]
dominant_left_singular_vectors - U, resonant-field patterns over the retained surfaces [n_retained × rank]
dominant_rational_index - rows of the rational_* arrays retained in the window [n_retained]
dominant_forcing_overlap - Vᴴ·b̃_x, the applied forcing's coefficient on each mode [rank]

Control-surface forcing/response spectra [numpert_total], in the three Pharr (2026) field representations (all tesla; no flux/weber is stored):

forcing_b/response_b - bare normal field b (Σ⁻¹·b̃)
forcing_b_rootarea/response_b_rootarea - root-area-weighted field b̃ (coordinate-invariant)
forcing_b_area/response_b_area - area-weighted field b̄ (= S·b̃; flux is Φ = A·b̄)

Control surface matrices [numperttotal × numperttotal], stored in the coordinate-invariant root-area-weighted field (b̃) space (issue #233 / Pharr 2026). Writing S ≡ rootarea_to_area_weight (b̃→b̄) and A ≡ surface_area, the brief internal flux-conform operator is R = S·A (Φ = R·b̃). Recover the area-weighted (b̄) forms by conforming with S (e.g. L_b̄ = S·L̃·S†); recover flux with Φ = A·b̄.

plasma_inductance - Λ̃ = R⁻¹·Λ·R⁻† (wt0-based plasma inductance, congruence)
surface_inductance - L̃ = R⁻¹·L·R⁻† (vacuum surface inductance, congruence)
permeability - P̃ = R⁻¹·P·R (plasma response operator P=Λ·L⁻¹, similarity)
reluctance - ϱ̃ = R†·ϱ·R (ϱ = L⁻¹·(Λ−L)·L⁻¹, congruence)
rootarea_to_area_weight - S = Σ/√A at psilim (b̃→b̄ recovery operator)
surface_area - scalar A = ∫J|∇ψ|dθ at psilim (flux recovery Φ = A·b̄)

Energies (Joules; Fortran gpout convention). Congruence-invariant scalars (energy = Φ†·G⁻¹·Φ = b̃†·G̃⁻¹·b̃), evaluated from the brief internal flux vectors Φx, Φtot with the well-conditioned flux-space inductances L, Λ:

vacuum_energy - Re( ⟨Φx, L⁻¹·Φx⟩ ) / 4 (energy to perturb the vacuum)
surface_energy - Re( ⟨Φtot, L⁻¹·Φtot⟩ ) / 4 (energy at the control surface)
plasma_energy - Σₙ Re( ⟨Φtot,n, Λₙₙ⁻¹·Φtot,n⟩ ) / 4 over the diagonal n blocks of Λ (energy to perturb the plasma; Fortran's "total energy")
toroidal_torque - Σₙ −2n·Im( ⟨Φtot,n, Λₙₙ⁻¹·Φtot,n⟩ / 4 ) [Park 2011 PoP 18 110702, eq. 1] — 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.

source

Functions

GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.compute_perturbed_equilibrium — Function

compute_perturbed_equilibrium(ffs, forcing, ctrl, intr; runtimes=nothing)::PerturbedEquilibriumState

Main entry point for perturbed equilibrium calculations.

Computes plasma response to external forcing and calculates singular layer coupling metrics. Every ForceFreeStates input — the equilibrium, the mode space, the metric and matrix fits, the free-boundary energies and the ξ solution — is read off ffs. Products the producing integrator could not supply gate the corresponding calculation: the step warns and is skipped instead of erroring.

Arguments

ffs: ForceFreeStates.ForceFreeStatesResult from the stability solve
forcing: the external-field description — a ForcingTermsControl (TOML path) or any ForcingTerms.RMPField
ctrl: Control parameters from [PerturbedEquilibrium] section
intr: Internal state variables
runtimes: optional collector; receives "forcing_terms" => dt for the forcing-mode materialization (coil Biot-Savart onto the plasma surface), when that step runs

Returns

PerturbedEquilibriumState: Calculation results

source

GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.write_outputs_to_HDF5 — Function

write_outputs_to_HDF5(
    state::PerturbedEquilibriumState,
    intr::PerturbedEquilibriumInternal,
    filename::String
)

Write perturbed equilibrium results to HDF5 file (appends to existing ForceFreeStates output).

Output Structure

PerturbedEquilibrium/
├── mode_m / mode_n  # (m, n) of each entry along every `mode` axis below [numpert_total]
├── ForcingModes/
│   ├── n              # Toroidal mode numbers
│   ├── m              # Poloidal mode numbers
│   └── amplitude      # ComplexF64 forcing amplitudes
├── forcing_b / forcing_b_root_area / forcing_b_area      # control-surface forcing spectrum (b, b̃, b̄) [numpert_total], tesla
├── response_b / response_b_root_area / response_b_area   # control-surface response spectrum (b, b̃, b̄) [numpert_total], tesla
├── Response/
│   ├── psi             # Radial abscissa ψ_N [npsi] shared by every response profile below
│   ├── xi_psi         # Radial displacement ξ^ψ = ξ·∇ψ (ComplexF64 [npsi, numpert_total])
│   ├── Jxi_psi        # J·ξ^ψ Jacobian-weighted (from gpeq_contra)
│   ├── b_psi_area_weighted       # b^ψ / ⟨J·|∇ψ|⟩_θ area-normalized (ComplexF64 [npsi, numpert_total])
│   ├── b_n            # Physical normal field b_n (ComplexF64 [npsi, numpert_total])
│   ├── xi_n           # Physical normal displacement xi_n (ComplexF64 [npsi, numpert_total])
│   ├── Jb_theta
│   └── Jb_zeta
├── ResponseMatrices/         # [numpert_total × numpert_total], root-area-weighted field (b̃) space; R = S·A
│   ├── plasma_inductance     # Λ̃ = R⁻¹·Λ·R⁻†
│   ├── surface_inductance    # L̃ = R⁻¹·L·R⁻†
│   ├── permeability          # P̃ = R⁻¹·P·R  (P = Λ·L⁻¹)
│   ├── reluctance            # ϱ̃ = R†·ϱ·R
│   ├── rootarea_to_area_weight_operator  # S = Σ/√A at psilim; recover area-weighted field b̄ = S·b̃
│   └── surface_area          # scalar A = ∫J|∇ψ|dθ; recover flux via Φ = A·b̄
├── SingularCoupling/
│   ├── C_resonant_area_weighted_field     # [n_rational × numpert_total] coupling matrix (b̃-space input, resonant area-weighted field b^r=Φ^r/A^r [T])
│   ├── C_resonant_current
│   ├── C_island_width_sq
│   ├── C_penetrated_area_weighted_field
│   ├── C_Delta_prime
│   ├── resonant_area_weighted_field       # [n_rational] applied vector = C̃ · b̃_x (resonant area-weighted field b^r [T])
│   ├── resonant_current
│   ├── island_width_sq
│   ├── penetrated_area_weighted_field
│   ├── Delta_prime
│   ├── island_half_width    # [n_rational] Float64
│   ├── chirikov_parameter
│   ├── rational_psi         # [n_rational] surface metadata
│   ├── rational_q
│   ├── rational_m
│   ├── rational_n
│   └── DominantMode/        # SVD U·diag(σ)·Vᴴ of C_resonant_area_weighted_field over the core-window rational surfaces (ψ_N ≤ 0.9)
│       ├── singular_values          # σ, descending [rank]
│       ├── right_singular_vectors   # V: applied-b̃ spectra ranked by resonant drive [numpert_total × rank]
│       ├── 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

source

Several toroidal modes

With nn_low < nn_high every mode axis in PerturbedEquilibrium/ runs over the full numpert_total = mpert·npert space, with m varying fastest and one block of mpert entries for each n. PerturbedEquilibrium/mode_m and mode_n label each entry. The equilibrium is axisymmetric, so different n never couple:

Λ, L, P and ϱ are block-diagonal in n, and each block equals the matrix of a single-n run.
Energies and torque are sums of the single-n values; the plasma energy Φn†·Λnn⁻¹·Φn/4 and the torque −2n·Im(Φn†·Λnn⁻¹·Φn)/4 use only the diagonal n blocks of Λ: off-n blocks vanish for an axisymmetric equilibrium, so numerical leakage there must not enter.
Profiles in Response/ hold one block of columns for each n.
SingularCoupling/ has one row for each resonant (surface, n) pair, sorted by ψ across all n. An integer-q surface appears once for each n that resonates there. Those rows are not neighbors of each other when the Chirikov parameter is computed.

set_psilim_via_dmlim is ignored for multi-n runs, so set the edge with qhigh or psihigh, away from any rational m/n in the range. KineticForces still takes a single n.

Dominant resonant-coupling mode

The singular-coupling matrix C_resonant_area_weighted_field maps an applied root-area-weighted field spectrum b̃ on the control surface to the resonant field at each rational surface. Its singular-value decomposition over a chosen set of rational surfaces ranks the applied spectra by how strongly they drive resonant field there: the first right singular vector is the dominant mode, the spectrum the plasma is most sensitive to, and the singular values are coordinate-invariant.

The run always writes the full coupling matrix, so the surface window is an analysis choice made afterwards — never a reason to re-run. A ResonantCoupling bundles the matrix with the labels and normalization needed to evaluate arbitrary applied spectra against it, and is built the same way from a finished run in memory or from its gpec.h5:

using GeneralizedPerturbedEquilibrium.PerturbedEquilibrium

rc = ResonantCoupling("gpec.h5")                 # post hoc; or ResonantCoupling(pe_state, ffs) in memory
dom = dominant_coupling(rc; psi_low=0.0, psi_high=0.9)   # SVD over the surfaces in the window

b̃ = rootarea_field(rc, coil_modes)               # unit-norm forcing modes → root-area-weighted field
c = coupling_overlap(dom, b̃)                     # Vᴴ·b̃: c[1] is the overlap with the dominant mode
dom.singular_values[1] * abs(c[1])               # resonant field the dominant mode drives

rootarea_field takes care of the mode ordering and the R⁻¹ conform between the unit-norm convention the forcing loaders and coil integration produce and the b̃ basis the matrix acts on; coupling_overlap takes care of the conjugation. Singular vectors carry an arbitrary global phase, so compare abs of overlaps across runs, not the complex value.

As a summary the run also stores the core-window decomposition (ψ_N ≤ CORE_PSI_HIGH = 0.9) under PerturbedEquilibrium/SingularCoupling/DominantMode/, with forcing_overlap holding the run's own forcing coefficients Vᴴ·b̃_x.

GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.ResonantCoupling — Type

ResonantCoupling

Everything needed to project an applied control-surface field onto the resonant responses of a perturbed-equilibrium solve, detached from the run that produced it. Build it in memory from a PerturbedEquilibriumState and its ForceFreeStatesResult, or post hoc from a gpec.h5; then window and decompose it with dominant_coupling and evaluate coil spectra with rootarea_field and coupling_overlap as often as needed without re-running.

Fields

C::Matrix{ComplexF64} - coupling from the applied root-area-weighted field b̃ to the resonant area-weighted field, one row per resonant (surface, n) pair [n_rational × numpert_total]
rational_psi, rational_q, rational_m, rational_n - row labels [n_rational]
m_modes, n_modes - column labels: the poloidal and toroidal mode number of each entry of an applied spectrum [numpert_total]
flux_conform::Matrix{ComplexF64} - R = S·A, mapping b̃ to the unit-norm flux the forcing loaders and coil integration produce, Φ_x = R·b̃ [numpert_total × numpert_total]

source

GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.DominantCoupling — Type

DominantCoupling

Singular-value decomposition U·diag(σ)·Vᴴ of a resonant coupling matrix over a set of retained rational surfaces, produced by dominant_coupling.

Fields

singular_values - σ, descending [rank]
right_singular_vectors - V: applied-b̃ spectra ranked by resonant drive [numpert_total × rank]
left_singular_vectors - U: resonant-field patterns over the retained surfaces [n_retained × rank]
rational_index - rows of the coupling matrix (and its rational_* labels) that were retained [n_retained]
m_modes, n_modes - the (m, n) basis V is expressed on, copied from the coupling it came from [numpert_total]; empty when built from a bare matrix, which carries no labels

source

GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.dominant_coupling — Function

dominant_coupling(rc::ResonantCoupling; psi_low=CORE_PSI_LOW, psi_high=CORE_PSI_HIGH) -> DominantCoupling
dominant_coupling(C, rational_psi; psi_low=CORE_PSI_LOW, psi_high=CORE_PSI_HIGH) -> DominantCoupling

Singular-value decomposition of the resonant coupling matrix restricted to the rational surfaces with psi_low ≤ ψ_N ≤ psi_high, by default the core window ψ_N ≤ 0.9 (CORE_PSI_HIGH). With C_w = U·diag(σ)·Vᴴ on the retained rows, the right singular vectors V[:, k] are the applied-field spectra ordered by how strongly they drive resonant field inside the window, and the singular values are coordinate-invariant. The overlap of an applied spectrum b̃ with mode k is dot(V[:, k], b̃) — see coupling_overlap — and the resonant field it drives is σ[k] times that coefficient; k = 1 is the dominant mode [Park 2007b]. The window is an analysis choice, so re-evaluate it freely on the same rc.

Throws ArgumentError when the window contains no rational surface.

source

GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.CORE_PSI_HIGH — Constant

Upper edge of the core window, ψ_N = 0.9, the default window of the dominant-mode SVD. Rational surfaces closer to the edge couple strongly to the applied field on any tokamak and would dominate the decomposition, while the locking physics the overlap feeds concerns the core.

source

GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.rootarea_field — Function

rootarea_field(rc::ResonantCoupling, modes) -> Vector{ComplexF64}

Root-area-weighted control-surface field b̃ of an applied spectrum given in the unit-norm (Φx) convention — what the forcing-file loaders and the coil integration produce. modes is a Vector{ForcingMode}, placed on rc's (m, n) column ordering (modes outside the basis are ignored), or an already-ordered amplitude vector. The result is conformed with R⁻¹ (`Φx = R·b̃`) and is the vector the coupling matrix and its singular vectors act on.

source

GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.coupling_overlap — Function

coupling_overlap(dom::DominantCoupling, b̃) -> Vector{ComplexF64}
coupling_overlap(dom, rc::ResonantCoupling, modes) -> Vector{ComplexF64}

Coefficients Vᴴ·b̃ of an applied root-area-weighted field on the singular modes of dom; the first entry is the overlap with the dominant mode.

These are unnormalized, in the units of b̃ itself. The dimensionless overlap δ that the error-field literature quotes divides by the axis toroidal field, and the resonant fraction of a coil's own spectrum divides by ‖b̃‖ instead. ErrorFields.CoilOverlap carries all three together so a caller never has to know which normalization a bare number was in. dom.singular_values .* coupling_overlap(dom, b̃) is the resonant field each mode drives, with dom.left_singular_vectors giving its pattern over the retained surfaces. The second form conforms modes through rootarea_field first.

source

GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.check_mode_basis — Function

check_mode_basis(dom::DominantCoupling, m_modes, n_modes, what::AbstractString)

Error when dom was decomposed on a different (m, n) ordering than the spectrum being projected.

V is only meaningful against the ordering it was decomposed on, and a decomposition from another equilibrium or another mlow can carry the same number of columns while meaning something else — so a length check passes and the projection returns a plausible number. A no-op when dom carries no basis, which is the bare-matrix construction.

source

GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.compute_dominant_coupling! — Function

compute_dominant_coupling!(state, ctrl)

Store the dominant resonant-coupling decomposition of state.C_resonant_area_weighted_field over the core window (ψ_N ≤ CORE_PSI_HIGH) on state as a run summary, with the applied forcing's coefficients on each mode when the response spectrum is available. Other windows are evaluated post hoc on a ResonantCoupling. Returns without change when no coupling matrix exists.

source

Plotting per-surface results against ψ or q

SingularCoupling/ quantities are indexed by rational-surface index, not by q: with multi-n runs a single q value can host several resonances, so the index is the only unambiguous axis. Both rational_psi and rational_q are attached to that axis as HDF5 dimension scales, so plotting against either is direct:

h5open("gpec.h5", "r") do f
    g = f["PerturbedEquilibrium/SingularCoupling"]
    q = read(g["rational_q"])
    b_res = abs.(read(g["resonant_area_weighted_field"]))
    scatter(q, b_res; xlabel="q", ylabel="|b^r| [T]")   # or read(g["rational_psi"]) for ψ_N
end

In Python the same scales are visible through h5py's dimension API (dset.dims[0]["psi_rational"], dset.dims[0]["q_rational"]), so xarray-style tooling can label the axis automatically.