ForceFreeStates - BUG FIX - Power-normalized Jacobian-invariant eigenmode energies via generalized eigenproblem

benchmarks/benchmark_solovev_kinetic_stability.jl

@@ -14,7 +14,7 @@ band (`delta_mlow=delta_mhigh=0`), same `nl`, collision operator, grid, and
 edge truncation (`sas_flag`/`dmlim`).
 
 The comparison is the least-stable total-energy eigenvalue:
-  - Julia:   `FreeBoundaryStability/XiNorm/eigenmode_energies[1]` in `gpec.h5`
+  - Julia:   `FreeBoundaryStability/eigenmode_energies[1]` in `gpec.h5`
   - Fortran: the "Energies: ... real = …, imaginary = …" line printed by `dcon`
              (== `W_t_eigenvalue[1]`; parsing the log avoids a NetCDF dependency).
 
@@ -157,7 +157,7 @@ function run_julia_reference()
     GPE.main([rundir])
     wall = time() - t0
     et = h5open(joinpath(rundir, "gpec.h5"), "r") do h5
-        read(h5["FreeBoundaryStability/XiNorm/eigenmode_energies"])
+        read(h5["FreeBoundaryStability/eigenmode_energies"])
     end
     return real(et[1]), imag(et[1]), wall
 end

benchmarks/compare_gal_vs_el.jl

@@ -3,7 +3,7 @@
 #
 # EL :  ξ_EL(ψ)  = U_EL(ψ) · (U_EL_edge \ w)          (fundamental matrix, integration/xi_psi)
 # gal:  ξ_gal(ψ) = U_gal(ψ) · w                       (identity-at-edge ⇒ coefficient is w itself)
-# w = eigenvector of the total energy operator W = W_plasma + W_vacuum (FreeBoundaryStability/XiNorm).
+# w = eigenvector of the total energy operator W = W_plasma + W_vacuum (FreeBoundaryStability).
 #
 # Needs ONE gpec.h5 from a run with populate_dense_xi=true (EL dense u_store), gal_match_flag=true,
 # gal_ideal_flag=true (gal ideal matched), vac_flag=true (energy operator).
@@ -19,8 +19,8 @@ ksel = length(ARGS) >= 3 ? ARGS[3] : "highest"
 to_c(a) = eltype(a) <: Complex ? ComplexF64.(a) : map(x -> ComplexF64(x.re, x.im), a)
 
 et, wt, u1, psiE, gxi, psiG, issing, mlow, sing_psi = h5open(h5path) do f
-    (to_c(read(f["FreeBoundaryStability/XiNorm/eigenmode_energies"])),
-        to_c(read(f["FreeBoundaryStability/XiNorm/W_freeboundary_eigenmodes"])),
+    (to_c(read(f["FreeBoundaryStability/eigenmode_energies"])),
+        to_c(read(f["FreeBoundaryStability/W_freeboundary_eigenmodes"])),
         to_c(read(f["integration/xi_psi"])), read(f["integration/psi"]),
         to_c(read(f["galerkin/match/xi"])), read(f["galerkin/solution/psi"]),
         Bool.(read(f["galerkin/solution/issing"])), read(f["info/mlow"]),

benchmarks/plot_xi_eigenmode.jl

@@ -17,8 +17,8 @@ outpng = length(ARGS) >= 2 ? ARGS[2] : joinpath(@__DIR__, "xi_eigenmode.png")
 to_c(a) = eltype(a) <: Complex ? ComplexF64.(a) : map(x -> ComplexF64(x.re, x.im), a)
 
 et, wt, u1, psi, mlow, sing_psi = h5open(h5path) do f
-    (to_c(read(f["FreeBoundaryStability/XiNorm/eigenmode_energies"])),
-        to_c(read(f["FreeBoundaryStability/XiNorm/W_freeboundary_eigenmodes"])),
+    (to_c(read(f["FreeBoundaryStability/eigenmode_energies"])),
+        to_c(read(f["FreeBoundaryStability/W_freeboundary_eigenmodes"])),
         to_c(read(f["integration/xi_psi"])),
         read(f["integration/psi"]),
         read(f["info/mlow"]),

docs/src/conventions.md

@@ -318,9 +318,16 @@ The internal flux-conform operator is just ``R = \Sigma\sqrt{A} = `` `rootarea_t
   `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``.
+  `surface_area` ``A`` recovers flux (``\Phi = A\,\bar b``).
+- **Eigenmode energies** (`FreeBoundaryStability/eigenmode_energies`) — the generalized eigenvalues of the
+  pencil ``W\,v = \lambda\,N\,v``, where ``N`` is the power-normalization (surface-norm) matrix built from
+  the Jacobian Fourier coefficients: ``\xi^\dagger N\,\xi = \oint J\,|\xi(\theta)|^2\,d\theta\,/\,V' =
+  \langle|\xi|^2\rangle``. The eigenvalues are stationary values of the power quotient
+  ``\xi^\dagger W \xi / \xi^\dagger N \xi`` — energies per unit surface-averaged squared displacement
+  (a normalization convention, not raw joules). Because ``W`` and ``N`` transform by the same congruence
+  under a change of working (Jacobian) coordinate, the spectrum is invariant to that coordinate choice.
+  ``N`` is not an inertia matrix (the stability stage carries no velocities); it is the DCON power
+  normalization promoted from a post-hoc per-eigenvector division into the metric of the eigenproblem.
 
 
 ## Rotation Velocity Conventions (KineticForces)

regression-harness/cases/diiid_n1.toml

@@ -8,44 +8,44 @@ 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 root-area-weighted (Φ-space, Jacobian-invariant); ξ-space counterparts live under FreeBoundaryStability/XiNorm/ and are not tracked.
+# Energies — leading generalized (W,N) pencil eigenvalues at the final truncation (psilim): power-normalized (⟨|ξ|²⟩ = 1 metric) and invariant to the working-coordinate Jacobian.
 [quantities.et_real]
 h5path = "FreeBoundaryStability/eigenmode_energies"
 type = "complex_vector"
 extract = "real_first"
-label = "root-area-weighted total energy Re(et[1])"
+label = "total energy Re(et[1])"
 noise_threshold = 1e-10
 order = 10
 
 [quantities.et_imag]
 h5path = "FreeBoundaryStability/eigenmode_energies"
 type = "complex_vector"
 extract = "imag_first"
-label = "root-area-weighted total energy Im(et[1])"
+label = "total energy Im(et[1])"
 noise_threshold = 1e-10
 order = 11
 
 [quantities.ep_real]
 h5path = "FreeBoundaryStability/eigenmode_plasma_energies"
 type = "complex_vector"
 extract = "real_first"
-label = "root-area-weighted plasma energy Re(ep[1])"
+label = "plasma energy Re(ep[1])"
 noise_threshold = 1e-10
 order = 12
 
 [quantities.ev_real]
 h5path = "FreeBoundaryStability/eigenmode_vacuum_energies"
 type = "complex_vector"
 extract = "real_first"
-label = "root-area-weighted vacuum energy Re(ev[1])"
+label = "vacuum energy Re(ev[1])"
 noise_threshold = 1e-10
 order = 13
 
 [quantities.vacuum_eigenvalue]
-h5path = "FreeBoundaryStability/XiNorm/vacuum_eigenvalue"
+h5path = "FreeBoundaryStability/vacuum_eigenvalue"
 type = "real_scalar"
 extract = "value"
-label = "vacuum matrix min eigenvalue (XiNorm)"
+label = "vacuum matrix min eigenvalue"
 noise_threshold = 1e-10
 order = 14
 
@@ -54,23 +54,23 @@ order = 14
 h5path = "FreeBoundaryStability/eigenmode_plasma_energies"
 type = "complex_vector"
 extract = "all_complex"
-label = "root-area-weighted plasma energy (all)"
+label = "plasma energy (all)"
 noise_threshold = 1e-10
 order = 20
 
 [quantities.ev_all]
 h5path = "FreeBoundaryStability/eigenmode_vacuum_energies"
 type = "complex_vector"
 extract = "all_complex"
-label = "root-area-weighted vacuum energy (all)"
+label = "vacuum energy (all)"
 noise_threshold = 1e-10
 order = 21
 
 [quantities.et_all]
 h5path = "FreeBoundaryStability/eigenmode_energies"
 type = "complex_vector"
 extract = "all_complex"
-label = "root-area-weighted total energy (all)"
+label = "total energy (all)"
 noise_threshold = 1e-10
 order = 22
 

regression-harness/cases/solovev_kinetic_calculated.toml

@@ -9,36 +9,36 @@ example_dir = "examples/Solovev_kinetic_calculated_example"
 
 # Energies — leading eigenvalues
 [quantities.et_real]
-h5path = "FreeBoundaryStability/XiNorm/eigenmode_energies"
+h5path = "FreeBoundaryStability/eigenmode_energies"
 type = "complex_vector"
 extract = "real_first"
 label = "total energy Re(et[1])"
 noise_threshold = 1e-10
 
 [quantities.et_imag]
-h5path = "FreeBoundaryStability/XiNorm/eigenmode_energies"
+h5path = "FreeBoundaryStability/eigenmode_energies"
 type = "complex_vector"
 extract = "imag_first"
 label = "total energy Im(et[1])"
 noise_threshold = 1e-10
 
 [quantities.ep_real]
-h5path = "FreeBoundaryStability/XiNorm/eigenmode_plasma_energies"
+h5path = "FreeBoundaryStability/eigenmode_plasma_energies"
 type = "complex_vector"
 extract = "real_first"
 label = "plasma energy Re(ep[1])"
 noise_threshold = 1e-10
 
 [quantities.ev_real]
-h5path = "FreeBoundaryStability/XiNorm/eigenmode_vacuum_energies"
+h5path = "FreeBoundaryStability/eigenmode_vacuum_energies"
 type = "complex_vector"
 extract = "real_first"
 label = "vacuum energy Re(ev[1])"
 noise_threshold = 1e-10
 
 # Energies — full eigenvalue arrays
 [quantities.et_all]
-h5path = "FreeBoundaryStability/XiNorm/eigenmode_energies"
+h5path = "FreeBoundaryStability/eigenmode_energies"
 type = "complex_vector"
 extract = "all_complex"
 label = "total energy (all)"

regression-harness/cases/solovev_kinetic_nuzero.toml

@@ -9,36 +9,36 @@ example_dir = "examples/Solovev_kinetic_calculated_example"
 
 # Energies — leading eigenvalues
 [quantities.et_real]
-h5path = "FreeBoundaryStability/XiNorm/eigenmode_energies"
+h5path = "FreeBoundaryStability/eigenmode_energies"
 type = "complex_vector"
 extract = "real_first"
 label = "total energy Re(et[1])"
 noise_threshold = 1e-10
 
 [quantities.et_imag]
-h5path = "FreeBoundaryStability/XiNorm/eigenmode_energies"
+h5path = "FreeBoundaryStability/eigenmode_energies"
 type = "complex_vector"
 extract = "imag_first"
 label = "total energy Im(et[1])"
 noise_threshold = 1e-10
 
 [quantities.ep_real]
-h5path = "FreeBoundaryStability/XiNorm/eigenmode_plasma_energies"
+h5path = "FreeBoundaryStability/eigenmode_plasma_energies"
 type = "complex_vector"
 extract = "real_first"
 label = "plasma energy Re(ep[1])"
 noise_threshold = 1e-10
 
 [quantities.ev_real]
-h5path = "FreeBoundaryStability/XiNorm/eigenmode_vacuum_energies"
+h5path = "FreeBoundaryStability/eigenmode_vacuum_energies"
 type = "complex_vector"
 extract = "real_first"
 label = "vacuum energy Re(ev[1])"
 noise_threshold = 1e-10
 
 # Energies — full eigenvalue arrays
 [quantities.et_all]
-h5path = "FreeBoundaryStability/XiNorm/eigenmode_energies"
+h5path = "FreeBoundaryStability/eigenmode_energies"
 type = "complex_vector"
 extract = "all_complex"
 label = "total energy (all)"

regression-harness/cases/solovev_multi_n.toml

@@ -7,28 +7,28 @@ 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 root-area-weighted (Φ-space, Jacobian-invariant); ξ-space counterparts live under FreeBoundaryStability/XiNorm/ and are not tracked.
+# Energies — leading generalized (W,N) pencil eigenvalues at the final truncation (psilim): power-normalized (⟨|ξ|²⟩ = 1 metric) and invariant to the working-coordinate Jacobian.
 [quantities.et_real]
 h5path = "FreeBoundaryStability/eigenmode_energies"
 type = "complex_vector"
 extract = "real_first"
-label = "root-area-weighted total energy Re(et[1])"
+label = "total energy Re(et[1])"
 noise_threshold = 1e-10
 order = 10
 
 [quantities.et_imag]
 h5path = "FreeBoundaryStability/eigenmode_energies"
 type = "complex_vector"
 extract = "imag_first"
-label = "root-area-weighted total energy Im(et[1])"
+label = "total energy Im(et[1])"
 noise_threshold = 1e-10
 order = 11
 
 [quantities.vacuum_eigenvalue]
-h5path = "FreeBoundaryStability/XiNorm/vacuum_eigenvalue"
+h5path = "FreeBoundaryStability/vacuum_eigenvalue"
 type = "real_scalar"
 extract = "value"
-label = "vacuum matrix min eigenvalue (XiNorm)"
+label = "vacuum matrix min eigenvalue"
 noise_threshold = 1e-10
 order = 14
 
@@ -37,23 +37,23 @@ order = 14
 h5path = "FreeBoundaryStability/eigenmode_plasma_energies"
 type = "complex_vector"
 extract = "all_complex"
-label = "root-area-weighted plasma energy (all)"
+label = "plasma energy (all)"
 noise_threshold = 1e-10
 order = 20
 
 [quantities.ev_all]
 h5path = "FreeBoundaryStability/eigenmode_vacuum_energies"
 type = "complex_vector"
 extract = "all_complex"
-label = "root-area-weighted vacuum energy (all)"
+label = "vacuum energy (all)"
 noise_threshold = 1e-10
 order = 21
 
 [quantities.et_all]
 h5path = "FreeBoundaryStability/eigenmode_energies"
 type = "complex_vector"
 extract = "all_complex"
-label = "root-area-weighted total energy (all)"
+label = "total energy (all)"
 noise_threshold = 1e-10
 order = 22
 

regression-harness/cases/solovev_n1.toml

@@ -7,44 +7,44 @@ 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 root-area-weighted (Φ-space, Jacobian-invariant); ξ-space counterparts live under FreeBoundaryStability/XiNorm/ and are not tracked.
+# Energies — leading generalized (W,N) pencil eigenvalues at the final truncation (psilim): power-normalized (⟨|ξ|²⟩ = 1 metric) and invariant to the working-coordinate Jacobian.
 [quantities.et_real]
 h5path = "FreeBoundaryStability/eigenmode_energies"
 type = "complex_vector"
 extract = "real_first"
-label = "root-area-weighted total energy Re(et[1])"
+label = "total energy Re(et[1])"
 noise_threshold = 1e-10
 order = 10
 
 [quantities.et_imag]
 h5path = "FreeBoundaryStability/eigenmode_energies"
 type = "complex_vector"
 extract = "imag_first"
-label = "root-area-weighted total energy Im(et[1])"
+label = "total energy Im(et[1])"
 noise_threshold = 1e-10
 order = 11
 
 [quantities.ep_real]
 h5path = "FreeBoundaryStability/eigenmode_plasma_energies"
 type = "complex_vector"
 extract = "real_first"
-label = "root-area-weighted plasma energy Re(ep[1])"
+label = "plasma energy Re(ep[1])"
 noise_threshold = 1e-10
 order = 12
 
 [quantities.ev_real]
 h5path = "FreeBoundaryStability/eigenmode_vacuum_energies"
 type = "complex_vector"
 extract = "real_first"
-label = "root-area-weighted vacuum energy Re(ev[1])"
+label = "vacuum energy Re(ev[1])"
 noise_threshold = 1e-10
 order = 13
 
 [quantities.vacuum_eigenvalue]
-h5path = "FreeBoundaryStability/XiNorm/vacuum_eigenvalue"
+h5path = "FreeBoundaryStability/vacuum_eigenvalue"
 type = "real_scalar"
 extract = "value"
-label = "vacuum matrix min eigenvalue (XiNorm)"
+label = "vacuum matrix min eigenvalue"
 noise_threshold = 1e-10
 order = 14
 
@@ -53,23 +53,23 @@ order = 14
 h5path = "FreeBoundaryStability/eigenmode_plasma_energies"
 type = "complex_vector"
 extract = "all_complex"
-label = "root-area-weighted plasma energy (all)"
+label = "plasma energy (all)"
 noise_threshold = 1e-10
 order = 20
 
 [quantities.ev_all]
 h5path = "FreeBoundaryStability/eigenmode_vacuum_energies"
 type = "complex_vector"
 extract = "all_complex"
-label = "root-area-weighted vacuum energy (all)"
+label = "vacuum energy (all)"
 noise_threshold = 1e-10
 order = 21
 
 [quantities.et_all]
 h5path = "FreeBoundaryStability/eigenmode_energies"
 type = "complex_vector"
 extract = "all_complex"
-label = "root-area-weighted total energy (all)"
+label = "total energy (all)"
 noise_threshold = 1e-10
 order = 22
 

src/ForceFreeStates/EulerLagrange.jl

@@ -939,9 +939,6 @@ 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 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.
     # Steps where free_compute_total hits a singular wp solve are left as NaN per the EdgeScanState contract.
     for j in 1:N_edge
@@ -954,10 +951,6 @@ 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.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