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

@@ -317,9 +317,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)

examples/Solovev_kinetic_calculated_example/gpec.toml

@@ -1,6 +1,5 @@
 [Equilibrium]
 eq_type = "sol"                         # Type of the input 2D equilibrium file
-eq_filename = "sol.toml"                # Path to equilibrium file
 jac_type = "pest"                       # Coordinate system (hamada, pest, boozer, equal_arc)
 power_bp = 0                            # Poloidal field power exponent for Jacobian
 power_b = 0                             # Toroidal field power exponent for Jacobian
@@ -48,3 +47,15 @@ singfac_min = 1e-4            # Fractional distance from rational q at which ide
 ucrit = 1e3                   # Maximum fraction of solutions allowed before re-normalized
 write_outputs_to_HDF5 = true
 verbose = false
+
+[SOL_INPUT]
+mr = 128       # Number of radial grid zones
+mz = 128       # Number of axial grid zones
+ma = 128       # Number of flux grid zones
+e = 1.6        # Elongation
+a = 0.33       # Minor radius
+r0 = 1.0       # Major radius
+q0 = 1.9       # Safety factor at the magnetic axis (O-point)
+p0fac = 1      # Scale on-axis pressure (β changes; Φ, q fixed)
+b0fac = 1      # Scale toroidal field at constant β (Bt changes; shape, β fixed)
+f0fac = 1      # Scale toroidal field at constant pressure (β, q change; Φ, p, Bp fixed)

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