Conform PerturbedEquilibrium response & coupling matrices to coordinate-invariant field space (#233, first slice)

CLAUDE.md

@@ -61,6 +61,11 @@ The PerturbedEquilibrium module implements GPEC-style perturbed equilibrium calc
   - Published: Physical Review Letters **99**, 195003 (2007)
   - Describes: Plasma response to resonant magnetic perturbations (RMP)
 
+- **Park et al. (2008)**: "Spectral asymmetry due to magnetic coordinates"
+  - Location: `docs/resources/2008-Park-Spectral_asymmetry_due_to_magnetic_coordinates.pdf`
+  - Published: Physics of Plasmas **15**, 064501 (2008)
+  - Describes: Coordinate-dependence of the perturbed-field Fourier spectrum, motivating area-normalization of the resonant harmonic to obtain the coordinate-invariant resonant field
+
 - **Park et al. (2009)**: "Importance of plasma response to nonaxisymmetric perturbations in tokamaks"
   - Location: `docs/resources/2009-Park-Importance_of_plasma_response_to_nonaxisymmetric_perturbations_in_tokamaks-compressed.pdf`
   - Published: Physics of Plasmas **16**, 056115 (2009)

benchmarks/benchmark_against_fortran_run.jl

@@ -327,17 +327,22 @@ function load_julia_outputs(h5_path::String)
         julia["rational_q"]         = haskey(f, "$sc/rational_q")         ? read(f, "$sc/rational_q")         : Float64[]
         julia["rational_n"]         = haskey(f, "$sc/rational_n")         ? read(f, "$sc/rational_n")         : Int[]
         julia["rational_m_res"]     = haskey(f, "$sc/rational_m_res")     ? read(f, "$sc/rational_m_res")     : Int[]
-        julia["resonant_flux"]      = haskey(f, "$sc/resonant_flux")      ? read(f, "$sc/resonant_flux")      : ComplexF64[]
+        julia["resonant_area_weighted_field"]     = haskey(f, "$sc/resonant_area_weighted_field")     ? read(f, "$sc/resonant_area_weighted_field")     : ComplexF64[]
         julia["resonant_current"]   = haskey(f, "$sc/resonant_current")   ? read(f, "$sc/resonant_current")   : ComplexF64[]
         julia["island_half_width"]  = haskey(f, "$sc/island_half_width")  ? read(f, "$sc/island_half_width")  : Float64[]
         julia["chirikov_parameter"] = haskey(f, "$sc/chirikov_parameter") ? read(f, "$sc/chirikov_parameter") : Float64[]
         julia["delta_prime"]        = haskey(f, "$sc/delta_prime")        ? read(f, "$sc/delta_prime")        : ComplexF64[]
 
         pe = "perturbed_equilibrium"
-        julia["forcing_vec"]  = haskey(f, "$pe/forcing_vec")  ? read(f, "$pe/forcing_vec")  : ComplexF64[]
-        julia["response_vec"] = haskey(f, "$pe/response_vec") ? read(f, "$pe/response_vec") : ComplexF64[]
+        # GPEC stores the control-surface spectra as fields (tesla), not flux. Reconstruct the
+        # poloidal flux for the Fortran Phi_x/Phi_tot comparison as Φ = A·b̄ (area-weighted field).
+        A_surf  = haskey(f, "$pe/response_matrices/surface_area") ? read(f, "$pe/response_matrices/surface_area") : 1.0
+        fb_area = haskey(f, "$pe/forcing_b_area")  ? read(f, "$pe/forcing_b_area")  : ComplexF64[]
+        rb_area = haskey(f, "$pe/response_b_area") ? read(f, "$pe/response_b_area") : ComplexF64[]
+        julia["forcing_vec"]  = isempty(fb_area) ? ComplexF64[] : A_surf .* fb_area   # Φ_x   = A·b̄_x
+        julia["response_vec"] = isempty(rb_area) ? ComplexF64[] : A_surf .* rb_area   # Φ_tot = A·b̄_tot
         julia["b_n"]          = haskey(f, "$pe/response/b_n")        ? read(f, "$pe/response/b_n")        : Matrix{ComplexF64}(undef, 0, 0)
-        julia["Jbgradpsi"]    = haskey(f, "$pe/response/psi_area")    ? read(f, "$pe/response/psi_area")    : Matrix{ComplexF64}(undef, 0, 0)
+        julia["Jbgradpsi"]    = haskey(f, "$pe/response/b_psi_area_weighted") ? read(f, "$pe/response/b_psi_area_weighted") : Matrix{ComplexF64}(undef, 0, 0)
         julia["xi_psi"]       = haskey(f, "$pe/response/xi_psi")     ? read(f, "$pe/response/xi_psi")     : Matrix{ComplexF64}(undef, 0, 0)
         julia["xi_n"]         = haskey(f, "$pe/response/xi_n")       ? read(f, "$pe/response/xi_n")       : Matrix{ComplexF64}(undef, 0, 0)
         julia["clebsch_psi1"] = haskey(f, "$pe/response/clebsch_psi1")  ? read(f, "$pe/response/clebsch_psi1")  : Matrix{ComplexF64}(undef, 0, 0)
@@ -754,7 +759,7 @@ function build_comparison_table(fort, julia, fortran_dir, bench_dir, nn)
         row   = find_julia_row(julia, q_int, nn)
 
         jp_r  = (row > 0 && !isempty(julia["rational_psi"]))     ? julia["rational_psi"][row]        : NaN
-        j_phi = (row > 0 && !isempty(julia["resonant_flux"]))    ? abs(julia["resonant_flux"][row])   : NaN
+        j_phi = (row > 0 && !isempty(julia["resonant_area_weighted_field"]))   ? abs(julia["resonant_area_weighted_field"][row])  : NaN
         j_wisl= (row > 0 && !isempty(julia["island_half_width"])) ? 2*julia["island_half_width"][row] : NaN
         j_kch = (row > 0 && !isempty(julia["chirikov_parameter"])) ? julia["chirikov_parameter"][row]  : NaN
 
@@ -974,7 +979,7 @@ function generate_plots(fort, julia, bench_dir, nn)
     # Phi_res (resonant flux / area)
     f_phires_vals = isempty(fort["Phi_res"]) ? Float64[] : abs.(fort["Phi_res"])
     j_phires_psi  = julia["rational_psi"]
-    j_phires_vals = isempty(julia["resonant_flux"]) ? Float64[] : abs.(julia["resonant_flux"])
+    j_phires_vals = isempty(julia["resonant_area_weighted_field"]) ? Float64[] : abs.(julia["resonant_area_weighted_field"])
     p9 = _singcoup_panel("|Φ_res| [T]", "|Φ_res| vs ψ  (n=$nn)",
                           f_psi_rat, f_phires_vals, f_q_labels,
                           j_phires_psi, j_phires_vals)

docs/src/citations.md

@@ -66,6 +66,14 @@ Establishes the self-consistent plasma response calculation. Derives the permeab
 
 ---
 
+> J.-K. Park, A. H. Boozer, and J. E. Menard, "Spectral asymmetry due to magnetic coordinates,"
+> *Physics of Plasmas* **15**, 064501 (2008).
+> DOI: [10.1063/1.2932110](https://doi.org/10.1063/1.2932110)
+
+Shows that the Fourier spectrum of a perturbed field on a flux surface is asymmetric and depends on the magnetic working coordinate, so the raw spectrum is not by itself a physical measure of the field. Motivates normalizing the resonant harmonic by the scalar surface area to obtain the coordinate-independent resonant field — see the [Conventions Reference](conventions.md).
+
+---
+
 > M. Pharr, E. Bursch, N. C. Logan, P. Lunia, J.-K. Park, and C. Paz-Soldan, "Coordinate-invariant flux-surface Fourier analysis in tokamaks,"
 > Preprint (2026).
 > [arXiv:2606.02901](https://arxiv.org/abs/2606.02901)
@@ -163,4 +171,4 @@ Provides the complete theory and implementation details, including NTV in the pr
 > *Physics of Plasmas* **24**, 032505 (2017).
 > DOI: [10.1063/1.4977898](https://doi.org/10.1063/1.4977898)
 
-Describes the fully self-consistent coupling between the perturbed equilibrium and neoclassical toroidal viscosity (NTV), providing the theoretical foundation for the the KineticForces functionality.
+Describes the fully self-consistent coupling between the perturbed equilibrium and neoclassical toroidal viscosity (NTV), providing the theoretical foundation for the KineticForces functionality.

docs/src/conventions.md

@@ -57,7 +57,7 @@ kernel ``\exp(-im\theta)``; the inverse transform reconstructs with ``\exp(+im\t
 
 ### Toroidal Coordinate ``\zeta`` and ``\phi``
 
-- The magnetic cooridante toroidal angle is ``\zeta = \phi/(2\pi) + \nu(\psi,\theta)``, where ``\nu`` is
+- The magnetic coordinate toroidal angle is ``\zeta = \phi/(2\pi) + \nu(\psi,\theta)``, where ``\nu`` is
   a single-valued straight-field-line offset that depends on the working coordinate. PEST
   coordinates have ``\nu = 0``.
 - The physical toroidal angle is reconstructed as
@@ -219,7 +219,11 @@ b^r = \frac{\Phi^r}{A^r} \qquad [\mathrm{T}].
 
 Dividing the resonant flux ``\Phi^r`` by the surface area turns it into a field amplitude in Tesla
 that is invariant under changes of the poloidal-angle (working) coordinate. This is the quantity
-reported as `resonant_flux`.
+reported as `resonant_flux`. The need for this normalization is the central point of Park, Boozer,
+and Menard (2008) (see [Citations](citations.md)): the raw Fourier spectrum of the perturbed field is
+*spectrally asymmetric* and changes with the magnetic working coordinate, so only the resonant
+harmonic scaled by the scalar surface area ``A^r`` is a coordinate-independent physical measure of
+the resonant field.
 
 Note that there is strictly NO resonant field in ideal MHD perturbed equilibrium calculations, and
 only the effective resonant field is non-zero in these cases. Effective resonant fields are computed
@@ -249,6 +253,74 @@ The ``\sqrt{A}`` weighting is the unique one for which the spectrum transforms u
 working coordinates, so the power-normalized amplitudes — and quantities derived from them, such as
 the singular values of the resonant coupling matrix — are coordinate-invariant.
 
+### The Three Field Amplitudes
+
+Following Pharr (2026), *Coordinate-invariant flux-surface Fourier analysis in tokamaks* (see
+[Citations](citations.md)), GPEC works with **three** Fourier decompositions of the normal field, **all
+in field units (tesla)**, distinguished only by the **area weight applied to the Fourier integrand**.
+These are the authoritative names used in prose, in code, and in the HDF5 output:
+
+| Pharr symbol | FT integrand weight ``W`` | English name | HDF5 token |
+|---|---|---|---|
+| ``b`` (bare) | ``1`` | **normal field** | `b_n` |
+| ``\bar b`` (bar) | ``\mathcal{J}\,\lvert\nabla\psi\rvert`` | **area-weighted field** | `area` |
+| ``\tilde b`` (tilde) | ``\sqrt{\mathcal{J}\,\lvert\nabla\psi\rvert}`` | **root-area-weighted field** | `root_area` |
+
+All three carry units of tesla because each integral is divided by the appropriate power of the
+coordinate-invariant scalar surface area ``A``:
+
+```math
+b_m = \oint (\mathbf{b}\cdot\hat{\mathbf n})\, e^{-i m\theta}\, d\theta, \qquad
+\bar b_m = \frac{1}{A}\oint \mathcal{J}\lvert\nabla\psi\rvert\,(\mathbf{b}\cdot\hat{\mathbf n})\, e^{-i m\theta}\, d\theta, \qquad
+\tilde b_m = \frac{1}{\sqrt{A}}\oint \sqrt{\mathcal{J}\lvert\nabla\psi\rvert}\,(\mathbf{b}\cdot\hat{\mathbf n})\, e^{-i m\theta}\, d\theta .
+```
+
+Only the **square-root-area weighting** ``\tilde b`` transforms unitarily between working coordinates, so
+its 2-norm (and quantities derived from it, such as the singular values of the resonant coupling matrix)
+are coordinate-invariant. Note that the "root-area-weighted field" sometimes appears as "power-normalized
+field" in early GPEC literature, as the latter names a Parseval *property* of ``\tilde b``. The **full-area weighting** ``\bar b`` is coordinate-invariant for the
+pitch-resonant ``m = nq`` harmonic on a rational surface — the **resonant area-weighted field**
+``\bar b^{\,r} = \Phi^r/A^r`` as described above.
+
+GPEC operates and outputs **only in these field representations — poloidal flux ``\Phi`` (weber) is never
+stored.** Flux appears, briefly, only internally when a user supplies flux-valued forcing, and is
+recovered on demand as the scalar product ``\Phi = A\,\bar b``.
+
+### Translation Operators
+
+With ``\Sigma`` ≡ `sqrtamat` (the mode-space ``\sqrt{}``weight convolution) and the scalar surface area
+``A`` ≡ `jarea`, the three fields are related by (`Equilibrium/CoordinateInvariant.jl`):
+
+```math
+\tilde b = \Sigma\, b, \qquad
+\bar b = (\Sigma/\sqrt{A})\,\tilde b, \qquad
+\Phi = A\,\bar b = \Sigma\sqrt{A}\;\tilde b .
+```
+
+- **`rootarea_to_area_weight`** ``= \Sigma/\sqrt{A}`` maps ``\tilde b \to \bar b``.
+- **`area_to_rootarea_weight`** ``= \sqrt{A}\,\Sigma^{-1}`` is its inverse (``\bar b \to \tilde b``).
+
+The internal flux-conform operator is just ``R = \Sigma\sqrt{A} = `` `rootarea_to_area_weight` ``\cdot A``.
+
+### Field Representations in GPEC Output
+
+- **`forcing_b` / `forcing_b_root_area` / `forcing_b_area`** (and the `response_*` triplet) — the
+  control-surface forcing and response spectra in the bare (``b``), root-area-weighted (``\tilde b``) and
+  area-weighted (``\bar b``) representations, under `perturbed_equilibrium/`.
+- **`b_n`** — the bare normal field ``\mathbf{b}\cdot\hat{\mathbf n}`` (and the area-weighted radial field
+  `b_psi_area_weighted`), under `perturbed_equilibrium/response/`.
+- **`resonant_area_weighted_field`** / **`C_resonant_area_weighted_field`** — the resonant area-weighted
+  field ``\bar b^{\,r} = \Phi^r/A^r`` and its coupling matrix, under
+  `perturbed_equilibrium/singular_coupling/`. The sibling `penetrated_area_weighted_field` follows the
+  same convention.
+- **Root-area-weighted (``\tilde b``) space** — the control-surface response matrices (`permeability`,
+  `reluctance`, `plasma_inductance`, `surface_inductance`) are stored in this coordinate-invariant space
+  under `perturbed_equilibrium/response_matrices/`. The stored `rootarea_to_area_weight_operator` ``S``
+  recovers the area-weighted field forms (``L_{\bar b} = S\,\tilde L\,S^\dagger``) and the scalar
+  `surface_area` ``A`` recovers flux (``\Phi = A\,\bar b``). The coordinate-invariant ideal-MHD energies
+  written by the stability stage (`FreeBoundaryStability/eigenmode_energies`) are the ``\tilde b``
+  quadratic form scaled by the scalar ``c = A``: ``\mathrm{d}W = c\,\tilde b^\dagger W_t\,\tilde b``.
+
 
 ## Rotation Velocity Conventions (KineticForces)
 

docs/src/equilibrium.md

@@ -12,7 +12,7 @@ Key responsibilities of the module:
 
 - Read equilibrium input files and TOML configuration (see `EquilibriumConfig`).
 - Provide convenient constructors for analytic / model equilibria (Large
-	Aspect Ratio, Solev'ev).
+	Aspect Ratio, Solovev).
 - Build spline representations used throughout the code (1D cubic and
 	2D bicubic splines).
 - Run the direct or inverse equilibrium solver and post-process results
@@ -89,7 +89,7 @@ Modules = [GeneralizedPerturbedEquilibrium.Equilibrium]
 - `PlasmaEquilibrium` — the runtime structure containing spline fields,
 	geometry, profiles, and computed diagnostics (q-profile, separatrix,
 	etc.).
-- `LargeAspectRatioConfig`, `SolevevConfig` — convenience structs to
+- `LargeAspectRatioConfig`, `SolovevConfig` — convenience structs to
 	construct analytic/model equilibria when using `eq_type = "lar"` or
 	`eq_type = "sol"`.
 
@@ -134,7 +134,7 @@ using GeneralizedPerturbedEquilibrium
 larcfg = GeneralizedPerturbedEquilibrium.Equilibrium.LargeAspectRatioConfig(lar_r0=10.0, lar_a=1.0, beta0=1e-3)
 pe = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(control=Dict("eq_filename"=>"unused","eq_type"=>"lar")), larcfg)
 
-println("Built LAR equilibrium with a = ", lorcfg.lar_a)
+println("Built LAR equilibrium with a = ", larcfg.lar_a)
 ```
 
 ## Notes and Caveats

docs/src/set_up.md

@@ -46,7 +46,7 @@
         sudo mv julia-1.11.3 /opt/
         ```
 
-        ☆ Ensure these commands match the tarball you installed. These commands match the above tarball and might need to be modified for you installation.
+        ☆ Ensure these commands match the tarball you installed. These commands match the above tarball and might need to be modified for your installation.
 
     3. Add Julia to PATH:
 

regression-harness/cases/diiid_n1.toml

@@ -8,36 +8,36 @@ name = "diiid_n1"
 description = "DIII-D-like equilibrium, n=1, ideal + perturbed equilibrium"
 example_dir = "examples/DIIID-like_ideal_example"
 
-# Energies — leading eigenvalues at the final truncation (psilim). eigenmode_* are now power-normalized (Φ-space, Jacobian-invariant); ξ-space counterparts live under FreeBoundaryStability/XiNorm/ and are not tracked.
+# Energies — leading eigenvalues at the final truncation (psilim). eigenmode_* are now root-area-weighted (Φ-space, Jacobian-invariant); ξ-space counterparts live under FreeBoundaryStability/XiNorm/ and are not tracked.
 [quantities.et_real]
 h5path = "FreeBoundaryStability/eigenmode_energies"
 type = "complex_vector"
 extract = "real_first"
-label = "Power Normalized total energy Re(et[1])"
+label = "root-area-weighted total energy Re(et[1])"
 noise_threshold = 1e-10
 order = 10
 
 [quantities.et_imag]
 h5path = "FreeBoundaryStability/eigenmode_energies"
 type = "complex_vector"
 extract = "imag_first"
-label = "Power Normalized total energy Im(et[1])"
+label = "root-area-weighted total energy Im(et[1])"
 noise_threshold = 1e-10
 order = 11
 
 [quantities.ep_real]
 h5path = "FreeBoundaryStability/eigenmode_plasma_energies"
 type = "complex_vector"
 extract = "real_first"
-label = "Power Normalized plasma energy Re(ep[1])"
+label = "root-area-weighted plasma energy Re(ep[1])"
 noise_threshold = 1e-10
 order = 12
 
 [quantities.ev_real]
 h5path = "FreeBoundaryStability/eigenmode_vacuum_energies"
 type = "complex_vector"
 extract = "real_first"
-label = "Power Normalized vacuum energy Re(ev[1])"
+label = "root-area-weighted vacuum energy Re(ev[1])"
 noise_threshold = 1e-10
 order = 13
 
@@ -54,23 +54,23 @@ order = 14
 h5path = "FreeBoundaryStability/eigenmode_plasma_energies"
 type = "complex_vector"
 extract = "all_complex"
-label = "Power Normalized plasma energy (all)"
+label = "root-area-weighted plasma energy (all)"
 noise_threshold = 1e-10
 order = 20
 
 [quantities.ev_all]
 h5path = "FreeBoundaryStability/eigenmode_vacuum_energies"
 type = "complex_vector"
 extract = "all_complex"
-label = "Power Normalized vacuum energy (all)"
+label = "root-area-weighted vacuum energy (all)"
 noise_threshold = 1e-10
 order = 21
 
 [quantities.et_all]
 h5path = "FreeBoundaryStability/eigenmode_energies"
 type = "complex_vector"
 extract = "all_complex"
-label = "Power Normalized total energy (all)"
+label = "root-area-weighted total energy (all)"
 noise_threshold = 1e-10
 order = 22
 
@@ -299,19 +299,19 @@ label = "Chirikov parameter"
 noise_threshold = 1e-8
 order = 82
 
-[quantities.resonant_flux]
-h5path = "perturbed_equilibrium/singular_coupling/resonant_flux"
+[quantities.resonant_area_weighted_field]
+h5path = "perturbed_equilibrium/singular_coupling/resonant_area_weighted_field"
 type = "complex_vector"
 extract = "norm"
-label = "||resonant flux||"
+label = "||resonant area-weighted field||"
 noise_threshold = 1e-8
 order = 83
 
-[quantities.resonant_flux_all]
-h5path = "perturbed_equilibrium/singular_coupling/resonant_flux"
+[quantities.resonant_area_weighted_field_all]
+h5path = "perturbed_equilibrium/singular_coupling/resonant_area_weighted_field"
 type = "complex_vector"
 extract = "all_complex"
-label = "resonant flux Phi_res"
+label = "resonant area-weighted field b^r"
 noise_threshold = 1e-8
 
 # Perturbed equilibrium: energies

regression-harness/cases/solovev_multi_n.toml

@@ -7,20 +7,20 @@ name = "solovev_multi_n"
 description = "Solovev analytical equilibrium, multi-n, ideal stability"
 example_dir = "examples/Solovev_ideal_example_multi_n"
 
-# Energies — leading eigenvalues at the final truncation (psilim). eigenmode_* are now power-normalized (Φ-space, Jacobian-invariant); ξ-space counterparts live under FreeBoundaryStability/XiNorm/ and are not tracked.
+# Energies — leading eigenvalues at the final truncation (psilim). eigenmode_* are now root-area-weighted (Φ-space, Jacobian-invariant); ξ-space counterparts live under FreeBoundaryStability/XiNorm/ and are not tracked.
 [quantities.et_real]
 h5path = "FreeBoundaryStability/eigenmode_energies"
 type = "complex_vector"
 extract = "real_first"
-label = "Power Normalized total energy Re(et[1])"
+label = "root-area-weighted total energy Re(et[1])"
 noise_threshold = 1e-10
 order = 10
 
 [quantities.et_imag]
 h5path = "FreeBoundaryStability/eigenmode_energies"
 type = "complex_vector"
 extract = "imag_first"
-label = "Power Normalized total energy Im(et[1])"
+label = "root-area-weighted total energy Im(et[1])"
 noise_threshold = 1e-10
 order = 11
 
@@ -37,23 +37,23 @@ order = 14
 h5path = "FreeBoundaryStability/eigenmode_plasma_energies"
 type = "complex_vector"
 extract = "all_complex"
-label = "Power Normalized plasma energy (all)"
+label = "root-area-weighted plasma energy (all)"
 noise_threshold = 1e-10
 order = 20
 
 [quantities.ev_all]
 h5path = "FreeBoundaryStability/eigenmode_vacuum_energies"
 type = "complex_vector"
 extract = "all_complex"
-label = "Power Normalized vacuum energy (all)"
+label = "root-area-weighted vacuum energy (all)"
 noise_threshold = 1e-10
 order = 21
 
 [quantities.et_all]
 h5path = "FreeBoundaryStability/eigenmode_energies"
 type = "complex_vector"
 extract = "all_complex"
-label = "Power Normalized total energy (all)"
+label = "root-area-weighted total energy (all)"
 noise_threshold = 1e-10
 order = 22