diff --git a/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl b/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl index d8251de38..70af8f5ac 100644 --- a/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl +++ b/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl @@ -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 @@ -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̄) @@ -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] @@ -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) diff --git a/src/PerturbedEquilibrium/Response.jl b/src/PerturbedEquilibrium/Response.jl index d305a079b..594f3b866 100644 --- a/src/PerturbedEquilibrium/Response.jl +++ b/src/PerturbedEquilibrium/Response.jl @@ -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( diff --git a/src/PerturbedEquilibrium/Utils.jl b/src/PerturbedEquilibrium/Utils.jl index fc1eb9c0a..8707206d1 100644 --- a/src/PerturbedEquilibrium/Utils.jl +++ b/src/PerturbedEquilibrium/Utils.jl @@ -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) @@ -110,47 +111,47 @@ 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 @@ -158,36 +159,36 @@ function write_outputs_to_HDF5( 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 @@ -195,34 +196,34 @@ function write_outputs_to_HDF5( 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) @@ -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" => @@ -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