Skip to content

PE - DOCS - Delineate the boundary-response torque from the KineticForces NTV torque #424

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
logan-nc wants to merge 1 commit into
base: develop
Choose a base branch
from
Draft

PE - DOCS - Delineate the boundary-response torque from the KineticForces NTV torque #424

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
29 changes: 19 additions & 10 deletions src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,8 @@ Medium Priority (defer for MWE):
- `singular_point_method::String` - Method for singular point treatment (default: "standard")

Regularization:
# High Priority (MWE)
# High Priority (MWE)

- `reg_spot::Float64` - Regularization width for singular surface smoothing (default: 0.05). Set to 0 to disable. Must be ≥ 0.
"""
@kwdef struct PerturbedEquilibriumControl
Expand Down Expand Up @@ -121,6 +122,7 @@ Metadata [n_rational] — identifies each (surface, n) row:

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̄)
Expand All @@ -145,7 +147,14 @@ 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, Λ⁻¹·Φ_tot⟩ ) / 4 (energy to perturb the plasma; Fortran's "total energy") # Response fields in mode space [npsi, mpert]
- `toroidal_torque` - -2·n·Im( ⟨Φ_tot, Λ⁻¹·Φ_tot⟩ / 4 )
- `toroidal_torque` - -2·n·Im( ⟨Φ_tot, Λ⁻¹·Φ_tot⟩ / 4 ) — the **boundary-response torque**:
the net toroidal torque implied by the anti-Hermitian part of the plasma inductance at the
control surface. For exact self-consistent solutions this equals the volume-integrated
kinetic torque of the same Euler-Lagrange model (the δW surface-term identity), so any
numerical difference from a ψ-resolved torque profile is a solution/grid-quality diagnostic,
not physics. It is a *distinct construction* from the NTV torque under `KineticForces/`
(independent energy-space bounce-average evaluation); agreement between the two tests whether
the EL matrices faithfully discretize the drift-kinetic operator. Zero for ideal (Hermitian) runs.
"""
@kwdef mutable struct PerturbedEquilibriumState
# Radial grid (FFS ODE integration ψ_n values) [npsi]
Expand Down Expand Up @@ -185,18 +194,18 @@ well-conditioned flux-space inductances L, Λ:
rational_surface_idx::Vector{Int} = Int[]

# Control-surface forcing/response spectra in the three weightings of field representations [numpert_total], tesla
forcing_b::Vector{ComplexF64} = ComplexF64[] # bare normal field b (forcing Φ_x)
forcing_b_rootarea::Vector{ComplexF64} = ComplexF64[] # root-area-weighted field b̃ (coordinate-invariant)
forcing_b_area::Vector{ComplexF64} = ComplexF64[] # area-weighted field b̄
response_b::Vector{ComplexF64} = ComplexF64[] # bare normal field b (response Φ_tot = P·Φ_x)
forcing_b::Vector{ComplexF64} = ComplexF64[] # bare normal field b (forcing Φ_x)
forcing_b_rootarea::Vector{ComplexF64} = ComplexF64[] # root-area-weighted field b̃ (coordinate-invariant)
forcing_b_area::Vector{ComplexF64} = ComplexF64[] # area-weighted field b̄
response_b::Vector{ComplexF64} = ComplexF64[] # bare normal field b (response Φ_tot = P·Φ_x)
response_b_rootarea::Vector{ComplexF64} = ComplexF64[] # root-area-weighted field b̃
response_b_area::Vector{ComplexF64} = ComplexF64[] # area-weighted field b̄
response_b_area::Vector{ComplexF64} = ComplexF64[] # area-weighted field b̄

# Control surface matrices [numpert_total × numpert_total], root-area-weighted field (b̃) space
plasma_inductance::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) # Λ̃ (field space)
plasma_inductance::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) # Λ̃ (field space)
surface_inductance::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) # L̃ (field space)
permeability::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) # P̃ = R⁻¹·Λ·L⁻¹·R
reluctance::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) # ϱ̃ = R†·L⁻¹·(Λ−L)·L⁻¹·R
permeability::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) # P̃ = R⁻¹·Λ·L⁻¹·R
reluctance::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) # ϱ̃ = R†·L⁻¹·(Λ−L)·L⁻¹·R
rootarea_to_area_weight::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) # S = Σ/√A at psilim: b̃→b̄ recovery operator
surface_area::Float64 = 0.0 # scalar control-surface area A = ∫J|∇ψ|dθ (flux: Φ = A·b̄; conform R = S·A)

Expand Down
2 changes: 2 additions & 0 deletions src/PerturbedEquilibrium/Response.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -95,6 +95,8 @@ function compute_plasma_response!(
state.vacuum_energy = real(vy)
state.surface_energy = real(sy)
state.plasma_energy = real(py) # Fortran's "total energy" is this pengy
# Boundary-response torque: distinct construction from the KineticForces NTV torque —
# see the PerturbedEquilibriumState docstring for the delineation of GPEC torque outputs.
state.toroidal_torque = -2 * nn * imag(py)

xi_modes, b_modes = reconstruct_physical_fields(
Expand Down
109 changes: 57 additions & 52 deletions src/PerturbedEquilibrium/Utils.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,8 +12,9 @@ avoiding repeated index arithmetic throughout the code.
## Mode Indexing Convention

For linear index i ∈ [1, numpert_total]:
- m_modes[i] = (i-1) % mpert + mlow
- n_modes[i] = (i-1) ÷ mpert + nlow

- m_modes[i] = (i-1) % mpert + mlow
- n_modes[i] = (i-1) ÷ mpert + nlow

This matches the convention used in ForceFreeStates where modes are ordered as:
(m1,n1), (m2,n1), ..., (mpert,n1), (m1,n2), (m2,n2), ..., (mpert,npert)
Expand Down Expand Up @@ -110,119 +111,119 @@ function write_outputs_to_HDF5(

# Forcing modes
forcing_group = haskey(pe_group, "ForcingModes") ? pe_group["ForcingModes"] : create_group(pe_group, "ForcingModes")
forcing_group["n"] = [mode.n for mode in intr.forcing_modes]
forcing_group["m"] = [mode.m for mode in intr.forcing_modes]
forcing_group["n"] = [mode.n for mode in intr.forcing_modes]
forcing_group["m"] = [mode.m for mode in intr.forcing_modes]
forcing_group["amplitude"] = [mode.amplitude for mode in intr.forcing_modes]

# Control-surface forcing/response spectra in the three Pharr field representations
# (all tesla; flux/weber is never stored). b̃ = root-area-weighted (coordinate-invariant).
!isempty(state.forcing_b) && (pe_group["forcing_b"] = state.forcing_b)
!isempty(state.forcing_b_rootarea) && (pe_group["forcing_b_root_area"] = state.forcing_b_rootarea)
!isempty(state.forcing_b_area) && (pe_group["forcing_b_area"] = state.forcing_b_area)
!isempty(state.response_b) && (pe_group["response_b"] = state.response_b)
!isempty(state.forcing_b) && (pe_group["forcing_b"] = state.forcing_b)
!isempty(state.forcing_b_rootarea) && (pe_group["forcing_b_root_area"] = state.forcing_b_rootarea)
!isempty(state.forcing_b_area) && (pe_group["forcing_b_area"] = state.forcing_b_area)
!isempty(state.response_b) && (pe_group["response_b"] = state.response_b)
!isempty(state.response_b_rootarea) && (pe_group["response_b_root_area"] = state.response_b_rootarea)
!isempty(state.response_b_area) && (pe_group["response_b_area"] = state.response_b_area)
!isempty(state.response_b_area) && (pe_group["response_b_area"] = state.response_b_area)

# Control surface matrices [numpert_total × numpert_total], in coordinate-invariant
# root-area-weighted field (b̃) space. Recover the area-weighted field b̄ with the stored
# operator S ≡ rootarea_to_area_weight (b̄ = S·b̃): e.g. L_b̄ = S·L̃·S†; recover flux with the
# scalar surface_area A: Φ = A·b̄ (internally R = S·A, Φ = R·b̃). [Pharr 2026]
mat_group = haskey(pe_group, "ResponseMatrices") ? pe_group["ResponseMatrices"] : create_group(pe_group, "ResponseMatrices")
!isempty(state.plasma_inductance) && (mat_group["plasma_inductance"] = state.plasma_inductance)
!isempty(state.plasma_inductance) && (mat_group["plasma_inductance"] = state.plasma_inductance)
!isempty(state.surface_inductance) && (mat_group["surface_inductance"] = state.surface_inductance)
!isempty(state.permeability) && (mat_group["permeability"] = state.permeability)
!isempty(state.reluctance) && (mat_group["reluctance"] = state.reluctance)
!isempty(state.permeability) && (mat_group["permeability"] = state.permeability)
!isempty(state.reluctance) && (mat_group["reluctance"] = state.reluctance)
!isempty(state.rootarea_to_area_weight) && (mat_group["rootarea_to_area_weight_operator"] = state.rootarea_to_area_weight)
(state.surface_area != 0.0) && (mat_group["surface_area"] = state.surface_area)
(state.surface_area != 0.0) && (mat_group["surface_area"] = state.surface_area)

# Response fields (ComplexF64 directly)
response_group = haskey(pe_group, "Response") ? pe_group["Response"] : create_group(pe_group, "Response")
!isempty(state.psi_grid) && (response_group["psi"] = state.psi_grid)
have_xi = !isnothing(state.xi_modes)
have_b = have_xi && !isnothing(state.b_modes)
response_group["xi_psi"] = have_xi ? state.xi_modes.psi : ComplexF64[]
response_group["b_psi_area_weighted"] = have_b ? state.b_modes.b_psi_area_weighted : ComplexF64[]
response_group["Jb_theta"] = have_b ? state.b_modes.theta : ComplexF64[]
response_group["Jb_zeta"] = have_b ? state.b_modes.zeta : ComplexF64[]
response_group["b_n"] = !isnothing(state.b_n_modes) ? state.b_n_modes : ComplexF64[]
response_group["xi_n"] = !isnothing(state.xi_n_modes) ? state.xi_n_modes : ComplexF64[]
have_b = have_xi && !isnothing(state.b_modes)
response_group["xi_psi"] = have_xi ? state.xi_modes.psi : ComplexF64[]
response_group["b_psi_area_weighted"] = have_b ? state.b_modes.b_psi_area_weighted : ComplexF64[]
response_group["Jb_theta"] = have_b ? state.b_modes.theta : ComplexF64[]
response_group["Jb_zeta"] = have_b ? state.b_modes.zeta : ComplexF64[]
response_group["b_n"] = !isnothing(state.b_n_modes) ? state.b_n_modes : ComplexF64[]
response_group["xi_n"] = !isnothing(state.xi_n_modes) ? state.xi_n_modes : ComplexF64[]

# Clebsch displacements for PENTRC (matches Fortran gpout_xclebsch)
if have_xi
response_group["xi_clebsch_psi"] = state.xi_modes.clebsch_psi
response_group["dxi_clebsch_psidpsi"] = state.xi_modes.clebsch_psi1
response_group["xi_clebsch_psi"] = state.xi_modes.clebsch_psi
response_group["dxi_clebsch_psidpsi"] = state.xi_modes.clebsch_psi1
response_group["xi_clebsch_alpha"] = state.xi_modes.clebsch_alpha
end

# Contravariant displacement (from gpeq_contra, all J-weighted)
if have_xi
response_group["Jxi_psi"] = state.xi_modes.psi_J
response_group["Jxi_theta"] = state.xi_modes.theta
response_group["Jxi_zeta"] = state.xi_modes.zeta
response_group["Jxi_zeta"] = state.xi_modes.zeta
end

# Covariant components (from gpeq_cova)
if have_xi
response_group["xi_cov_psi"] = state.xi_modes.cova_psi
response_group["xi_cov_psi"] = state.xi_modes.cova_psi
response_group["xi_cov_theta"] = state.xi_modes.cova_theta
response_group["xi_cov_zeta"] = state.xi_modes.cova_zeta
response_group["xi_cov_zeta"] = state.xi_modes.cova_zeta
end
if have_xi
response_group["Jxi_theta_reg"] = state.xi_modes.theta_reg
response_group["Jxi_zeta_reg"] = state.xi_modes.zeta_reg
response_group["Jxi_zeta_reg"] = state.xi_modes.zeta_reg
end
if have_b
response_group["Jb_theta_reg"] = state.b_modes.theta_reg
response_group["Jb_zeta_reg"] = state.b_modes.zeta_reg
response_group["b_cov_psi"] = state.b_modes.cova_psi
response_group["Jb_zeta_reg"] = state.b_modes.zeta_reg
response_group["b_cov_psi"] = state.b_modes.cova_psi
response_group["b_cov_theta"] = state.b_modes.cova_theta
response_group["b_cov_zeta"] = state.b_modes.cova_zeta
response_group["b_cov_zeta"] = state.b_modes.cova_zeta
end

# R,Z,φ cylindrical components in mode-space (from gpeq_rzphi)
if have_xi
response_group["xi_R"] = state.xi_modes.R
response_group["xi_Z"] = state.xi_modes.Z
response_group["xi_R"] = state.xi_modes.R
response_group["xi_Z"] = state.xi_modes.Z
response_group["xi_phi"] = state.xi_modes.phi
end
if have_b
response_group["b_R"] = state.b_modes.R
response_group["b_Z"] = state.b_modes.Z
response_group["b_R"] = state.b_modes.R
response_group["b_Z"] = state.b_modes.Z
response_group["b_phi"] = state.b_modes.phi
end

# Singular coupling
coupling_group = haskey(pe_group, "SingularCoupling") ? pe_group["SingularCoupling"] : create_group(pe_group, "SingularCoupling")

# Coupling matrices [n_rational × numpert_total]
!isempty(state.C_resonant_area_weighted_field) && (coupling_group["C_resonant_area_weighted_field"] = state.C_resonant_area_weighted_field)
!isempty(state.C_resonant_area_weighted_field) && (coupling_group["C_resonant_area_weighted_field"] = state.C_resonant_area_weighted_field)
!isempty(state.C_resonant_current) && (coupling_group["C_resonant_current"] = state.C_resonant_current)
!isempty(state.C_island_width_sq) && (coupling_group["C_island_width_sq"] = state.C_island_width_sq)
!isempty(state.C_island_width_sq) && (coupling_group["C_island_width_sq"] = state.C_island_width_sq)
!isempty(state.C_penetrated_area_weighted_field) && (coupling_group["C_penetrated_area_weighted_field"] = state.C_penetrated_area_weighted_field)
!isempty(state.C_delta_prime) && (coupling_group["C_Delta_prime"] = state.C_delta_prime)
!isempty(state.C_delta_prime) && (coupling_group["C_Delta_prime"] = state.C_delta_prime)

# Applied resonant vectors [n_rational]
!isempty(state.resonant_area_weighted_field) && (coupling_group["resonant_area_weighted_field"] = state.resonant_area_weighted_field)
!isempty(state.resonant_current) && (coupling_group["resonant_current"] = state.resonant_current)
!isempty(state.island_width_sq) && (coupling_group["island_width_sq"] = state.island_width_sq)
!isempty(state.penetrated_area_weighted_field) && (coupling_group["penetrated_area_weighted_field"] = state.penetrated_area_weighted_field)
!isempty(state.delta_prime) && (coupling_group["Delta_prime"] = state.delta_prime)
!isempty(state.forcing_solution_weights) && (coupling_group["forcing_solution_weights"] = state.forcing_solution_weights)
!isempty(state.resonant_area_weighted_field) && (coupling_group["resonant_area_weighted_field"] = state.resonant_area_weighted_field)
!isempty(state.resonant_current) && (coupling_group["resonant_current"] = state.resonant_current)
!isempty(state.island_width_sq) && (coupling_group["island_width_sq"] = state.island_width_sq)
!isempty(state.penetrated_area_weighted_field) && (coupling_group["penetrated_area_weighted_field"] = state.penetrated_area_weighted_field)
!isempty(state.delta_prime) && (coupling_group["Delta_prime"] = state.delta_prime)
!isempty(state.forcing_solution_weights) && (coupling_group["forcing_solution_weights"] = state.forcing_solution_weights)
!isempty(state.rational_area) && (coupling_group["rational_area"] = state.rational_area)
!isempty(state.island_half_width) && (coupling_group["island_half_width"] = state.island_half_width)
!isempty(state.island_half_width) && (coupling_group["island_half_width"] = state.island_half_width)
!isempty(state.chirikov_parameter) && (coupling_group["chirikov_parameter"] = state.chirikov_parameter)

# Metadata [n_rational]
!isempty(state.rational_psi) && (coupling_group["rational_psi"] = state.rational_psi)
!isempty(state.rational_q) && (coupling_group["rational_q"] = state.rational_q)
!isempty(state.rational_m_res) && (coupling_group["rational_m"] = state.rational_m_res)
!isempty(state.rational_n) && (coupling_group["rational_n"] = state.rational_n)
!isempty(state.rational_psi) && (coupling_group["rational_psi"] = state.rational_psi)
!isempty(state.rational_q) && (coupling_group["rational_q"] = state.rational_q)
!isempty(state.rational_m_res) && (coupling_group["rational_m"] = state.rational_m_res)
!isempty(state.rational_n) && (coupling_group["rational_n"] = state.rational_n)

# Energies
energy_group = haskey(pe_group, "Energies") ? pe_group["Energies"] : create_group(pe_group, "Energies")
energy_group["vacuum_energy"] = state.vacuum_energy
energy_group["surface_energy"] = state.surface_energy
energy_group["plasma_energy"] = state.plasma_energy
energy_group["vacuum_energy"] = state.vacuum_energy
energy_group["surface_energy"] = state.surface_energy
energy_group["plasma_energy"] = state.plasma_energy
energy_group["toroidal_torque"] = state.toroidal_torque

annotate_pe!(pe_group)
Expand Down Expand Up @@ -261,7 +262,8 @@ const PE_H5_ANNOTATIONS = [
"Response/xi_cov_psi" => (; long_name="covariant radial displacement ξ_ψ = ξ·e_ψ", units="m^2", dims=("psi", "mode"), attach=(1 => "Response/psi",)),
"Response/xi_cov_theta" => (; long_name="covariant poloidal displacement ξ_θ = ξ·e_θ", units="m^2", dims=("psi", "mode"), attach=(1 => "Response/psi",)),
"Response/xi_cov_zeta" => (; long_name="covariant toroidal displacement ξ_ζ = ξ·e_ζ", units="m^2", dims=("psi", "mode"), attach=(1 => "Response/psi",)),
"Response/xi_clebsch_psi" => (; long_name="Clebsch displacement component ξ^ψ (PENTRC input, gpout_xclebsch convention)", dims=("psi", "mode"), attach=(1 => "Response/psi",)),
"Response/xi_clebsch_psi" =>
(; long_name="Clebsch displacement component ξ^ψ (PENTRC input, gpout_xclebsch convention)", dims=("psi", "mode"), attach=(1 => "Response/psi",)),
"Response/dxi_clebsch_psidpsi" =>
(; long_name="regularized ψ_N derivative of ξ^ψ (× singfac²/(singfac²+reg_spot²))", dims=("psi", "mode"), attach=(1 => "Response/psi",)),
"Response/xi_clebsch_alpha" =>
Expand Down Expand Up @@ -345,7 +347,10 @@ const PE_H5_ANNOTATIONS = [
"Energies/vacuum_energy" => (; long_name="perturbed vacuum energy", units="J"),
"Energies/surface_energy" => (; long_name="perturbed surface energy", units="J"),
"Energies/plasma_energy" => (; long_name="perturbed plasma energy", units="J"),
"Energies/toroidal_torque" => (; long_name="net toroidal torque on the plasma", units="N*m")
"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"
)
]

# Attach long_name/units/dims + dimension scales (declared in-table) to the
Expand Down
Loading