diff --git a/CLAUDE.md b/CLAUDE.md index b4878bfc0..bb7a65689 100644 --- a/CLAUDE.md +++ b/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) diff --git a/benchmarks/benchmark_against_fortran_run.jl b/benchmarks/benchmark_against_fortran_run.jl index d5125af13..4af6e58eb 100644 --- a/benchmarks/benchmark_against_fortran_run.jl +++ b/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) diff --git a/docs/resources/2008-Park-Spectral_asymmetry_due_to_magnetic_coordinates.pdf b/docs/resources/2008-Park-Spectral_asymmetry_due_to_magnetic_coordinates.pdf new file mode 100644 index 000000000..1351fb014 Binary files /dev/null and b/docs/resources/2008-Park-Spectral_asymmetry_due_to_magnetic_coordinates.pdf differ diff --git a/docs/src/citations.md b/docs/src/citations.md index 57017a172..dfb65fb65 100644 --- a/docs/src/citations.md +++ b/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. diff --git a/docs/src/conventions.md b/docs/src/conventions.md index 8886aee2f..268fde014 100644 --- a/docs/src/conventions.md +++ b/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) diff --git a/docs/src/equilibrium.md b/docs/src/equilibrium.md index 76f4cfc00..eefd99835 100644 --- a/docs/src/equilibrium.md +++ b/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 diff --git a/docs/src/set_up.md b/docs/src/set_up.md index 4d7614ae7..38c9dcebf 100644 --- a/docs/src/set_up.md +++ b/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: diff --git a/regression-harness/cases/diiid_n1.toml b/regression-harness/cases/diiid_n1.toml index 0b822d249..1a69b7279 100644 --- a/regression-harness/cases/diiid_n1.toml +++ b/regression-harness/cases/diiid_n1.toml @@ -8,12 +8,12 @@ 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 @@ -21,7 +21,7 @@ order = 10 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 @@ -29,7 +29,7 @@ order = 11 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 @@ -37,7 +37,7 @@ order = 12 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,7 +54,7 @@ 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 @@ -62,7 +62,7 @@ order = 20 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 @@ -70,7 +70,7 @@ order = 21 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 diff --git a/regression-harness/cases/solovev_multi_n.toml b/regression-harness/cases/solovev_multi_n.toml index 553c71f40..c67882ca6 100644 --- a/regression-harness/cases/solovev_multi_n.toml +++ b/regression-harness/cases/solovev_multi_n.toml @@ -7,12 +7,12 @@ 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 @@ -20,7 +20,7 @@ order = 10 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,7 +37,7 @@ 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 @@ -45,7 +45,7 @@ order = 20 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 @@ -53,7 +53,7 @@ order = 21 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 diff --git a/regression-harness/cases/solovev_n1.toml b/regression-harness/cases/solovev_n1.toml index 4ff03feeb..eecbf8ccc 100644 --- a/regression-harness/cases/solovev_n1.toml +++ b/regression-harness/cases/solovev_n1.toml @@ -7,12 +7,12 @@ name = "solovev_n1" description = "Solovev analytical equilibrium, n=1, ideal stability" example_dir = "examples/Solovev_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 @@ -20,7 +20,7 @@ order = 10 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 @@ -28,7 +28,7 @@ order = 11 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 @@ -36,7 +36,7 @@ order = 12 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 @@ -53,7 +53,7 @@ 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 @@ -61,7 +61,7 @@ order = 20 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 @@ -69,7 +69,7 @@ order = 21 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 diff --git a/src/Analysis/PerturbedEquilibrium.jl b/src/Analysis/PerturbedEquilibrium.jl index 3c738ddb3..f1219a2e9 100644 --- a/src/Analysis/PerturbedEquilibrium.jl +++ b/src/Analysis/PerturbedEquilibrium.jl @@ -18,13 +18,14 @@ function _has_pe_data(h5path, key) end """ - plot_resonant_flux(h5path; save_path=nothing) + plot_resonant_area_weighted_field_amplitude(h5path; save_path=nothing) -Scatter plot of `|Φ_res|` per rational surface vs ψ_N. One marker series per toroidal -mode n. Integer-valued q rational surfaces are annotated. +Scatter plot of the resonant area-weighted field `|b^r|` per rational surface vs ψ_N. One marker series +per toroidal mode n. Integer-valued q rational surfaces are annotated. The resonant area-weighted field +`b^r = Φ^r/A^r` is the resonant flux normalized by the scalar surface area, in tesla [Pharr 2026]. Requires the perturbed equilibrium module to have been run and -`singular_coupling/resonant_flux` to be present in the HDF5 file. +`singular_coupling/resonant_area_weighted_field` to be present in the HDF5 file. ### Arguments @@ -38,29 +39,29 @@ Requires the perturbed equilibrium module to have been run and A `Plots.jl` plot object. """ -function plot_resonant_flux(h5path; save_path=nothing) +function plot_resonant_area_weighted_field_amplitude(h5path; save_path=nothing) base = "perturbed_equilibrium/singular_coupling/" - _has_pe_data(h5path, base * "resonant_flux") || - return plot(; title="No resonant flux data — run with perturbed equilibrium enabled", legend=false) + _has_pe_data(h5path, base * "resonant_area_weighted_field") || + return plot(; title="No resonant area-weighted field data — run with perturbed equilibrium enabled", legend=false) - resonant_flux, rational_psi, rational_q, rational_n = h5open(h5path, "r") do fid - read(fid[base * "resonant_flux"]), + resonant_area_weighted_field, rational_psi, rational_q, rational_n = h5open(h5path, "r") do fid + read(fid[base * "resonant_area_weighted_field"]), read(fid[base * "rational_psi"]), read(fid[base * "rational_q"]), read(fid[base * "rational_n"]) end - p = plot(; xlabel="Norm. Poloidal Flux", ylabel="|Φ_res|", - title="Resonant flux |Φ_res| per surface", legend=:outertopright, + p = plot(; xlabel="Norm. Poloidal Flux", ylabel="|b^r| [T]", + title="Resonant area-weighted field |b^r| per surface", legend=:outertopright, left_margin=10Plots.mm, bottom_margin=5Plots.mm) for nn in unique(rational_n) mask = rational_n .== nn - scatter!(p, rational_psi[mask], abs.(resonant_flux[mask]); + scatter!(p, rational_psi[mask], abs.(resonant_area_weighted_field[mask]); label="n=$nn", markersize=7, markerstrokewidth=0) for k in findall(mask) abs(rational_q[k] - round(rational_q[k])) < 0.05 || continue - annotate!(p, rational_psi[k], abs(resonant_flux[k]), + annotate!(p, rational_psi[k], abs(resonant_area_weighted_field[k]), text(" q=$(round(Int, rational_q[k]))", 8, :left, :black)) end end @@ -240,17 +241,17 @@ function plot_driven_delta_prime(h5path; save_path=nothing) end """ - plot_resonant_field(h5path; save_path=nothing) + plot_resonant_area_weighted_field(h5path; save_path=nothing) Five-panel summary of resonant coupling quantities at each rational surface vs ψ_N: - - `|Φ_res|`: resonant flux (`plot_resonant_flux`) + - `|b^r|`: resonant area-weighted field (`plot_resonant_area_weighted_field_amplitude`) - `Re(Δ')`: tearing stability parameter (`plot_driven_delta_prime`) - `|I_res|`: resonant current - `w/2`: island half-width (`plot_island_widths`) - `K`: Chirikov overlap parameter (`plot_chirikov_parameter`) -Inspired by `plot_resonant_field.py` from OMFIT GPEC. +Inspired by `plot_resonant_flux.py` from OMFIT GPEC. ### Arguments @@ -264,8 +265,8 @@ Inspired by `plot_resonant_field.py` from OMFIT GPEC. A `Plots.jl` plot object. """ -function plot_resonant_field(h5path; save_path=nothing) - p1 = plot_resonant_flux(h5path) +function plot_resonant_area_weighted_field(h5path; save_path=nothing) + p1 = plot_resonant_area_weighted_field_amplitude(h5path) p2 = plot_driven_delta_prime(h5path) p3 = _plot_resonant_current(h5path) p4 = plot_island_widths(h5path) @@ -342,7 +343,7 @@ Inspired by `plot_spectrograms.py` from OMFIT GPEC. - `component`: Response field component to plot; one of `:xi_psi`, `:b_psi`, `:b_theta`, `:b_zeta` (default: `:xi_psi`). `:b_psi` reads the area-normalized - `psi_area` dataset. + `b_psi_area_weighted` dataset. - `save_path`: If provided, save the figure to this path (default: `nothing`) ### Returns @@ -352,7 +353,7 @@ A `Plots.jl` plot object. function plot_mode_spectrogram(h5path; component=:xi_psi, save_path=nothing) comp_map = Dict( :xi_psi => "xi_psi", - :b_psi => "psi_area", # area-normalized b^ψ + :b_psi => "b_psi_area_weighted", # area-normalized b^ψ :b_theta => "b_theta", :b_zeta => "b_zeta", ) @@ -421,7 +422,7 @@ end Three-panel composite summary of perturbed equilibrium results: - - Top-left: Resonant flux per surface (`plot_resonant_flux`) + - Top-left: Resonant area-weighted field per surface (`plot_resonant_area_weighted_field_amplitude`) - Top-right: Edge |b_ψ| spectrum - Bottom: ξ_ψ mode spectrogram (`plot_mode_spectrogram`) @@ -438,7 +439,7 @@ Three-panel composite summary of perturbed equilibrium results: A `Plots.jl` plot object. """ function plot_perturbed_equilibrium_summary(h5path; save_path=nothing) - p_islands = plot_resonant_flux(h5path) + p_islands = plot_resonant_area_weighted_field_amplitude(h5path) p_bpsi = _plot_bpsi_edge_spectrum(h5path) p_spectro = plot_mode_spectrogram(h5path; component=:xi_psi) @@ -452,11 +453,11 @@ end # Internal helper — |b_psi(m)| spectrum at the outermost psi surface function _plot_bpsi_edge_spectrum(h5path) base = "perturbed_equilibrium/response/" - _has_pe_data(h5path, base * "psi_area") || + _has_pe_data(h5path, base * "b_psi_area_weighted") || return plot(; title="No b_psi data — run with perturbed equilibrium enabled", legend=false) data, mlow, mhigh = h5open(h5path, "r") do fid - read(fid[base * "psi_area"]), + read(fid[base * "b_psi_area_weighted"]), read(fid["info/mlow"]), read(fid["info/mhigh"]) end diff --git a/src/Equilibrium/CoordinateInvariant.jl b/src/Equilibrium/CoordinateInvariant.jl new file mode 100644 index 000000000..fb5c2492d --- /dev/null +++ b/src/Equilibrium/CoordinateInvariant.jl @@ -0,0 +1,136 @@ +""" +Coordinate-invariant (root-area-weighted field) flux-surface operators. + +These building blocks implement the √area weighting of Pharr (2026), +"Coordinate-invariant flux-surface Fourier analysis in tokamaks". They are the +single source of truth for translating a flux-surface field component between the +three field representations (all in tesla), shared by the ForceFreeStates and +PerturbedEquilibrium modules: + + - `b` bare / normal field (Fourier weight W = 1) + - `b̃` root-area-weighted field (W = √(J|∇ψ|)) — the coordinate-invariant field + - `b̄` area-weighted field (W = J|∇ψ|) + +With `Σ ≡ sqrtamat` (the mpert×mpert √weight convolution) and the scalar surface area +`A ≡ jarea = ∫ J|∇ψ| dθ`, the mode-space vectors are related by + + b̃ = Σ·b, b̄ = (Σ/√A)·b̃, Φ = A·b̄ = Σ·√A·b̃ (poloidal flux, weber) + +The √-weight (`b̃`) basis is the one in which an operator's singular values / spectra are +independent of the straight-field-line (working) coordinate; `b` and `b̄` are not +coordinate-invariant and are provided only as field/flux recovery views. Poloidal flux is +recovered, when ever needed, as the scalar product `Φ = A·b̄` — it is never stored. + +See `scripts/test_power_norm_invariance.jl` for the numerical invariance proof of +the underlying √weight identity and angle map. +""" + +""" + compute_sqrt_jac_delpsi(equil, psi, mtheta) -> Vector{Float64} + +Compute √(J·|∇ψ|) at `mtheta` equally-spaced θ points on the flux surface at `psi`. +This is the √weight function that maps a field component `b(θ)` to its +root-area-weighted form `√(J|∇ψ|)·b(θ)` in θ-space. +""" +function compute_sqrt_jac_delpsi(equil::PlasmaEquilibrium, psi::Float64, mtheta::Int) + sqrt_jac_delpsi = Vector{Float64}(undef, mtheta) + + hint2d = (Ref(1), Ref(1)) + for itheta in 0:(mtheta-1) + theta = itheta / mtheta # normalized to [0, 1) + m = flux_surface_metric(equil, psi, theta; hint=hint2d) + sqrt_jac_delpsi[itheta+1] = sqrt(abs(m.jac * m.delpsi)) + end + + return sqrt_jac_delpsi +end + +""" + compute_sqrtamat(equil, psi, ft) -> Matrix{ComplexF64} + +Build the √A convolution matrix `sqrtamat`. Produces a Hermitian Toeplitz matrix +with entries sqrtamat[m',k] = ŵ_{m_k − m'} where ŵ_n = (1/N)Σ w_j exp(+inθ_j) and +w(θ) = √(J·|∇ψ|). + +Operationally, sqrtamat is the mode-space √weight operator: for a field b with +Fourier coefficients b_fft, it satisfies the identity + `‖sqrtamat·b_fft‖² = N² · ∫ |b|² · J|∇ψ| dθ` +which is Jacobian-invariant on a given flux surface (see +`scripts/test_power_norm_invariance.jl`). + +The backward step uses exp(−imθ)/(1/N) normalization paired with the Julia +forward FT exp(+imθ) so that round-trip = identity and the convolution +structure is correct. +""" +function compute_sqrtamat( + equil::PlasmaEquilibrium, + psi::Float64, + ft::Utilities.FourierTransforms.FourierTransform +) + mpert = ft.mpert + mtheta = ft.mtheta + + sqrt_jdp = compute_sqrt_jac_delpsi(equil, psi, mtheta) + sqrtamat = zeros(ComplexF64, mpert, mpert) + + e_k = zeros(ComplexF64, mpert) + for k in 1:mpert + e_k .= 0.0 + e_k[k] = 1.0 + 0.0im + + # Standard backward FT: f(θ_j) = (1/N) Σ_m c_m exp(-imθ_j) + # exp(-imθ) = cos(mθ) - i·sin(mθ), so: + # Re(f) = (1/N)(cslth·Re(c) + snlth·Im(c)) + # Im(f) = (1/N)(cslth·Im(c) - snlth·Re(c)) + real_part = (ft.cslth * real.(e_k) .+ ft.snlth * imag.(e_k)) ./ mtheta + imag_part = (ft.cslth * imag.(e_k) .- ft.snlth * real.(e_k)) ./ mtheta + theta_vec = complex.(real_part, imag_part) + + # Multiply pointwise by √(J·|∇ψ|) in theta-space + theta_vec .*= sqrt_jdp + + # Forward FT: theta-space → mode-space (exp(+imθ), Julia convention) + sqrtamat[:, k] .= ft(theta_vec) + end + + return sqrtamat +end + +""" + rootarea_to_area_weight(equil, psi, ft) -> Matrix{ComplexF64} + +Build the root-area-weighted → area-weighted field operator `Σ/√A = sqrtamat ./ √jarea` at the +flux surface `psi`, where `jarea = ∫ J|∇ψ| dθ` is the scalar flux-surface area. It maps the +coordinate-invariant root-area-weighted field `b̃` to the area-weighted field `b̄`: `b̄ = (Σ/√A)·b̃` +(both in tesla). Poloidal flux is the scalar product `Φ = A·b̄` (so `Φ = Σ·√A·b̃`), recovered only +when a user supplies/requests flux — it is never stored. + +The √-weight (`b̃`) basis is the one in which operator singular values / spectra are independent of +the straight-field-line (working) coordinate — see `field_space_response_matrices`. `b̄` is *not* +coordinate-invariant; it is only a field/flux recovery view. [Pharr 2026] +""" +function rootarea_to_area_weight( + equil::PlasmaEquilibrium, + psi::Float64, + ft::Utilities.FourierTransforms.FourierTransform +) + sqrtamat = compute_sqrtamat(equil, psi, ft) + jarea = flux_surface_area(equil, psi, ft.mtheta) + return sqrtamat ./ sqrt(jarea) +end + +""" + area_to_rootarea_weight(equil, psi, ft) -> Matrix{ComplexF64} + +Inverse of [`rootarea_to_area_weight`](@ref): the area-weighted → root-area-weighted field operator +`√A·Σ⁻¹ = √jarea · inv(sqrtamat)`, mapping `b̄ → b̃` (`b̃ = √A·Σ⁻¹·b̄`). [Pharr 2026] +""" +function area_to_rootarea_weight( + equil::PlasmaEquilibrium, + psi::Float64, + ft::Utilities.FourierTransforms.FourierTransform +) + sqrtamat = compute_sqrtamat(equil, psi, ft) + jarea = flux_surface_area(equil, psi, ft.mtheta) + return sqrt(jarea) .* inv(sqrtamat) +end diff --git a/src/Equilibrium/Equilibrium.jl b/src/Equilibrium/Equilibrium.jl index 3cefd43cd..d11e1d649 100644 --- a/src/Equilibrium/Equilibrium.jl +++ b/src/Equilibrium/Equilibrium.jl @@ -8,10 +8,12 @@ using TOML using FastInterpolations using AdaptiveArrayPools import StaticArrays: @MMatrix, SVector +import ..Utilities # --- Internal Module Structure --- include("EquilibriumTypes.jl") include("FluxSurfaceMetrics.jl") +include("CoordinateInvariant.jl") include("ReadEquilibrium.jl") include("DirectEquilibrium.jl") include("DirectEquilibriumArcLength.jl") @@ -26,6 +28,7 @@ export setup_equilibrium, EquilibriumConfig, PlasmaEquilibrium, EquilibriumParam ProfileSplines, GeometryProfileSplines, compute_geometry_profiles, KineticProfileSplines, load_kinetic_profiles export flux_surface_metric, flux_surface_area +export compute_sqrt_jac_delpsi, compute_sqrtamat, rootarea_to_area_weight, area_to_rootarea_weight # --- Constants --- const mu0 = 4π * 1e-7 diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl index c124479d7..daa373376 100644 --- a/src/ForceFreeStates/EulerLagrange.jl +++ b/src/ForceFreeStates/EulerLagrange.jl @@ -939,7 +939,7 @@ function findmax_dW_edge!(odet::OdeState, ctrl::ForceFreeStatesControl, equil::E # Create a rough spline for wv matrix between psiedge -> psilim so we can approximate dW es.wvmat = free_compute_wv_spline(ctrl, equil, intr) - # Pre-compute sqrtamat + jarea spline for power-normalized eigenvalues + # Pre-compute sqrtamat + jarea spline for root-area-weighted eigenvalues es.sqrtamat_spline = free_compute_sqrtamat_spline(ctrl, equil, intr) # Loop with compact index j into EdgeScanState; ODE index is edge_start + j - 1. @@ -954,10 +954,10 @@ function findmax_dW_edge!(odet::OdeState, ctrl::ForceFreeStatesControl, equil::E es.plasma_energy[j] = result.plasma_energy es.vacuum_energy[j] = result.vacuum_energy es.vacuum_eigenvalue[j] = result.vacuum_eigenvalue - es.pn_total_eigenvalue[j] = result.pn_total_eigenvalue - es.pn_plasma_energy[j] = result.pn_plasma_energy - es.pn_vacuum_energy[j] = result.pn_vacuum_energy - es.pn_vacuum_eigenvalue[j] = result.pn_vacuum_eigenvalue + es.rootA_total_eigenvalue[j] = result.rootA_total_eigenvalue + es.rootA_plasma_energy[j] = result.rootA_plasma_energy + es.rootA_vacuum_energy[j] = result.rootA_vacuum_energy + es.rootA_vacuum_eigenvalue[j] = result.rootA_vacuum_eigenvalue catch e e isa LinearAlgebra.SingularException || rethrow() end diff --git a/src/ForceFreeStates/ForceFreeStates.jl b/src/ForceFreeStates/ForceFreeStates.jl index 3eb3d2c77..cc144731f 100644 --- a/src/ForceFreeStates/ForceFreeStates.jl +++ b/src/ForceFreeStates/ForceFreeStates.jl @@ -29,7 +29,7 @@ include("FixedKineticMatrices.jl") include("Kinetic.jl") include("FixedBoundaryStability.jl") include("Utils.jl") -include("PowerNorm.jl") +include("RootAreaWeighted.jl") include("Free.jl") include("Riccati.jl") diff --git a/src/ForceFreeStates/ForceFreeStatesStructs.jl b/src/ForceFreeStates/ForceFreeStatesStructs.jl index 67398f3ac..337eb158b 100644 --- a/src/ForceFreeStates/ForceFreeStatesStructs.jl +++ b/src/ForceFreeStates/ForceFreeStatesStructs.jl @@ -381,8 +381,8 @@ Populated in `Free.jl`. - `wt0::Array{ComplexF64, 2}` - Free-boundary total-energy matrix W = wp + wv before diagonalisation (numpert_total × numpert_total). **ξ-space.** - `wp::Array{ComplexF64, 2}` - Plasma energy matrix (numpert_total × numpert_total). **ξ-space.** - `wv::Array{ComplexF64, 2}` - Vacuum energy matrix (numpert_total × numpert_total). **ξ-space.** - - `pn_wt0, pn_wp, pn_wv::Array{ComplexF64,2}` - Power-normalized flux (Φ-space) counterparts of `wt0`, `wp`, `wv` at the plasma edge (see PowerNorm.jl). Jacobian-invariant up to the M†·W·M transform. - - `pn_wt::Array{ComplexF64,2}` - Φ-space eigenvector matrix of `pn_wt0` (columns sorted most-unstable first, phase-normalized so each column's largest-magnitude entry is real-positive). + - `rootA_wt0, rootA_wp, rootA_wv::Array{ComplexF64,2}` - Root-area-weighted flux (Φ-space) counterparts of `wt0`, `wp`, `wv` at the plasma edge (see RootAreaWeighted.jl). Jacobian-invariant up to the M†·W·M transform. + - `rootA_wt::Array{ComplexF64,2}` - Φ-space eigenvector matrix of `rootA_wt0` (columns sorted most-unstable first, phase-normalized so each column's largest-magnitude entry is real-positive). - `ep::Vector{ComplexF64}` - Plasma eigenvalues - `ev::Vector{ComplexF64}` - Vacuum eigenvalues - `et::Vector{ComplexF64}` - Total eigenvalues of plasma + vacuum @@ -408,16 +408,16 @@ Populated in `Free.jl`. n_tor_idx::Vector{Int} = zeros(Int, numpert_total) vacuum_eigenvalue::Float64 = NaN - # Power-normalized flux eigenvalues (Jacobian-invariant) - pn_et::Vector{ComplexF64} = fill(complex(NaN), numpert_total) - pn_ep::Vector{ComplexF64} = fill(complex(NaN), numpert_total) - pn_ev::Vector{ComplexF64} = fill(complex(NaN), numpert_total) + # Root-area-weighted flux eigenvalues (Jacobian-invariant) + rootA_et::Vector{ComplexF64} = fill(complex(NaN), numpert_total) + rootA_ep::Vector{ComplexF64} = fill(complex(NaN), numpert_total) + rootA_ev::Vector{ComplexF64} = fill(complex(NaN), numpert_total) - # Power-normalized flux matrices at the plasma edge (Jacobian-invariant). - pn_wt0::Array{ComplexF64,2} = fill(complex(NaN), numpert_total, numpert_total) - pn_wp::Array{ComplexF64,2} = fill(complex(NaN), numpert_total, numpert_total) - pn_wv::Array{ComplexF64,2} = fill(complex(NaN), numpert_total, numpert_total) - pn_wt::Array{ComplexF64,2} = fill(complex(NaN), numpert_total, numpert_total) + # Root-area-weighted flux matrices at the plasma edge (Jacobian-invariant). + rootA_wt0::Array{ComplexF64,2} = fill(complex(NaN), numpert_total, numpert_total) + rootA_wp::Array{ComplexF64,2} = fill(complex(NaN), numpert_total, numpert_total) + rootA_wv::Array{ComplexF64,2} = fill(complex(NaN), numpert_total, numpert_total) + rootA_wt::Array{ComplexF64,2} = fill(complex(NaN), numpert_total, numpert_total) grri::Array{Float64,2} = Array{Float64}(undef, 2 * numpoints, 2 * numpert_total) grre::Array{Float64,2} = Array{Float64}(undef, 2 * numpoints, 2 * numpert_total) @@ -432,17 +432,17 @@ EdgeScanState Holds the state and results for the edge dW stability scan over ψ ∈ [psiedge, psilim]. Initialized and populated by `findmax_dW_edge!`; results written to HDF5 under `EdgeScan/` -(power-norm values) and `EdgeScan/XiNorm/` (ξ-space values, for benchmarking). +(root-area-weighted values) and `EdgeScan/XiNorm/` (ξ-space values, for benchmarking). ## Fields - `wvmat` - Precomputed wv matrix spline (raw, no singfac); singfac applied analytically in `free_compute_total`. - `wv_hint::Base.RefValue{Int}` - Search hint for wvmat spline (different grid from equilibrium profiles). - - `sqrtamat_spline` - Precomputed √A convolution matrix + surface area spline for power-norm eigenvalues. + - `sqrtamat_spline` - Precomputed √A convolution matrix + surface area spline for root-area-weighted eigenvalues. - `sqrtamat_hint::Base.RefValue{Int}` - Search hint for sqrtamat spline. - `psi, q` - ψ and q values at each edge scan step. - `total_eigenvalue, plasma_energy, vacuum_energy, vacuum_eigenvalue` - ξ-space energy components at each step (NaN for failed steps). Kept for Fortran benchmarking. - - `pn_total_eigenvalue, pn_plasma_energy, pn_vacuum_energy, pn_vacuum_eigenvalue` - Power-normalized flux (Φ-space) energies (Jacobian-invariant; NaN for failed/singular steps). These are the values written to `EdgeScan/` by default. + - `rootA_total_eigenvalue, rootA_plasma_energy, rootA_vacuum_energy, rootA_vacuum_eigenvalue` - Root-area-weighted flux (Φ-space) energies (Jacobian-invariant; NaN for failed/singular steps). These are the values written to `EdgeScan/` by default. """ @kwdef mutable struct EdgeScanState numpert_total::Int @@ -452,7 +452,7 @@ Initialized and populated by `findmax_dW_edge!`; results written to HDF5 under ` wvmat::CubicSeriesInterpolant{Float64,ComplexF64} = _empty_series_interp_complex(numpert_total^2) wv_hint::Base.RefValue{Int} = Ref(1) - # Power-norm sqrtamat + jarea spline (mpert^2 + 1 series: flattened sqrtamat then jarea) + # Root-area-weighted sqrtamat + jarea spline (mpert^2 + 1 series: flattened sqrtamat then jarea) sqrtamat_spline::CubicSeriesInterpolant{Float64,ComplexF64} = _empty_series_interp_complex(numpert_total^2 + 1) sqrtamat_hint::Base.RefValue{Int} = Ref(1) @@ -464,11 +464,11 @@ Initialized and populated by `findmax_dW_edge!`; results written to HDF5 under ` vacuum_energy::Vector{ComplexF64} = fill(complex(NaN), N_edge) vacuum_eigenvalue::Vector{Float64} = fill(NaN, N_edge) - # Power-normalized flux eigenvalues (Jacobian-invariant) - pn_total_eigenvalue::Vector{ComplexF64} = fill(complex(NaN), N_edge) - pn_plasma_energy::Vector{ComplexF64} = fill(complex(NaN), N_edge) - pn_vacuum_energy::Vector{ComplexF64} = fill(complex(NaN), N_edge) - pn_vacuum_eigenvalue::Vector{Float64} = fill(NaN, N_edge) + # Root-area-weighted flux eigenvalues (Jacobian-invariant) + rootA_total_eigenvalue::Vector{ComplexF64} = fill(complex(NaN), N_edge) + rootA_plasma_energy::Vector{ComplexF64} = fill(complex(NaN), N_edge) + rootA_vacuum_energy::Vector{ComplexF64} = fill(complex(NaN), N_edge) + rootA_vacuum_eigenvalue::Vector{Float64} = fill(NaN, N_edge) end EdgeScanState(numpert_total::Int, N_edge::Int) = EdgeScanState(; numpert_total, N_edge) diff --git a/src/ForceFreeStates/Free.jl b/src/ForceFreeStates/Free.jl index 2faaf04de..31cec0f50 100644 --- a/src/ForceFreeStates/Free.jl +++ b/src/ForceFreeStates/Free.jl @@ -94,20 +94,20 @@ and data dumping. vac_data.ev[ipert] = wvt[ipert, ipert] end - # Eigenspectrum of W_Φ at psilim — Jacobian-invariant energy values; see PowerNorm.jl. + # Eigenspectrum of W_Φ at psilim — Jacobian-invariant energy values; see RootAreaWeighted.jl. # Computed directly here (no spline); spline is used by the edge scan. mtheta_eq = length(equil.rzphi_ys) - ft_pn = Utilities.FourierTransforms.FourierTransform(mtheta_eq, mpert, mlow) - sqrtamat_pn = compute_sqrtamat(equil, psilim, ft_pn) - jarea_pn = Equilibrium.flux_surface_area(equil, psilim, mtheta_eq) - pn_result = compute_power_norm_eigenvalues(vac_data.wt0, wp, vac_data.wv, sqrtamat_pn, jarea_pn, equil, psilim, intr; all_eigenvalues=true) - vac_data.pn_et .= pn_result.pn_et_all - vac_data.pn_ep .= pn_result.pn_ep_all - vac_data.pn_ev .= pn_result.pn_ev_all - vac_data.pn_wt0 .= pn_result.wt_pn - vac_data.pn_wp .= pn_result.wp_pn - vac_data.pn_wv .= pn_result.wv_pn - vac_data.pn_wt .= pn_result.pn_eigenvectors + ft_rootA = Utilities.FourierTransforms.FourierTransform(mtheta_eq, mpert, mlow) + sqrtamat_rootA = Equilibrium.compute_sqrtamat(equil, psilim, ft_rootA) + jarea_rootA = Equilibrium.flux_surface_area(equil, psilim, mtheta_eq) + rootA_result = compute_rootarea_eigenvalues(vac_data.wt0, wp, vac_data.wv, sqrtamat_rootA, jarea_rootA, equil, psilim, intr; all_eigenvalues=true) + vac_data.rootA_et .= rootA_result.rootA_et_all + vac_data.rootA_ep .= rootA_result.rootA_ep_all + vac_data.rootA_ev .= rootA_result.rootA_ev_all + vac_data.rootA_wt0 .= rootA_result.wt_rootA + vac_data.rootA_wp .= rootA_result.wp_rootA + vac_data.rootA_wv .= rootA_result.wv_rootA + vac_data.rootA_wt .= rootA_result.rootA_eigenvectors # Normalize eigenvectors based on scaled wt coeffs = odet.u[:, :, 1, end] \ (vac_data.wt .* (2π * equil.psio * 1e-3)) @@ -220,7 +220,7 @@ wv matrix spline to `free_compute_wv_spline` and pass it in `odet.edge_scan.wvma # Compute total energy matrix and eigen-decomposition wt .= wp .+ wv - # Save wt before eigen (which overwrites the input) for power-norm computation + # Save wt before eigen (which overwrites the input) for root-area-weighted computation wt_saved = copy(wt) Ev = eigen(wt) @@ -261,7 +261,7 @@ wv matrix spline to `free_compute_wv_spline` and pass it in `odet.edge_scan.wvma # negative. Clamp to zero to enforce the physical constraint. vacuum_eigenvalue = real(max(0.0, minimum(real.(eigvals(Hermitian(wv))))) / norm) - # Eigenspectrum of W_Φ — Jacobian-invariant energy values; see PowerNorm.jl. + # Eigenspectrum of W_Φ — Jacobian-invariant energy values; see RootAreaWeighted.jl. # Evaluate sqrtamat + jarea from the pre-computed spline. mpert = intr.mpert sqrtamat_flat = Vector{ComplexF64}(undef, mpert^2 + 1) @@ -269,9 +269,9 @@ wv matrix spline to `free_compute_wv_spline` and pass it in `odet.edge_scan.wvma sqrtamat_local = reshape(@view(sqrtamat_flat[1:(mpert^2)]), mpert, mpert) jarea_local = real(sqrtamat_flat[end]) - pn_result = compute_power_norm_eigenvalues(wt_saved, wp, wv, sqrtamat_local, jarea_local, equil, odet.psifac, intr) + rootA_result = compute_rootarea_eigenvalues(wt_saved, wp, wv, sqrtamat_local, jarea_local, equil, odet.psifac, intr) return (total_eigenvalue=eigenvalues[1], plasma_energy=plasma_energy, vacuum_energy=vacuum_energy, vacuum_eigenvalue=vacuum_eigenvalue, - pn_total_eigenvalue=pn_result.pn_total_eigenvalue, pn_plasma_energy=pn_result.pn_plasma_energy, pn_vacuum_energy=pn_result.pn_vacuum_energy, - pn_vacuum_eigenvalue=pn_result.pn_vacuum_eigenvalue) + rootA_total_eigenvalue=rootA_result.rootA_total_eigenvalue, rootA_plasma_energy=rootA_result.rootA_plasma_energy, rootA_vacuum_energy=rootA_result.rootA_vacuum_energy, + rootA_vacuum_eigenvalue=rootA_result.rootA_vacuum_eigenvalue) end diff --git a/src/ForceFreeStates/PowerNorm.jl b/src/ForceFreeStates/PowerNorm.jl deleted file mode 100644 index 0b1d503b1..000000000 --- a/src/ForceFreeStates/PowerNorm.jl +++ /dev/null @@ -1,289 +0,0 @@ -""" -Power-normalized flux eigenvalue computation for the FFS edge scan. - -Transforms the energy matrix W from ξ-space into power-normalized flux Φ-space and -returns the eigenvalues of `W_Φ = M†·W·M`. The **eigenvalues** (energies) are -coordinate-invariant: energy is a physically meaningful scalar and does not depend -on the choice of straight-field-line coordinate. - - dW = ξ†·W·ξ = Φ†·M†·W·M·Φ - -Transformation chain: - - T = diag(i·2π·χ₁·singfac) converts ξ → flux (χ₁ = 2π·ψ₀, singfac = m − n·q). - - ptof = sqrtamat·√jarea converts the power-normalized field → flux. - - M = T⁻¹·ptof maps Φ → ξ. - -What is **NOT** invariant across Jacobian choices (see -`scripts/test_power_norm_invariance.jl` for numerical proof): - - ξ itself (W is Jacobian-dependent). - - Φ itself — T depends on the m-labelling of the chosen Jacobian, so - ‖Φ‖ = ‖M⁻¹ξ‖ drifts with Jacobian. - - The θ-space field reconstructed from Φ. - -What **IS** invariant (verified numerically in the test script, within the -area-integral precision floor ~1e-6): - - Flux-surface area A = ∫ J|∇ψ| dθ. - - The √weight operator identity ‖sqrtamat·b_fft‖² = N²·∫|b|²·J|∇ψ| dθ. - - The angle-map convmat (b_jac2 = convmat·b_jac1) to machine precision. - - The eigenspectrum of `W_Φ`. -""" - -""" - compute_sqrt_jac_delpsi(equil, psi, mtheta) -> Vector{Float64} - -Compute √(J·|∇ψ|) at `mtheta` equally-spaced θ points on the flux surface at `psi`. -This is the √weight function that maps a field component `b(θ)` to its -power-normalized form `√(J|∇ψ|)·b(θ)` in θ-space. -""" -function compute_sqrt_jac_delpsi(equil::Equilibrium.PlasmaEquilibrium, psi::Float64, mtheta::Int) - sqrt_jac_delpsi = Vector{Float64}(undef, mtheta) - - hint2d = (Ref(1), Ref(1)) - for itheta in 0:(mtheta-1) - theta = itheta / mtheta # normalized to [0, 1) - m = Equilibrium.flux_surface_metric(equil, psi, theta; hint=hint2d) - sqrt_jac_delpsi[itheta+1] = sqrt(abs(m.jac * m.delpsi)) - end - - return sqrt_jac_delpsi -end - -""" - compute_sqrtamat(equil, psi, ft) -> Matrix{ComplexF64} - -Build the √A convolution matrix `sqrtamat`. Produces a Hermitian Toeplitz matrix -with entries sqrtamat[m',k] = ŵ_{m_k − m'} where ŵ_n = (1/N)Σ w_j exp(+inθ_j) and -w(θ) = √(J·|∇ψ|). - -Operationally, sqrtamat is the mode-space √weight operator: for a field b with -Fourier coefficients b_fft, it satisfies the identity - `‖sqrtamat·b_fft‖² = N² · ∫ |b|² · J|∇ψ| dθ` -which is Jacobian-invariant on a given flux surface (see -`scripts/test_power_norm_invariance.jl`). - -The backward step uses exp(−imθ)/(1/N) normalization paired with the Julia -forward FT exp(+imθ) so that round-trip = identity and the convolution -structure is correct. -""" -function compute_sqrtamat( - equil::Equilibrium.PlasmaEquilibrium, - psi::Float64, - ft::Utilities.FourierTransforms.FourierTransform -) - mpert = ft.mpert - mtheta = ft.mtheta - - sqrt_jdp = compute_sqrt_jac_delpsi(equil, psi, mtheta) - sqrtamat = zeros(ComplexF64, mpert, mpert) - - e_k = zeros(ComplexF64, mpert) - for k in 1:mpert - e_k .= 0.0 - e_k[k] = 1.0 + 0.0im - - # Standard backward FT: f(θ_j) = (1/N) Σ_m c_m exp(-imθ_j) - # exp(-imθ) = cos(mθ) - i·sin(mθ), so: - # Re(f) = (1/N)(cslth·Re(c) + snlth·Im(c)) - # Im(f) = (1/N)(cslth·Im(c) - snlth·Re(c)) - real_part = (ft.cslth * real.(e_k) .+ ft.snlth * imag.(e_k)) ./ mtheta - imag_part = (ft.cslth * imag.(e_k) .- ft.snlth * real.(e_k)) ./ mtheta - theta_vec = complex.(real_part, imag_part) - - # Multiply pointwise by √(J·|∇ψ|) in theta-space - theta_vec .*= sqrt_jdp - - # Forward FT: theta-space → mode-space (exp(+imθ), Julia convention) - sqrtamat[:, k] .= ft(theta_vec) - end - - return sqrtamat -end - -""" - compute_power_norm_eigenvalues(wt, wp, wv, sqrtamat, jarea, equil, psi, intr; all_eigenvalues=false) - -Transform W from ξ-space to power-normalized flux Φ-space and eigendecompose. - -The transformation chain: - 1. T = diag(i·2π·χ₁·singfac) converts ξ → flux (singfac = m − n·q). - 2. ptof = sqrtamat·√jarea converts the power-normalized field → flux. - 3. M = T⁻¹·ptof maps Φ → ξ (power-norm to displacement). - 4. W_Φ = M†·W·M = ptof†·T⁻†·W·T⁻¹·ptof. - -Only the eigenspectrum of W_Φ is physically invariant across Jacobian choices — -it is the coordinate-independent energy. The eigenvectors in this basis are -normalized under Julia's `eigen` (‖v‖₂ = 1) rather than the ∫|b|²dA = 1 -convention, so `v` should not be reused as a physical Φ-space mode shape — here -it is only used to split λ = v†Wt v into ⟨v,Wp v⟩ + ⟨v,Wv v⟩, both of which are -spectrally invariant because Wp + Wv = Wt. - -Returns the standard eigenvalues of W_Φ sorted ascending (most negative/unstable -first), matching the ξ-space convention in `free_run!` and `free_compute_total`. -Also returns `pn_vacuum_eigenvalue = max(0, λ_min(Hermitian(wv_Φ)))` — the Φ-space -counterpart of the ξ-space `vacuum_eigenvalue` diagnostic. -Returns NaN when any singfac ≈ 0 (rational surface crossing makes T singular). -""" -function compute_power_norm_eigenvalues( - wt::AbstractMatrix{ComplexF64}, - wp::AbstractMatrix{ComplexF64}, - wv::AbstractMatrix{ComplexF64}, - sqrtamat::AbstractMatrix{ComplexF64}, - jarea::Float64, - equil::Equilibrium.PlasmaEquilibrium, - psi::Float64, - intr::ForceFreeStatesInternal; - all_eigenvalues::Bool=false -) - Npert = intr.numpert_total - mpert = intr.mpert - npert = intr.npert - chi1 = 2π * equil.psio - - # Compute singfac = m - n·q for each mode - q_at_psi = equil.profiles.q_spline(psi) - singfac = vec((intr.mlow:intr.mhigh) .- q_at_psi .* (intr.nlow:intr.nhigh)') - - # Check for rational surface proximity — T is singular there - if any(abs.(singfac) .< 1e-6) - if all_eigenvalues - nan_vec = fill(complex(NaN), Npert) - nan_mat = fill(complex(NaN), Npert, Npert) - return (pn_total_eigenvalue=complex(NaN), pn_plasma_energy=complex(NaN), pn_vacuum_energy=complex(NaN), - pn_vacuum_eigenvalue=NaN, pn_et_all=nan_vec, pn_ep_all=nan_vec, pn_ev_all=nan_vec, - wt_pn=nan_mat, wp_pn=nan_mat, wv_pn=nan_mat, pn_eigenvectors=nan_mat) - else - return (pn_total_eigenvalue=complex(NaN), pn_plasma_energy=complex(NaN), pn_vacuum_energy=complex(NaN), - pn_vacuum_eigenvalue=NaN) - end - end - - # T = diag(i·2π·chi1·singfac): ξ → flux - T_diag_inv = Vector{ComplexF64}(undef, Npert) - for i in 1:Npert - T_diag_inv[i] = 1.0 / (im * 2π * chi1 * singfac[i]) - end - - # ptof = sqrtamat * sqrt(jarea): power-norm field → flux - # For multi-n, expand mpert×mpert sqrtamat to Npert×Npert block-diagonal - ptof_block = sqrtamat .* sqrt(jarea) - if npert == 1 - ptof_full = ptof_block - else - ptof_full = zeros(ComplexF64, Npert, Npert) - for in in 1:npert - r = ((in - 1) * mpert + 1):(in * mpert) - ptof_full[r, r] .= ptof_block - end - end - - # M = diag(T_inv) · ptof: row-scale ptof by T_inv (maps Φ → ξ) - M = similar(ptof_full) - for i in 1:Npert - M[i, :] .= T_diag_inv[i] .* ptof_full[i, :] - end - - # Transform energy matrices: W_Φ = M† · W · M - wt_pn = M' * wt * M - wp_pn = M' * wp * M - wv_pn = M' * wv * M - - # Smallest eigenvalue of the vacuum matrix alone in Φ-space, clamped to zero. - # wv_pn should be PSD by congruence of the PSD ξ-space wv; numerical noise - # can make the smallest eigenvalue slightly negative. - pn_vacuum_eigenvalue = real(max(0.0, minimum(real.(eigvals(Hermitian(wv_pn)))))) - - # Eigendecompose the total energy in power-norm space. Only the eigenspectrum - # is Jacobian-invariant; the eigenvectors are coordinate-dependent. - Ev = eigen(wt_pn) - eigenvalues = Ev.values - vectors = Ev.vectors - - # Sort ascending by real part (most negative/most unstable first) — matches ξ-space convention in free_run!/free_compute_total - eindex = sortperm(real.(eigenvalues); rev=true) - - # λ = v†Wt v = v†Wp v + v†Wv v ; the plasma/vacuum split is the projection - # onto the same eigenvector and is therefore also Jacobian-invariant. - if all_eigenvalues - pn_et_all = Vector{ComplexF64}(undef, Npert) - pn_ep_all = Vector{ComplexF64}(undef, Npert) - pn_ev_all = Vector{ComplexF64}(undef, Npert) - pn_eigenvectors = Matrix{ComplexF64}(undef, Npert, Npert) - - for i in 1:Npert - v = vectors[:, eindex[Npert+1-i]] - pn_eigenvectors[:, i] .= v - pn_et_all[i] = eigenvalues[eindex[Npert+1-i]] - pn_ep_all[i] = dot(v, wp_pn * v) - pn_ev_all[i] = dot(v, wv_pn * v) - end - - # Phase convention: rotate each column so its largest-magnitude entry is real-positive - # (matches the ξ-space convention in free_run!). Magnitudes are preserved since the - # eigenvectors come from Julia's `eigen` with ‖v‖₂ = 1, the natural Φ-space norm. - for isol in 1:Npert - imax = argmax(abs.(@view pn_eigenvectors[:, isol])) - phase = abs(pn_eigenvectors[imax, isol]) / pn_eigenvectors[imax, isol] - @view(pn_eigenvectors[:, isol]) .*= phase - end - - return (pn_total_eigenvalue=pn_et_all[1], pn_plasma_energy=pn_ep_all[1], pn_vacuum_energy=pn_ev_all[1], - pn_vacuum_eigenvalue=pn_vacuum_eigenvalue, - pn_et_all=pn_et_all, pn_ep_all=pn_ep_all, pn_ev_all=pn_ev_all, - wt_pn=wt_pn, wp_pn=wp_pn, wv_pn=wv_pn, pn_eigenvectors=pn_eigenvectors) - else - idx = eindex[Npert] - v = vectors[:, idx] - pn_total = eigenvalues[idx] - pn_plasma = ComplexF64(dot(v, wp_pn * v)) - pn_vacuum = ComplexF64(dot(v, wv_pn * v)) - - return (pn_total_eigenvalue=pn_total, pn_plasma_energy=pn_plasma, pn_vacuum_energy=pn_vacuum, - pn_vacuum_eigenvalue=pn_vacuum_eigenvalue) - end -end - -""" - free_compute_sqrtamat_spline(ctrl, equil, intr) -> CubicSeriesInterpolant - -Pre-compute sqrtamat and jarea over the edge scan ψ range and return a cubic series -interpolant. Uses the same q-evenly-spaced grid as `free_compute_wv_spline`. - -The interpolant stores `mpert^2 + 1` complex series per grid point: - - First mpert^2 values: flattened sqrtamat matrix - - Last value: jarea (stored as complex with zero imaginary part) -""" -function free_compute_sqrtamat_spline(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal) - profiles = equil.profiles - mtheta_eq = length(equil.rzphi_ys) - mpert = intr.mpert - - # Same grid as wv spline: evenly spaced in q - qedge = profiles.q_spline(ctrl.psiedge) - npsi = max(4, ceil(Int, (intr.qlim - qedge) * intr.nhigh * 4)) - psi_array = zeros(Float64, npsi + 1) - - n_series = mpert^2 + 1 # sqrtamat (flattened) + jarea - data_array = zeros(ComplexF64, npsi + 1, n_series) - - # Create FourierTransform once (same for all psi since mtheta/mpert/mlow don't change) - ft = Utilities.FourierTransforms.FourierTransform(mtheta_eq, mpert, intr.mlow) - - for i in 1:(npsi+1) - qi = qedge + (intr.qlim - qedge) * ((i - 1) / npsi) - psii = ctrl.psiedge + (intr.psilim - ctrl.psiedge) * ((i - 1) / npsi) - psi_array[i] = find_zero( - (psi -> profiles.q_spline(psi) - qi, - psi -> profiles.q_deriv(psi)), - psii, Roots.Newton() - ) - - sqrtamat = compute_sqrtamat(equil, psi_array[i], ft) - jarea = Equilibrium.flux_surface_area(equil, psi_array[i], mtheta_eq) - - # Flatten sqrtamat into first mpert^2 entries, jarea as last entry - data_array[i, 1:mpert^2] .= vec(sqrtamat) - data_array[i, end] = complex(jarea, 0.0) - end - - return cubic_interp(psi_array, Series(data_array); extrap=ExtendExtrap()) -end diff --git a/src/ForceFreeStates/RootAreaWeighted.jl b/src/ForceFreeStates/RootAreaWeighted.jl new file mode 100644 index 000000000..bbd451539 --- /dev/null +++ b/src/ForceFreeStates/RootAreaWeighted.jl @@ -0,0 +1,221 @@ +""" +Root-area-weighted field eigenvalue computation for the FFS edge scan. + +Transforms the energy matrix W from ξ-space into the coordinate-invariant root-area-weighted +field (b̃) space and returns the eigenvalues of the scaled quadratic form `c·W_t` (energies). +The **eigenvalues** (energies) are coordinate-invariant: energy is a physically meaningful +scalar and does not depend on the choice of straight-field-line coordinate. + + dW = ξ†·W·ξ = c·b̃†·(M†·W·M)·b̃ = c·b̃†·W_t·b̃ + +Transformation chain (no poloidal flux Φ appears — it would only be the scalar product `Φ = A·b̄`): + - T = diag(i·2π·χ₁·singfac) converts ξ → flux-like amplitude (χ₁ = 2π·ψ₀, singfac = m − n·q). + - Σ = sqrtamat is the √weight (the bare → root-area-weighted field map). + - M = T⁻¹·Σ maps the root-area-weighted field b̃ → ξ. + - c = jarea = ∫ J|∇ψ| dθ is the scalar surface area carrying the energy's area normalization. + +What is **NOT** invariant across Jacobian choices (see +`scripts/test_power_norm_invariance.jl` for numerical proof): + - ξ itself (W is Jacobian-dependent). + - The θ-space field reconstructed from the eigenvector. + +What **IS** invariant (verified numerically in the test script, within the +area-integral precision floor ~1e-6): + - Flux-surface area A = ∫ J|∇ψ| dθ. + - The √weight operator identity ‖sqrtamat·b_fft‖² = N²·∫|b|²·J|∇ψ| dθ. + - The angle-map convmat (b_jac2 = convmat·b_jac1) to machine precision. + - The eigenspectrum of `c·W_t`. +""" + +# The √weight building blocks `compute_sqrt_jac_delpsi`, `compute_sqrtamat`, and the field-translation +# operators now live in `Equilibrium/CoordinateInvariant.jl` so the ForceFreeStates and +# PerturbedEquilibrium modules share one coordinate-invariant definition. + +""" + compute_rootarea_eigenvalues(wt, wp, wv, sqrtamat, jarea, equil, psi, intr; all_eigenvalues=false) + +Transform W from ξ-space to the root-area-weighted field (b̃) space and eigendecompose. + +The transformation chain: + 1. T = diag(i·2π·χ₁·singfac) converts ξ → flux-like amplitude (singfac = m − n·q). + 2. Σ = sqrtamat is the √weight (bare → root-area-weighted field). + 3. M = T⁻¹·Σ maps the root-area-weighted field b̃ → ξ. + 4. Energy matrix = c·W_t = c·M†·W·M with the scalar c = jarea (the surface area). + +Only the eigenspectrum of `c·W_t` is physically invariant across Jacobian choices — +it is the coordinate-independent energy. The eigenvectors in this basis are +normalized under Julia's `eigen` (‖v‖₂ = 1) rather than the ∫|b|²dA = 1 +convention, so `v` should not be reused as a physical b̃-space mode shape — here +it is only used to split λ = v†Wt v into ⟨v,Wp v⟩ + ⟨v,Wv v⟩, both of which are +spectrally invariant because Wp + Wv = Wt. + +Returns the standard eigenvalues of `c·W_t` sorted ascending (most negative/unstable +first), matching the ξ-space convention in `free_run!` and `free_compute_total`. +Also returns `rootA_vacuum_eigenvalue = max(0, λ_min(Hermitian(c·wv_t)))` — the b̃-space +counterpart of the ξ-space `vacuum_eigenvalue` diagnostic. +Returns NaN when any singfac ≈ 0 (rational surface crossing makes T singular). +""" +function compute_rootarea_eigenvalues( + wt::AbstractMatrix{ComplexF64}, + wp::AbstractMatrix{ComplexF64}, + wv::AbstractMatrix{ComplexF64}, + sqrtamat::AbstractMatrix{ComplexF64}, + jarea::Float64, + equil::Equilibrium.PlasmaEquilibrium, + psi::Float64, + intr::ForceFreeStatesInternal; + all_eigenvalues::Bool=false +) + Npert = intr.numpert_total + mpert = intr.mpert + npert = intr.npert + chi1 = 2π * equil.psio + + # Compute singfac = m - n·q for each mode + q_at_psi = equil.profiles.q_spline(psi) + singfac = vec((intr.mlow:intr.mhigh) .- q_at_psi .* (intr.nlow:intr.nhigh)') + + # Check for rational surface proximity — T is singular there + if any(abs.(singfac) .< 1e-6) + if all_eigenvalues + nan_vec = fill(complex(NaN), Npert) + nan_mat = fill(complex(NaN), Npert, Npert) + return (rootA_total_eigenvalue=complex(NaN), rootA_plasma_energy=complex(NaN), rootA_vacuum_energy=complex(NaN), + rootA_vacuum_eigenvalue=NaN, rootA_et_all=nan_vec, rootA_ep_all=nan_vec, rootA_ev_all=nan_vec, + wt_rootA=nan_mat, wp_rootA=nan_mat, wv_rootA=nan_mat, rootA_eigenvectors=nan_mat) + else + return (rootA_total_eigenvalue=complex(NaN), rootA_plasma_energy=complex(NaN), rootA_vacuum_energy=complex(NaN), + rootA_vacuum_eigenvalue=NaN) + end + end + + # T = diag(i·2π·chi1·singfac): ξ → flux-like amplitude + T_diag_inv = Vector{ComplexF64}(undef, Npert) + for i in 1:Npert + T_diag_inv[i] = 1.0 / (im * 2π * chi1 * singfac[i]) + end + + # √weight operator Σ = sqrtamat maps the bare field → root-area-weighted field b̃. + # For multi-n, expand the mpert×mpert block to Npert×Npert block-diagonal. + sqrtamat_full = sqrtamat + if npert != 1 + sqrtamat_full = zeros(ComplexF64, Npert, Npert) + for in in 1:npert + r = ((in - 1) * mpert + 1):(in * mpert) + sqrtamat_full[r, r] .= sqrtamat + end + end + + # M = diag(T_inv) · Σ: row-scale Σ by T_inv (maps the root-area-weighted field b̃ → ξ) + M = similar(sqrtamat_full) + for i in 1:Npert + M[i, :] .= T_diag_inv[i] .* sqrtamat_full[i, :] + end + + # Energy is the b̃ quadratic form scaled by the scalar surface area c = jarea: + # dW = c·b̃†·W_t·b̃, W_t = M†·W·M, c = jarea [Pharr 2026, Eq. for coordinate-invariant energy] + # Folding the scalar c here keeps the eigenspectrum (the coordinate-invariant energy) physical [J]. + wt_rootA = jarea .* (M' * wt * M) + wp_rootA = jarea .* (M' * wp * M) + wv_rootA = jarea .* (M' * wv * M) + + # Smallest eigenvalue of the vacuum matrix alone in b̃-space, clamped to zero. + # wv_rootA should be PSD by congruence of the PSD ξ-space wv; numerical noise + # can make the smallest eigenvalue slightly negative. + rootA_vacuum_eigenvalue = real(max(0.0, minimum(real.(eigvals(Hermitian(wv_rootA)))))) + + # Eigendecompose the total energy in root-area-weighted space. Only the eigenspectrum + # is Jacobian-invariant; the eigenvectors are coordinate-dependent. + Ev = eigen(wt_rootA) + eigenvalues = Ev.values + vectors = Ev.vectors + + # Sort ascending by real part (most negative/most unstable first) — matches ξ-space convention in free_run!/free_compute_total + eindex = sortperm(real.(eigenvalues); rev=true) + + # λ = v†Wt v = v†Wp v + v†Wv v ; the plasma/vacuum split is the projection + # onto the same eigenvector and is therefore also Jacobian-invariant. + if all_eigenvalues + rootA_et_all = Vector{ComplexF64}(undef, Npert) + rootA_ep_all = Vector{ComplexF64}(undef, Npert) + rootA_ev_all = Vector{ComplexF64}(undef, Npert) + rootA_eigenvectors = Matrix{ComplexF64}(undef, Npert, Npert) + + for i in 1:Npert + v = vectors[:, eindex[Npert+1-i]] + rootA_eigenvectors[:, i] .= v + rootA_et_all[i] = eigenvalues[eindex[Npert+1-i]] + rootA_ep_all[i] = dot(v, wp_rootA * v) + rootA_ev_all[i] = dot(v, wv_rootA * v) + end + + # Phase convention: rotate each column so its largest-magnitude entry is real-positive + # (matches the ξ-space convention in free_run!). Magnitudes are preserved since the + # eigenvectors come from Julia's `eigen` with ‖v‖₂ = 1, the natural b̃-space norm. + for isol in 1:Npert + imax = argmax(abs.(@view rootA_eigenvectors[:, isol])) + phase = abs(rootA_eigenvectors[imax, isol]) / rootA_eigenvectors[imax, isol] + @view(rootA_eigenvectors[:, isol]) .*= phase + end + + return (rootA_total_eigenvalue=rootA_et_all[1], rootA_plasma_energy=rootA_ep_all[1], rootA_vacuum_energy=rootA_ev_all[1], + rootA_vacuum_eigenvalue=rootA_vacuum_eigenvalue, + rootA_et_all=rootA_et_all, rootA_ep_all=rootA_ep_all, rootA_ev_all=rootA_ev_all, + wt_rootA=wt_rootA, wp_rootA=wp_rootA, wv_rootA=wv_rootA, rootA_eigenvectors=rootA_eigenvectors) + else + idx = eindex[Npert] + v = vectors[:, idx] + rootA_total = eigenvalues[idx] + rootA_plasma = ComplexF64(dot(v, wp_rootA * v)) + rootA_vacuum = ComplexF64(dot(v, wv_rootA * v)) + + return (rootA_total_eigenvalue=rootA_total, rootA_plasma_energy=rootA_plasma, rootA_vacuum_energy=rootA_vacuum, + rootA_vacuum_eigenvalue=rootA_vacuum_eigenvalue) + end +end + +""" + free_compute_sqrtamat_spline(ctrl, equil, intr) -> CubicSeriesInterpolant + +Pre-compute sqrtamat and jarea over the edge scan ψ range and return a cubic series +interpolant. Uses the same q-evenly-spaced grid as `free_compute_wv_spline`. + +The interpolant stores `mpert^2 + 1` complex series per grid point: + - First mpert^2 values: flattened sqrtamat matrix + - Last value: jarea (stored as complex with zero imaginary part) +""" +function free_compute_sqrtamat_spline(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal) + profiles = equil.profiles + mtheta_eq = length(equil.rzphi_ys) + mpert = intr.mpert + + # Same grid as wv spline: evenly spaced in q + qedge = profiles.q_spline(ctrl.psiedge) + npsi = max(4, ceil(Int, (intr.qlim - qedge) * intr.nhigh * 4)) + psi_array = zeros(Float64, npsi + 1) + + n_series = mpert^2 + 1 # sqrtamat (flattened) + jarea + data_array = zeros(ComplexF64, npsi + 1, n_series) + + # Create FourierTransform once (same for all psi since mtheta/mpert/mlow don't change) + ft = Utilities.FourierTransforms.FourierTransform(mtheta_eq, mpert, intr.mlow) + + for i in 1:(npsi+1) + qi = qedge + (intr.qlim - qedge) * ((i - 1) / npsi) + psii = ctrl.psiedge + (intr.psilim - ctrl.psiedge) * ((i - 1) / npsi) + psi_array[i] = find_zero( + (psi -> profiles.q_spline(psi) - qi, + psi -> profiles.q_deriv(psi)), + psii, Roots.Newton() + ) + + sqrtamat = Equilibrium.compute_sqrtamat(equil, psi_array[i], ft) + jarea = Equilibrium.flux_surface_area(equil, psi_array[i], mtheta_eq) + + # Flatten sqrtamat into first mpert^2 entries, jarea as last entry + data_array[i, 1:mpert^2] .= vec(sqrtamat) + data_array[i, end] = complex(jarea, 0.0) + end + + return cubic_interp(psi_array, Series(data_array); extrap=ExtendExtrap()) +end diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 6829acc1c..72b26283d 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -587,17 +587,17 @@ function write_outputs_to_HDF5( out_h5["integration/crit"] = odet.crit_store # Write edge stability scan data (only present when psiedge < psilim). - # Power-normalized flux (Φ-space) energies are the default — they are Jacobian- + # Root-area-weighted (Φ-space) energies are the default — they are Jacobian- # invariant. The ξ-space values sit under EdgeScan/XiNorm/ and are retained for # benchmarking against the Fortran GPEC lineage. if !isempty(odet.edge_scan.psi) es = odet.edge_scan out_h5["EdgeScan/psi"] = es.psi out_h5["EdgeScan/q"] = es.q - out_h5["EdgeScan/total_energy"] = es.pn_total_eigenvalue - out_h5["EdgeScan/plasma_energy"] = es.pn_plasma_energy - out_h5["EdgeScan/vacuum_energy"] = es.pn_vacuum_energy - out_h5["EdgeScan/vacuum_eigenvalue"] = es.pn_vacuum_eigenvalue + out_h5["EdgeScan/total_energy"] = es.rootA_total_eigenvalue + out_h5["EdgeScan/plasma_energy"] = es.rootA_plasma_energy + out_h5["EdgeScan/vacuum_energy"] = es.rootA_vacuum_energy + out_h5["EdgeScan/vacuum_eigenvalue"] = es.rootA_vacuum_eigenvalue out_h5["EdgeScan/XiNorm/total_energy"] = es.total_eigenvalue out_h5["EdgeScan/XiNorm/plasma_energy"] = es.plasma_energy out_h5["EdgeScan/XiNorm/vacuum_energy"] = es.vacuum_energy @@ -645,19 +645,19 @@ function write_outputs_to_HDF5( out_h5["singular/kinetic/scan_cond"] = intr.kinsing_scan_cond out_h5["singular/kinetic/scan_threshold"] = intr.kinsing_scan_threshold - # Write free-boundary stability data. Power-normalized flux (Φ-space) is the + # Write free-boundary stability data. Root-area-weighted (Φ-space) is the # default — Jacobian-invariant. ξ-space counterparts sit under # FreeBoundaryStability/XiNorm/ for Fortran benchmarking. # W_freeboundary_eigenmodes holds the eigenvector matrix of W_freeboundary with # columns sorted most-unstable first; the same phase normalization is applied in # both spaces (largest-magnitude entry made real-positive). - out_h5["FreeBoundaryStability/W_freeboundary"] = ctrl.vac_flag ? vac_data.pn_wt0 : ComplexF64[] - out_h5["FreeBoundaryStability/W_plasma"] = ctrl.vac_flag ? vac_data.pn_wp : ComplexF64[] - out_h5["FreeBoundaryStability/W_vacuum"] = ctrl.vac_flag ? vac_data.pn_wv : ComplexF64[] - out_h5["FreeBoundaryStability/W_freeboundary_eigenmodes"] = ctrl.vac_flag ? vac_data.pn_wt : ComplexF64[] - out_h5["FreeBoundaryStability/eigenmode_energies"] = ctrl.vac_flag ? vac_data.pn_et : ComplexF64[] - out_h5["FreeBoundaryStability/eigenmode_plasma_energies"] = ctrl.vac_flag ? vac_data.pn_ep : ComplexF64[] - out_h5["FreeBoundaryStability/eigenmode_vacuum_energies"] = ctrl.vac_flag ? vac_data.pn_ev : ComplexF64[] + out_h5["FreeBoundaryStability/W_freeboundary"] = ctrl.vac_flag ? vac_data.rootA_wt0 : ComplexF64[] + out_h5["FreeBoundaryStability/W_plasma"] = ctrl.vac_flag ? vac_data.rootA_wp : ComplexF64[] + out_h5["FreeBoundaryStability/W_vacuum"] = ctrl.vac_flag ? vac_data.rootA_wv : ComplexF64[] + out_h5["FreeBoundaryStability/W_freeboundary_eigenmodes"] = ctrl.vac_flag ? vac_data.rootA_wt : ComplexF64[] + out_h5["FreeBoundaryStability/eigenmode_energies"] = ctrl.vac_flag ? vac_data.rootA_et : ComplexF64[] + out_h5["FreeBoundaryStability/eigenmode_plasma_energies"] = ctrl.vac_flag ? vac_data.rootA_ep : ComplexF64[] + out_h5["FreeBoundaryStability/eigenmode_vacuum_energies"] = ctrl.vac_flag ? vac_data.rootA_ev : ComplexF64[] out_h5["FreeBoundaryStability/XiNorm/W_freeboundary"] = ctrl.vac_flag ? vac_data.wt0 : ComplexF64[] out_h5["FreeBoundaryStability/XiNorm/W_plasma"] = ctrl.vac_flag ? vac_data.wp : ComplexF64[] out_h5["FreeBoundaryStability/XiNorm/W_vacuum"] = ctrl.vac_flag ? vac_data.wv : ComplexF64[] diff --git a/src/PerturbedEquilibrium/FieldReconstruction.jl b/src/PerturbedEquilibrium/FieldReconstruction.jl index 541537b30..b599098fe 100644 --- a/src/PerturbedEquilibrium/FieldReconstruction.jl +++ b/src/PerturbedEquilibrium/FieldReconstruction.jl @@ -60,7 +60,7 @@ Tuple of (xi_modes, b_modes) NamedTuples: - `xi_modes.zeta_reg`: ξ^ζ regularized (= xmz, from gpeq_contra with reg_spot smoothing) - `xi_modes.cova_psi/theta/zeta`: covariant displacement (from gpeq_cova) - `b_modes.psi`: b^ψ [npsi, mpert] - - `b_modes.psi_area`: b^ψ / ⟨J·|∇ψ|⟩_θ (area-normalized, for b_n computation) + - `b_modes.b_psi_area_weighted`: b^ψ / ⟨J·|∇ψ|⟩_θ (area-normalized, for b_n computation) - `b_modes.theta`: b^θ [npsi, mpert] - `b_modes.zeta`: b^ζ [npsi, mpert] - `b_modes.theta_reg/zeta_reg`: regularized b^θ, b^ζ (from gpeq_sol with reg_spot smoothing) @@ -181,7 +181,7 @@ function reconstruct_physical_fields( ) b_modes = ( psi=b_psi_modes, # b^ψ (no Jacobian) — used for b_n normal projection - psi_area=Jb_psi_modes, # b^ψ / ⟨J·|∇ψ|⟩_θ (area-normalized, for b_n) + b_psi_area_weighted=Jb_psi_modes, # b^ψ / ⟨J·|∇ψ|⟩_θ (area-normalized, for b_n) theta=b_theta_modes, # b^θ unregularized zeta=b_zeta_modes, # b^ζ unregularized theta_reg=b_theta_reg, # b^θ regularized (from gpeq_sol with reg_spot) diff --git a/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl b/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl index b065294c8..9474abf0f 100644 --- a/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl +++ b/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl @@ -78,22 +78,22 @@ Results from perturbed equilibrium calculations. Response fields (mode space): - `xi_modes::Union{Nothing, NamedTuple}` - Displacement (psi, theta, zeta) [npsi, mpert] - - `b_modes::Union{Nothing, NamedTuple}` - Magnetic field; psi=b^ψ, psi_area=b^ψ/⟨J·|∇ψ|⟩_θ, theta/zeta=unregularized, theta_reg/zeta_reg=regularized [npsi, mpert] + - `b_modes::Union{Nothing, NamedTuple}` - Magnetic field; psi=b^ψ, b_psi_area_weighted=b^ψ/⟨J·|∇ψ|⟩_θ, theta/zeta=unregularized, theta_reg/zeta_reg=regularized [npsi, mpert] - `b_n_modes::Union{Nothing, Matrix{ComplexF64}}` - Physical normal field b_n [npsi, mpert] - `xi_n_modes::Union{Nothing, Matrix{ComplexF64}}` - Physical normal displacement xi_n [npsi, mpert] Coupling matrices [n_rational × numpert_total] — 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 [mode_C, m_out]). - - `C_resonant_flux` - Φ_r/A coupling (singcoup row 1) + - `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_field` - Penetrated field coupling (singcoup row 4) + - `C_penetrated_area_weighted_field` - Penetrated area-weighted field coupling (singcoup row 4) - `C_delta_prime` - Δ' coupling (singcoup row 5) Applied resonant vectors [n_rational] = C · forcing_amplitudes. Matches Fortran `Phi_res`, `w_isl`, `K_isl`, `Delta`. - - `resonant_flux`, `resonant_current`, `island_width_sq`, `penetrated_field`, `delta_prime` + - `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(|island_width_sq|) from applied forcing @@ -102,13 +102,27 @@ Diagnostics [n_rational]: Metadata [n_rational] — identifies each (surface, n) row: - `rational_psi`, `rational_q`, `rational_m_res`, `rational_n`, `rational_surface_idx` -Control surface matrices [numpert_total × numpert_total]: - - `plasma_inductance` - Lambda (wt0-based plasma inductance) - - `surface_inductance` - L (vacuum surface inductance from Green's functions) - - `permeability` - P = Lambda * L^{-1} (plasma response matrix, Phi_tot = P * Phi_x) - - `reluctance` - Rho = L^{-1} * (Lambda - L) * L^{-1} - -Energies (Fortran gpout convention; Φ_x external flux, Φ_tot total flux, L/Λ inductances): +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 [numpert_total × numpert_total], 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, Λ⁻¹·Φ_tot⟩ ) / 4 (energy to perturb the plasma; Fortran's "total energy") @@ -125,17 +139,17 @@ Energies (Fortran gpout convention; Φ_x external flux, Φ_tot total flux, L/Λ xi_n_modes::Union{Nothing, Matrix{ComplexF64}} = nothing # physical normal displacement xi_n [npsi, mpert] # Coupling matrices [n_rational × numpert_total] - C_resonant_flux::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) + C_resonant_area_weighted_field::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) C_resonant_current::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) C_island_width_sq::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) - C_penetrated_field::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) + C_penetrated_area_weighted_field::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) C_delta_prime::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) # Applied resonant vectors [n_rational] = C · amp_vec - resonant_flux::Vector{ComplexF64} = ComplexF64[] + resonant_area_weighted_field::Vector{ComplexF64} = ComplexF64[] resonant_current::Vector{ComplexF64} = ComplexF64[] island_width_sq::Vector{ComplexF64} = ComplexF64[] - penetrated_field::Vector{ComplexF64} = ComplexF64[] + penetrated_area_weighted_field::Vector{ComplexF64} = ComplexF64[] delta_prime::Vector{ComplexF64} = ComplexF64[] # Diagnostics [n_rational] @@ -149,15 +163,21 @@ Energies (Fortran gpout convention; Φ_x external flux, Φ_tot total flux, L/Λ rational_n::Vector{Int} = Int[] rational_surface_idx::Vector{Int} = Int[] - # Control surface perturbation vectors [numpert_total] - forcing_vec::Vector{ComplexF64} = ComplexF64[] # Phi_x: external forcing in eigenmode basis - response_vec::Vector{ComplexF64} = ComplexF64[] # Phi_tot = P * Phi_x: total plasma response - - # Control surface matrices [numpert_total × numpert_total] - plasma_inductance::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) # Lambda - surface_inductance::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) # L - permeability::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) # P = Lambda * L^{-1} - reluctance::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) # Rho = L^{-1}*(Lambda-L)*L^{-1} + # Control-surface forcing/response spectra in the three Pharr 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) + response_b_rootarea::Vector{ComplexF64} = ComplexF64[] # root-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) + 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 + 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) # Energies — see the struct docstring for formulas vacuum_energy::Float64 = 0.0 diff --git a/src/PerturbedEquilibrium/Response.jl b/src/PerturbedEquilibrium/Response.jl index 5cc4ec327..df0682256 100644 --- a/src/PerturbedEquilibrium/Response.jl +++ b/src/PerturbedEquilibrium/Response.jl @@ -54,35 +54,59 @@ function compute_plasma_response!( L_inv = inv(surface_inductance) reluctance = L_inv * (plasma_inductance - surface_inductance) * L_inv - # Store control surface matrices in state for HDF5 output - state.plasma_inductance = plasma_inductance - state.surface_inductance = surface_inductance - state.permeability = permeability - state.reluctance = reluctance - - # Store permeability in internal state for singular coupling use + # Store permeability in internal state for singular coupling / field reconstruction. + # These consumers operate on the physical control-surface flux Φ_x, so the internal + # copy stays in flux space; only the stored/output quantities are conformed to fields below. intr.plasma_response = permeability - forcing_vector = map_forcing_to_eigenmodes(intr.forcing_modes, ffs_intr) - response_vector = compute_plasma_response_vector(permeability, forcing_vector) + # Conform the control-surface matrices to the coordinate-invariant root-area-weighted + # field (b̃) space for output (issue #233 / Pharr 2026). Store the b̃→b̄ operator S = Σ/√A + # and the scalar surface area A so users can recover the area-weighted field (b̄ = S·b̃) or + # flux (Φ = A·b̄) — see Utils.jl output docs. + rootarea_to_area_weight, surface_area = build_control_surface_rootarea_to_area_weight(equil, ffs_intr) + field_mats = field_space_response_matrices(plasma_inductance, surface_inductance, permeability, reluctance, rootarea_to_area_weight, surface_area) + state.plasma_inductance = field_mats.plasma_inductance + state.surface_inductance = field_mats.surface_inductance + state.permeability = field_mats.permeability + state.reluctance = field_mats.reluctance + state.rootarea_to_area_weight = rootarea_to_area_weight + state.surface_area = surface_area + + # Forcing and response on the control surface. Flux Φ appears only as a brief internal bridge: + # forcing arrives as Φ_x, the field reconstruction below consumes Φ_tot, and the b̃ spectra are + # formed via the conform operator R = S·A (Φ = R·b̃). + forcing_flux = map_forcing_to_eigenmodes(intr.forcing_modes, ffs_intr) + response_flux = compute_plasma_response_vector(permeability, forcing_flux) - state.forcing_vec = forcing_vector - state.response_vec = response_vector + # Output forcing/response in the three Pharr field representations (all tesla): + # b̃ (root-area-weighted) = R⁻¹·Φ, b (bare) = Σ⁻¹·b̃, b̄ (area-weighted) = S·b̃. + flux_conform = rootarea_to_area_weight .* surface_area # R = Σ·√A (b̃ → flux) + bare_to_rootarea = rootarea_to_area_weight .* sqrt(surface_area) # Σ (bare → b̃) + forcing_b_rootarea = flux_conform \ forcing_flux + response_b_rootarea = flux_conform \ response_flux + state.forcing_b_rootarea = forcing_b_rootarea + state.response_b_rootarea = response_b_rootarea + state.forcing_b = bare_to_rootarea \ forcing_b_rootarea + state.response_b = bare_to_rootarea \ response_b_rootarea + state.forcing_b_area = rootarea_to_area_weight * forcing_b_rootarea + state.response_b_area = rootarea_to_area_weight * response_b_rootarea - # Scalar energies and torque, ported from Fortran gpout. forcing_vector and - # response_vector are Phi_x and Phi_tot in the gpeq_weight=1 basis. + # Scalar energies and torque (Joules), ported from Fortran gpout. These are congruence-invariant + # (energy = Φ†·G⁻¹·Φ = b̃†·G̃⁻¹·b̃), so they are evaluated from the brief internal flux vectors with + # the well-conditioned flux-space inductances — the b̃ form routes through inv(R⁻¹LR⁻†) and is + # needlessly ill-conditioned. The result is a physical scalar, not a stored flux quantity. L_surf_inv = inv(surface_inductance) L_plas_inv = inv(plasma_inductance) - vy = dot(forcing_vector, L_surf_inv * forcing_vector) / 4 - sy = dot(response_vector, L_surf_inv * response_vector) / 4 - py = dot(response_vector, L_plas_inv * response_vector) / 4 + vy = dot(forcing_flux, L_surf_inv * forcing_flux) / 4 + sy = dot(response_flux, L_surf_inv * response_flux) / 4 + py = dot(response_flux, L_plas_inv * response_flux) / 4 state.vacuum_energy = real(vy) state.surface_energy = real(sy) state.plasma_energy = real(py) # Fortran's "total energy" is this pengy state.toroidal_torque = -2 * nn * imag(py) xi_modes, b_modes = reconstruct_physical_fields( - response_vector, flux_matrix, ForceFreeStates_results, equil, ffs_intr, intr, + response_flux, flux_matrix, ForceFreeStates_results, equil, ffs_intr, intr, metric, ffit, ctrl ) @@ -98,6 +122,6 @@ function compute_plasma_response!( state.xi_n_modes = xi_n_modes if ctrl.verbose - @info "Response complete: $(length(intr.forcing_modes)) forcing modes, max amplitude = $(@sprintf("%.3e", maximum(abs.(response_vector))))" + @info "Response complete: $(length(intr.forcing_modes)) forcing modes, max amplitude = $(@sprintf("%.3e", maximum(abs.(response_flux))))" end end diff --git a/src/PerturbedEquilibrium/ResponseMatrices.jl b/src/PerturbedEquilibrium/ResponseMatrices.jl index d6173d0a0..d929e9d33 100644 --- a/src/PerturbedEquilibrium/ResponseMatrices.jl +++ b/src/PerturbedEquilibrium/ResponseMatrices.jl @@ -298,6 +298,87 @@ function calc_permeability( return plasma_inductance / surface_inductance end +""" + build_control_surface_rootarea_to_area_weight( + equil::Equilibrium.PlasmaEquilibrium, + ffs_intr::ForceFreeStatesInternal + )::Tuple{Matrix{ComplexF64}, Float64} + +Build the numpert_total × numpert_total root-area-weighted → area-weighted field operator +`S = Σ/√A` at the control surface (psilim), and return it together with the scalar surface area +`A = jarea`. The mpert × mpert single-n block (Equilibrium.rootarea_to_area_weight) is repeated +block-diagonally over the `npert` toroidal harmonics, matching the numpert_total mode ordering used +by the response matrices. + +`S` maps a root-area-weighted control-surface field `b̃` to the area-weighted field `b̄` +(`b̄ = S·b̃`); poloidal flux is the scalar product `Φ = A·b̄` (so the b̃→flux conform operator is +`R = S·A`, used only internally). `S` and `A` are the recovery aids users need to express the stored +coordinate-invariant (b̃) matrices in the area-weighted field or recover flux — see +`field_space_response_matrices`. [Pharr 2026] +""" +function build_control_surface_rootarea_to_area_weight( + equil::Equilibrium.PlasmaEquilibrium, + ffs_intr::ForceFreeStatesInternal +)::Tuple{Matrix{ComplexF64}, Float64} + mpert = ffs_intr.mpert + npert = ffs_intr.npert + Npert = ffs_intr.numpert_total + + mtheta_eq = length(equil.rzphi_ys) + ft = Utilities.FourierTransforms.FourierTransform(mtheta_eq, mpert, ffs_intr.mlow) + S_block = Equilibrium.rootarea_to_area_weight(equil, ffs_intr.psilim, ft) + jarea = Equilibrium.flux_surface_area(equil, ffs_intr.psilim, mtheta_eq) + + npert == 1 && return (Matrix{ComplexF64}(S_block), jarea) + + S_full = zeros(ComplexF64, Npert, Npert) + for in in 1:npert + r = ((in - 1) * mpert + 1):(in * mpert) + S_full[r, r] .= S_block + end + return (S_full, jarea) +end + +""" + field_space_response_matrices( + plasma_inductance, surface_inductance, permeability, reluctance, + rootarea_to_area_weight, jarea + )::NamedTuple + +Express the control-surface response matrices in the coordinate-invariant root-area-weighted +field (b̃) space, given the flux-space matrices, the b̃→b̄ operator `S = rootarea_to_area_weight`, +and the scalar surface area `jarea`. The brief internal flux-conform operator is `R = S·jarea` +(`Φ = R·b̃`); poloidal flux never leaves this function. + +The matrices fall into two algebraic classes: + - **Operators** (map flux → flux): permeability `P` (Φ_tot = P·Φ_x) transforms by similarity + `P̃ = R⁻¹·P·R`. Its singular values are coordinate-invariant. + - **Quadratic generators** (energy = Φ†·G⁻¹·Φ): inductances `Λ`, `L` transform by congruence + `G̃ = R⁻¹·G·R⁻†`; the inverse-inductance-like reluctance `ϱ` (energy = Φ†·ϱ·Φ) transforms as + `ϱ̃ = R†·ϱ·R`. Their spectra are coordinate-invariant. + +These rules are mutually consistent: `P̃ = Λ̃·L̃⁻¹ = R⁻¹·Λ·L⁻¹·R = R⁻¹·P·R`, and +`ϱ̃ = L̃⁻¹·(Λ̃−L̃)·L̃⁻¹`. To recover the area-weighted (`b̄`) forms, conform with `S` instead of `R` +(e.g. `L_b̄ = S·L̃·S†`); recover flux with the scalar `A`: `Φ = A·b̄`. [Pharr 2026] +""" +function field_space_response_matrices( + plasma_inductance::Matrix{ComplexF64}, + surface_inductance::Matrix{ComplexF64}, + permeability::Matrix{ComplexF64}, + reluctance::Matrix{ComplexF64}, + rootarea_to_area_weight::Matrix{ComplexF64}, + jarea::Float64 +)::NamedTuple + R = rootarea_to_area_weight .* jarea # b̃→flux conform operator Σ·√A = (Σ/√A)·A + R_inv = inv(R) + return ( + plasma_inductance=R_inv * plasma_inductance * R_inv', + surface_inductance=R_inv * surface_inductance * R_inv', + permeability=R_inv * permeability * R, + reluctance=R' * reluctance * R + ) +end + """ map_forcing_to_eigenmodes( forcing_modes::Vector{ForcingMode}, diff --git a/src/PerturbedEquilibrium/SingularCoupling.jl b/src/PerturbedEquilibrium/SingularCoupling.jl index f3eca3f61..fd34d5f25 100644 --- a/src/PerturbedEquilibrium/SingularCoupling.jl +++ b/src/PerturbedEquilibrium/SingularCoupling.jl @@ -62,10 +62,10 @@ Matching Fortran GPEC output: `C_f_x_out` (coupling matrix) and `Phi_res`, `w_is Populates `state` with: Coupling matrices `[n_rational × numpert_total]` — C[row, j] = coupling when forcing mode j has unit amplitude: - - `C_resonant_flux`, `C_resonant_current`, `C_island_width_sq`, `C_penetrated_field`, `C_delta_prime` + - `C_resonant_area_weighted_field`, `C_resonant_current`, `C_island_width_sq`, `C_penetrated_area_weighted_field`, `C_delta_prime` Applied resonant vectors `[n_rational]` = C · amp_vec: - - `resonant_flux`, `resonant_current`, `island_width_sq`, `penetrated_field`, `delta_prime` + - `resonant_area_weighted_field`, `resonant_current`, `island_width_sq`, `penetrated_area_weighted_field`, `delta_prime` Diagnostics `[n_rational]`: `island_half_width`, `chirikov_parameter` @@ -125,10 +125,10 @@ function compute_singular_coupling_metrics!( ctrl.verbose && @info "Found $n_rational resonant (surface, n) pairs" # Phase 2: Allocate output arrays - state.C_resonant_flux = zeros(ComplexF64, n_rational, numpert_total) + state.C_resonant_area_weighted_field = zeros(ComplexF64, n_rational, numpert_total) state.C_resonant_current = zeros(ComplexF64, n_rational, numpert_total) state.C_island_width_sq = zeros(ComplexF64, n_rational, numpert_total) - state.C_penetrated_field = zeros(ComplexF64, n_rational, numpert_total) + state.C_penetrated_area_weighted_field = zeros(ComplexF64, n_rational, numpert_total) state.C_delta_prime = zeros(ComplexF64, n_rational, numpert_total) state.rational_psi = zeros(Float64, n_rational) state.rational_q = zeros(Float64, n_rational) @@ -229,19 +229,28 @@ function compute_singular_coupling_metrics!( bwp1_l = 2π * im * chi1 * (singfac_l * xsp1_l - nn * q1_l * xsp_l) bwp1_r = 2π * im * chi1 * (singfac_r * xsp1_r - nn * q1_r * xsp_r) jump_vec[k] = bwp1_r - bwp1_l - # C_penetrated_field: midpoint of b^ψ at lpsi/rpsi divided by area. - # Matches Fortran gpout_resp: gpeq_interp_singsurf evaluates bwp_mn at respsi + # C_penetrated_area_weighted_field: midpoint of b^ψ at lpsi/rpsi divided by the scalar surface area. + # Matches Fortran gpout_resp: gpeq_interp_singsurf evaluates bwp_mn at respsi. + # LHS normalization audit (#233): the resonant flux Φ^r divided by the scalar area A^r + # is a genuine field amplitude in tesla and is coordinate-invariant [Pharr 2026; cf. + # the resonant-field definition in the Conventions Reference]. b_l = chi1 * singfac_l * 2π * im * xsp_l b_r = chi1 * singfac_r * 2π * im * xsp_r - state.C_penetrated_field[row, k] = (b_l + b_r) / 2 / area + state.C_penetrated_area_weighted_field[row, k] = (b_l + b_r) / 2 / area end + # LHS normalization audit (#233) — output scalar coordinate-invariance per row: + # - Δ' (1/length): the resonant-surface jump in ∂b^ψ/∂ψ over 2π·χ₁; the tearing index is + # coordinate-invariant (its sign/zero-crossing set the stability boundary) [Glasser 2016]. + # - resonant (shielding) current: j_c already integrates jac·|∇ψ| over the surface, so the + # Jacobian weighting is carried inside j_c — no separate area factor needed. + # - resonant flux → field: Φ^r/A^r [T], invariant [Park 2008; Pharr 2026]. state.C_delta_prime[row, :] = jump_vec ./ (twopi * chi1) state.C_resonant_current[row, :] = jump_vec .* (-j_c / (twopi * m_res)) - # Matches Fortran gpout_resp: singflx = L·fkaxmn, singbnoflxs = singflx/area, + # Matches Fortran gpout_resp: singflx = L·fkaxmn, resonant area-weighted field = singflx/area, # islandhwids = 4·singflx/(2π·shear·q·chi1) singflx_pre = (L_mm / (twopi * nn)) .* state.C_resonant_current[row, :] - state.C_resonant_flux[row, :] = singflx_pre ./ area + state.C_resonant_area_weighted_field[row, :] = singflx_pre ./ area if abs(shear) > 1e-10 state.C_island_width_sq[row, :] = (4.0 / (twopi * shear * sing_surf.q * chi1)) .* singflx_pre end @@ -258,18 +267,34 @@ function compute_singular_coupling_metrics!( end end - # Phase 4: Apply forcing amplitudes → R = C · amp_vec - amp_vec = zeros(ComplexF64, numpert_total) + # Phase 4: Apply forcing amplitudes → R = C · Φ_x. The applied resonant scalars are + # physical, coordinate-invariant quantities, so evaluate them from the flux-space C and + # the physical forcing Φ_x (bit-identical to pre-conform values). + forcing_flux = zeros(ComplexF64, numpert_total) for mode in intr.forcing_modes j = findfirst(k -> intr.m_modes[k] == mode.m && intr.n_modes[k] == mode.n, 1:numpert_total) - isnothing(j) || (amp_vec[j] = mode.amplitude) + isnothing(j) || (forcing_flux[j] = mode.amplitude) end - state.resonant_flux = state.C_resonant_flux * amp_vec - state.resonant_current = state.C_resonant_current * amp_vec - state.island_width_sq = state.C_island_width_sq * amp_vec - state.penetrated_field = state.C_penetrated_field * amp_vec - state.delta_prime = state.C_delta_prime * amp_vec + state.resonant_area_weighted_field = state.C_resonant_area_weighted_field * forcing_flux + state.resonant_current = state.C_resonant_current * forcing_flux + state.island_width_sq = state.C_island_width_sq * forcing_flux + state.penetrated_area_weighted_field = state.C_penetrated_area_weighted_field * forcing_flux + state.delta_prime = state.C_delta_prime * forcing_flux + + # Conform the stored coupling-matrix input basis to the coordinate-invariant root-area-weighted + # field (b̃) space (#233 / Pharr 2026): C̃ = C·R, so each stored row acts on the applied field + # b̃_x (Φ_x = R·b̃_x) and its singular values become coordinate-invariant. The conform operator + # R = S·A (Σ·√A) is the only place flux briefly appears; it is built from the b̃→b̄ operator S and + # the scalar surface area A. Done after the applied-vector evaluation above so those physical + # scalars carry no round-trip noise. + rootarea_to_area_weight, surface_area = build_control_surface_rootarea_to_area_weight(equil, ffs_intr) + flux_conform = rootarea_to_area_weight .* surface_area + state.C_resonant_area_weighted_field = state.C_resonant_area_weighted_field * flux_conform + state.C_resonant_current = state.C_resonant_current * flux_conform + state.C_island_width_sq = state.C_island_width_sq * flux_conform + state.C_penetrated_area_weighted_field = state.C_penetrated_area_weighted_field * flux_conform + state.C_delta_prime = state.C_delta_prime * flux_conform # Phase 5: Island diagnostics from applied resonant vectors compute_island_diagnostics!(state, n_rational) diff --git a/src/PerturbedEquilibrium/Utils.jl b/src/PerturbedEquilibrium/Utils.jl index 4c4f6a00c..9957fe7cd 100644 --- a/src/PerturbedEquilibrium/Utils.jl +++ b/src/PerturbedEquilibrium/Utils.jl @@ -58,24 +58,33 @@ perturbed_equilibrium/ │ ├── 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/ │ ├── xi_psi # Radial displacement ξ^ψ = ξ·∇ψ (ComplexF64 [npsi, mpert]) │ ├── xi_psi_J # J·ξ^ψ Jacobian-weighted (from gpeq_contra) -│ ├── psi_area # b^ψ / ⟨J·|∇ψ|⟩_θ area-normalized (ComplexF64 [npsi, mpert]) +│ ├── b_psi_area_weighted # b^ψ / ⟨J·|∇ψ|⟩_θ area-normalized (ComplexF64 [npsi, mpert]) │ ├── b_n # Physical normal field b_n (ComplexF64 [npsi, mpert]) │ ├── xi_n # Physical normal displacement xi_n (ComplexF64 [npsi, mpert]) │ ├── b_theta │ └── b_zeta +├── response_matrices/ # [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̄ ├── singular_coupling/ -│ ├── C_resonant_flux # [n_rational × numpert_total] coupling matrix +│ ├── 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_field +│ ├── C_penetrated_area_weighted_field │ ├── C_delta_prime -│ ├── resonant_flux # [n_rational] applied vector = C · amp_vec +│ ├── resonant_area_weighted_field # [n_rational] applied vector = C̃ · b̃_x (resonant area-weighted field b^r [T]) │ ├── resonant_current │ ├── island_width_sq -│ ├── penetrated_field +│ ├── penetrated_area_weighted_field │ ├── delta_prime │ ├── island_half_width # [n_rational] Float64 │ ├── chirikov_parameter @@ -104,23 +113,33 @@ function write_outputs_to_HDF5( forcing_group["m"] = [mode.m for mode in intr.forcing_modes] forcing_group["amplitude"] = [mode.amplitude for mode in intr.forcing_modes] - # Control surface perturbation vectors (Phi_x and Phi_tot = P*Phi_x) - !isempty(state.forcing_vec) && (pe_group["forcing_vec"] = state.forcing_vec) - !isempty(state.response_vec) && (pe_group["response_vec"] = state.response_vec) - - # Control surface matrices [numpert_total × numpert_total] + # 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.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) + + # 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, "response_matrices") ? pe_group["response_matrices"] : create_group(pe_group, "response_matrices") !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.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) # Response fields (ComplexF64 directly) response_group = haskey(pe_group, "response") ? pe_group["response"] : create_group(pe_group, "response") 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["psi_area"] = have_b ? state.b_modes.psi_area : ComplexF64[] + response_group["b_psi_area_weighted"] = have_b ? state.b_modes.b_psi_area_weighted : ComplexF64[] response_group["b_theta"] = have_b ? state.b_modes.theta : ComplexF64[] response_group["b_zeta"] = have_b ? state.b_modes.zeta : ComplexF64[] response_group["b_n"] = !isnothing(state.b_n_modes) ? state.b_n_modes : ComplexF64[] @@ -174,17 +193,17 @@ function write_outputs_to_HDF5( coupling_group = haskey(pe_group, "singular_coupling") ? pe_group["singular_coupling"] : create_group(pe_group, "singular_coupling") # Coupling matrices [n_rational × numpert_total] - !isempty(state.C_resonant_flux) && (coupling_group["C_resonant_flux"] = state.C_resonant_flux) + !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_penetrated_field) && (coupling_group["C_penetrated_field"] = state.C_penetrated_field) + !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) # Applied resonant vectors [n_rational] - !isempty(state.resonant_flux) && (coupling_group["resonant_flux"] = state.resonant_flux) + !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_field) && (coupling_group["penetrated_field"] = state.penetrated_field) + !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.island_half_width) && (coupling_group["island_half_width"] = state.island_half_width) !isempty(state.chirikov_parameter) && (coupling_group["chirikov_parameter"] = state.chirikov_parameter) diff --git a/test/runtests.jl b/test/runtests.jl index f424a7d23..8cecea160 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -24,6 +24,7 @@ else include("./runtests_fouriertransforms.jl") include("./runtests_vacuum.jl") include("./runtests_equil.jl") + include("./runtests_coordinate_invariant.jl") include("./runtests_eulerlagrange.jl") include("./runtests_riccati.jl") include("./runtests_parallel_integration.jl") diff --git a/test/runtests_coordinate_invariant.jl b/test/runtests_coordinate_invariant.jl new file mode 100644 index 000000000..39b20330a --- /dev/null +++ b/test/runtests_coordinate_invariant.jl @@ -0,0 +1,99 @@ +# Tests for the coordinate-invariant field-space (b, b̃, b̄) treatment of the +# PerturbedEquilibrium control-surface response matrices and the ForceFreeStates energy +# (issue #233 / Pharr 2026). +# +# The control-surface matrices are stored in the coordinate-invariant root-area-weighted +# field (b̃) space. With Σ ≡ sqrtamat (the √weight) and the scalar surface area A ≡ jarea: +# - b̃→b̄ operator S = Σ/√A (area-weighted field; b̄ = S·b̃) +# - bare → b̃ Σ = S·√A (b = Σ⁻¹·b̃) +# - b̃→flux conform R = S·A (Φ = R·b̃; poloidal flux is only the scalar Φ = A·b̄) +# These tests verify the transform rules in `field_space_response_matrices`, the field/flux +# recovery contracts, and the ForceFreeStates `c = jarea` energy refactor. + +using Test +using LinearAlgebra +using GeneralizedPerturbedEquilibrium + +const PE = GeneralizedPerturbedEquilibrium.PerturbedEquilibrium + +# Deterministic, well-conditioned synthetic operators (no RNG → reproducible). +function _synthetic_inputs(n::Int) + # Σ ≡ sqrtamat: a generic invertible complex operator standing in for the √weight convolution. + Σ = Matrix{ComplexF64}(undef, n, n) + for i in 1:n, j in 1:n + Σ[i, j] = (1.0 + 0.3 * cospi((i + 2j) / n)) * cis(0.2 * (i - j)) * (i == j ? 1.0 : 0.15) + end + jarea = 2.37 # scalar flux-surface area A + S = Σ ./ sqrt(jarea) # b̃→b̄ operator + # Hermitian positive-definite inductances (energy generators, flux space). + A = [cis(0.1 * (i - j)) / (1 + abs(i - j)) for i in 1:n, j in 1:n] + L = A * A' + n * I # Hermitian PD surface inductance + B = [cis(-0.07 * (i + j)) / (1 + abs(i - j)) for i in 1:n, j in 1:n] + Λ = B * B' + 2n * I # Hermitian PD plasma inductance + return Matrix{ComplexF64}(Σ), Matrix{ComplexF64}(S), jarea, Matrix{ComplexF64}(L), Matrix{ComplexF64}(Λ) +end + +@testset "Coordinate-invariant field-space response matrices" begin + n = 6 + Σ, S, jarea, L, Λ = _synthetic_inputs(n) + R = S .* jarea # b̃→flux conform Σ·√A + P = Λ / L # permeability Φ_tot = P·Φ_x + ϱ = inv(L) * (Λ - L) * inv(L) # reluctance + + fm = PE.field_space_response_matrices(Λ, L, P, ϱ, S, jarea) + + @testset "Flux recovery contract (round-trip via R = S·A)" begin + @test R * fm.permeability / R ≈ P rtol = 1e-10 + @test R * fm.surface_inductance * R' ≈ L rtol = 1e-10 + @test R * fm.plasma_inductance * R' ≈ Λ rtol = 1e-10 + @test (R') \ fm.reluctance / R ≈ ϱ rtol = 1e-10 + end + + @testset "Area-weighted (b̄) recovery via S (= flux/A²)" begin + # b̄-space inductance L_b̄ = S·L̃·S† = L/A² since S·R⁻¹ = A⁻¹·I. + @test S * fm.surface_inductance * S' ≈ L ./ jarea^2 rtol = 1e-10 + @test S * fm.plasma_inductance * S' ≈ Λ ./ jarea^2 rtol = 1e-10 + @test S * fm.permeability / S ≈ P rtol = 1e-10 # similarity: A cancels + end + + @testset "Internal consistency of the b̃ transform rules" begin + @test fm.permeability ≈ fm.plasma_inductance / fm.surface_inductance rtol = 1e-10 + L̃inv = inv(fm.surface_inductance) + @test fm.reluctance ≈ L̃inv * (fm.plasma_inductance - fm.surface_inductance) * L̃inv rtol = 1e-10 + end + + @testset "Energy-scalar invariance (flux ↔ b̃)" begin + Φ = ComplexF64[cis(0.3k) / k for k in 1:n] + b̃ = R \ Φ # root-area-weighted field + @test dot(Φ, L \ Φ) ≈ dot(b̃, fm.surface_inductance \ b̃) rtol = 1e-10 + @test dot(Φ, Λ \ Φ) ≈ dot(b̃, fm.plasma_inductance \ b̃) rtol = 1e-10 + @test dot(Φ, ϱ * Φ) ≈ dot(b̃, fm.reluctance * b̃) rtol = 1e-10 + end + + @testset "Three-field vector relations (b, b̃, b̄; flux = A·b̄)" begin + b̃ = ComplexF64[cis(0.21k) / (1 + k) for k in 1:n] + b̄ = S * b̃ # area-weighted field + b = (S .* sqrt(jarea)) \ b̃ # bare field b = Σ⁻¹·b̃ + Φ = R * b̃ # poloidal flux + @test Φ ≈ jarea .* b̄ rtol = 1e-12 # Φ = A·b̄ + @test b̄ ≈ Φ ./ jarea rtol = 1e-12 + @test (S .* sqrt(jarea)) * b ≈ b̃ rtol = 1e-12 # Σ·b = b̃ + end +end + +@testset "ForceFreeStates c = jarea energy refactor (net-zero)" begin + # The FFS energy is dW = c·b̃†·W_t·b̃ with W_t = (T⁻¹·Σ)†·W·(T⁻¹·Σ) and c = jarea. + # This must reproduce the eigenspectrum of the old flux operator M_old = T⁻¹·(Σ·√jarea). + n = 6 + Σ, _, jarea, _, _ = _synthetic_inputs(n) + A = [cis(0.13 * (i - j)) / (1 + abs(i - j)) for i in 1:n, j in 1:n] + W = A * A' + n * I # Hermitian PD energy matrix (ξ-space stand-in) + Tinv = ComplexF64[1.0 / (im * 2π * 1.7 * (i - 0.5)) for i in 1:n] # diag(T⁻¹) + + M = Tinv .* Σ # M = T⁻¹·Σ (new: pure √weight operator) + Mold = Tinv .* (Σ .* sqrt(jarea)) # M_old = T⁻¹·Σ·√jarea (old flux operator) + + e_new = sort(real.(eigvals(jarea .* (M' * W * M)))) + e_old = sort(real.(eigvals(Mold' * W * Mold))) + @test e_new ≈ e_old rtol = 1e-12 +end diff --git a/test/runtests_fullruns.jl b/test/runtests_fullruns.jl index 4eb166014..ac11df528 100644 --- a/test/runtests_fullruns.jl +++ b/test/runtests_fullruns.jl @@ -25,6 +25,8 @@ using HDF5 # not numeric regression — the mode-converged physical value is pinned in the # regression harness (examples/Solovev_kinetic_calculated_example). Assert only # nerfed-grid-robust facts: finite and positive (Solovev is stable). + # ξ-space (XiNorm) eigenvalues are read here because root-area weighting collapses + # these stable modes to ~0, making them unusable as a regression anchor. et = read(h5["FreeBoundaryStability/XiNorm/eigenmode_energies"]) @test isfinite(real(et[1])) @test real(et[1]) > 0 @@ -43,6 +45,8 @@ using HDF5 # harness on the mode-converged deck, not here — et[1] is a near-marginal, # ill-conditioned, FP-reassociation-sensitive quantity on this grid (sign not even # robust across platforms), so no value is pinned. + # ξ-space (XiNorm) eigenvalues are read here because root-area weighting collapses + # these modes to ~0, making them unusable as a regression anchor. et = read(h5["FreeBoundaryStability/XiNorm/eigenmode_energies"]) @test isfinite(real(et[1])) end