PerturbedEquilibrium - FEATURE - Decompose the plasma energy of free-boundary eigenmodes on the (psi,theta) grid

docs/development/architecture.md

@@ -92,6 +92,7 @@ Splines are provided by the external `FastInterpolations` package rather than by
      - `ResponseMatrices.jl` - Permeability matrix calculation
      - `FieldReconstruction.jl` - Mode-space field reconstruction
      - `Response.jl` - Plasma response computation
+     - `EnergyDecomposition.jl` - δW_p of free-boundary eigenmodes on the (ψ,θ) grid in the effective-field and standard forms, run as its own `[EnergyDecomposition]` stage
      - `SingularCoupling.jl` - **Singular surface analysis** including:
        - Delta prime (Δ') tearing stability parameter
        - Resonant flux and currents at rational surfaces

docs/development/hdf5-conventions.md

@@ -43,7 +43,7 @@ Top level (11 groups):
 | `ForceFreeStates/` | `Solutions/ForwardIntegration/` (u-solutions), `Solutions/GalerkinIntegration/` (closed ξ profiles in the shared layout, `Match/` diagnostics, the gal surface list, debug-gated `Basis/`), `EulerLagrangeMatrices/{Ideal,Kinetic}`, `FreeBoundaryStability/`, `EdgeScan/` |
 | `LocalStability/` | Mercier `D_I`, resistive interchange `D_R`, `ballooning_Delta_prime` on `psi`; the ballooning α boundary on `ballooning_psi` |
 | `SingularSurfaces/` | Per-rational-surface data: `rational_psi`/`rational_q`/`rational_m`/`rational_n`, GGJ coefficients, `Delta_prime_matrix`/`Delta_prime_raw`/`Delta_coil`/`pest3_A`/`pest3_B`/`pest3_Gamma` (Riccati or Galerkin alike), `Kinetic/` |
-| `PerturbedEquilibrium/` | `ForcingModes/`, `Response/`, `ResponseMatrices/`, `SingularCoupling/`, `Energies/`, control-surface spectra |
+| `PerturbedEquilibrium/` | `ForcingModes/`, `Response/`, `ResponseMatrices/`, `SingularCoupling/`, `Energies/`, `EnergyDecomposition/` (δW_p of free-boundary eigenmodes split into its terms on an `eigenmode_index` coordinate, optional `Densities/`), control-surface spectra |
 | `KineticForces/` | `<method>/` (torque/energy profiles, `EnergyIntegrals/`, `KineticMatrices/`); multi-ion runs add `PerSpecies/<species>/<method>/` with the same per-method layout, summing to the top-level total |
 | `ErrorFields/` | `CoilSensitivities/` (per-coil-set control-surface spectra and their rigid shift/tilt derivatives, `DominantMode/` full-window projection); `MonteCarlo/` (intrinsic and corrected `\|δ\|` histograms over the sampled tolerances, per batch and averaged); `Risk/` (threshold density, `P(lock\|δ)`, locking probabilities, `ToleranceScan/`); `NTV/` (correction-coil overlap and NTV torque couplings per kAt) |
 | `Tearing/` | `PerSurface/` (+ `DpMatrix/`), `Roots/`, `LayerWidths/`, `Diagnostics/{ValidRoots,Poles,FilteredRoots}`, `Scan/Surface_<k>/` |

docs/development/references.md

@@ -68,6 +68,30 @@ The PerturbedEquilibrium module implements GPEC-style perturbed equilibrium calc
   - Published: Physics of Plasmas **24**, 032505 (2017)
   - Describes: Self-consistent coupling with neoclassical effects
 
+The energy decomposition of free-boundary eigenmodes (`EnergyDecomposition.jl`) follows:
+
+- **Bernstein et al. (1958)**: "An energy principle for hydromagnetic stability problems"
+  - Location: not in `docs/resources/`
+  - Published: Proceedings of the Royal Society A **244**, 17 (1958)
+  - Link: https://doi.org/10.1098/rspa.1958.0023
+  - Describes: The ideal-MHD energy principle; the standard form of δW_p evaluated by the energy decomposition
+
+- **Bernstein (1983)**: in *Basic Plasma Physics I*, Handbook of Plasma Physics Vol. 1, edited by A. A. Galeev and R. N. Sudan (North-Holland, Amsterdam), p. 421
+  - Location: not in `docs/resources/`
+  - Describes: Derivation of the effective-field (Bernstein) form |δB_eff|²/μ₀ − K|ξ_n|² of δW_p
+
+- **Chance et al. (1992)**: "MHD stability of tokamak plasmas"
+  - Location: not in `docs/resources/`
+  - Published: PPPL-CFP-2687, 2nd Symposium on Plasma Dynamics, Trieste (1992)
+  - Link: https://www.osti.gov/biblio/7234275
+  - Describes: Both forms of δW_p (Eqs. 12a–12c) and the kernel K = K₁ + K₂ + K₃ (Eqs. 15–17): shear–current |∇ψ|²σS, parallel current μ₀B²σ², pressure–curvature 2p'κ_ψ
+
+- **Boozer (2004)**: "Physics of magnetically confined plasmas"
+  - Location: not in `docs/resources/`
+  - Published: Reviews of Modern Physics **76**, 1071 (2004)
+  - Link: https://doi.org/10.1103/RevModPhys.76.1071
+  - Describes: Sec. V.C, the energy principle, the effective field C = δB + ξ_n(μ₀ j × n̂) (Eq. 98), the (∇×C)·∇ψ = 0 constraint (Eq. 99) and w = C²/2μ₀ − w_d (Eq. 100)
+
 ## InnerLayer Module (SLAYER)
 
 The SLAYER slab inner-layer solver and the resistive-layer width diagnostics are based on:

docs/make.jl

@@ -37,6 +37,7 @@ makedocs(;
             "KineticForces" => "kinetic_forces.md",
             "Forcing Terms" => "forcing_terms.md",
             "Perturbed Equilibrium" => "perturbed_equilibrium.md",
+            "Energy Decomposition" => "energy_decomposition.md",
             "Error Fields" => "error_fields.md",
             "Tearing" => "inner_layer.md",
             "Analysis" => "analysis.md",

docs/src/api.md

@@ -97,6 +97,7 @@ remain TOML-driven.
 GeneralizedPerturbedEquilibrium.EulerLagrangeProblem
 GeneralizedPerturbedEquilibrium.solve
 GeneralizedPerturbedEquilibrium.perturbed_equilibrium
+GeneralizedPerturbedEquilibrium.energy_decomposition
 ```
 
 ## Integrator selectors

docs/src/citations.md

@@ -82,6 +82,39 @@ Establishes the area weightings that make flux-surface Fourier spectra of pertur
 
 ---
 
+### Energy decomposition of free-boundary eigenmodes
+
+> I. B. Bernstein, E. A. Frieman, M. D. Kruskal, and R. M. Kulsrud, "An energy principle for hydromagnetic stability problems,"
+> *Proceedings of the Royal Society A* **244**, 17 (1958).
+> DOI: [10.1098/rspa.1958.0023](https://doi.org/10.1098/rspa.1958.0023)
+
+The ideal-MHD energy principle: the standard form of the plasma potential energy δW_p that the `EnergyDecomposition` stage evaluates on the (ψ,θ) grid — see [Energy Decomposition](energy_decomposition.md).
+
+---
+
+> I. B. Bernstein, in *Basic Plasma Physics I*, Handbook of Plasma Physics Vol. 1, edited by A. A. Galeev and R. N. Sudan,
+> North-Holland, Amsterdam (1983), p. 421.
+
+Derives the effective-field (Bernstein) form of δW_p, in which the current and pressure terms are completed into |δB_eff|²/μ₀ − K|ξ_n|².
+
+---
+
+> M. S. Chance, Y.-C. Sun, S. C. Jardin, C. E. Kessel, and M. Okabayashi, "MHD stability of tokamak plasmas,"
+> PPPL-CFP-2687, 2nd Symposium on Plasma Dynamics, Trieste (1992).
+> Link: [OSTI 7234275](https://www.osti.gov/biblio/7234275)
+
+Writes both forms of δW_p side by side (Eqs. 12a–12c) and defines the destabilizing kernel K = K₁ + K₂ + K₃ with its shear–current, parallel-current and pressure–curvature pieces (Eqs. 15–17), the split reported as `shear_current`, `parallel_current_squared` and `pressure_curvature`.
+
+---
+
+> A. H. Boozer, "Physics of magnetically confined plasmas,"
+> *Reviews of Modern Physics* **76**, 1071 (2004).
+> DOI: [10.1103/RevModPhys.76.1071](https://doi.org/10.1103/RevModPhys.76.1071)
+
+Sec. V.C derives the energy principle in the form used here, introduces the effective field C = δB + ξ_n(μ₀ j × n̂) (Eq. 98), the constraint (∇×C)·∇ψ = 0 that the perturbed current lie in the perturbed flux surfaces (Eq. 99) and the split w = C²/2μ₀ − w_d (Eq. 100): the two consistency checks of the stage.
+
+---
+
 ## InnerLayer Module
 
 > A. H. Glasser, Z. R. Wang, and J.-K. Park, "Computation of resistive instabilities by matched asymptotic expansions,"

docs/src/energy_decomposition.md

@@ -0,0 +1,147 @@
+# Energy Decomposition of Free-Boundary Eigenmodes
+
+`EnergyDecomposition.jl` evaluates the ideal-MHD plasma potential energy ``\delta W_p`` of a
+free-boundary eigenmode of the ForceFreeStates solve on the ``(\psi,\theta)`` grid and splits
+it into its physical terms. It is the Julia counterpart of the Fortran GPEC energy-decomposition
+diagnostics (`effective_field_flag` and `standard_form_flag`, routines `gpout_recon` and
+`gpout_recon2`; their `docs/tex/recon/main.tex` note carries the full derivation); "Fortran" on
+this page and in the code comments refers to them. This page
+keeps only what is needed to read the output. The stage covers the free-boundary DCON
+eigenmodes only: the edge condition is the free-boundary eigenvector and the reference energy is
+its ``e_{p,k}``, which fixed-boundary or forced-response perturbations do not supply.
+
+## The two forms of ``\delta W_p``
+
+The ideal-MHD energy principle (Bernstein et al. 1958; Boozer 2004, Sec. V.C) writes the fluid
+potential energy of a displacement ``\boldsymbol\xi`` with perturbed field
+``\delta\mathbf B = \nabla\times(\boldsymbol\xi\times\mathbf B)`` in the *standard form*
+
+```math
+2\,\delta W_p = \int d\psi \oint d\theta\, \mathcal J \left[ \frac{|\delta\mathbf B|^2}{\mu_0}
+- \mathbf j\cdot(\delta\mathbf B\times\boldsymbol\xi^*) + (\nabla\cdot\boldsymbol\xi)^*\,(\boldsymbol\xi\cdot\nabla p) \right],
+```
+
+and, after completing the square with the equilibrium current (Bernstein 1983; Chance et al.
+1992, Eq. 12c; Boozer 2004, Eq. 100), in the *effective-field form*
+
+```math
+2\,\delta W_p = \int d\psi \oint d\theta\, \mathcal J \left[ \frac{|\delta\mathbf B_{\rm eff}|^2}{\mu_0}
+- K\,|\xi_n|^2 \right],
+\qquad
+\delta\mathbf B_{\rm eff} = \delta\mathbf B + \xi_n\,(\mu_0 \mathbf j\times\hat n),
+\quad \xi_n = \boldsymbol\xi\cdot\hat n .
+```
+
+Both drop the compressional ``\gamma p|\nabla\cdot\boldsymbol\xi|^2`` term, a reconstruction
+assumption shared with the Fortran. The destabilizing kernel
+``K = 2\,(\mathbf j\times\hat n)\cdot(\mathbf B\cdot\nabla\hat n)`` is split as in Chance et al.
+(1992, Eqs. 15–17) into
+
+```math
+K = K_1 + K_2 + K_3,\qquad
+K_1 = |\nabla\psi|^2\,\sigma S,\quad
+K_2 = \mu_0 B^2 \sigma^2,\quad
+K_3 = 2 p'\,\kappa_\psi,\qquad
+\sigma = \frac{\mathbf j\cdot\mathbf B}{B^2},
+```
+
+with ``S`` the local magnetic shear (Chance et al. 1992, Eq. 17; its surface average is the
+global shear ``q'``), ``\kappa_\psi = \boldsymbol\kappa\cdot\nabla\psi`` the curvature projection
+and ``p'`` the pressure gradient. The three ``K`` terms are the shear–current coupling, the
+parallel-current drive and the pressure–curvature (ballooning) drive.
+
+The effective field is the vector ``\mathbf C`` of Boozer (2004, Eq. 98): its curl is the current
+of the perturbed equilibrium, which has to lie in the perturbed flux surfaces, so
+``\nabla\psi\cdot\nabla\times\delta\mathbf B_{\rm eff} = \nabla\cdot(\delta\mathbf B_{\rm eff}\times\nabla\psi) = 0``
+(Boozer 2004, Eq. 99; Park 2009, §2.2). The decomposition evaluates this identity on every
+equilibrium knot and reports the residual; a non-zero value means the reconstructed field
+components are inconsistent. It is an algebraic identity there, so the residual sits at roundoff
+(1e-13). It validates that ``\delta\mathbf B_{\rm eff}`` was assembled consistently from
+``\delta\mathbf B``, ``\xi_n`` and the metric on each surface; it does not test the radial
+``\xi`` solution, the curvature kernel ``K``, the radial integral or the convergence of the totals.
+
+## Normalization and the reference energy
+
+The eigenmode is built from the free-boundary eigenvector ``\mathbf w_k`` of the
+``(W,N)`` pencil with unit power norm, exactly as `FreeBoundaryStability/eigenmode_plasma_energies`
+is. Every stored energy is scaled by ``2\mu_0/\psi_0^2`` so that `dW_plasma` reads in the same
+normalization, and `dW_plasma_reference` holds the corresponding entry ``e_{p,k}``. The relative
+mismatch of each form against the reference is stored and, above 5 %, reported with a warning.
+On the Solovev test fixture (`mpsi = 16`, `mpert = 8`) both forms reproduce ``e_{p,k}`` to ≲2 %
+(effective / standard form +1.7 % / +1.4 % for the least stable mode, −2.0 % / −2.0 % for the
+second), converging with resolution (+5.0e-3 / +1.7e-4 and +5.6e-4 / −7.8e-5 at `mpsi = 128`).
+The effective-field form needs the metric-consistent
+``B^2 = (\chi'/\mathcal J)^2(g_{\theta\theta} + 2q g_{\theta\zeta} + q^2 g_{\zeta\zeta})`` for
+``\sigma``, ``K_2`` and ``\kappa_\psi`` rather than the equilibrium's separate ``|B|`` spline.
+The standard-form remainder is pure ``\psi``-quadrature error of the DCON Lagrangian (converging
+roughly as ``h^2``–``h^3``). The remaining offset between the two forms is equilibrium noise from
+the fixture's loose `etol = 1e-7`: with `etol = 1e-10` it drops 50–300× to ≲1e-4 for both modes at
+`mpsi` = 64 and 128. The pressure-curvature term ``K_3`` (through ``\partial_\psi B^2``) carries
+most of it for the least stable mode but not all of it for the second, and at `etol = 1e-7` every
+term is rough near the edge.
+
+## Running it
+
+```toml
+[EnergyDecomposition]
+eigenmodes = [1]              # Free-boundary eigenmodes to decompose (1 = least stable)
+effective_field_form = true   # Decompose δW_p in the effective-field form ½∫[|δB_eff|²/μ₀ − K|ξ_n|²]
+standard_form = true          # Decompose δW_p in the standard form ½∫[|δB|²/μ₀ − j·(δB×ξ*) + (∇·ξ)*(ξ·∇p)]
+write_densities = false       # Also store the per-(ψ,θ) energy densities with R,Z for 2-D maps
+```
+
+The stage runs right after ForceFreeStates, also with `force_termination = true`, and needs the
+dense ``\xi`` solution (`integrator = "forward"`), the free-boundary eigenmodes
+(`vac_flag = true`), a single toroidal mode number and an ideal solve; otherwise it warns and is
+skipped. From a script:
+
+```julia
+result = GeneralizedPerturbedEquilibrium.main(["examples/Solovev_ideal_example"])
+energy = GeneralizedPerturbedEquilibrium.energy_decomposition(result.ffs; eigenmodes=[1, 2])
+```
+
+The radial nodes are the equilibrium ``\psi`` knots inside the solution domain plus the solution
+edge: on the knots the metric, the Euler–Lagrange matrices and the geometry are all exact, which
+is what keeps the curl identity at roundoff. ``\xi`` and ``\xi'`` at each node come from
+integrating the Euler–Lagrange system from the nearest saved ODE step that no singular surface
+separates from it, so they do not depend on how densely the forward sweep saved its solution, and
+``\xi_s`` is recomputed from the matrices at each node. The Fortran evaluates
+its profiles on the same knots but integrates its totals over the ODE `psifac` grid.
+
+## Output
+
+Everything lands in `PerturbedEquilibrium/EnergyDecomposition/` of `gpec.h5`, profiles on the
+decomposition ``\psi`` nodes with one column per decomposed eigenmode; `eigenmode_index` is the
+coordinate of that axis (the free-boundary eigenmode number, 1 = least stable):
+
+| Dataset | Content |
+|---|---|
+| `effective_b_squared` (+ `_psi`, `_theta`, `_zeta`) | ``\oint \mathcal J\,\lvert\delta\mathbf B_{\rm eff}\rvert^2/\mu_0\,d\theta`` and its component split |
+| `shear_current`, `parallel_current_squared`, `pressure_curvature` | ``\oint \mathcal J\,K_i\lvert\xi_n\rvert^2\,d\theta`` for ``K_1, K_2, K_3`` |
+| `b_squared`, `current_coupling`, `pressure_compression` | the three standard-form pieces |
+| `dW_cumulative`, `dW_cumulative_standard_form` | running ``\tfrac12\int d\psi`` of each form |
+| `dW_plasma`, `dW_plasma_standard_form`, `dW_plasma_reference` | the totals and ``e_{p,k}`` |
+| `dW_plasma_relative_error`, `dW_plasma_standard_form_relative_error` | mismatch against the reference |
+| `effective_b_curl_residual_max`, `_rms`, `_relative` | the ``\nabla\psi\cdot\nabla\times\delta\mathbf B_{\rm eff} = 0`` check |
+| `Densities/…` | per-``(\psi,\theta)`` integrands with `R`, `Z`, only with `write_densities` |
+
+A disabled form leaves its datasets zero-extent. The dataset names and the API docstrings below keep the
+code's shorthand `b` for the perturbed field ``\delta\mathbf B``.
+
+## References
+
+- I. B. Bernstein, E. A. Frieman, M. D. Kruskal, R. M. Kulsrud, "An energy principle for hydromagnetic stability problems", *Proc. R. Soc. London A* **244**, 17 (1958), [doi:10.1098/rspa.1958.0023](https://doi.org/10.1098/rspa.1958.0023) — the energy principle.
+- I. B. Bernstein, in *Basic Plasma Physics I*, Handbook of Plasma Physics Vol. 1, edited by A. A. Galeev and R. N. Sudan (North-Holland, Amsterdam, 1983), p. 421 — the derivation of the effective-field form.
+- M. S. Chance, Y.-C. Sun, S. C. Jardin, C. E. Kessel, M. Okabayashi, "MHD stability of tokamak plasmas", PPPL-CFP-2687, 2nd Symposium on Plasma Dynamics, Trieste (1992), [OSTI 7234275](https://www.osti.gov/biblio/7234275) — Eqs. (12a)–(12c) for the forms of ``\delta W_p``, Eqs. (15)–(17) for ``K = K_1 + K_2 + K_3`` and the local shear ``S``.
+- A. H. Boozer, "Physics of magnetically confined plasmas", *Rev. Mod. Phys.* **76**, 1071 (2004), Sec. V.C, [doi:10.1103/RevModPhys.76.1071](https://doi.org/10.1103/RevModPhys.76.1071) — the energy principle, the effective field ``\mathbf C`` (Eq. 98), the ``(\nabla\times\mathbf C)\cdot\nabla\psi = 0`` constraint (Eq. 99) and the split ``w = C^2/2\mu_0 - w_d`` (Eq. 100).
+- J.-K. Park, *Ideal Perturbed Equilibria in Tokamaks*, PhD thesis, Princeton University (2009), §2.2 — the effective field as the surface-current field.
+- A. H. Glasser, "The direct criterion of Newcomb for the ideal MHD stability of an axisymmetric toroidal plasma", *Phys. Plasmas* **23**, 072505 (2016), [doi:10.1063/1.4958328](https://doi.org/10.1063/1.4958328) — the DCON coordinates, matrices and shear used here.
+
+## API
+
+```@docs
+GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.EnergyDecompositionControl
+GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.EnergyDecompositionResult
+GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.decompose_energy
+GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.write_energy_decomposition!
+```

examples/DIIID-like_ideal_example/gpec.toml

@@ -54,6 +54,12 @@ integrator            = "forward" # Integration formalism: "forward" (dense ξ,
 set_psilim_via_dmlim  = true   # Truncate at (last_rational_q + dmlim)/n — TRUE for diverted equilibria (q → ∞ at separatrix)
 dmlim                 = 0.2    # Truncate integration at (last_rational_q + dmlim)/n (used when set_psilim_via_dmlim = true)
 
+[EnergyDecomposition]
+eigenmodes = [1]              # Free-boundary eigenmodes to decompose (1 = least stable)
+effective_field_form = true   # Decompose δW_p in the effective-field form ½∫[|b_eff|²/μ₀ − K|ξ_n|²]
+standard_form = true          # Decompose δW_p in the standard form ½∫[|b|²/μ₀ − j·(b×ξ*) + (∇·ξ)*(ξ·∇p)]
+write_densities = false       # Also store the per-(ψ,θ) energy densities with R,Z for 2-D maps
+
 [ForcingTerms]
 forcing_data_format = "coil"            # Format: "ascii", "hdf5", or "coil" (Biot-Savart from 3D wires)
 machine = "d3d"                         # Geometry prefix; resolves to bundled coil_geometries/d3d_*.dat

examples/Solovev_ideal_example/gpec.toml

@@ -75,6 +75,13 @@ save_interval = 1             # Save every Nth ODE step (1=all). Always saves ne
 integrator            = "forward" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic), "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix), or "galerkin" (RDCON outer-region singular Galerkin Δ′, with the optional RPEC inner-layer match)
 set_psilim_via_dmlim  = false  # FALSE for limited/analytical equilibria — rationals sparse, dmlim would chop too much edge
 dmlim                 = 0.2    # Truncate integration at (last_rational_q + dmlim)/n (used when set_psilim_via_dmlim = true)
+
+[EnergyDecomposition]
+eigenmodes = [1]              # Free-boundary eigenmodes to decompose (1 = least stable)
+effective_field_form = true   # Decompose δW_p in the effective-field form ½∫[|b_eff|²/μ₀ − K|ξ_n|²]
+standard_form = true          # Decompose δW_p in the standard form ½∫[|b|²/μ₀ − j·(b×ξ*) + (∇·ξ)*(ξ·∇p)]
+write_densities = false       # Also store the per-(ψ,θ) energy densities with R,Z for 2-D maps
+
 # Solovev analytic equilibrium parameters (eq_type = "sol"); see SolovevConfig in src/Equilibrium.
 [SOL_INPUT]
 mr = 128       # Number of radial grid zones