ForceFreeStates - FEATURE - Add perpendicular inertia and an ELITE-style growth-rate solve

docs/development/references.md

@@ -34,6 +34,25 @@ The ForceFreeStates module (ideal MHD stability analysis) implements methods fro
   - Published: Physics of Plasmas **25**, 032507 (2018)
   - Describes: Riccati method for ideal MHD eigenvalue problem
 
+- **Arnold (1967)**: "On a characteristic class entering into conditions of quantization"
+  - Published: Functional Analysis and Its Applications **1**, 1 (1967)
+  - Describes: the Maslov index of a path of Lagrangian subspaces, counted through the eigenvalues of a unitary (Cayley) image; basis of the exact conjugate-point count `nzero_conj` in `FixedBoundaryStability.jl`
+
+- **Robbin & Salamon (1993)**: "The Maslov index for paths"
+  - Published: Topology **32**, 827 (1993)
+  - Describes: the Maslov index as a signed count of crossings of a Lagrangian path with a reference subspace, which is what the Cayley eigenphases crossing π compute step by step
+
+- **Wilson, Snyder, Huysmans & Miller (2002)**: "Numerical studies of edge localized instabilities in tokamaks"
+  - Published: Physics of Plasmas **9**, 1277 (2002)
+  - Describes: the ELITE formulation; Eq. 9 is the perpendicular (incompressible) inertia that `PerpendicularInertia` implements; the growth rate is found by shooting with a secant iteration on an edge matching eigenvalue
+
+- **Snyder, Wilson et al. (2002)**: "Edge localized modes and the pedestal: A model based on coupled peeling–ballooning modes"
+  - Published: Physics of Plasmas **9**, 2037 (2002)
+  - Describes: ELITE's peeling–ballooning application and the γ > ω_*pi/2 diamagnetic-stabilization test that `IdealGrowthRateResult.omega_star_pi` serves
+
+- **Lunniss (2016)**: "Arbitrary-n ELITE" (PhD thesis, University of York; etheses.whiterose.ac.uk/16638)
+  - Describes: ELITE's radial finite-difference sweep (App. A.1), the edge matching matrix T(γ²) and the fictitious eigenvalue λ(γ²) driven to zero (§2.5.8), and the arbitrary-n inertia options; the comparison point for the count-based bracket of `growth_rate`, which ELITE does not have
+
 ## PerturbedEquilibrium Module
 
 The PerturbedEquilibrium module implements GPEC-style perturbed equilibrium calculations from:

docs/src/api.md

@@ -21,20 +21,24 @@ foreach(r -> save_to_hdf5("gpec.h5", r), (pe, kf, ef, tear))
 ```
 
 `solve(eq, alg; kwargs...)` is sugar for the problem form `solve(ForceBalanceProblem(eq; kwargs...), alg)`: the problem names what is solved (the force-balance equation on this equilibrium, with this mode range, wall and closure) and the algorithm names how.
-`main` translates each deck section into the problem, algorithm or control it names, calls these functions for every section present, saves each result into `[Run] output_file`, and returns `(; equil, local_stability, ffs, match, resp, pe, kinetic, error_fields, tearing)`, each `nothing` when its section is absent.
+`main` translates each deck section into the problem, algorithm or control it names, calls these functions for every section present, saves each result into `[Run] output_file`, and returns `(; equil, local_stability, ffs, match, growth_rate, resp, pe, kinetic, error_fields, tearing)`, each `nothing` when its section is absent.
 
 | Stage | Deck section | Call | Result |
 |---|---|---|---|
 | Equilibrium | `[Equilibrium]` (or `[Solovev]`, `[LargeAspectRatio]`, `[TJAnalytic]`) | `PlasmaEquilibrium(file; kwargs...)`, `PlasmaEquilibrium(config)`, `PlasmaEquilibrium(dd)` | `PlasmaEquilibrium` |
 | ForceFreeStates | `[ForceFreeStates]` + `[Forward]`, `[Riccati]` or `[Galerkin]` | `solve(ForceBalanceProblem(eq; ...), Forward() \| Riccati() \| Galerkin())` | `ForceFreeStatesResult` |
 | LocalStability | `[LocalStability]` + `[Ballooning]` | `solve(LocalStabilityProblem(eq; ...), Ballooning(; ...))` | `LocalStabilityResult` |
 | Match | `[Match]` + `[GGJ]` | `solve(MatchProblem(ffs; ...), GGJModel())` | `ForceFreeStatesResult` with `closure == :matched` |
+| GrowthRate | `[GrowthRate]` | `growth_rate(ForceBalanceProblem(eq; ...), Forward(), GrowthRateConfig(; ...); marginal=ffs)` | `IdealGrowthRateResult` |
 | Drive | `[ForcingTerms]` | `RMPField(ft_ctrl)` (a `ForcingTermsControl`), `RMPField(file)`, `RMPField(modes)` | `RMPField` |
 | PerturbedEquilibrium | `[PerturbedEquilibrium]` | `plasma_response(ffs; ...)`, then `perturbed_equilibrium(resp, rmp; ...)` | `PlasmaResponse`, `PerturbedEquilibriumState` |
 | KineticForces | `[KineticForces]` | `kinetic_torque(pe, KineticForcesControl(; ...))` | `KineticForcesResult` |
 | ErrorFields | `[ErrorFields]` | `error_fields(resp, pe.coil_sets, CoilConfig(...); ...)` | `ErrorFieldsResult` |
 | Tearing | `[Tearing]` + `[SLAYER]` or `[GGJ]` | `solve(TearingProblem(ffs; ...), SLAYERModel() \| GGJModel(); config=TearingConfig(...))` | `TearingResult` |
 
+`[GrowthRate]` runs the ideal growth-rate solve after the force-free-states stage; its keys are the `GrowthRateConfig` fields, and its `rho`, left at NaN, derives the mass density from the equilibrium's kinetic profiles.
+It needs `[ForceFreeStates] vac_flag = true` and `integrator = "Forward"`, and cannot share a deck with `[Match]`, since its marginal solve is the plain ideal one.
+
 The [Workflow](workflow.md) page describes what each stage computes; the module pages document the keywords, which are the deck keys of the matching section.
 
 ## Inspecting a result

docs/src/stability.md

@@ -141,6 +141,15 @@ The perpendicular diffusivities that set the layer widths are described with the
 - Both integrators use the adaptive ninth-order `Vern9` solver at `eulerlagrange_tolerance`; the Riccati basis is recovered by renormalizing whenever the solution exceeds `renorm_threshold`.
 - The two crossing conventions differ slightly, so Riccati and Forward free-boundary energies agree to within a few percent, not to machine precision.
 
+### Counting fixed-boundary modes
+
+Two counts of the conjugate points (fixed-boundary modes fitting inside the ψ interval) come out of the forward integration.
+
+- `nzero` is DCON's criterion: a sign change of the smallest-magnitude eigenvalue of W⁻¹ = U₁U₂⁻¹ between stored steps, confirmed at the interpolated zero. It is cheap and is what the `crit` scan in the output records, but it misses a crossing whenever a stiff harmonic holds a tiny eigenvalue of W⁻¹ near zero at every stored step, and it can report crossings where the frame has lost its Hermitian structure.
+- `nzero_conj` is the Maslov index of the solution frame (Arnold 1967; Robbin & Salamon 1993), computed at every accepted step and exact for a Lagrangian frame. The Cayley transform Ψ = (U₁ − i s U₂)(U₁ + i s U₂)⁻¹ of the column-normalized frame is unitary when U₁†U₂ is Hermitian, has an eigenvalue −1 exactly where det U₁ = 0 (a conjugate point) and +1 where U₂ is singular (a pole of W⁻¹, harmless), and is invariant under U → U·T, so the Gaussian reductions do not disturb it. Each eigenphase is followed through its eigenvector from one step to the next and counts one conjugate point when its shortest arc crosses π. The phases are seeded one accepted step into each chunk, because at the axis start and after an ideal crossing U₁ has vanishing columns whose phases sit exactly at π and belong to the interval's end, not its interior; a step in which a phase rotates by more than π/4 is bisected on the integrator's dense output until the arcs are short, since an eigenvalue of W⁻¹ can pass zero within one step near a conjugate point. The scale s is set from the frame one step into the chunk and never changed within it: a change of s moves every phase but none through π, yet a large s squashes the phases toward π until a pole of W⁻¹ becomes indistinguishable from a zero crossing, so rescaling is not a substitute for bisection.
+
+`nconj_flags` counts the steps where the count cannot be trusted: a rotation above π/2 after bisection, or a unitarity defect max||λ|−1| of Ψ above 1e-3, which means the integrated frame is no longer Lagrangian to working precision (the Solovev regression fixture shows this after its q = 2 crossing, where `nzero` also misfires). `maslov_defect` records the largest defect seen. The count costs two matrix inversions and one eigendecomposition per accepted step, roughly 10–15 % of a solve, so it is off by default; `growth_rate` switches it on for all of its solves and uses it, with the sign of the least free-boundary energy, as the number of modes faster than the trial γ, which is what makes its bracket monotone. Checked on the TJ analytic equilibrium at `qa = 2.9, p0 = 0.8`, where a fixed-boundary mode at γ² ≈ 1.2e8 is invisible to `nzero`: `nzero_conj` is 1 below that eigenvalue and 0 above, and the Morse index `nzero_conj` + (negative free-boundary eigenvalues) is non-increasing in γ² and continuous across the pole, as it must be.
+
 Local stability (Mercier, resistive interchange and ballooning) is a separate stage on the same equilibrium; see [Local Stability](ballooning.md).
 
 ## API reference
@@ -149,5 +158,5 @@ The Galerkin solver's functions are documented on the [Galerkin Solver](galerkin
 
 ```@autodocs
 Modules = [GeneralizedPerturbedEquilibrium.ForceFreeStates]
-Pages = ["ForceFreeStates.jl", "CoreTypes.jl", "Surfaces/Types.jl", "Riccati/Types.jl", "Matching/DeltaPrime.jl", "Result.jl", "Surfaces/ResistEval.jl", "Matching/ResonantMatch.jl", "Matching/MatchProblem.jl", "EulerLagrange.jl", "Surfaces/Finding.jl", "Surfaces/Asymptotics.jl", "Fourfit.jl", "Kinetic.jl", "FixedKineticMatrices.jl", "FixedBoundaryStability.jl", "Utils.jl", "Free.jl", "Riccati/Propagators.jl", "Riccati/Crossings.jl", "Riccati/DeltaPrimeBVP.jl", "Riccati/Driver.jl"]
+Pages = ["ForceFreeStates.jl", "CoreTypes.jl", "Surfaces/Types.jl", "Riccati/Types.jl", "Matching/DeltaPrime.jl", "Result.jl", "Surfaces/ResistEval.jl", "Matching/ResonantMatch.jl", "Matching/MatchProblem.jl", "EulerLagrange.jl", "Surfaces/Finding.jl", "Surfaces/Asymptotics.jl", "Fourfit.jl", "Kinetic.jl", "Inertia.jl", "FixedKineticMatrices.jl", "FixedBoundaryStability.jl", "Utils.jl", "Free.jl", "Riccati/Propagators.jl", "Riccati/Crossings.jl", "Riccati/DeltaPrimeBVP.jl", "Riccati/Driver.jl"]
 ```

docs/src/workflow.md

@@ -21,6 +21,7 @@ Every stage after `Equilibrium` reads the force-free states, and `KineticForces`
 |---|---|---|
 | Is the equilibrium read and mapped correctly; what are q and β_N? | `[Equilibrium]` (+ `[ForceFreeStates]`, which triggers the write) | `Equilibrium/q_95`, `Equilibrium/beta_N`, `Equilibrium/Profiles/` |
 | Is the plasma ideal-kink stable, with or without a wall? | `[ForceFreeStates]` (+ `[Wall]`) | `ForceFreeStates/FreeBoundaryStability/W_freeboundary`, `…/eigenmode_energies` |
+| How fast does the least-stable free-boundary ideal mode grow? | `[GrowthRate]` (+ `[ForceFreeStates] vac_flag = true`, Forward) | `ForceFreeStates/GrowthRate/` |
 | Is it Mercier, resistive-interchange or ballooning unstable? | `[LocalStability]` (+ `[Ballooning]`) | `LocalStability/D_I`, `LocalStability/D_R`, `LocalStability/ballooning_Delta_prime` |
 | What is the tearing stability index Δ′ at each rational surface? | `[ForceFreeStates] integrator = "Riccati"` or `"Galerkin"` (+ `[Riccati]` or `[Galerkin]`) | `SingularSurfaces/Delta_prime_matrix` |
 | How does the plasma respond to an RMP or error field; do islands overlap? | `[ForcingTerms]` + `[PerturbedEquilibrium]` | `PerturbedEquilibrium/SingularCoupling/resonant_area_weighted_field`, `…/island_half_width`, `…/chirikov_parameter` |

regression-harness/cases/lar_reversed_shear_growth_rate_n1.toml

@@ -0,0 +1,79 @@
+# Regression case: ELITE-style ideal growth rate of the large-aspect-ratio reversed-shear n=1 example.
+# The example is free-boundary unstable at n=1, so [GrowthRate] root-finds a nonzero γ with a uniform
+# deuterium mass density (1e20 m⁻³). Tracks γ, γ², the root-find history and the ideal energy it starts from.
+# Each [quantities.*] block names an HDF5 path, how to extract it, and its noise floor.
+[case]
+name = "lar_reversed_shear_growth_rate_n1"
+description = "Large-aspect-ratio reversed-shear equilibrium, n=1, ELITE-style ideal growth rate with uniform deuterium density"
+example_dir = "examples/LAR_reversed_shear_example"
+
+[overrides]
+"GrowthRate.rho" = 3.3e-7                # Uniform mass density [kg/m³] (deuterium at 1e20 m⁻³); NaN derives ρ(ψ) from the kinetic profiles
+
+# Growth rate (γ ≈ 3.0e3 1/s); absolute floors match a 1e-10 shift in λ_min near the root
+[quantities.gamma]
+h5path = "ForceFreeStates/GrowthRate/gamma"
+type = "real_scalar"
+extract = "value"
+label = "growth rate gamma (1/s)"
+noise_threshold = 1e-5
+order = 10
+
+[quantities.gamma_sq]
+h5path = "ForceFreeStates/GrowthRate/gamma_sq"
+type = "real_scalar"
+extract = "value"
+label = "growth rate squared gamma^2 (1/s^2)"
+noise_threshold = 1e-1
+order = 11
+
+# Root-find history: every γ² evaluated, the least-stable free-boundary energy there, and the fixed-boundary zero count
+[quantities.history_gamma_sq]
+h5path = "ForceFreeStates/GrowthRate/History/gamma_sq"
+type = "real_vector"
+extract = "all_real"
+label = "root-find gamma^2 history"
+noise_threshold = 1e-1
+order = 20
+
+[quantities.history_lambda_min]
+h5path = "ForceFreeStates/GrowthRate/History/lambda_min"
+type = "real_vector"
+extract = "all_real"
+label = "root-find lambda_min history"
+noise_threshold = 1e-10
+order = 21
+
+[quantities.history_nzero]
+h5path = "ForceFreeStates/GrowthRate/History/nzero"
+type = "real_vector"
+extract = "all_real"
+label = "root-find nzero history (integer counts)"
+noise_threshold = 0
+order = 22
+
+# Ideal marginal solve the root find starts from
+[quantities.et_real]
+h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_energies"
+type = "complex_vector"
+extract = "real_first"
+label = "total energy Re(et[1])"
+noise_threshold = 1e-10
+order = 30
+
+[quantities.msing]
+h5path = "SingularSurfaces/rational_count"
+type = "int_scalar"
+extract = "value"
+label = "# singular surfaces"
+noise_threshold = 0
+order = 31
+
+# Runtime (special: not from H5)
+[quantities.runtime]
+h5path = ""
+type = "runtime"
+extract = "value"
+label = "Runtime (s)"
+noise_threshold = 0.0
+order = 999

src/Deck.jl

@@ -78,7 +78,7 @@ const _RETIRED_SECTIONS = Dict("DEBUG" => "[DEBUG] is retired: output_benchmark_
 
 # Every section a deck may hold, besides the analytic-equilibrium sections.
 const _DECK_SECTIONS = ("Run", "Equilibrium", "Wall", "LocalStability", "Ballooning", "ForceFreeStates", "Forward", "Riccati",
-    "Galerkin", "Match", "GGJ", "Tearing", "SLAYER", "ForcingTerms", "PerturbedEquilibrium", "KineticForces", "ErrorFields")
+    "Galerkin", "Match", "GGJ", "Tearing", "SLAYER", "GrowthRate", "ForcingTerms", "PerturbedEquilibrium", "KineticForces", "ErrorFields")
 const _INTEGRATORS = ("Forward" => Forward, "Riccati" => Riccati, "Galerkin" => Galerkin)
 
 # Stop the run when a parsed gpec.toml section sets keys that no longer exist, so the deck is updated.
@@ -377,6 +377,16 @@ function tearing_from_deck(inputs::AbstractDict)
     return (; problem, model=inner_model == "SLAYER" ? slayer_model(inputs) : ggj_model(inputs), config)
 end
 
+"""
+    growth_rate_from_deck(inputs) -> Union{Nothing,GrowthRateConfig}
+
+The `GrowthRateConfig` of `[GrowthRate]`; `nothing` without that section.
+"""
+function growth_rate_from_deck(inputs::AbstractDict)
+    haskey(inputs, "GrowthRate") || return nothing
+    return GrowthRateConfig(; section_keywords(inputs["GrowthRate"], "GrowthRate", _struct_fields(GrowthRateConfig))...)
+end
+
 """
     perturbed_equilibrium_from_deck(inputs) -> Union{Nothing,NamedTuple}
 
@@ -447,7 +457,7 @@ function parse_deck(inputs::AbstractDict, path::AbstractString)
     end
     has_ffs = haskey(inputs, "ForceFreeStates")
     if !has_ffs
-        for section in ("Wall", "Forward", "Riccati", "Galerkin", "Match", "GGJ", "Tearing", "SLAYER", "LocalStability", "Ballooning",
+        for section in ("Wall", "Forward", "Riccati", "Galerkin", "Match", "GGJ", "Tearing", "SLAYER", "GrowthRate", "LocalStability", "Ballooning",
             "ForcingTerms", "PerturbedEquilibrium", "KineticForces", "ErrorFields")
             haskey(inputs, section) && error("[$section] is set but there is no [ForceFreeStates] section, so the run stops after the equilibrium")
         end
@@ -462,6 +472,15 @@ function parse_deck(inputs::AbstractDict, path::AbstractString)
         error("[KineticForces] is set but neither [PerturbedEquilibrium] (the torque) nor [ForceFreeStates] kinetic = \"KineticForces\" uses it")
     (haskey(inputs, "KineticForces") || kinetic_source) && isempty(get(get(inputs, "Equilibrium", Dict{String,Any}()), "kinetic_file", "")) &&
         error("this run needs kinetic profiles ([KineticForces], or kinetic = \"KineticForces\"): set kinetic_file in [Equilibrium]")
+    growth_rate = growth_rate_from_deck(inputs)
+    # The growth-rate solve needs the plain ideal solve as its marginal one and the free-boundary energy of a Forward solve.
+    if growth_rate !== nothing
+        (; keywords, algorithm) = force_balance
+        haskey(inputs, "Match") && error("[GrowthRate] needs the ideal solve as its marginal one: remove the [Match] section")
+        get(keywords, :kinetic, nothing) === nothing || error("[GrowthRate] needs an ideal force-balance solve: remove kinetic from [ForceFreeStates]")
+        get(keywords, :vac_flag, false) || error("[GrowthRate] needs the free-boundary energy: set [ForceFreeStates] vac_flag = true")
+        algorithm isa Forward || error("[GrowthRate] needs integrator = \"Forward\" in [ForceFreeStates]")
+    end
     match = match_from_deck(inputs)
     # The MatchProblem guards, checked before the force-free-states solve rather than after it.
     if match !== nothing
@@ -483,7 +502,7 @@ function parse_deck(inputs::AbstractDict, path::AbstractString)
         error("[GGJ] is set but neither [Match] nor [Tearing] inner_model = \"GGJ\" uses it")
     haskey(inputs, "SLAYER") && tearing_model != "SLAYER" && !(has_ffs && get(inputs["ForceFreeStates"], "psilim_from_layer_overlap", false) === true) &&
         error("[SLAYER] is set but neither [Tearing] inner_model = \"SLAYER\" nor [ForceFreeStates] psilim_from_layer_overlap = true uses it")
-    return (; run=run_settings(inputs, path), force_balance, local_stability=local_stability_from_deck(inputs), match, tearing,
+    return (; run=run_settings(inputs, path), force_balance, local_stability=local_stability_from_deck(inputs), match, tearing, growth_rate,
         kinetic_forces=haskey(inputs, "KineticForces") ? kinetic_forces_control(inputs) : nothing, perturbed_equilibrium, forcing,
         error_fields=error_fields_from_deck(inputs, path))
 end
@@ -578,6 +597,7 @@ function deck_inputs(deck::NamedTuple, config::Equilibrium.EquilibriumConfig, an
         merge!(inputs, _tearing_sections(filter(p -> p.first != :critical_resonant_field, problem), get(problem, :critical_resonant_field, nothing),
             model, config))
     end
+    deck.growth_rate === nothing || (inputs["GrowthRate"] = _deck_table(deck.growth_rate))
     deck.kinetic_forces === nothing || merge!(inputs, deck_sections(deck.kinetic_forces))
     pe = deck.perturbed_equilibrium
     pe === nothing || (inputs["PerturbedEquilibrium"] = _deck_table(merge(pe.response, pe.drive)))
@@ -625,6 +645,7 @@ deck_sections(ls::LocalStabilityResult) = Dict{String,Any}("LocalStability" => _
 deck_sections(m::ForceFreeStates.MatchResult) =
     merge(deck_sections(m.model), Dict{String,Any}("Match" => merge(_deck_table(m.problem, (:ffs,)), Dict{String,Any}("inner_model" => "GGJ"))))
 deck_sections(kf::KineticForcesResult) = deck_sections(kf.control)
+deck_sections(g::IdealGrowthRateResult) = g.config === nothing ? Dict{String,Any}() : Dict{String,Any}("GrowthRate" => _deck_table(g.config))
 deck_sections(t::TearingResult) = t.model === nothing ? Dict{String,Any}() : _tearing_sections(t.problem, t.problem.critical_resonant_field, t.model, t.config)
 deck_sections(ef::ErrorFieldsResult) = _error_fields_sections(ef.settings, ef.tolerances === nothing ? "" : _RECORDED_TOLERANCES)