diff --git a/docs/development/architecture.md b/docs/development/architecture.md index 37aad79b3..402d0152f 100644 --- a/docs/development/architecture.md +++ b/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 diff --git a/docs/development/hdf5-conventions.md b/docs/development/hdf5-conventions.md index c4631c984..52725701f 100644 --- a/docs/development/hdf5-conventions.md +++ b/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/` | `/` (torque/energy profiles, `EnergyIntegrals/`, `KineticMatrices/`); multi-ion runs add `PerSpecies///` 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_/` | diff --git a/docs/development/references.md b/docs/development/references.md index 43f3daf74..d97bfdfb0 100644 --- a/docs/development/references.md +++ b/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: diff --git a/docs/make.jl b/docs/make.jl index 6fb5a3ed5..44daab3e6 100644 --- a/docs/make.jl +++ b/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", diff --git a/docs/src/api.md b/docs/src/api.md index b2406c53d..6033d3c58 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -97,6 +97,7 @@ remain TOML-driven. GeneralizedPerturbedEquilibrium.EulerLagrangeProblem GeneralizedPerturbedEquilibrium.solve GeneralizedPerturbedEquilibrium.perturbed_equilibrium +GeneralizedPerturbedEquilibrium.energy_decomposition ``` ## Integrator selectors diff --git a/docs/src/citations.md b/docs/src/citations.md index dd8ef7546..ebffffb3f 100644 --- a/docs/src/citations.md +++ b/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," diff --git a/docs/src/energy_decomposition.md b/docs/src/energy_decomposition.md new file mode 100644 index 000000000..97d52f0c5 --- /dev/null +++ b/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! +``` diff --git a/examples/DIIID-like_ideal_example/gpec.toml b/examples/DIIID-like_ideal_example/gpec.toml index 2aca73d8a..e74bbd4c4 100644 --- a/examples/DIIID-like_ideal_example/gpec.toml +++ b/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 diff --git a/examples/Solovev_ideal_example/gpec.toml b/examples/Solovev_ideal_example/gpec.toml index f74d26fa6..d3d6b4e9a 100644 --- a/examples/Solovev_ideal_example/gpec.toml +++ b/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 diff --git a/regression-harness/cases/diiid_n1.toml b/regression-harness/cases/diiid_n1.toml index 1e033e5ca..2f9315f9d 100644 --- a/regression-harness/cases/diiid_n1.toml +++ b/regression-harness/cases/diiid_n1.toml @@ -422,6 +422,49 @@ label = "ballooning Delta' profile (checksum)" noise_threshold = 0 order = 74 +# Energy decomposition of the least-stable free-boundary eigenmode — the two δW_p route totals +# (power-normalized like eigenmode_plasma_energies), their mismatch against ep[1], and the +# (∇×b_eff)·∇ψ consistency residual. +[quantities.energy_dW_plasma] +h5path = "PerturbedEquilibrium/EnergyDecomposition/dW_plasma" +type = "real_vector" +extract = "all_real" +label = "energy decomposition dW_p (effective field)" +noise_threshold = 1e-10 +order = 75 + +[quantities.energy_dW_plasma_standard_form] +h5path = "PerturbedEquilibrium/EnergyDecomposition/dW_plasma_standard_form" +type = "complex_vector" +extract = "all_complex" +label = "energy decomposition dW_p (standard form)" +noise_threshold = 1e-10 +order = 76 + +[quantities.energy_dW_relative_error] +h5path = "PerturbedEquilibrium/EnergyDecomposition/dW_plasma_relative_error" +type = "real_vector" +extract = "all_real" +label = "energy decomposition rel. error vs ep" +noise_threshold = 1e-10 +order = 77 + +[quantities.energy_dW_standard_form_relative_error] +h5path = "PerturbedEquilibrium/EnergyDecomposition/dW_plasma_standard_form_relative_error" +type = "real_vector" +extract = "all_real" +label = "energy decomposition std-form rel. error vs ep" +noise_threshold = 1e-10 +order = 78 + +[quantities.energy_curl_residual] +h5path = "PerturbedEquilibrium/EnergyDecomposition/effective_b_curl_residual_relative" +type = "real_vector" +extract = "all_real" +label = "energy decomposition curl residual (relative)" +noise_threshold = 1e-9 +order = 79 + # Error fields — per-coil-set overlap with the full-window dominant mode and its rigid shift/tilt # sensitivities. Singular vectors carry an arbitrary phase, so only magnitudes are tracked; the # in-plane components of a whole n=1-phased array vanish by symmetry, so the norms track the diff --git a/regression-harness/cases/solovev_n1.toml b/regression-harness/cases/solovev_n1.toml index d7f4224e5..80a2d2a1f 100644 --- a/regression-harness/cases/solovev_n1.toml +++ b/regression-harness/cases/solovev_n1.toml @@ -192,6 +192,49 @@ label = "ca_left (checksum)" noise_threshold = 0 order = 72 +# Energy decomposition of the least-stable free-boundary eigenmode — the two δW_p route totals +# (power-normalized like eigenmode_plasma_energies), their mismatch against ep[1], and the +# (∇×b_eff)·∇ψ consistency residual. +[quantities.energy_dW_plasma] +h5path = "PerturbedEquilibrium/EnergyDecomposition/dW_plasma" +type = "real_vector" +extract = "all_real" +label = "energy decomposition dW_p (effective field)" +noise_threshold = 1e-10 +order = 80 + +[quantities.energy_dW_plasma_standard_form] +h5path = "PerturbedEquilibrium/EnergyDecomposition/dW_plasma_standard_form" +type = "complex_vector" +extract = "all_complex" +label = "energy decomposition dW_p (standard form)" +noise_threshold = 1e-10 +order = 81 + +[quantities.energy_dW_relative_error] +h5path = "PerturbedEquilibrium/EnergyDecomposition/dW_plasma_relative_error" +type = "real_vector" +extract = "all_real" +label = "energy decomposition rel. error vs ep" +noise_threshold = 1e-10 +order = 82 + +[quantities.energy_dW_standard_form_relative_error] +h5path = "PerturbedEquilibrium/EnergyDecomposition/dW_plasma_standard_form_relative_error" +type = "real_vector" +extract = "all_real" +label = "energy decomposition std-form rel. error vs ep" +noise_threshold = 1e-10 +order = 83 + +[quantities.energy_curl_residual] +h5path = "PerturbedEquilibrium/EnergyDecomposition/effective_b_curl_residual_relative" +type = "real_vector" +extract = "all_real" +label = "energy decomposition curl residual (relative)" +noise_threshold = 1e-9 +order = 84 + # Runtime (special: not from H5) [quantities.runtime] h5path = "" diff --git a/src/Equilibrium/Equilibrium.jl b/src/Equilibrium/Equilibrium.jl index 25fdcc46d..d8f0eb0ba 100644 --- a/src/Equilibrium/Equilibrium.jl +++ b/src/Equilibrium/Equilibrium.jl @@ -32,7 +32,7 @@ export setup_equilibrium, EquilibriumConfig, PlasmaEquilibrium, EquilibriumParam ProfileSplines, GeometryProfileSplines, compute_geometry_profiles, KineticProfileSplines, load_kinetic_profiles, shift_exb_rotation, KineticProfileData, read_kinetic_file, write_kinetic_h5 -export flux_surface_metric, flux_surface_area +export flux_surface_metric, flux_surface_area, flux_surface_basis, metric_bsq_interpolant export wants_two_pass, refined_psi_grid, merge_mandatory_nodes, bracket_mandatory_nodes, enforce_min_spacing, implied_knot_count export compute_sqrt_jac_delpsi, compute_sqrtamat, rootarea_to_area_weight, area_to_rootarea_weight diff --git a/src/Equilibrium/FluxSurfaceMetrics.jl b/src/Equilibrium/FluxSurfaceMetrics.jl index cce67db2d..f565d47f3 100644 --- a/src/Equilibrium/FluxSurfaceMetrics.jl +++ b/src/Equilibrium/FluxSurfaceMetrics.jl @@ -17,14 +17,101 @@ function flux_surface_metric(equil::PlasmaEquilibrium, psi::Float64, theta::Floa eta = 2π * (theta + deta) r = equil.ro + rfac * cos(eta) - # Metric components for |∇ψ| in flux coordinates (GPEC gpout.f) - w11 = (1.0 + deta_y) * (2π)^2 * rfac * r / jac - w12 = -r2_y * π * r / (rfac * jac) + w11, w12 = _grad_psi_poloidal(rfac, r, jac, r2_y, deta_y) delpsi = sqrt(w11^2 + w12^2) return (; r, jac, delpsi) end +# Poloidal-plane components (w11, w12) of ∇ψ in the rzphi straight-field-line coordinates (GPEC gpout.f). +@inline _grad_psi_poloidal(rfac, r, jac, r2_y, deta_y) = ((1.0 + deta_y) * (2π)^2 * rfac * r / jac, -r2_y * π * r / (rfac * jac)) + +""" + flux_surface_basis(equil, psi, theta; hint=(Ref(1), Ref(1))) -> NamedTuple + flux_surface_basis(equil, ipsi, itheta) -> NamedTuple + +Full straight-field-line geometry at one point: `r`, `z`, Jacobian `jac` and its ψ-derivative +`jac_psi`, `delpsi2` = |∇ψ|², `dpdt` = ∇ψ·∇θ, `dpdz` = ∇ψ·∇ζ, and the covariant metric `g11`, +`g12`, `g13`, `g22`, `g23`, `g33` of the tangent basis e_i = ∂x/∂q^i. The gradient and tangent bases +are the `w` and `v` matrices of the Fortran `gpeq_c`; |∇ψ| and `jac` agree with +`flux_surface_metric`. Derivatives are with respect to ψ_N. The `(ipsi, itheta)` method reads the +values and first derivatives stored at equilibrium grid point `(rzphi_xs[ipsi], rzphi_ys[itheta])`, +as `make_metric` does. +""" +function flux_surface_basis(equil::PlasmaEquilibrium, psi::Float64, theta::Float64; hint=(Ref(1), Ref(1))) + pt = (psi, theta) + dx = DerivOp(1, 0) + dy = DerivOp(0, 1) + return _flux_surface_basis(equil, theta, + equil.rzphi_rsquared(pt; hint=hint), equil.rzphi_rsquared(pt; deriv=dx, hint=hint), equil.rzphi_rsquared(pt; deriv=dy, hint=hint), + equil.rzphi_offset(pt; hint=hint), equil.rzphi_offset(pt; deriv=dx, hint=hint), equil.rzphi_offset(pt; deriv=dy, hint=hint), + equil.rzphi_nu(pt; deriv=dx, hint=hint), equil.rzphi_nu(pt; deriv=dy, hint=hint), + equil.rzphi_jac(pt; hint=hint), equil.rzphi_jac(pt; deriv=dx, hint=hint)) +end + +function flux_surface_basis(equil::PlasmaEquilibrium, ipsi::Int, itheta::Int) + r2 = equil.rzphi_rsquared.nodal_derivs.partials + deta = equil.rzphi_offset.nodal_derivs.partials + nu = equil.rzphi_nu.nodal_derivs.partials + jac = equil.rzphi_jac.nodal_derivs.partials + return _flux_surface_basis(equil, equil.rzphi_ys[itheta], + r2[1, ipsi, itheta], r2[2, ipsi, itheta], r2[3, ipsi, itheta], + deta[1, ipsi, itheta], deta[2, ipsi, itheta], deta[3, ipsi, itheta], + nu[2, ipsi, itheta], nu[3, ipsi, itheta], + jac[1, ipsi, itheta], jac[2, ipsi, itheta]) +end + +function _flux_surface_basis(equil::PlasmaEquilibrium, theta, r2, r2_x, r2_y, deta, deta_x, deta_y, nu_x, nu_y, jac, jac_psi) + rfac = sqrt(abs(r2)) + eta = 2π * (theta + deta) + r = equil.ro + rfac * cos(eta) + z = equil.zo + rfac * sin(eta) + + # Gradient basis in the poloidal plane: rows ∇ψ, ∇θ, ∇ζ (Fortran gpeq_c w-matrix). + w11, w12 = _grad_psi_poloidal(rfac, r, jac, r2_y, deta_y) + w21 = -(2π)^2 * rfac * r * deta_x / jac + w22 = π * r * r2_x / (rfac * jac) + w31 = (2π * r * rfac / jac) * (deta_x * nu_y - nu_x * (1.0 + deta_y)) + w32 = (r / (2.0 * rfac * jac)) * (nu_x * r2_y - r2_x * nu_y) + + # Tangent basis e_ψ, e_θ, e_ζ without the Jacobian factor (Fortran gpeq_c v-matrix). + v11 = r2_x / (2.0 * rfac) + v12 = deta_x * 2π * rfac + v13 = nu_x * r + v21 = r2_y / (2.0 * rfac) + v22 = (1.0 + deta_y) * 2π * rfac + v23 = nu_y * r + v33 = 2π * r + + return (; r, z, jac, jac_psi, delpsi2=w11^2 + w12^2, dpdt=w11 * w21 + w12 * w22, dpdz=w11 * w31 + w12 * w32, + g11=v11^2 + v12^2 + v13^2, g12=v11 * v21 + v12 * v22 + v13 * v23, g13=v13 * v33, + g22=v21^2 + v22^2 + v23^2, g23=v23 * v33, g33=v33^2) +end + +""" + metric_bsq_interpolant(equil) -> CubicInterpolantND + +Bicubic interpolant of B² = (χ'/J)² (g_θθ + 2q g_θζ + q² g_ζζ) on the equilibrium (ψ,θ) grid, the +field strength of the straight-field-line metric that the Euler-Lagrange matrices and the perturbed +field are built from. The equilibrium's `eqfun_B` (from ψ₀|∇ψ| and F) agrees with it only to the +accuracy of the coordinate transformation. +""" +function metric_bsq_interpolant(equil::PlasmaEquilibrium) + xs = equil.rzphi_xs + ys = equil.rzphi_ys + chi1 = 2π * equil.psio + vals = zeros(length(xs), length(ys)) + for i in eachindex(xs) + q = equil.profiles.q_spline(xs[i]) + for k in 1:(length(ys)-1) + b = flux_surface_basis(equil, i, k) + vals[i, k] = (chi1 / b.jac)^2 * (b.g22 + 2.0 * q * b.g23 + q^2 * b.g33) + end + vals[i, end] = vals[i, 1] + end + return cubic_interp((xs, ys), vals; bc=(CubicFit(), PeriodicBC()), extrap=(ExtendExtrap(), WrapExtrap())) +end + """ flux_surface_area(equil, psi, mtheta) -> Float64 diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl index 2224c53e2..c5ed03595 100644 --- a/src/ForceFreeStates/EulerLagrange.jl +++ b/src/ForceFreeStates/EulerLagrange.jl @@ -1438,20 +1438,20 @@ end Evaluate Ξ_s = -A⁻¹(B·Ξ′_ψ + C·Ξ_ψ) [Glasser Phys. Plasmas 2016 112506 eq. 18] at `psieval`, writing into `xi_s`. `du1` and `u1` are the Ξ′_ψ and Ξ_ψ blocks at the same ψ, i.e. slices of a -`el_derivatives!` result and its input state. +`el_derivatives!` result and its input state, or a single displacement vector and its derivative. Split out of the derivative kernel because Ξ_s is needed only at saved nodes, not at every Runge-Kutta stage. Ideal runs factor the Hermitian A by Cholesky; with `kinetic=true` A picks up non-Hermitian contributions and needs an LU. """ -@with_pool pool function compute_node_xi_s!(xi_s::AbstractMatrix{ComplexF64}, du1::AbstractMatrix{ComplexF64}, - u1::AbstractMatrix{ComplexF64}, mats::MatrixSplines, psieval::Float64; kinetic::Bool=false, hint::Base.RefValue{Int}=Ref(1)) +@with_pool pool function compute_node_xi_s!(xi_s::AbstractVecOrMat{ComplexF64}, du1::AbstractVecOrMat{ComplexF64}, + u1::AbstractVecOrMat{ComplexF64}, mats::MatrixSplines, psieval::Float64; kinetic::Bool=false, hint::Base.RefValue{Int}=Ref(1)) Npert = size(u1, 1) amat = acquire!(pool, ComplexF64, Npert, Npert) bmat = similar!(pool, amat) cmat = similar!(pool, amat) - tmp_mat = similar!(pool, amat) + tmp_mat = acquire!(pool, ComplexF64, size(u1)...) # A/B/C of the active model: the kinetic A is non-Hermitian, hence the factorization split below active_mats = kinetic ? mats.kinetic : mats.ideal diff --git a/src/ForceFreeStates/Result.jl b/src/ForceFreeStates/Result.jl index da8fe3648..0ed7d5f42 100644 --- a/src/ForceFreeStates/Result.jl +++ b/src/ForceFreeStates/Result.jl @@ -254,3 +254,29 @@ function build_result( solution, odet, wp, free_energies, delta_prime, gal_data ) end + +""" + solution_state_at(ffs, psi; rtol=ffs.control.eulerlagrange_tolerance) -> Array{ComplexF64,3} + +State `(u₁, u₂)` of a forward run's stored Euler-Lagrange solution at `psi`, in the basis of +`ffs.solution.u_store`. A saved step is returned as stored; otherwise the system is integrated to +`psi` from the nearest saved step that no singular surface separates from it. The stored basis is a +constant right-multiple of the integrated one and the system is linear in `u`, so this is the state +the forward sweep would have saved at `psi`, without interpolating between saved steps. +""" +function solution_state_at(ffs::ForceFreeStatesResult, psi::Float64; rtol::Float64=ffs.control.eulerlagrange_tolerance) + sol = ffs.solution + psis = @view sol.psi_store[1:sol.step] + i = searchsortedlast(psis, psi) + (i >= 1 && psis[i] == psi) && return sol.u_store[:, :, :, i] + sing = [s.psifac for s in ffs.surfaces] + reachable(j) = 1 <= j <= sol.step && !any(s -> min(psis[j], psi) < s < max(psis[j], psi), sing) + from = filter(reachable, [i, i + 1]) + isempty(from) && throw(ArgumentError("ψ = $psi is cut off from both neighbouring saved steps by a singular surface")) + j = argmin(j -> abs(psis[j] - psi), from) + kinetic = ffs.mats.kinetic !== nothing + hints = (Ref(1), Ref(1)) + rhs!(du, u, _, t) = (el_derivatives!(du, u, kinetic, ffs.equil, ffs.mats, ffs, t, hints...); nothing) + prob = ODEProblem(rhs!, sol.u_store[:, :, :, j], (psis[j], psi)) + return solve_el(prob, rtol).u[end] +end diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 6322d2c79..c3c20f480 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -276,6 +276,8 @@ function main_from_inputs( push!(runtimes, "force_free_states" => ffs_dt) @info "Force-Free States completed in $(@sprintf("%.3f", ffs_dt)) s" + run_energy_decomposition(ffs_result, inputs) + # Early exit if user only requested force-free states (SLAYER still runs). if ctrl.force_termination slayer_result = run_slayer_stage(ffs_result, inputs, nothing; runtimes=runtimes) @@ -973,6 +975,65 @@ function run_perturbed_equilibrium( return pe_state end +""" + run_energy_decomposition(result, inputs) -> EnergyDecompositionResult or nothing + +Decompose the plasma energy of the free-boundary eigenmodes named in the `[EnergyDecomposition]` +section against a published force-free-states `result` and write the group into the run's +`gpec.h5`. Returns `nothing` when the deck has no such section, or when the solve cannot supply +what the decomposition needs: the forward integrator's dense ξ solution, the free-boundary +eigenmodes, a single toroidal mode number and ideal matrices (warned and skipped, like every +optional stage). +""" +function run_energy_decomposition(result::ForceFreeStatesResult, inputs::Dict{String,Any}) + "EnergyDecomposition" in keys(inputs) || return nothing + @info "\n Energy Decomposition\n$_SECTION" + ed_start = time() + ed_ctrl = PerturbedEquilibrium.EnergyDecompositionControl(; (Symbol(k) => v for (k, v) in inputs["EnergyDecomposition"])...) + ForceFreeStates.require_solution(result, "energy decomposition") || return nothing + ForceFreeStates.require(result, :free_boundary, "energy decomposition") || return nothing + reason = if result.integrator !== :forward + "it needs the forward integrator's ξ solution" + elseif result.nlow != result.nhigh + "it handles a single toroidal mode number (nn_low == nn_high)" + elseif result.mats.kinetic !== nothing + "it decomposes ideal eigenmodes only (kinetic_factor = 0)" + else + "" + end + if !isempty(reason) + @warn "Skipping energy decomposition: $reason" + return nothing + end + energy = energy_decomposition(result; eigenmodes=ed_ctrl.eigenmodes, effective_field_form=ed_ctrl.effective_field_form, + standard_form=ed_ctrl.standard_form, write_densities=ed_ctrl.write_densities) + @info "Energy Decomposition completed in $(@sprintf("%.3f", time() - ed_start)) s" + return energy +end + +""" + energy_decomposition(ffs; kwargs...) -> EnergyDecompositionResult + +Decompose the plasma energy of free-boundary eigenmodes of a forward solve `ffs` and write the +result to the run's `gpec.h5` when the solve wrote its outputs. Keyword arguments are the +`PerturbedEquilibrium.EnergyDecompositionControl` fields plus `verbose`. + +```julia +energy = energy_decomposition(ffs; eigenmodes=[1, 2]) +``` +""" +function energy_decomposition(ffs::ForceFreeStatesResult; kwargs...) + ctrl = ffs.control + energy = PerturbedEquilibrium.decompose_energy(ffs; kwargs...) + if ctrl.write_outputs_to_HDF5 + h5open(joinpath(ffs.dir_path, ctrl.HDF5_filename), "cw") do h5 + PerturbedEquilibrium.write_energy_decomposition!(h5, energy) + end + @info "Energy decomposition written to $(ctrl.HDF5_filename)" + end + return energy +end + """ perturbed_equilibrium(ffs, rmp; forcing_modes=nothing, coil_sets=nothing, runtimes=nothing, kwargs...) -> PerturbedEquilibriumState @@ -1886,7 +1947,7 @@ function write_imas(dd, result) end export main, write_imas -export solve, perturbed_equilibrium +export solve, perturbed_equilibrium, energy_decomposition export PlasmaEquilibrium, EulerLagrangeProblem, Forward, Riccati, Galerkin, ResistiveMatch, ForceFreeStatesResult, RMPField include("Precompile.jl") diff --git a/src/PerturbedEquilibrium/EnergyDecomposition.jl b/src/PerturbedEquilibrium/EnergyDecomposition.jl new file mode 100644 index 000000000..a4d764ad6 --- /dev/null +++ b/src/PerturbedEquilibrium/EnergyDecomposition.jl @@ -0,0 +1,946 @@ +""" +Energy decomposition of free-boundary eigenmodes on the (ψ,θ) grid. + +Evaluates the ideal-MHD plasma potential energy δW_p of a ForceFreeStates eigenmode in two +forms and splits it into its physical terms (Fortran GPEC `effective_field_flag`/`standard_form_flag`, +routines `gpout_recon`/`gpout_recon2`, theory note `docs/tex/recon/main.tex` there, Bernstein et al. 1958): + + - effective-field form 2δW_p = ∫dψ ∮dθ J [ |b_eff|²/μ₀ − K|ξ_n|² ], b_eff = b + ξ_n (μ₀ j × n̂) + - standard form 2δW_p = ∫dψ ∮dθ J [ |b|²/μ₀ − j·(b×ξ*) + (∇·ξ)*(ξ·∇p) ] + +both without the compressional γp|∇·ξ|² term (a reconstruction assumption shared with the +Fortran). K = K1 + K2 + K3 with K1 = |∇ψ|²σS (shear-current), K2 = μ₀B²σ² (parallel current), +K3 = 2p'κ_ψ (pressure-curvature) and σ = j·B/B². Every stored energy is scaled by 2μ₀/ψ₀² so the +totals read in the normalization of `eigenmode_plasma_energies`. + +Array conventions follow the gpeq arrays: `Jxi_psi` = Jξ^ψ, `Jb_*` = Jb^i (Jacobian-weighted +contravariant), `b_cov_*` = b_i (covariant), `xi_n` = ξ·n̂. +""" + +using ..Utilities.PhysicalConstants: MU_0 + +""" + EnergyDecompositionControl + +User-facing controls of the `[EnergyDecomposition]` section of gpec.toml. + +## Fields + + - `eigenmodes::Vector{Int}` - Free-boundary eigenmodes to decompose (1 = least stable) + - `effective_field_form::Bool` - Decompose δW_p in the effective-field form ½∫[|b_eff|²/μ₀ − K|ξ_n|²] + - `standard_form::Bool` - Decompose δW_p in the standard form ½∫[|b|²/μ₀ − j·(b×ξ*) + (∇·ξ)*(ξ·∇p)] + - `write_densities::Bool` - Also store the per-(ψ,θ) energy densities with R,Z for 2-D maps +""" +@kwdef struct EnergyDecompositionControl + eigenmodes::Vector{Int} = [1] + effective_field_form::Bool = true + standard_form::Bool = true + write_densities::Bool = false +end + +""" + EnergyDecompositionResult + +Energy decomposition of the requested free-boundary eigenmodes. Field names are the dataset +names under `PerturbedEquilibrium/EnergyDecomposition/`; profiles are `(psi × eigenmode)` matrices, +totals `(eigenmode)` vectors, densities `(psi × theta × eigenmode)` arrays. Energies are scaled +by 2μ₀/ψ₀² to the normalization of `eigenmode_plasma_energies`. Fields of a disabled form or of +unrequested densities stay empty. + +## Fields + + - `psi` - ψ_N of the ξ solution grid + - `eigenmode_index` - free-boundary eigenmode index of each decomposed mode + - `effective_b_squared`, `effective_b_squared_psi/theta/zeta` - ∮J|b_eff|²/μ₀ dθ and its component split + - `shear_current`, `parallel_current_squared`, `pressure_curvature` - ∮J K_i|ξ_n|² dθ for K1, K2, K3 + - `b_squared`, `current_coupling`, `pressure_compression` - standard-form pieces ∮J|b|²/μ₀, ∮J j·(b×ξ*), ∮J(∇·ξ)*(ξ·∇p) dθ + - `dW_cumulative`, `dW_cumulative_standard_form` - running ½∫dψ of each form + - `dW_plasma`, `dW_plasma_standard_form` - the two totals + - `dW_plasma_reference` - `ep[k]` of the same eigenmode from FreeBoundaryStability + - `dW_plasma_relative_error`, `dW_plasma_standard_form_relative_error` - (total − reference)/|reference| + - `effective_b_curl_residual_max/rms/relative` - the (∇×b_eff)·∇ψ = 0 check + - `theta`, `R`, `Z` - the θ grid and cylindrical coordinates of the densities + - `*_density` - per-point integrands of the profiles above (integral = Σ density·J/mtheta) +""" +@kwdef mutable struct EnergyDecompositionResult + psi::Vector{Float64} = Float64[] + eigenmode_index::Vector{Int} = Int[] + effective_b_squared::Matrix{Float64} = zeros(0, 0) + effective_b_squared_psi::Matrix{Float64} = zeros(0, 0) + effective_b_squared_theta::Matrix{Float64} = zeros(0, 0) + effective_b_squared_zeta::Matrix{Float64} = zeros(0, 0) + shear_current::Matrix{Float64} = zeros(0, 0) + parallel_current_squared::Matrix{Float64} = zeros(0, 0) + pressure_curvature::Matrix{Float64} = zeros(0, 0) + b_squared::Matrix{Float64} = zeros(0, 0) + current_coupling::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) + pressure_compression::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) + dW_cumulative::Matrix{Float64} = zeros(0, 0) + dW_cumulative_standard_form::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) + dW_plasma::Vector{Float64} = Float64[] + dW_plasma_standard_form::Vector{ComplexF64} = ComplexF64[] + dW_plasma_reference::Vector{Float64} = Float64[] + dW_plasma_relative_error::Vector{Float64} = Float64[] + dW_plasma_standard_form_relative_error::Vector{Float64} = Float64[] + effective_b_curl_residual_max::Vector{Float64} = Float64[] + effective_b_curl_residual_rms::Vector{Float64} = Float64[] + effective_b_curl_residual_relative::Vector{Float64} = Float64[] + theta::Vector{Float64} = Float64[] + R::Matrix{Float64} = zeros(0, 0) + Z::Matrix{Float64} = zeros(0, 0) + effective_b_squared_density::Array{Float64,3} = zeros(0, 0, 0) + effective_b_squared_psi_density::Array{Float64,3} = zeros(0, 0, 0) + effective_b_squared_theta_density::Array{Float64,3} = zeros(0, 0, 0) + effective_b_squared_zeta_density::Array{Float64,3} = zeros(0, 0, 0) + shear_current_density::Array{Float64,3} = zeros(0, 0, 0) + parallel_current_squared_density::Array{Float64,3} = zeros(0, 0, 0) + pressure_curvature_density::Array{Float64,3} = zeros(0, 0, 0) + effective_b_curl_residual_density::Array{ComplexF64,3} = zeros(ComplexF64, 0, 0, 0) + b_squared_density::Array{Float64,3} = zeros(0, 0, 0) + current_coupling_density::Array{ComplexF64,3} = zeros(ComplexF64, 0, 0, 0) + pressure_compression_density::Array{ComplexF64,3} = zeros(ComplexF64, 0, 0, 0) +end + +""" + SurfaceGeometry(mtheta) + +Equilibrium geometry of one flux surface on the periodic θ grid, filled by `surface_geometry!`. + +## Fields + + - `R`, `Z` - cylindrical coordinates of the θ points + - `jac`, `jac_psi` - Jacobian J and its ψ_N-derivative + - `delpsi2`, `dpdt`, `dpdz` - |∇ψ|², ∇ψ·∇θ and ∇ψ·∇ζ + - `g11`, `g12`, `g13`, `g22`, `g23`, `g33` - covariant metric of the tangent basis e_i = ∂x/∂q^i + - `bsq`, `bsq_psi`, `bsq_theta` - metric B² and its ψ_N- and θ-derivatives + - `mu0_j_theta`, `mu0_j_zeta` - equilibrium current μ₀j^θ = −F'/J and μ₀j^ζ = −(μ₀p)'/χ' − qF'/J +""" +struct SurfaceGeometry + R::Vector{Float64} + Z::Vector{Float64} + jac::Vector{Float64} + jac_psi::Vector{Float64} + delpsi2::Vector{Float64} + dpdt::Vector{Float64} + dpdz::Vector{Float64} + g11::Vector{Float64} + g12::Vector{Float64} + g13::Vector{Float64} + g22::Vector{Float64} + g23::Vector{Float64} + g33::Vector{Float64} + bsq::Vector{Float64} + bsq_psi::Vector{Float64} + bsq_theta::Vector{Float64} + mu0_j_theta::Vector{Float64} + mu0_j_zeta::Vector{Float64} +end + +# Per-surface scratch struct with one zeroed length-`mtheta` vector per field. +_zeros_struct(::Type{T}, mtheta::Int) where {T} = T((zeros(eltype(fieldtype(T, i)), mtheta) for i in 1:fieldcount(T))...) + +SurfaceGeometry(mtheta::Int) = _zeros_struct(SurfaceGeometry, mtheta) + +""" + surface_geometry!(geom, equil, psi, thetas; bsq, knot=0, hint=(Ref(1), Ref(1))) -> geom + +Fill `geom` at flux surface `psi` on the θ points `thetas` from `Equilibrium.flux_surface_basis`, +with B² and its derivatives from the metric interpolant `bsq` (`Equilibrium.metric_bsq_interpolant`). +When `psi` is the equilibrium knot `equil.rzphi_xs[knot]` and `thetas` are the θ knots, the geometry +and B² are read from the values stored at the knots rather than evaluated through the interpolants. +""" +function surface_geometry!(geom::SurfaceGeometry, equil::Equilibrium.PlasmaEquilibrium, psi::Float64, thetas::AbstractVector{Float64}; + bsq::FastInterpolations.CubicInterpolantND, knot::Int=0, hint=(Ref(1), Ref(1))) + profiles = equil.profiles + F1 = profiles.F_deriv(psi) + mu0p1 = profiles.P_deriv(psi) + q = profiles.q_spline(psi) + chi1 = 2π * equil.psio + bsq_knots = bsq.nodal_derivs.partials + for (k, theta) in enumerate(thetas) + if knot > 0 + b = Equilibrium.flux_surface_basis(equil, knot, k) + geom.bsq[k] = bsq_knots[1, knot, k] + geom.bsq_psi[k] = bsq_knots[2, knot, k] + geom.bsq_theta[k] = bsq_knots[3, knot, k] + else + b = Equilibrium.flux_surface_basis(equil, psi, theta; hint=hint) + geom.bsq[k] = bsq((psi, theta); hint=hint) + geom.bsq_psi[k] = bsq((psi, theta); deriv=DerivOp(1, 0), hint=hint) + geom.bsq_theta[k] = bsq((psi, theta); deriv=DerivOp(0, 1), hint=hint) + end + geom.R[k] = b.r + geom.Z[k] = b.z + geom.jac[k] = b.jac + geom.jac_psi[k] = b.jac_psi + geom.delpsi2[k] = b.delpsi2 + geom.dpdt[k] = b.dpdt + geom.dpdz[k] = b.dpdz + geom.g11[k] = b.g11 + geom.g12[k] = b.g12 + geom.g13[k] = b.g13 + geom.g22[k] = b.g22 + geom.g23[k] = b.g23 + geom.g33[k] = b.g33 + geom.mu0_j_theta[k] = -F1 / b.jac + geom.mu0_j_zeta[k] = -mu0p1 / chi1 - F1 * q / b.jac + end + return geom +end + +""" + CurvatureKernel(mtheta) + +Per-θ pieces of the destabilizing kernel on one surface, filled by `curvature_kernel!`. + +## Fields + + - `shear` - DCON local shear S + - `curvature` - curvature projection κ_ψ = κ·∇ψ + - `sigma` - parallel current density ratio σ = j·B/B² + - `K1`, `K2`, `K3` - K1 = |∇ψ|²σS, K2 = μ₀B²σ², K3 = 2p'κ_ψ +""" +struct CurvatureKernel + shear::Vector{Float64} + curvature::Vector{Float64} + sigma::Vector{Float64} + K1::Vector{Float64} + K2::Vector{Float64} + K3::Vector{Float64} +end + +CurvatureKernel(mtheta::Int) = _zeros_struct(CurvatureKernel, mtheta) + +""" + curvature_kernel!(kern, geom, equil, psi, thetas_ext; hint=Ref(1)) -> kern + +Fill the curvature kernel at `psi` from the geometry `geom` of the periodic θ points +`thetas_ext[1:end-1]` (`thetas_ext` closes the period at 1.0). Fortran `gpeq_shear`, +`gpeq_curvature`, `gpeq_K`; main.tex §"Expression of K in DCON coordinates": + + S = (χ'²/J) [ q' + ∂_θ( (q ∇ψ·∇θ − ∇ψ·∇ζ)/|∇ψ|² ) ] (periodic cubic-spline derivative) + κ_ψ = (|∇ψ|²/B²) [ (μ₀p)' + ½∂_ψB² + ½∂_θB² (∇ψ·∇θ)/|∇ψ|² ] + σ = j·B/B² with j^θ = −F'/(μ₀J), j^ζ = q j^θ − p'/χ', B^θ = χ'/J, B^ζ = qχ'/J + +where F' = d(2πF)/dψ and p' = (μ₀p)'/μ₀ come from the equilibrium profile splines. +""" +function curvature_kernel!(kern::CurvatureKernel, geom::SurfaceGeometry, equil::Equilibrium.PlasmaEquilibrium, psi::Float64, + thetas_ext::Vector{Float64}; hint=Ref(1)) + profiles = equil.profiles + q = profiles.q_spline(psi; hint=hint) + q1 = profiles.q_deriv(psi; hint=hint) + mu0p1 = profiles.P_deriv(psi; hint=hint) + p1 = mu0p1 / MU_0 + chi1 = 2π * equil.psio + mtheta = length(geom.jac) + + # Shear geometry term on the closed θ grid; its θ-derivative comes from a periodic cubic spline. + ratio = Vector{Float64}(undef, mtheta + 1) + for k in 1:mtheta + ratio[k] = (q * geom.dpdt[k] - geom.dpdz[k]) / geom.delpsi2[k] + end + ratio[mtheta+1] = ratio[1] + dratio = deriv1(cubic_interp(thetas_ext, ratio; bc=PeriodicBC())) + + for k in 1:mtheta + jac = geom.jac[k] + kern.shear[k] = (chi1^2 / jac) * (q1 + dratio(thetas_ext[k])) + # Deviation from Fortran gpeq_curvature/gpeq_K: κ_ψ, σ and K2 use the metric B², not eqfun B², so K2|ξ_n|² matches |V|². + kern.curvature[k] = (geom.delpsi2[k] / geom.bsq[k]) * (mu0p1 + 0.5 * geom.bsq_psi[k] + 0.5 * geom.bsq_theta[k] * geom.dpdt[k] / geom.delpsi2[k]) + jt = geom.mu0_j_theta[k] / MU_0 + jz = geom.mu0_j_zeta[k] / MU_0 + bt = chi1 / jac + bz = q * chi1 / jac + jdotb = geom.g22[k] * bt * jt + geom.g33[k] * bz * jz + geom.g23[k] * (bt * jz + bz * jt) + sigma = jdotb / geom.bsq[k] + kern.sigma[k] = sigma + kern.K1[k] = geom.delpsi2[k] * sigma * kern.shear[k] + kern.K2[k] = MU_0 * geom.bsq[k] * sigma^2 + kern.K3[k] = 2.0 * p1 * kern.curvature[k] + end + return kern +end + +# ψ_N gap below which the last equilibrium knot counts as the solution edge. +const _EDGE_KNOT_TOL = 1e-10 + +""" + decomposition_grid(ffs) -> (psi, is_knot) + +Radial nodes of the decomposition: the equilibrium ψ knots inside the ξ solution domain, where +the metric coefficients, the Euler-Lagrange matrices and the geometry are all exact and the +standard-form integrand equals the DCON Lagrangian density (1e-11 pointwise), plus the solution +edge appended when the last knot falls short of it so the radial integral reaches the plasma edge. +The Fortran `gpout_recon` uses these knots for its profiles but integrates its energy check on +the DCON `psifac` grid. `is_knot` marks the exact +nodes; the (∇×b_eff)·∇ψ residual is an algebraic identity there and is reported on them only. +""" +function decomposition_grid(ffs::ForceFreeStatesResult) + solution = ffs.solution + psi_lo = solution.psi_store[1] + psi_hi = solution.psi_store[solution.step] + psi = [x for x in ffs.equil.rzphi_xs if psi_lo <= x <= psi_hi] + is_knot = trues(length(psi)) + if isempty(psi) || psi[end] < psi_hi - _EDGE_KNOT_TOL + push!(psi, psi_hi) + push!(is_knot, false) + end + return psi, is_knot +end + +""" + node_solution(ffs, psi) -> (u1, du1) + +Euler-Lagrange basis Ξ_ψ and dΞ_ψ/dψ at the radial nodes `psi`, `(numpert_total × numpert_total × +length(psi))` in the basis of the stored solution. Each node state comes from +`ForceFreeStates.solution_state_at`, so ξ is the ODE's own value there rather than a spline through +the saved steps, whose spacing follows the integrator's tolerance control. +""" +function node_solution(ffs::ForceFreeStatesResult, psi::Vector{Float64}) + npert = ffs.numpert_total + u1 = zeros(ComplexF64, npert, npert, length(psi)) + du1 = zeros(ComplexF64, npert, npert, length(psi)) + du = zeros(ComplexF64, npert, npert, 2) + hints = (Ref(1), Ref(1)) + for (i, p) in enumerate(psi) + u = ForceFreeStates.solution_state_at(ffs, p) + ForceFreeStates.el_derivatives!(du, u, false, ffs.equil, ffs.mats, ffs, p, hints...) + u1[:, :, i] .= @view u[:, :, 1] + du1[:, :, i] .= @view du[:, :, 1] + end + return u1, du1 +end + +""" + eigenmode_modes(ffs, k, psi, u1, du1) -> NamedTuple + +Mode-space displacement and perturbed field of free-boundary eigenmode `k` on the radial nodes +`psi`, `(length(psi) × mpert)` matrices named after the gpeq arrays: `xi_psi` = ξ^ψ, +`xi_psi1` = ∂ξ^ψ/∂ψ, `Jxi_psi/theta/zeta` = Jξ^i, `Jb_psi/theta/zeta` = Jb^i. The edge boundary +condition is the eigenvector `wt[:, k]` (unit power norm), which fixes the combination of the node +basis `u1`, `du1` from `node_solution`. ξ_s = −A⁻¹(Bξ' + Cξ) is recomputed at every node +from the ideal Euler-Lagrange matrices there (`ForceFreeStates.compute_node_xi_s!`, Fortran +`gpeq_sol`), which keeps the fields algebraically consistent with the metric at the equilibrium +knots. The rest of the chain is the unregularized one `reconstruct_physical_fields` runs for a +forced response. +""" +function eigenmode_modes(ffs::ForceFreeStatesResult, k::Int, psi::Vector{Float64}, u1::Array{ComplexF64,3}, du1::Array{ComplexF64,3}) + solution = ffs.solution + equil = ffs.equil + u1_edge = @view solution.u_store[:, :, 1, solution.step] + alpha = u1_edge \ Vector{ComplexF64}(@view ffs.free_boundary.wt[:, k]) + + nnode = length(psi) + mpert = ffs.mpert + xi_psi = zeros(ComplexF64, nnode, mpert) + xi_psi1 = zeros(ComplexF64, nnode, mpert) + xi_s = zeros(ComplexF64, nnode, mpert) + hint = Ref(1) + for (i, p) in enumerate(psi) + mul!(view(xi_psi, i, :), view(u1, :, :, i), alpha) + mul!(view(xi_psi1, i, :), view(du1, :, :, i), alpha) + ForceFreeStates.compute_node_xi_s!(view(xi_s, i, :), view(xi_psi1, i, :), view(xi_psi, i, :), ffs.mats, p; hint=hint) + end + + chi1 = 2π * equil.psio + Jb_psi, Jb_theta, Jb_zeta = compute_perturbed_field_modes(xi_psi, xi_psi1, xi_s, psi, equil, ffs) + # Deviation from Fortran gpeq_sol/gpeq_contra: always unregularized (reg_spot = 0), the Fortran regularizes when reg_flag is set. + Jxi_psi, Jxi_theta, Jxi_zeta, _, _ = compute_contra_displacements(xi_psi, xi_psi1, xi_s ./ chi1, psi, equil, ffs, ffs.metric; reg_spot=0.0) + return (; xi_psi, xi_psi1, Jxi_psi, Jxi_theta, Jxi_zeta, Jb_psi, Jb_theta, Jb_zeta) +end + +""" + SurfaceFields(mtheta) + +θ-space eigenmode fields on one surface, filled by `surface_fields!`. + +## Fields + + - `xi_psi`, `xi_psi1` - ξ^ψ and ∂ξ^ψ/∂ψ + - `Jxi_psi`, `Jxi_theta`, `Jxi_zeta` - Jacobian-weighted contravariant displacement Jξ^i + - `dJxi_theta` - ∂_θ(Jξ^θ) + - `Jb_psi`, `Jb_theta`, `Jb_zeta` - Jacobian-weighted contravariant perturbed field Jb^i + - `b_cov_psi`, `b_cov_theta`, `b_cov_zeta` - covariant perturbed field b_i = g_ij Jb^j/J + - `xi_n` - normal displacement ξ_n = Jξ^ψ/(J|∇ψ|) +""" +struct SurfaceFields + xi_psi::Vector{ComplexF64} + xi_psi1::Vector{ComplexF64} + Jxi_psi::Vector{ComplexF64} + Jxi_theta::Vector{ComplexF64} + Jxi_zeta::Vector{ComplexF64} + dJxi_theta::Vector{ComplexF64} + Jb_psi::Vector{ComplexF64} + Jb_theta::Vector{ComplexF64} + Jb_zeta::Vector{ComplexF64} + b_cov_psi::Vector{ComplexF64} + b_cov_theta::Vector{ComplexF64} + b_cov_zeta::Vector{ComplexF64} + xi_n::Vector{ComplexF64} +end + +SurfaceFields(mtheta::Int) = _zeros_struct(SurfaceFields, mtheta) + +""" + surface_fields!(sf, ft, modes, ipsi, mvals, geom, work_modes) -> sf + +Inverse-transform row `ipsi` of every matrix in `modes` onto the θ grid of `ft`, take the +θ-derivative of Jξ^θ spectrally, lower the field with the local metric of `geom` +(b_i = g_ij Jb^j/J), and form Jξ^ψ = J·ξ^ψ and ξ_n = ξ^ψ/|∇ψ| pointwise. The pointwise products +keep every θ harmonic, so |V|² and K₂|ξ_n|² cancel exactly in the effective-field form; truncating +them to the mode range leaves a ψ-resolution-independent, mpert-dependent deficit. Jξ^θ and Jξ^ζ +still come from the mode-truncated `jmat` convolution of `compute_contra_displacements`. +`work_modes` is an `mpert` scratch vector. +""" +function surface_fields!(sf::SurfaceFields, ft::Utilities.FourierTransforms.FourierTransform, modes::NamedTuple, ipsi::Int, + mvals::Vector{Int}, geom::SurfaceGeometry, work_modes::Vector{ComplexF64}) + for name in (:xi_psi, :xi_psi1, :Jxi_theta, :Jxi_zeta, :Jb_psi, :Jb_theta, :Jb_zeta) + Utilities.FourierTransforms.inverse_transform!(getfield(sf, name), ft, view(getfield(modes, name), ipsi, :)) + end + # Deviation from Fortran gpeq_firstform: ∂_θ(Jξ^θ) is spectral (2πi m per mode), not a periodic-spline derivative. + for i in eachindex(mvals) + work_modes[i] = 2π * im * mvals[i] * modes.Jxi_theta[ipsi, i] + end + Utilities.FourierTransforms.inverse_transform!(sf.dJxi_theta, ft, work_modes) + # Deviation from Fortran gpeq_contra/gpeq_cova/gpeq_normal: Jξ^ψ, b_i and ξ_n are pointwise, not band-truncated. + for k in eachindex(geom.jac) + jac = geom.jac[k] + sf.Jxi_psi[k] = jac * sf.xi_psi[k] + sf.b_cov_psi[k] = (geom.g11[k] * sf.Jb_psi[k] + geom.g12[k] * sf.Jb_theta[k] + geom.g13[k] * sf.Jb_zeta[k]) / jac + sf.b_cov_theta[k] = (geom.g12[k] * sf.Jb_psi[k] + geom.g22[k] * sf.Jb_theta[k] + geom.g23[k] * sf.Jb_zeta[k]) / jac + sf.b_cov_zeta[k] = (geom.g13[k] * sf.Jb_psi[k] + geom.g23[k] * sf.Jb_theta[k] + geom.g33[k] * sf.Jb_zeta[k]) / jac + sf.xi_n[k] = sf.Jxi_psi[k] / (jac * sqrt(geom.delpsi2[k])) + end + return sf +end + +""" + EffectiveField(mtheta) + +θ-space components of the effective field b_eff = b + ξ_n(μ₀ j × n̂) on one surface, filled by +`effective_field!`. + +## Fields + + - `Jc_psi`, `Jc_theta`, `Jc_zeta` - Jacobian-weighted contravariant components J b_eff^i + - `c_cov_psi`, `c_cov_theta`, `c_cov_zeta` - covariant components b_eff,i +""" +struct EffectiveField + Jc_psi::Vector{ComplexF64} + Jc_theta::Vector{ComplexF64} + Jc_zeta::Vector{ComplexF64} + c_cov_psi::Vector{ComplexF64} + c_cov_theta::Vector{ComplexF64} + c_cov_zeta::Vector{ComplexF64} +end + +EffectiveField(mtheta::Int) = _zeros_struct(EffectiveField, mtheta) + +""" + effective_field!(eff, geom, sf) -> eff + +Build b_eff = b + ξ_n(μ₀ j × n̂) from the perturbed field and Jξ^ψ of `sf` and the equilibrium +current of `geom` (Fortran `gpeq_c`, main.tex Eqs. 148–151): + + J b_eff^ψ = J b^ψ + J b_eff^θ = J b^θ + (Jξ^ψ)/(J|∇ψ|²) [ μ₀j^θ g_θζ + μ₀j^ζ g_ζζ ] + J b_eff^ζ = J b^ζ − (Jξ^ψ)/(J|∇ψ|²) [ μ₀j^θ g_θθ + μ₀j^ζ g_θζ ] + b_eff,ψ = b_ψ + (Jξ^ψ)/|∇ψ|² [ μ₀j^θ (∇ψ·∇ζ) − μ₀j^ζ (∇ψ·∇θ) ] + b_eff,θ = b_θ + (Jξ^ψ) μ₀j^ζ + b_eff,ζ = b_ζ − (Jξ^ψ) μ₀j^θ +""" +function effective_field!(eff::EffectiveField, geom::SurfaceGeometry, sf::SurfaceFields) + # Deviation from Fortran gpeq_c: b_eff^i and b_eff,i are formed pointwise in θ, not band-truncated. + for k in eachindex(geom.jac) + jac = geom.jac[k] + mu0jt = geom.mu0_j_theta[k] + mu0jz = geom.mu0_j_zeta[k] + d2 = geom.delpsi2[k] + xw = sf.Jxi_psi[k] + eff.Jc_psi[k] = sf.Jb_psi[k] + eff.Jc_theta[k] = sf.Jb_theta[k] + xw / (d2 * jac) * (mu0jt * geom.g23[k] + mu0jz * geom.g33[k]) + eff.Jc_zeta[k] = sf.Jb_zeta[k] - xw / (d2 * jac) * (mu0jt * geom.g22[k] + mu0jz * geom.g23[k]) + eff.c_cov_psi[k] = sf.b_cov_psi[k] + xw / d2 * (mu0jt * geom.dpdz[k] - mu0jz * geom.dpdt[k]) + eff.c_cov_theta[k] = sf.b_cov_theta[k] + xw * mu0jz + eff.c_cov_zeta[k] = sf.b_cov_zeta[k] - xw * mu0jt + end + return eff +end + +""" + curl_residual!(res, ft, eff, geom, mvals, n, theta_modes, zeta_modes, work_modes, work_fun) -> scale + +Evaluate (∇×b_eff)·∇ψ = (∂_θ b_eff,ζ − ∂_ζ b_eff,θ)/J on the θ grid into `res` (Fortran +`gpeq_cveri`): the covariant θ and ζ components are transformed to mode space, differentiated +exactly (∂_θ → 2πi m, ∂_ζ → −2πi n), transformed back and divided by J. Returns +max|∂_θ b_eff,ζ/J|, the size of the cancelling terms, for the relative residual. The four +trailing arguments are scratch buffers (`mpert`, `mpert`, `mpert`, `mtheta`). +""" +function curl_residual!(res::Vector{ComplexF64}, ft::Utilities.FourierTransforms.FourierTransform, eff::EffectiveField, geom::SurfaceGeometry, + mvals::Vector{Int}, n::Int, theta_modes::Vector{ComplexF64}, zeta_modes::Vector{ComplexF64}, work_modes::Vector{ComplexF64}, + work_fun::Vector{ComplexF64}) + Utilities.FourierTransforms.transform!(theta_modes, ft, eff.c_cov_theta) + Utilities.FourierTransforms.transform!(zeta_modes, ft, eff.c_cov_zeta) + for i in eachindex(mvals) + work_modes[i] = 2π * im * (mvals[i] * zeta_modes[i] + n * theta_modes[i]) + end + Utilities.FourierTransforms.inverse_transform!(res, ft, work_modes) + res ./= geom.jac + for i in eachindex(mvals) + work_modes[i] = 2π * im * mvals[i] * zeta_modes[i] + end + Utilities.FourierTransforms.inverse_transform!(work_fun, ft, work_modes) + work_fun ./= geom.jac + return maximum(abs, work_fun) +end + +""" + SurfaceDensities(mtheta) + +Per-θ energy densities on one surface (integral = Σ density·J/mtheta, see `surface_integral`), +filled by `surface_densities!`. Field names match the profiles of `EnergyDecompositionResult`. + +## Fields + + - `effective_b_squared` - |b_eff|²/μ₀ + - `effective_b_squared_psi`, `effective_b_squared_theta`, `effective_b_squared_zeta` - its component contributions + - `shear_current`, `parallel_current_squared`, `pressure_curvature` - K1|ξ_n|², K2|ξ_n|², K3|ξ_n|² + - `b_squared` - standard form |b|²/μ₀ + - `current_coupling` - standard form j·(b×ξ*) + - `pressure_compression` - standard form (∇·ξ)*(ξ·∇p) +""" +struct SurfaceDensities + effective_b_squared::Vector{Float64} + effective_b_squared_psi::Vector{Float64} + effective_b_squared_theta::Vector{Float64} + effective_b_squared_zeta::Vector{Float64} + shear_current::Vector{Float64} + parallel_current_squared::Vector{Float64} + pressure_curvature::Vector{Float64} + b_squared::Vector{Float64} + current_coupling::Vector{ComplexF64} + pressure_compression::Vector{ComplexF64} +end + +SurfaceDensities(mtheta::Int) = _zeros_struct(SurfaceDensities, mtheta) + +# Per-surface profiles of each form; the names are shared by SurfaceDensities and EnergyDecompositionResult. +const _ED_EFFECTIVE_PROFILES = (:effective_b_squared, :effective_b_squared_psi, :effective_b_squared_theta, :effective_b_squared_zeta, + :shear_current, :parallel_current_squared, :pressure_curvature) +const _ED_STANDARD_PROFILES = (:b_squared, :current_coupling, :pressure_compression) + +# Result field holding the per-(ψ,θ) density of a profile. +_density_field(name::Symbol) = Symbol(name, :_density) + +""" + surface_densities!(den, geom, kern, sf, eff, equil, psi, n; effective_field_form, standard_form, hint=Ref(1)) -> den + +|b_eff|² = Re[(J b_eff^i)* b_eff,i]/J and |b|² = Re[(J b^i)* b_i]/J; the current term is +j·(b×ξ*) = [ −μ₀j^θ (Jb^ψ (Jξ^ζ)* − Jb^ζ (Jξ^ψ)*) + μ₀j^ζ (Jb^ψ (Jξ^θ)* − Jb^θ (Jξ^ψ)*) ]/(μ₀J), and +∇·ξ = ∂_ψξ^ψ + (∂_ψJ/J)ξ^ψ + [∂_θ(Jξ^θ) − 2πi n Jξ^ζ]/J with ξ·∇p = p'ξ^ψ. +""" +function surface_densities!(den::SurfaceDensities, geom::SurfaceGeometry, kern::CurvatureKernel, sf::SurfaceFields, eff::EffectiveField, + equil::Equilibrium.PlasmaEquilibrium, psi::Float64, n::Int; effective_field_form::Bool, standard_form::Bool, hint=Ref(1)) + p1 = equil.profiles.P_deriv(psi; hint=hint) / MU_0 + for k in eachindex(geom.jac) + jac = geom.jac[k] + if effective_field_form + dpsi = real(conj(eff.Jc_psi[k]) * eff.c_cov_psi[k]) / (MU_0 * jac) + dthe = real(conj(eff.Jc_theta[k]) * eff.c_cov_theta[k]) / (MU_0 * jac) + dzet = real(conj(eff.Jc_zeta[k]) * eff.c_cov_zeta[k]) / (MU_0 * jac) + den.effective_b_squared_psi[k] = dpsi + den.effective_b_squared_theta[k] = dthe + den.effective_b_squared_zeta[k] = dzet + den.effective_b_squared[k] = dpsi + dthe + dzet + xin2 = abs2(sf.xi_n[k]) + den.shear_current[k] = kern.K1[k] * xin2 + den.parallel_current_squared[k] = kern.K2[k] * xin2 + den.pressure_curvature[k] = kern.K3[k] * xin2 + end + if standard_form + mu0jt = geom.mu0_j_theta[k] + mu0jz = geom.mu0_j_zeta[k] + den.b_squared[k] = real(conj(sf.Jb_psi[k]) * sf.b_cov_psi[k] + conj(sf.Jb_theta[k]) * sf.b_cov_theta[k] + conj(sf.Jb_zeta[k]) * sf.b_cov_zeta[k]) / (MU_0 * jac) + det_term = + -mu0jt * (sf.Jb_psi[k] * conj(sf.Jxi_zeta[k]) - sf.Jb_zeta[k] * conj(sf.Jxi_psi[k])) + + mu0jz * (sf.Jb_psi[k] * conj(sf.Jxi_theta[k]) - sf.Jb_theta[k] * conj(sf.Jxi_psi[k])) + den.current_coupling[k] = det_term / (MU_0 * jac) + divxi = sf.xi_psi1[k] + (geom.jac_psi[k] / jac) * sf.xi_psi[k] + sf.dJxi_theta[k] / jac - (2π * im * n) * sf.Jxi_zeta[k] / jac + den.pressure_compression[k] = conj(divxi) * (p1 * sf.xi_psi[k]) + end + end + return den +end + +""" + surface_integral(density, jac) -> Number + +∮ J·density dθ over the periodic θ grid as the Riemann sum Σ density·J/mtheta, exact for the +band-limited integrands of the decomposition and identical to the Fortran sums. +""" +function surface_integral(density::AbstractVector, jac::AbstractVector{Float64}) + acc = zero(eltype(density)) + for k in eachindex(density) + acc += density[k] * jac[k] + end + return acc / length(density) +end + +# Relative (∇×b_eff)·∇ψ residual above which the effective field is reported as inconsistent. +const _CURL_RESIDUAL_WARN_TOL = 1e-8 +# Relative mismatch between a reconstructed δW_p and eigenmode_plasma_energies that triggers a warning. +const _ENERGY_MISMATCH_WARN_TOL = 0.05 + +# Scale from the Bernstein integrals (2μ₀δW form) to the normalization of eigenmode_plasma_energies. +_energy_scale(equil::Equilibrium.PlasmaEquilibrium) = 2 * MU_0 / equil.psio^2 + +# Cubic-spline integral over ψ with the running value at every node (Fortran recon_int = "spline"). +# Repeated ψ nodes are collapsed for the fit and the running value is carried across them. +function _radial_integral(psi::Vector{Float64}, y::AbstractVector{<:Real}) + keep = trues(length(psi)) + for i in 2:length(psi) + keep[i] = psi[i] > psi[i-1] + end + itp = cubic_interp(psi[keep], Float64.(y[keep]); bc=CubicFit()) + total = FastInterpolations.integrate(itp) + cum_keep = vec(FastInterpolations.cumulative_integrate(itp)) + cum = Vector{Float64}(undef, length(psi)) + j = 0 + for i in eachindex(psi) + j += keep[i] + cum[i] = cum_keep[j] + end + return total, cum +end + +function _radial_integral(psi::Vector{Float64}, y::AbstractVector{<:Complex}) + tr, cr = _radial_integral(psi, real.(y)) + ti, ci = _radial_integral(psi, imag.(y)) + return complex(tr, ti), complex.(cr, ci) +end + +# Replace result field `f` by zeros of its element type and shape `dims`. +_zero_field!(res::EnergyDecompositionResult, f::Symbol, dims::Int...) = setfield!(res, f, zeros(eltype(fieldtype(EnergyDecompositionResult, f)), dims...)) + +# Allocate the result arrays of the enabled forms; disabled forms keep their empty defaults. +function _allocate_result!(res::EnergyDecompositionResult, npsi::Int, mtheta::Int, nk::Int, thetas::Vector{Float64}, + effective_field_form::Bool, standard_form::Bool, write_densities::Bool) + res.dW_plasma_reference = zeros(nk) + if effective_field_form + for f in (_ED_EFFECTIVE_PROFILES..., :dW_cumulative) + _zero_field!(res, f, npsi, nk) + end + for f in (:dW_plasma, :dW_plasma_relative_error, :effective_b_curl_residual_max, :effective_b_curl_residual_rms, :effective_b_curl_residual_relative) + _zero_field!(res, f, nk) + end + end + if standard_form + for f in (_ED_STANDARD_PROFILES..., :dW_cumulative_standard_form) + _zero_field!(res, f, npsi, nk) + end + for f in (:dW_plasma_standard_form, :dW_plasma_standard_form_relative_error) + _zero_field!(res, f, nk) + end + end + if write_densities + res.theta = copy(thetas) + res.R = zeros(npsi, mtheta) + res.Z = zeros(npsi, mtheta) + effective_field_form && _zero_field!(res, :effective_b_curl_residual_density, npsi, mtheta, nk) + for f in _enabled_profiles(effective_field_form, standard_form) + _zero_field!(res, _density_field(f), npsi, mtheta, nk) + end + end + return res +end + +_enabled_profiles(effective_field_form::Bool, standard_form::Bool) = + ((effective_field_form ? _ED_EFFECTIVE_PROFILES : ())..., (standard_form ? _ED_STANDARD_PROFILES : ())...) + +""" + decompose_energy(ffs; eigenmodes=[1], effective_field_form=true, standard_form=true, write_densities=false, verbose=ffs.control.verbose) + -> EnergyDecompositionResult + +Decompose the plasma potential energy of the free-boundary eigenmodes `eigenmodes` of a forward +ForceFreeStates solve on the (ψ,θ) grid (Fortran `gpout_recon` and `gpout_recon2`). For each +eigenmode the unregularized gpeq chain gives ξ and b in mode space on the radial nodes of +`decomposition_grid`; every node is transformed to θ space, the effective field +b_eff = b + ξ_n(μ₀ j × n̂) and the curvature kernel K are built, and the per-surface integrals of +the enabled forms are integrated in ψ with a cubic spline. All energies are scaled by 2μ₀/ψ₀² so +`dW_plasma` reads in the normalization of `eigenmode_plasma_energies`, whose entry `ep[k]` is +stored as `dW_plasma_reference`. Two checks run every time: the relative (∇×b_eff)·∇ψ residual +over the equilibrium knots and the relative mismatch of each total against the reference; both +warn above their tolerance and are stored as data. + +It covers the free-boundary DCON eigenmodes only: the edge condition is the eigenvector +`wt[:, k]` and the reference is `ep[k]`, which fixed-boundary or forced-response perturbations +do not supply. Requires the dense ξ solution (`integrator = "forward"`), the free-boundary +eigenmodes (`vac_flag = true`), an ideal solve, a single toroidal mode number, and at least one +enabled form; anything else throws an `ArgumentError`. +""" +function decompose_energy(ffs::ForceFreeStatesResult; eigenmodes::Vector{Int}=[1], effective_field_form::Bool=true, standard_form::Bool=true, + write_densities::Bool=false, verbose::Bool=ffs.control.verbose) + (effective_field_form || standard_form) || throw(ArgumentError("EnergyDecomposition: enable at least one of effective_field_form and standard_form")) + ffs.solution === nothing && throw(ArgumentError("EnergyDecomposition needs the dense ξ solution of a forward run (integrator = \"forward\")")) + ffs.free_boundary === nothing && throw(ArgumentError("EnergyDecomposition needs the free-boundary eigenmodes (vac_flag = true)")) + ffs.nlow == ffs.nhigh || throw(ArgumentError("EnergyDecomposition handles a single toroidal mode number (nn_low == nn_high)")) + ffs.mats.kinetic === nothing || throw(ArgumentError("EnergyDecomposition decomposes ideal eigenmodes only (kinetic_factor = 0)")) + free = ffs.free_boundary + equil = ffs.equil + nmodes = size(free.wt, 2) + for k in eigenmodes + 1 <= k <= nmodes || throw(ArgumentError("EnergyDecomposition: eigenmode $k is outside 1:$nmodes")) + end + + # Deviation from Fortran gpout_recon: the totals integrate over the equilibrium knots, the Fortran ones over the ODE psifac grid. + psi, is_knot = decomposition_grid(ffs) + knots = [on_knot ? searchsortedfirst(equil.rzphi_xs, p) : 0 for (p, on_knot) in zip(psi, is_knot)] + npsi = length(psi) + thetas_ext = collect(equil.rzphi_ys) + mtheta = length(thetas_ext) - 1 + thetas = thetas_ext[1:mtheta] + mpert = ffs.mpert + mvals = collect(ffs.mlow:ffs.mhigh) + nn = ffs.nlow + scale = _energy_scale(equil) + nk = length(eigenmodes) + profiles = _enabled_profiles(effective_field_form, standard_form) + + res = EnergyDecompositionResult(; psi=copy(psi), eigenmode_index=copy(eigenmodes)) + _allocate_result!(res, npsi, mtheta, nk, thetas, effective_field_form, standard_form, write_densities) + + # Per-thread scratch; the ψ loops are threaded :static so threadid() indexes it safely. + ft = Utilities.FourierTransforms.FourierTransform(mtheta, mpert, ffs.mlow) + nt = Threads.maxthreadid() + sfs = [SurfaceFields(mtheta) for _ in 1:nt] + effs = [EffectiveField(mtheta) for _ in 1:nt] + dens = [SurfaceDensities(mtheta) for _ in 1:nt] + curls = [zeros(ComplexF64, mtheta) for _ in 1:nt] + mode_bufs = [(zeros(ComplexF64, mpert), zeros(ComplexF64, mpert), zeros(ComplexF64, mpert)) for _ in 1:nt] + fun_bufs = [zeros(ComplexF64, mtheta) for _ in 1:nt] + hints2d = [(Ref(1), Ref(1)) for _ in 1:nt] + hints1d = [Ref(1) for _ in 1:nt] + + # Equilibrium geometry and curvature kernel of every surface, shared by all eigenmodes. + bsq = Equilibrium.metric_bsq_interpolant(equil) + geoms = [SurfaceGeometry(mtheta) for _ in 1:npsi] + kerns = [CurvatureKernel(mtheta) for _ in 1:npsi] + Threads.@threads :static for ipsi in 1:npsi + tid = Threads.threadid() + surface_geometry!(geoms[ipsi], equil, psi[ipsi], thetas; bsq, knot=knots[ipsi], hint=hints2d[tid]) + effective_field_form && curvature_kernel!(kerns[ipsi], geoms[ipsi], equil, psi[ipsi], thetas_ext; hint=hints1d[tid]) + end + if write_densities + for ipsi in 1:npsi + res.R[ipsi, :] .= geoms[ipsi].R + res.Z[ipsi, :] .= geoms[ipsi].Z + end + end + + # Per-surface (∇×b_eff)·∇ψ statistics of the current eigenmode, filled by the threaded loop. + curl_max = zeros(npsi) + curl_sumsq = zeros(npsi) + curl_scale = zeros(npsi) + + # The θ transforms are small BLAS matvecs; one BLAS thread per Julia thread avoids oversubscription. + blas_threads = BLAS.get_num_threads() + BLAS.set_num_threads(1) + try + u1, du1 = node_solution(ffs, psi) + for (ik, k) in enumerate(eigenmodes) + modes = eigenmode_modes(ffs, k, psi, u1, du1) + Threads.@threads :static for ipsi in 1:npsi + tid = Threads.threadid() + geom = geoms[ipsi] + sf = sfs[tid] + eff = effs[tid] + den = dens[tid] + theta_modes, zeta_modes, work_modes = mode_bufs[tid] + surface_fields!(sf, ft, modes, ipsi, mvals, geom, work_modes) + if effective_field_form + effective_field!(eff, geom, sf) + curl_scale[ipsi] = curl_residual!(curls[tid], ft, eff, geom, mvals, nn, theta_modes, zeta_modes, work_modes, fun_bufs[tid]) + curl_max[ipsi] = maximum(abs, curls[tid]) + curl_sumsq[ipsi] = sum(abs2, curls[tid]) + end + surface_densities!(den, geom, kerns[ipsi], sf, eff, equil, psi[ipsi], nn; effective_field_form, standard_form, hint=hints1d[tid]) + for name in profiles + getfield(res, name)[ipsi, ik] = scale * surface_integral(getfield(den, name), geom.jac) + end + if write_densities + for name in profiles + getfield(res, _density_field(name))[ipsi, :, ik] .= scale .* getfield(den, name) + end + effective_field_form && (res.effective_b_curl_residual_density[ipsi, :, ik] .= curls[tid]) + end + end + + reference = real(free.ep[k]) + res.dW_plasma_reference[ik] = reference + if effective_field_form + tot_eb, cum_eb = _radial_integral(psi, @view res.effective_b_squared[:, ik]) + tot_k1, cum_k1 = _radial_integral(psi, @view res.shear_current[:, ik]) + tot_k2, cum_k2 = _radial_integral(psi, @view res.parallel_current_squared[:, ik]) + tot_k3, cum_k3 = _radial_integral(psi, @view res.pressure_curvature[:, ik]) + res.dW_cumulative[:, ik] .= 0.5 .* (cum_eb .- cum_k1 .- cum_k2 .- cum_k3) + res.dW_plasma[ik] = 0.5 * (tot_eb - tot_k1 - tot_k2 - tot_k3) + res.dW_plasma_relative_error[ik] = (res.dW_plasma[ik] - reference) / abs(reference) + res.effective_b_curl_residual_max[ik] = maximum(curl_max[is_knot]) + res.effective_b_curl_residual_rms[ik] = sqrt(sum(curl_sumsq[is_knot]) / (count(is_knot) * mtheta)) + res.effective_b_curl_residual_relative[ik] = res.effective_b_curl_residual_max[ik] / maximum(curl_scale[is_knot]) + end + if standard_form + tot_b2, cum_b2 = _radial_integral(psi, @view res.b_squared[:, ik]) + tot_cc, cum_cc = _radial_integral(psi, @view res.current_coupling[:, ik]) + tot_pc, cum_pc = _radial_integral(psi, @view res.pressure_compression[:, ik]) + res.dW_cumulative_standard_form[:, ik] .= 0.5 .* (cum_b2 .- cum_cc .+ cum_pc) + res.dW_plasma_standard_form[ik] = 0.5 * (tot_b2 - tot_cc + tot_pc) + res.dW_plasma_standard_form_relative_error[ik] = (real(res.dW_plasma_standard_form[ik]) - reference) / abs(reference) + end + + verbose && print_energy_decomposition_summary(res, ik) + if effective_field_form + rel = res.effective_b_curl_residual_relative[ik] + rel > _CURL_RESIDUAL_WARN_TOL && + @warn "EnergyDecomposition: (∇×b_eff)·∇ψ of eigenmode $k does not vanish (relative residual $(@sprintf("%.2e", rel)) > $(_CURL_RESIDUAL_WARN_TOL))" + err = res.dW_plasma_relative_error[ik] + abs(err) > _ENERGY_MISMATCH_WARN_TOL && + @warn "EnergyDecomposition: effective-field δW_p of eigenmode $k differs from eigenmode_plasma_energies by $(@sprintf("%.2f", 100 * err)) %" + end + if standard_form + err = res.dW_plasma_standard_form_relative_error[ik] + abs(err) > _ENERGY_MISMATCH_WARN_TOL && + @warn "EnergyDecomposition: standard-form δW_p of eigenmode $k differs from eigenmode_plasma_energies by $(@sprintf("%.2f", 100 * err)) %" + end + end + finally + BLAS.set_num_threads(blas_threads) + end + return res +end + +# Total of a stored profile over ψ (the profiles are already scaled). +_profile_total(res::EnergyDecompositionResult, profile::AbstractMatrix, ik::Int) = first(_radial_integral(res.psi, @view profile[:, ik])) + +""" + print_energy_decomposition_summary(res, ik) + +Log the totals of decomposed eigenmode number `ik` of `res` and the two consistency checks. +""" +function print_energy_decomposition_summary(res::EnergyDecompositionResult, ik::Int) + k = res.eigenmode_index[ik] + lines = ["Energy decomposition of free-boundary eigenmode $k (power-normalized, per unit ⟨|ξ|²⟩)"] + if !isempty(res.dW_plasma) + eb = _profile_total(res, res.effective_b_squared, ik) + kk = _profile_total(res, res.shear_current, ik) + _profile_total(res, res.parallel_current_squared, ik) + _profile_total(res, res.pressure_curvature, ik) + push!(lines, @sprintf(" |b_eff|²/μ₀ = %+.4e K|ξ_n|² = %+.4e δW_p = %+.4e", eb, kk, res.dW_plasma[ik])) + push!( + lines, + @sprintf(" shear = %+.4e parallel = %+.4e curvature = %+.4e", + _profile_total(res, res.shear_current, ik), _profile_total(res, res.parallel_current_squared, ik), _profile_total(res, res.pressure_curvature, ik)) + ) + end + if !isempty(res.dW_plasma_standard_form) + push!( + lines, + @sprintf(" standard form: |b|²/μ₀ = %+.4e current coupling = %+.4e pressure = %+.4e δW_p = %+.4e", + _profile_total(res, res.b_squared, ik), real(_profile_total(res, res.current_coupling, ik)), + real(_profile_total(res, res.pressure_compression, ik)), real(res.dW_plasma_standard_form[ik])) + ) + end + ref_line = @sprintf(" ForceFreeStates ep[%d] = %+.4e", k, res.dW_plasma_reference[ik]) + isempty(res.dW_plasma) || (ref_line *= @sprintf(" relative difference = %.2e", res.dW_plasma_relative_error[ik])) + isempty(res.dW_plasma_standard_form) || (ref_line *= @sprintf(" (standard form %.2e)", res.dW_plasma_standard_form_relative_error[ik])) + push!(lines, ref_line) + if !isempty(res.effective_b_curl_residual_max) + push!( + lines, + @sprintf(" (∇×b_eff)·∇ψ residual: max %.2e, rms %.2e, relative %.2e", + res.effective_b_curl_residual_max[ik], res.effective_b_curl_residual_rms[ik], res.effective_b_curl_residual_relative[ik]) + ) + end + @info join(lines, "\n") + return nothing +end + +# Datasets of the group, in write order; names are the struct fields. +const _ED_MAIN_FIELDS = (:psi, :eigenmode_index, :effective_b_squared, :effective_b_squared_psi, :effective_b_squared_theta, + :effective_b_squared_zeta, :shear_current, :parallel_current_squared, :pressure_curvature, :b_squared, :current_coupling, + :pressure_compression, :dW_cumulative, :dW_cumulative_standard_form, :dW_plasma, :dW_plasma_standard_form, :dW_plasma_reference, + :dW_plasma_relative_error, :dW_plasma_standard_form_relative_error, :effective_b_curl_residual_max, :effective_b_curl_residual_rms, + :effective_b_curl_residual_relative) +const _ED_DENSITY_FIELDS = (map(_density_field, _ED_EFFECTIVE_PROFILES)..., :effective_b_curl_residual_density, map(_density_field, _ED_STANDARD_PROFILES)...) + +const _ED_NORM = "power-normalized like eigenmode_plasma_energies (per unit ⟨|ξ|²⟩, scaled by 2μ₀/ψ₀²)" +const _ED_TOTAL = (; dims=("eigenmode",), attach=(1 => "eigenmode_index",)) +const _ED_PROFILE = (; dims=("psi", "eigenmode"), attach=(1 => "psi", 2 => "eigenmode_index")) +const _ED_DENSITY = (; dims=("psi", "theta", "eigenmode"), attach=(1 => "psi", 2 => "Densities/theta", 3 => "eigenmode_index")) + +# Metadata table for PerturbedEquilibrium/EnergyDecomposition/ (paths relative to the group). +const ENERGY_DECOMPOSITION_H5_ANNOTATIONS = [ + "psi" => (; long_name="normalized poloidal flux ψ_N of the decomposition nodes (equilibrium knots plus the solution edge)", scale="psi"), + "eigenmode_index" => (; long_name="free-boundary eigenmode index (1 = least stable) of each decomposed mode, the eigenmode coordinate", scale="eigenmode"), + "effective_b_squared" => (; long_name="∮J|b_eff|²/μ₀ dθ with b_eff = b + ξ_n μ₀ j×n̂, $_ED_NORM", _ED_PROFILE...), + "effective_b_squared_psi" => (; long_name="ψ-component contribution Re[(J b_eff^ψ)* b_eff,ψ]/(μ₀J) to effective_b_squared", _ED_PROFILE...), + "effective_b_squared_theta" => (; long_name="θ-component contribution Re[(J b_eff^θ)* b_eff,θ]/(μ₀J) to effective_b_squared", _ED_PROFILE...), + "effective_b_squared_zeta" => (; long_name="ζ-component contribution Re[(J b_eff^ζ)* b_eff,ζ]/(μ₀J) to effective_b_squared", _ED_PROFILE...), + "shear_current" => (; long_name="∮J K1|ξ_n|² dθ with K1 = |∇ψ|²σS (shear-current coupling), $_ED_NORM", _ED_PROFILE...), + "parallel_current_squared" => (; long_name="∮J K2|ξ_n|² dθ with K2 = μ₀B²σ² (parallel current squared), $_ED_NORM", _ED_PROFILE...), + "pressure_curvature" => (; long_name="∮J K3|ξ_n|² dθ with K3 = 2p'κ_ψ (pressure-curvature drive), $_ED_NORM", _ED_PROFILE...), + "b_squared" => (; long_name="standard form: ∮J|b|²/μ₀ dθ, $_ED_NORM", _ED_PROFILE...), + "current_coupling" => (; long_name="standard form: ∮J j·(b×ξ*) dθ (complex per surface; its ψ-integral is real to roundoff), $_ED_NORM", _ED_PROFILE...), + "pressure_compression" => (; long_name="standard form: ∮J (∇·ξ)*(ξ·∇p) dθ (complex per surface; its ψ-integral is real to roundoff), $_ED_NORM", _ED_PROFILE...), + "dW_cumulative" => (; long_name="running ½∫dψ [effective_b_squared − shear_current − parallel_current_squared − pressure_curvature]", _ED_PROFILE...), + "dW_cumulative_standard_form" => (; long_name="running ½∫dψ [b_squared − current_coupling + pressure_compression]", _ED_PROFILE...), + "dW_plasma" => (; long_name="plasma potential energy δW_p in the effective-field form, $_ED_NORM", _ED_TOTAL...), + "dW_plasma_standard_form" => (; long_name="plasma potential energy δW_p in the standard form (imaginary part is roundoff), $_ED_NORM", _ED_TOTAL...), + "dW_plasma_reference" => (; long_name="eigenmode_plasma_energies entry of the same eigenmode, the reference for both forms", _ED_TOTAL...), + "dW_plasma_relative_error" => (; long_name="(dW_plasma − dW_plasma_reference)/|dW_plasma_reference|", _ED_TOTAL...), + "dW_plasma_standard_form_relative_error" => (; long_name="(Re dW_plasma_standard_form − dW_plasma_reference)/|dW_plasma_reference|", _ED_TOTAL...), + "effective_b_curl_residual_max" => + (; long_name="max over the equilibrium knots and θ of |(∇×b_eff)·∇ψ|, which vanishes for a consistent effective field", _ED_TOTAL...), + "effective_b_curl_residual_rms" => (; long_name="rms over the equilibrium knots and θ of |(∇×b_eff)·∇ψ|", _ED_TOTAL...), + "effective_b_curl_residual_relative" => (; long_name="effective_b_curl_residual_max divided by max|∂_θ b_eff,ζ/J|, the size of the cancelling terms", _ED_TOTAL...), + "Densities/theta" => (; long_name="poloidal angle θ ∈ [0,1) of the density grid", scale="theta"), + "Densities/R" => (; long_name="major radius of the (ψ,θ) grid points", units="m", dims=("psi", "theta"), attach=(1 => "psi", 2 => "Densities/theta")), + "Densities/Z" => (; long_name="vertical position of the (ψ,θ) grid points", units="m", dims=("psi", "theta"), attach=(1 => "psi", 2 => "Densities/theta")), + "Densities/effective_b_squared_density" => (; long_name="|b_eff|²/μ₀ per point (profile = Σ density·J/mtheta), $_ED_NORM", _ED_DENSITY...), + "Densities/effective_b_squared_psi_density" => (; long_name="ψ-component of |b_eff|²/μ₀ per point", _ED_DENSITY...), + "Densities/effective_b_squared_theta_density" => (; long_name="θ-component of |b_eff|²/μ₀ per point", _ED_DENSITY...), + "Densities/effective_b_squared_zeta_density" => (; long_name="ζ-component of |b_eff|²/μ₀ per point", _ED_DENSITY...), + "Densities/shear_current_density" => (; long_name="K1|ξ_n|² per point, K1 = |∇ψ|²σS", _ED_DENSITY...), + "Densities/parallel_current_squared_density" => (; long_name="K2|ξ_n|² per point, K2 = μ₀B²σ²", _ED_DENSITY...), + "Densities/pressure_curvature_density" => (; long_name="K3|ξ_n|² per point, K3 = 2p'κ_ψ", _ED_DENSITY...), + "Densities/effective_b_curl_residual_density" => (; long_name="(∇×b_eff)·∇ψ per point", _ED_DENSITY...), + "Densities/b_squared_density" => (; long_name="standard form: |b|²/μ₀ per point", _ED_DENSITY...), + "Densities/current_coupling_density" => (; long_name="standard form: j·(b×ξ*) per point", _ED_DENSITY...), + "Densities/pressure_compression_density" => (; long_name="standard form: (∇·ξ)*(ξ·∇p) per point", _ED_DENSITY...) +] + +""" + write_energy_decomposition!(h5, res) -> nothing + +Write `res` to `PerturbedEquilibrium/EnergyDecomposition/` of the open file `h5`, replacing the group +if it exists, and apply the metadata table. Fields of a disabled form are written zero-extent; +`Densities/` appears only when the densities were requested. +""" +function write_energy_decomposition!(h5::HDF5.File, res::EnergyDecompositionResult) + pe_group = haskey(h5, "PerturbedEquilibrium") ? h5["PerturbedEquilibrium"] : create_group(h5, "PerturbedEquilibrium") + haskey(pe_group, "EnergyDecomposition") && delete_object(pe_group, "EnergyDecomposition") + g = create_group(pe_group, "EnergyDecomposition") + # Disabled forms hold zero-extent arrays of the documented rank, the file's not-computed sentinel. + for name in _ED_MAIN_FIELDS + g[String(name)] = getfield(res, name) + end + if !isempty(res.theta) + d = create_group(g, "Densities") + d["theta"] = res.theta + d["R"] = res.R + d["Z"] = res.Z + for name in _ED_DENSITY_FIELDS + d[String(name)] = getfield(res, name) + end + end + Utilities.HDF5Annotations.annotate!(g, ENERGY_DECOMPOSITION_H5_ANNOTATIONS) + return nothing +end diff --git a/src/PerturbedEquilibrium/FieldReconstruction.jl b/src/PerturbedEquilibrium/FieldReconstruction.jl index 56d74b012..6ef432ea6 100644 --- a/src/PerturbedEquilibrium/FieldReconstruction.jl +++ b/src/PerturbedEquilibrium/FieldReconstruction.jl @@ -125,7 +125,7 @@ function reconstruct_physical_fields( # Compute contravariant displacement via Jacobian convolution (matches Fortran gpeq_contra) xwp_modes, xwt_modes, xwz_modes, xmt_modes, xmz_modes = compute_contra_displacements( xi_psi_modes, clebsch_psi1, clebsch_alpha, - psi_grid, equil, ffs, metric, ctrl + psi_grid, equil, ffs, metric; reg_spot=ctrl.reg_spot ) # Compute covariant components via metric tensor contraction (matches Fortran gpeq_cova) @@ -483,7 +483,7 @@ end """ compute_contra_displacements( xi_psi_modes, clebsch_psi1, clebsch_alpha, - psi_grid, equil, ffs, metric, ctrl + psi_grid, equil, ffs, metric; reg_spot ) -> (xwp_modes, xwt_modes, xwz_modes, xmt_modes, xmz_modes) Compute contravariant displacement via Jacobian mode coupling convolution. @@ -496,7 +496,8 @@ Matches Fortran `gpeq_contra`: Uses regularized xmp1 and xms for the theta/zeta components, then optionally applies additional regularization to get xmt/xmz. -Returns both unregularized (xwt, xwz) and regularized (xmt, xmz) versions. +Returns both unregularized (xwt, xwz) and regularized (xmt, xmz) versions. `reg_spot` is the +regularization width applied to xmt/xmz (0 disables it). """ function compute_contra_displacements( xi_psi_modes::Matrix{ComplexF64}, @@ -505,13 +506,12 @@ function compute_contra_displacements( psi_grid::Vector{Float64}, equil::Equilibrium.PlasmaEquilibrium, ffs::ForceFreeStatesResult, - metric::MetricData, - ctrl::PerturbedEquilibriumControl + metric::MetricData; + reg_spot::Float64 ) npsi, numpert_total = size(xi_psi_modes) (; mpert, mlow, nlow, npert) = ffs chi1 = 2π * equil.psio - reg_spot = ctrl.reg_spot fc = metric.fourier_coeffs xwp_modes = zeros(ComplexF64, npsi, numpert_total) diff --git a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl index 1151e5425..dc60986f7 100644 --- a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl +++ b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl @@ -23,6 +23,7 @@ import DelimitedFiles: readdlm include("PerturbedEquilibriumStructs.jl") include("ResponseMatrices.jl") include("FieldReconstruction.jl") +include("EnergyDecomposition.jl") include("Response.jl") include("SingularCoupling.jl") include("ResonantCoupling.jl") @@ -33,6 +34,8 @@ export PerturbedEquilibriumControl export PerturbedEquilibriumInternal export PerturbedEquilibriumState export ForcingMode +export EnergyDecompositionControl, EnergyDecompositionResult +export decompose_energy, write_energy_decomposition! # Export main functions export compute_perturbed_equilibrium diff --git a/test/runtests.jl b/test/runtests.jl index e18762e18..1df433c8b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -60,6 +60,7 @@ else include("./runtests_kinetic.jl") include("./runtests_multiion.jl") include("./runtests_fullruns.jl") + include("./runtests_energy_decomposition.jl") include("./runtests_coils.jl") include("./runtests_imas.jl") include("./runtests_rerun_from_h5.jl") diff --git a/test/runtests_energy_decomposition.jl b/test/runtests_energy_decomposition.jl new file mode 100644 index 000000000..ba1c2b039 --- /dev/null +++ b/test/runtests_energy_decomposition.jl @@ -0,0 +1,284 @@ +using Test +using HDF5 +using GeneralizedPerturbedEquilibrium +using GeneralizedPerturbedEquilibrium.PerturbedEquilibrium +using GeneralizedPerturbedEquilibrium.ForceFreeStates +using GeneralizedPerturbedEquilibrium.Equilibrium + +include("h5_metadata_check.jl") + +# Solovev forward fixture (mpsi = 16, vac_flag = true) with the [EnergyDecomposition] section +# appended, run once in a scratch directory and shared by every testset below. +const _ED_FIXTURE = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example") +const _ED_SECTION = """ + +[EnergyDecomposition] +eigenmodes = [1] # Free-boundary eigenmodes to decompose (1 = least stable) +effective_field_form = true # Decompose δW_p in the effective-field form +standard_form = true # Decompose δW_p in the standard form +write_densities = true # Store the per-(ψ,θ) densities so their shapes can be checked +""" + +function _run_energy_fixture() + dir = mktempdir() + # The fixture keeps HDF5 output off for speed; the stage tests need the file. + deck = replace(read(joinpath(_ED_FIXTURE, "gpec.toml"), String), "write_outputs_to_HDF5 = false" => "write_outputs_to_HDF5 = true") + write(joinpath(dir, "gpec.toml"), deck * _ED_SECTION) + return dir, GeneralizedPerturbedEquilibrium.main([dir]) +end + +# Copy of a force-free-states result with some fields replaced. +function _with(ffs, changes::Pair{Symbol}...) + names = fieldnames(typeof(ffs)) + vals = Any[getfield(ffs, f) for f in names] + for (f, v) in changes + vals[findfirst(==(f), names)] = v + end + return ForceFreeStates.ForceFreeStatesResult(vals...) +end + +const _ED_DIR, _ED_RES = _run_energy_fixture() +const _ED_FFS = _ED_RES.ffs +const _ED_BSQ = Equilibrium.metric_bsq_interpolant(_ED_FFS.equil) + +# Decompositions shared by the testsets below. +const _ED_FULL = PerturbedEquilibrium.decompose_energy(_ED_FFS; eigenmodes=[1], write_densities=true, verbose=false) +const _ED_PLAIN = PerturbedEquilibrium.decompose_energy(_ED_FFS; verbose=false) +const _ED_ONLY_EFF = PerturbedEquilibrium.decompose_energy(_ED_FFS; standard_form=false, verbose=false) +const _ED_ONLY_STD = PerturbedEquilibrium.decompose_energy(_ED_FFS; effective_field_form=false, verbose=false) +const _ED_TWO = PerturbedEquilibrium.decompose_energy(_ED_FFS; eigenmodes=[1, 2], verbose=false) + +@testset "EnergyDecomposition: surface geometry" begin + equil = _ED_FFS.equil + mtheta = length(equil.rzphi_ys) - 1 + thetas = [(k - 1) / mtheta for k in 1:mtheta] + geom = PerturbedEquilibrium.SurfaceGeometry(mtheta) + psi = 0.5 + PerturbedEquilibrium.surface_geometry!(geom, equil, psi, thetas; bsq=_ED_BSQ) + hint = (Ref(1), Ref(1)) + for k in 1:mtheta + m = Equilibrium.flux_surface_metric(equil, psi, thetas[k]; hint=hint) + @test geom.jac[k] ≈ m.jac rtol = 1e-12 + @test sqrt(geom.delpsi2[k]) ≈ m.delpsi rtol = 1e-12 + @test hypot(geom.R[k] - equil.ro, geom.Z[k] - equil.zo)^2 ≈ equil.rzphi_rsquared((psi, thetas[k]); hint=hint) rtol = 1e-10 + end + @test all(geom.bsq .> 0) + @test all(geom.g22 .> 0) && all(geom.g33 .> 0) + @test all(geom.g22 .* geom.g33 .- geom.g23 .^ 2 .> 0) + @test all(isfinite, geom.jac_psi) && all(isfinite, geom.bsq_psi) && all(isfinite, geom.bsq_theta) + # B² is required: the curvature kernel divides by it. + @test_throws UndefKeywordError PerturbedEquilibrium.surface_geometry!(geom, equil, psi, thetas) + # On an equilibrium knot the stored knot values reproduce the interpolant evaluation. + knot = length(equil.rzphi_xs) ÷ 2 + on_knot = PerturbedEquilibrium.SurfaceGeometry(mtheta) + PerturbedEquilibrium.surface_geometry!(on_knot, equil, equil.rzphi_xs[knot], thetas; bsq=_ED_BSQ, knot) + PerturbedEquilibrium.surface_geometry!(geom, equil, equil.rzphi_xs[knot], thetas; bsq=_ED_BSQ) + for name in fieldnames(PerturbedEquilibrium.SurfaceGeometry) + @test getfield(on_knot, name) ≈ getfield(geom, name) rtol = 1e-10 + end +end + +@testset "EnergyDecomposition: curvature kernel" begin + equil = _ED_FFS.equil + thetas_ext = collect(equil.rzphi_ys) + mtheta = length(thetas_ext) - 1 + geom = PerturbedEquilibrium.SurfaceGeometry(mtheta) + kern = PerturbedEquilibrium.CurvatureKernel(mtheta) + psi = 0.5 + PerturbedEquilibrium.surface_geometry!(geom, equil, psi, thetas_ext[1:mtheta]; bsq=_ED_BSQ) + PerturbedEquilibrium.curvature_kernel!(kern, geom, equil, psi, thetas_ext) + @test all(isfinite, kern.shear) && all(isfinite, kern.curvature) && all(isfinite, kern.sigma) + # The Solovev equilibrium carries a finite parallel current: σ must not vanish identically. + @test maximum(abs, kern.sigma) > 0 + # Flux-surface average of the shear geometry term vanishes, so ⟨J S⟩/χ'² reduces to q' ⟨1⟩. + chi1 = 2π * equil.psio + q1 = equil.profiles.q_deriv(psi) + @test sum(kern.shear .* geom.jac) / mtheta ≈ chi1^2 * q1 rtol = 1e-6 +end + +@testset "EnergyDecomposition: solution state between saved steps" begin + ffs = _ED_FFS + sol = ffs.solution + # A saved step comes back as stored. + @test ForceFreeStates.solution_state_at(ffs, sol.psi_store[5]) == sol.u_store[:, :, :, 5] + # Just left of a mid-interval point is integrated forward from the left step, just right of it + # backward from the right step; the two must meet. + sing = [s.psifac for s in ffs.surfaces] + i = 9 + findfirst(j -> !any(s -> sol.psi_store[j] < s < sol.psi_store[j+1], sing), 10:sol.step-1) + a, b = sol.psi_store[i], sol.psi_store[i+1] + delta = 1e-9 * (b - a) + from_left = ForceFreeStates.solution_state_at(ffs, (a + b) / 2 - delta) + from_right = ForceFreeStates.solution_state_at(ffs, (a + b) / 2 + delta) + @test maximum(abs, from_left - from_right) / maximum(abs, from_left) < 1e-6 +end + +@testset "EnergyDecomposition: surface densities" begin + ffs = _ED_FFS + equil = ffs.equil + psi_grid, _ = PerturbedEquilibrium.decomposition_grid(ffs) + modes = PerturbedEquilibrium.eigenmode_modes(ffs, 1, psi_grid, PerturbedEquilibrium.node_solution(ffs, psi_grid)...) + @test size(modes.Jxi_psi) == (length(psi_grid), ffs.mpert) + thetas_ext = collect(equil.rzphi_ys) + mtheta = length(thetas_ext) - 1 + mvals = collect(ffs.mlow:ffs.mhigh) + ft = GeneralizedPerturbedEquilibrium.Utilities.FourierTransforms.FourierTransform(mtheta, ffs.mpert, ffs.mlow) + geom = PerturbedEquilibrium.SurfaceGeometry(mtheta) + kern = PerturbedEquilibrium.CurvatureKernel(mtheta) + sf = PerturbedEquilibrium.SurfaceFields(mtheta) + eff = PerturbedEquilibrium.EffectiveField(mtheta) + den = PerturbedEquilibrium.SurfaceDensities(mtheta) + work = zeros(ComplexF64, ffs.mpert) + ipsi = length(psi_grid) ÷ 2 + psi = psi_grid[ipsi] + PerturbedEquilibrium.surface_geometry!(geom, equil, psi, thetas_ext[1:mtheta]; bsq=_ED_BSQ) + PerturbedEquilibrium.curvature_kernel!(kern, geom, equil, psi, thetas_ext) + PerturbedEquilibrium.surface_fields!(sf, ft, modes, ipsi, mvals, geom, work) + PerturbedEquilibrium.effective_field!(eff, geom, sf) + PerturbedEquilibrium.surface_densities!(den, geom, kern, sf, eff, equil, psi, ffs.nlow; effective_field_form=true, standard_form=true) + @test den.effective_b_squared ≈ den.effective_b_squared_psi .+ den.effective_b_squared_theta .+ den.effective_b_squared_zeta + @test all(den.effective_b_squared .>= 0) + @test all(den.b_squared .>= 0) + @test all(den.parallel_current_squared .>= 0) + @test PerturbedEquilibrium.surface_integral(den.effective_b_squared, geom.jac) ≈ sum(den.effective_b_squared .* geom.jac) / mtheta + # The standard-form current term is real up to roundoff at the surface-integral level. + cc = PerturbedEquilibrium.surface_integral(den.current_coupling, geom.jac) + @test abs(imag(cc)) < 1e-8 * max(abs(real(cc)), 1e-300) + # Disabling a form leaves its densities untouched (zero). + fill!(den.b_squared, 0.0) + PerturbedEquilibrium.surface_densities!(den, geom, kern, sf, eff, equil, psi, ffs.nlow; effective_field_form=true, standard_form=false) + @test all(iszero, den.b_squared) + # The profile name tuples cover every density. + @test Set((PerturbedEquilibrium._ED_EFFECTIVE_PROFILES..., PerturbedEquilibrium._ED_STANDARD_PROFILES...)) == + Set(fieldnames(PerturbedEquilibrium.SurfaceDensities)) +end + +@testset "EnergyDecomposition: decompose_energy totals, checks, and toggles" begin + ffs = _ED_FFS + res = _ED_FULL + npsi = length(res.psi) + mtheta = length(ffs.equil.rzphi_ys) - 1 + @test res.eigenmode_index == [1] + @test res.psi == first(PerturbedEquilibrium.decomposition_grid(ffs)) + @test all(diff(res.psi) .> 0) + @test res.psi[end] == ffs.solution.psi_store[ffs.solution.step] + @test res.effective_b_curl_residual_relative[1] < 1e-8 + @test res.effective_b_curl_residual_rms[1] <= res.effective_b_curl_residual_max[1] + @test res.dW_plasma_reference[1] ≈ real(ffs.free_boundary.ep[1]) + # Tolerances cover the coarse mpsi = 16 fixture; the docs list the observed mismatch and its convergence. + @test abs(res.dW_plasma_relative_error[1]) < 0.03 + @test abs(res.dW_plasma_standard_form_relative_error[1]) < 0.03 + @test res.dW_cumulative[1, 1] == 0.0 + @test res.dW_cumulative[end, 1] ≈ res.dW_plasma[1] + @test res.dW_cumulative_standard_form[end, 1] ≈ res.dW_plasma_standard_form[1] + @test res.effective_b_squared ≈ res.effective_b_squared_psi .+ res.effective_b_squared_theta .+ res.effective_b_squared_zeta + @test size(res.effective_b_squared_density) == (npsi, mtheta, 1) + @test size(res.current_coupling_density) == (npsi, mtheta, 1) + @test size(res.R) == (npsi, mtheta) + @test res.theta == [(k - 1) / mtheta for k in 1:mtheta] + ipsi = npsi ÷ 2 + jac = [ffs.equil.rzphi_jac((res.psi[ipsi], t)) for t in res.theta] + @test sum(res.effective_b_squared_density[ipsi, :, 1] .* jac) / mtheta ≈ res.effective_b_squared[ipsi, 1] rtol = 1e-10 + @test sum(res.b_squared_density[ipsi, :, 1] .* jac) / mtheta ≈ res.b_squared[ipsi, 1] rtol = 1e-10 + + plain = _ED_PLAIN + @test isempty(plain.theta) && isempty(plain.R) && isempty(plain.effective_b_squared_density) + @test plain.dW_plasma == res.dW_plasma + + only_eff = _ED_ONLY_EFF + @test isempty(only_eff.b_squared) && isempty(only_eff.dW_plasma_standard_form) && isempty(only_eff.dW_cumulative_standard_form) + @test only_eff.dW_plasma == res.dW_plasma + + only_std = _ED_ONLY_STD + @test isempty(only_std.dW_plasma) && isempty(only_std.effective_b_squared) && isempty(only_std.effective_b_curl_residual_max) + @test only_std.dW_plasma_standard_form == res.dW_plasma_standard_form + + two = _ED_TWO + @test size(two.effective_b_squared) == (npsi, 2) + @test two.dW_plasma[1] == res.dW_plasma[1] + @test two.dW_plasma_reference[2] ≈ real(ffs.free_boundary.ep[2]) + @test abs(two.dW_plasma_relative_error[2]) < 0.07 + + @test_throws ArgumentError PerturbedEquilibrium.decompose_energy(ffs; effective_field_form=false, standard_form=false) + @test_throws ArgumentError PerturbedEquilibrium.decompose_energy(ffs; eigenmodes=[0]) + @test_throws ArgumentError PerturbedEquilibrium.decompose_energy(ffs; eigenmodes=[size(ffs.free_boundary.wt, 2) + 1]) +end + +@testset "EnergyDecomposition: HDF5 writer" begin + res = _ED_FULL + path = joinpath(mktempdir(), "energy.h5") + h5open(path, "w") do h5 + PerturbedEquilibrium.write_energy_decomposition!(h5, res) + # A second write replaces the group instead of failing on existing datasets. + PerturbedEquilibrium.write_energy_decomposition!(h5, res) + end + h5open(path, "r") do h5 + g = h5["PerturbedEquilibrium/EnergyDecomposition"] + for name in ("psi", "eigenmode_index", "effective_b_squared", "shear_current", "b_squared", "current_coupling", + "dW_cumulative", "dW_cumulative_standard_form", "dW_plasma", "dW_plasma_standard_form", "dW_plasma_reference", + "dW_plasma_relative_error", "dW_plasma_standard_form_relative_error", "effective_b_curl_residual_max", + "effective_b_curl_residual_rms", "effective_b_curl_residual_relative", "Densities/theta", "Densities/R", + "Densities/effective_b_squared_density", "Densities/current_coupling_density") + @test haskey(g, name) + end + @test read(g["dW_plasma"]) == res.dW_plasma + @test read(g["current_coupling"]) == res.current_coupling + @test isempty(_collect_metadata_violations(g)) + @test HDF5.API.h5ds_is_scale(g["psi"]) + @test HDF5.API.h5ds_is_attached(g["effective_b_squared"], g["psi"], 1) + @test HDF5.API.h5ds_is_attached(g["Densities/effective_b_squared_density"], g["Densities/theta"], 1) + # eigenmode_index is the coordinate of the eigenmode axis (C index 0 is the last Julia axis). + @test HDF5.API.h5ds_is_scale(g["eigenmode_index"]) + @test HDF5.API.h5ds_is_attached(g["dW_plasma"], g["eigenmode_index"], 0) + @test HDF5.API.h5ds_is_attached(g["effective_b_squared"], g["eigenmode_index"], 0) + @test HDF5.API.h5ds_is_attached(g["Densities/effective_b_squared_density"], g["eigenmode_index"], 0) + @test attrs(g["Densities/R"])["units"] == "m" + end + # A disabled form writes zero-extent datasets, never missing ones. + path2 = joinpath(mktempdir(), "energy2.h5") + h5open(path2, "w") do h5 + PerturbedEquilibrium.write_energy_decomposition!(h5, _ED_ONLY_EFF) + end + h5open(path2, "r") do h5 + g = h5["PerturbedEquilibrium/EnergyDecomposition"] + @test length(read(g["dW_plasma_standard_form"])) == 0 + @test length(read(g["b_squared"])) == 0 + # Zero-extent sentinels keep the documented rank so the dimension metadata stays consistent. + @test size(read(g["b_squared"])) == (0, 0) + @test size(read(g["dW_cumulative_standard_form"])) == (0, 0) + @test !haskey(g, "Densities") + end +end + +@testset "EnergyDecomposition: main stage from the TOML section" begin + # The fixture deck carries [EnergyDecomposition] with write_densities = true. + # The stage adds nothing to main's return; its result lives in gpec.h5 only. + @test keys(_ED_RES)[1:3] == (:ffs, :pe, :slayer) + @test !any(k -> occursin("energy", String(k)), keys(_ED_RES)) + h5open(joinpath(_ED_DIR, "gpec.h5"), "r") do h5 + g = h5["PerturbedEquilibrium/EnergyDecomposition"] + @test read(g["eigenmode_index"]) == [1] + @test length(read(g["dW_plasma"])) == 1 + @test haskey(g, "Densities/R") + @test isempty(_collect_metadata_violations(g)) + end + # No section, no stage. + @test GeneralizedPerturbedEquilibrium.run_energy_decomposition(_ED_FFS, Dict{String,Any}()) === nothing + # Runs the decomposition cannot handle are skipped with a warning instead of aborting main. + inputs = Dict{String,Any}("EnergyDecomposition" => Dict{String,Any}("eigenmodes" => [1])) + for unsupported in (_with(_ED_FFS, :nhigh => _ED_FFS.nhigh + 1), _with(_ED_FFS, :integrator => :riccati)) + @test_logs (:warn, r"Skipping energy decomposition") match_mode = :any begin + @test GeneralizedPerturbedEquilibrium.run_energy_decomposition(unsupported, inputs) === nothing + end + end + # The scripting entry writes into the run's gpec.h5 (here a scratch copy) and returns the result. + scratch = mktempdir() + cp(joinpath(_ED_DIR, "gpec.h5"), joinpath(scratch, "gpec.h5")) + scripted = GeneralizedPerturbedEquilibrium.energy_decomposition(_with(_ED_FFS, :dir_path => scratch); eigenmodes=[2], standard_form=false, verbose=false) + @test scripted.eigenmode_index == [2] + h5open(joinpath(scratch, "gpec.h5"), "r") do h5 + @test read(h5["PerturbedEquilibrium/EnergyDecomposition/eigenmode_index"]) == [2] + end + h5open(joinpath(_ED_DIR, "gpec.h5"), "r") do h5 + @test read(h5["PerturbedEquilibrium/EnergyDecomposition/eigenmode_index"]) == [1] + end +end