Jacobian-invariant power-norm flux eigenvalues

[logan-nc]

Summary

• Adds a power-normalized flux (Φ-space) eigenvalue channel to the FFS edge scan. The energy matrix W is transformed via W_Φ = M† W M where M = T⁻¹ · ptof, with T = diag(i·2π·χ₁·(m−n·q)) and ptof = sqrtamat · √jarea. The eigenvalues of W_Φ are the coordinate-independent energy — they do not depend on Jacobian choice.
• Writes edge_scan/pn_total_energy, pn_plasma_energy, pn_vacuum_energy alongside the existing et/ep/ev ξ-space channels.
• Accompanying unit-level test scripts/test_power_norm_invariance.jl establishes a numerical floor of ~6.9e-7 for Jacobian invariance of the underlying operators (√weight identity and convmat angle-map) at psihigh=0.994, mtheta=1024.
• scripts/compare_jacobians_power_norm.jl runs the full stack for {hamada, pest, park, boozer} and reports residual spreads at matched ψ against the unit-test floor.

What is invariant, what isn't

Only the eigenspectrum of W_Φ is Jacobian-invariant. The Φ-space eigenvectors and the θ-space field reconstructed from Φ still carry Jacobian dependence (T depends on m-labelling; ‖Φ‖ = ‖M⁻¹ξ‖ drifts). The module docstring in PowerNorm.jl spells this out explicitly.

Observed behaviour

• n=1 (full stack, mtheta=1024): all four jacobians agree to ~15 % at edge probe surfaces.
• n=3: hamada/pest/boozer cluster to ~20 %; park is anomalous (sign-flipping at 3/5 surfaces) — consistent with θ-grid aliasing of high-m edge modes under park's J ∝ B weighting, which concentrates grid density on the high-field side. Likely fixable with larger mtheta; not pursued in this PR.

Test plan

• regress --cases diiid_n1 --refs develop,feature/power-norm-eigenvalues — confirm no baseline drift in existing scalars.
• julia --project=. scripts/test_power_norm_invariance.jl — confirm operator floor ~7e-7.
• julia --project=. scripts/compare_jacobians_power_norm.jl (with const NN = 1) — confirm n=1 full-stack residual spread ~15 %.
• Inspect examples/DIIID-like_ideal_example/jacobian_comparison_n1.png and gpec_edge_eigenvalues_n1.png.

🤖 Generated with Claude Code

[logan-nc]

This branch makes the power normalized flux the default space for reported free boundary eigenmodes and eigenvalues. This removes the jacobian dependence of the ξ based decomposition used in DCON. While coordinate dependence was a problem for running in different jacobians, that wasn't the real issue. The real issue is that the actual magnetic coordinate system of even the same jacobian evolves sharply in the far edge of the plasma - so using a coordinate dependent normalization results in wildly unintuitive variations of the result with small differences in the truncation. When we fix this by going into a coordinate invariant normalization (i.e. dW = Φ•W•Φ where √(J∇ψ)b = ∑Φexp(i(mθ - nφ)) then the eigenvalue energies suddenly don't depend on the truncation so much! This makes more sense and makes the choice of truncation much less stressful (we still need to avoid false external kink instabilities from truncating just inside a rational though).

Here is a n=1 edge truncation scan comparison:

Here is n=3:

Ignore the park coordinate result in the n=3 one. It's clearly a numerical instability from my choice of m and/or the equilibrium reconstruction. It isn't worth fixing - the main point is clear.

[logan-nc]

HDF5 schema updates in the latest commits

The edge_scan/ and vacuum/ groups have been restructured into a CamelCase layout with power-normalized flux (Φ-space) values as the default. ξ-space values are preserved under XiNorm/ subgroups for Fortran benchmarking.

Key changes

• edge_scan/ → EdgeScan/ — top-level fields are Φ-space; EdgeScan/XiNorm/ holds ξ-space counterparts.
• vacuum/ → FreeBoundaryStability/ + SurfaceGeometries/ — energy matrices and eigenmodes move to FreeBoundaryStability/; Cartesian point clouds move to SurfaceGeometries/{Plasma,Wall}/{x,y,z}.
• W_freeboundary, W_plasma, W_vacuum — now stored in Φ-space (M†·W·M, with M = T⁻¹·sqrtamat·√jarea). ξ-space versions live under XiNorm/.
• W_freeboundary_eigenmodes — now the Φ-space eigenvector matrix of W_freeboundary, with columns sorted most-unstable-first and phase-normalized so each column's largest-magnitude entry is real-positive. Columns satisfy W_freeboundary · v[i] = eigenmode_energies[i] · v[i] to ~1e-13. The ξ-space eigenvector matrix now lives at XiNorm/W_freeboundary_eigenmodes.
• New diagnosticEdgeScan/vacuum_eigenvalue — smallest eigenvalue of W_vacuum in Φ-space, clamped at zero.

Full gpec.h5 tree (Solovev n=1 example)

gpec.h5
├── EdgeScan/ (power-norm default)
│ ├── psi, q
│ ├── total_energy, plasma_energy, vacuum_energy, vacuum_eigenvalue
│ └── XiNorm/
│ └── total_energy, plasma_energy, vacuum_energy, vacuum_eigenvalue
├── FreeBoundaryStability/ (power-norm default)
│ ├── W_freeboundary (Φ: M†·wt0·M)
│ ├── W_plasma (Φ: M†·wp·M)
│ ├── W_vacuum (Φ: M†·wv·M)
│ ├── W_freeboundary_eigenmodes (Φ: eigvecs of W_freeboundary,
│ │ sorted+phased)
│ ├── eigenmode_energies
│ ├── eigenmode_plasma_energies
│ ├── eigenmode_vacuum_energies
│ └── XiNorm/
│ ├── W_freeboundary (ξ: wt0)
│ ├── W_plasma (ξ: wp)
│ ├── W_vacuum (ξ: wv)
│ ├── W_freeboundary_eigenmodes (ξ: wt, kinetic-normalized)
│ ├── eigenmode_energies
│ ├── eigenmode_plasma_energies
│ └── eigenmode_vacuum_energies
├── SurfaceGeometries/
│ ├── Plasma/{x, y, z}
│ └── Wall/{x, y, z}
│
├── equil/ # → follow-on rename to Equilibrium/ (issue #226)
├── info/ # → Info/ (issue #226)
├── input/ # → Input/ (holds EQUIL_CONTROL/, ForceFreeStates/)
├── integration/ # → EulerLagrangeSolutions/ (issue #226)
├── locstab/ # → LocalStability/ (issue #226)
├── matrices/ # → Matrices/ (holds ideal/) (issue #226)
├── singular/ # → SingularSurfaces/ (issue #226)
└── splines/ # → Splines/ (holds profiles/, rzphi/) (issue #226)

Consumers updated to match: src/Analysis/ForceFreeStates.jl, the benchmarks, the fullruns smoke test, the regression-harness/cases/*.toml files, and all scripts/. Kinetic Solovev test thresholds rebaselined to Φ-space values. Issue #226 tracks the remaining lowercase-group renames and the CLAUDE.md naming-convention section.

[logan-nc]

Development scripts used during this PR (not being merged)

During development I wrote a few exploratory scripts to prove out Jacobian invariance of the power-norm eigenspectrum and to pin down numerical settings that actually exercise invariance. Documenting them here for posterity before removing them — they were scaffolding, not product.

scripts/test_jac_invariance.jl — the decisive first test (pivotal)

Small harness that ran the full DIIID-like example once per Jacobian (hamada, pest) through GeneralizedPerturbedEquilibrium.main, then read edge_scan/total_energy (ξ-space) and edge_scan/pn_total_energy (Φ-space / √A-normalized) from each gpec.h5 and compared them at the edge. This is the test that first confirmed the √A-normalized (power-norm) eigenspectrum is Jacobian-invariant in the full stack, not just in isolated operators.

Representative result at psilim (DIIID-like, n=1):

quantity	hamada	pest	rel. diff across jac
ξ-space et[1]	+1.5047e+00	≈ +1.24e+00	8–154% across the edge scan
power-norm pn_et[1]	+1.5787e+00	≈ +1.579e+00	0.1–0.5% across the edge scan

The ξ-space spread was O(10%) at best and grew to >100% near rational surfaces, while the power-norm spread sat at ~10⁻³ across the smooth interior of the edge scan. That disparity is what motivated the rest of the PR — all the downstream schema work (CamelCase restructure, XiNorm/ subgroup, W_freeboundary* in Φ-space by default) exists because this script showed Φ-space is the invariant you want to publish.

scripts/test_power_norm_invariance.jl — unit-level invariance floor

Isolated-operator test of the √weight transform and convmat angle-map across Jacobians, short-circuiting the full Euler–Lagrange pipeline. Established the numerical precision floor at:

• PSIHIGH=0.994, N_TOR=3, MLOW=-128, MHIGH=128, MTHETA=1024
• Floor: ~6.9e-7 relative area residual
• Both the √weight operator and the convmat angle-map agreed across Jacobians to within ~2× the floor

This set the bar for what "invariant" means numerically in this codebase and told us which knobs (mtheta, psihigh, absolute m-band) were actually binding on the full-stack comparison.

scripts/compare_jacobians_power_norm.jl — full-stack 4-way comparison

Runs the full pipeline on the DIIID-like example four times (hamada / pest / boozer / park), all with apples-to-apples overrides:

psihigh = 0.994 # inside the proven-invariant regime
mtheta = 1024 # matches unit-test θ-grid
cyl_flag = true # absolute mlow/mhigh across jacobians
delta_mlow = -48
delta_mhigh = 48 # mpert = 97 on all four

Before these overrides, the four runs silently disagreed because cyl_flag=false makes mlow/mhigh depend on each Jacobian's qmin/qmax, and mtheta=256 left aliasing as a free parameter. With the overrides, the Φ-space eigenvalue spread across Jacobians collapsed toward the unit-test floor.

scripts/compare_eigenvalues_edge_scan.jl — ξ-space vs Φ-space visualization

Plots the edge-scan eigenvalues from a gpec.h5 in both representations side-by-side. The punchline:

• ξ-space (FreeBoundaryStability/XiNorm/*): curves disagree across Jacobians by O(1).
• Φ-space (FreeBoundaryStability/*, the new default): curves overlap across Jacobians to the floor.

This is the visual proof of the thing the schema restructure is built around — only the Φ-space spectrum is invariant, so that's what ships as the default output.

Throwaway diagnostic

• scripts/diag_pn.jl — ad-hoc prints while wiring pn_wt0/pn_wp/pn_wv into VacuumData (sqrtamat Hermiticity, Toeplitz structure, M-matrix operator norm, singfac values).

Key results these scripts produced

• First-of-a-kind confirmation (test_jac_invariance.jl): √A-normalized (power-norm) eigenvalues are Jacobian-invariant at ~0.1–0.5% across the edge scan, vs ξ-space eigenvalues at 8–154% — a 2–3 order-of-magnitude improvement.
• Numerical floor on invariance at the operator level: ~6.9e-7.
• Park n=3 anomaly (seen early on) was confirmed to be θ-aliasing, not a correctness bug — n=1 with mtheta=1024 is well-behaved and agrees with the other Jacobians at the floor.
• W_plasma + W_vacuum = W_freeboundary in Φ-space holds to 8.9e-16 (additivity preserved by the M†·W·M transform).
• W_freeboundary · v = λ · v verified to 2.6e-13 on Solovev n=1 for the stored Φ-space eigenmodes.

Being removed

These scripts live in scripts/ and were tracked in commits ffad3994 and ab14b09f. Removing them (and the untracked diagnostic diag_pn.jl) in the next commit — they've served their purpose.

[logan-nc]

@ebursch@matt-pharr this is ready to review - you have both thought about these things a lot. It's cool to see them proven and working to our benefit here.

[logan-nc]

@ebursch and @matt-pharr please mark formal reviews here!

Something to think about as you look at it and think about your paper: Should we use sqrt(A) weighted fields (Phi-tilde) or sqrt(A) normalized displacements (xi-tilde? More consistent with DCON? Still invariant?). Note the T matrix above translates xi to flux and has a m-nq in it, which could do funny things near rationals.

[logan-nc]

@ebursch and @matt-pharr please submit reviews today - it's been 1 week since our big meeting where we agreed to close PRs. 90% of the line changes are just h5 variable name improvements. There is really just the one new PowerNorm.jl file and corresponding renormalization of the energy matrices in free.jl to double check.

[ebursch]

This all looks good to me. Made a few minor suggestions to improve codebase consistency. I will leave approval to @matt-pharr since I am less familiar with the broader code generally.

[logan-nc]

@matt-pharr I toggled on auto-merge based on @ebursch's positive review... If the checks fail, that gives you a chance to speak now or forever hold your peace. Otherwise, I'll fix and merge tomorrow