ForceFreeStates - REFACTOR - Reduce whole-struct passing and in-place struct mutation

docs/src/stability.md

@@ -295,8 +295,8 @@ equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(
 intr  = FFS.ForceFreeStatesInternal(; dir_path=ex)
 intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(;
     (Symbol(k) => v for (k, v) in inputs["Wall"])...)
-FFS.sing_lim!(intr, ctrl, equil)
 intr.nlow = ctrl.nn_low; intr.nhigh = ctrl.nn_high; intr.npert = 1
+FFS.sing_lim!(intr, ctrl, equil)
 FFS.sing_find!(intr, equil)
 intr.mlow  = min(intr.nlow * equil.params.qmin, 0) - 4 - ctrl.delta_mlow
 intr.mhigh = trunc(Int, intr.nhigh * equil.params.qmax) + ctrl.delta_mhigh

src/ForceFreeStates/Ballooning.jl

@@ -78,6 +78,19 @@ function compute_ballooning_stability!(
 
 end
 
+"""
+    compute_local_stability(ctrl, plasma_eq) -> CubicSeriesInterpolant
+
+Local stability profile spline over `plasma_eq.profiles.xs`, with the columns filled by
+[`compute_ballooning_stability!`](@ref): 1 = `D_I·ψ`, 2 = `D_R·ψ`, 4 = ballooning `Δ'`.
+"""
+function compute_local_stability(ctrl::ForceFreeStatesControl, plasma_eq::Equilibrium.PlasmaEquilibrium)
+    xs = plasma_eq.profiles.xs
+    locstab_fs = zeros(Float64, length(xs), 5)
+    compute_ballooning_stability!(ctrl, locstab_fs, plasma_eq)
+    return cubic_interp(xs, Series(locstab_fs); extrap=ExtendExtrap())
+end
+
 """
     resistive_interchange_h(flux_surface_index, plasma_eq)
 

src/ForceFreeStates/EulerLagrange.jl

@@ -137,25 +137,16 @@ function balance_integration_chunks(chunks::Vector{IntegrationChunk}, ctrl::Forc
 end
 
 """
-    eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal)
-
-Main driver for integrating the Euler-Lagrange equations across the plasma and detecting singular surfaces.
-Formerly `ode_run`. Has the same functionality as `ode_run` in the Fortran code, with the addition of
-a single dump to the `euler.h5` file at the end of integration instead of multiple dumps
-to `euler.bin` throughout the integration. We have made the control logic more clear
-by pre-computing all integration chunks upfront and using a for loop to iterate through them,
-eliminating the while-loop logic and making integration bounds explicit at each step.
-We now perform significant post-processing after integration including finding the peak dW
-in the edge region and evaluating the stability criterion over the entire integration,
-which were previously done during integration in the Fortran code.
+    eulerlagrange_integration(ctrl, equil, ffit, intr) -> (odet, propagators, chunks, S_left)
 
-### TODOs
-
-restype functionality if we decide to do this
+Integrate the Euler-Lagrange equations from the axis to `intr.psilim`, crossing each singular
+surface on the way (Fortran `ode_run`). Dispatches on `ctrl` to the parallel propagator BVP
+(`use_parallel`), the dual Riccati formulation (`use_riccati`), or
+[`serial_eulerlagrange_integration`](@ref).
 
-### Returns
-
-An OdeState struct containing the final state of the ODE solver after integration is complete.
+Only the parallel branch populates `propagators` / `chunks` / `S_left`, which
+`compute_delta_prime_matrix!` consumes for the Δ' BVP; the other two return `nothing` for all
+three.
 """
 function eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal)
 
@@ -166,6 +157,20 @@ function eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil::Equilibr
     elseif ctrl.use_riccati
         return (riccati_eulerlagrange_integration(ctrl, equil, ffit, intr), nothing, nothing, nothing)
     end
+    return serial_eulerlagrange_integration(ctrl, equil, ffit, intr)
+end
+
+"""
+    serial_eulerlagrange_integration(ctrl, equil, ffit, intr; verbose=ctrl.verbose) -> (odet, nothing, nothing, nothing)
+
+Serial shooting branch of [`eulerlagrange_integration`](@ref): integrates chunk by chunk,
+applying Gaussian reduction whenever a solution norm ratio exceeds `ctrl.ucrit` and undoing it
+via `transform_u!` at the end, so `odet.u_store` comes back dense in the axis basis. Call
+directly to force this branch regardless of `ctrl.use_parallel` / `ctrl.use_riccati`; `verbose`
+overrides `ctrl.verbose` for progress logging.
+"""
+function serial_eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal;
+    verbose::Bool=ctrl.verbose)
 
     # Initialization
     odet = OdeState(intr.numpert_total, ctrl.numsteps_init, ctrl.numunorms_init, intr.msing)
@@ -183,15 +188,15 @@ function eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil::Equilibr
     chunks = chunk_el_integration_bounds(odet, ctrl, intr)
 
     # Print initial integration condition
-    if ctrl.verbose
+    if verbose
         @info "   ψ = $((@sprintf "%.3f" odet.psifac)),  q = $((@sprintf "%.3f" equil.profiles.q_spline(odet.psifac)))"
     end
 
     # Iterate through each integration chunk
     for chunk in chunks
         # Integrate this region and display progress
         integrate_el_region!(odet, ctrl, equil, ffit, intr, chunk)
-        if ctrl.verbose
+        if verbose
             @info "   ψ = $((@sprintf "%.3f" odet.psifac)),  q = $((@sprintf "%.3f" odet.q)),  steps = $(odet.total_steps)"
         end
 
@@ -230,20 +235,20 @@ function eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil::Equilibr
             intr.psilim = odet.psi_store[end]
             intr.qlim = odet.q_store[end]
             odet.u .= odet.u_store[:, :, :, end]
-            if ctrl.verbose
+            if verbose
                 @info "Truncating integration at peak edge dW (LEGACY — Δ'/δW unreliable): ψ = $((@sprintf "%.3f" odet.psi_store[odet.step])),  q = $((@sprintf "%.3f" odet.q_store[odet.step]))"
             end
         else
             odet.psifac = saved_psifac
             odet.u .= saved_u
-            if ctrl.verbose
+            if verbose
                 @info "Edge-dW peak (diagnostic): ψ = $((@sprintf "%.3f" odet.psi_store[peak_step])),  q = $((@sprintf "%.3f" odet.q_store[peak_step])); integration domain unchanged"
             end
         end
     end
 
     # Evaluate stability criterion (critical determinant) of saved solutions
-    if ctrl.verbose
+    if verbose
         @info "Evaluating fixed-boundary stability criterion"
     end
     odet.nzero = evaluate_stability_criterion!(odet, equil.profiles)

src/ForceFreeStates/ForceFreeStatesStructs.jl

@@ -140,14 +140,10 @@ A mutable struct holding internal state variables for stability calculations.
   - `mlow::Int` - Lowest poloidal mode number
   - `mhigh::Int` - Highest poloidal mode number
   - `mpert::Int` - Number of poloidal modes (mhigh - mlow + 1)
-  - `nlow::Int` - Lowest toroidal mode number
-  - `nhigh::Int` - Highest toroidal mode number
+  - `nlow::Int` - Lowest toroidal mode number, resolved from `ctrl.nn_low`/`nn_high`
+  - `nhigh::Int` - Highest toroidal mode number, resolved from `ctrl.nn_low`/`nn_high`
   - `npert::Int` - Number of toroidal modes (nhigh - nlow + 1)
   - `numpert_total::Int` - Total number of modes (mpert × npert)
-  - `keq_out::Bool` - Flag to output equilibrium quantities (not yet implemented)
-  - `theta_out::Bool` - Flag to output theta coordinate data (not yet implemented)
-  - `xlmda_out::Bool` - Flag to output eigenvalue data (not yet implemented)
-  - `sol_base::Int` - Base index for solution vectors (not yet implemented)
   - `msing::Int` - Number of ideal singular surfaces
   - `kmsing::Int` - Number of kinetic singular surfaces (det(F̄) near-zeros)
   - `sing::Vector{SingType}` - Vector of ideal singular surface data
@@ -158,7 +154,6 @@ A mutable struct holding internal state variables for stability calculations.
   - `psilim::Float64` - Flux limit for integration
   - `qlim::Float64` - Safety factor at psilim
   - `q1lim::Float64` - Safety factor derivative at psilim
-  - `locstab::CubicSeriesInterpolant` - Spline for local stability analysis
   - `wall_settings::Vacuum.WallShapeSettings` - Wall shape settings for vacuum calculations
 """
 @kwdef mutable struct ForceFreeStatesInternal
@@ -170,10 +165,6 @@ A mutable struct holding internal state variables for stability calculations.
     nhigh::Int = 0
     npert::Int = 0
     numpert_total::Int = 0
-    keq_out::Bool = false
-    theta_out::Bool = false
-    xlmda_out::Bool = false
-    sol_base::Int = 50
     msing::Int = 0
     kmsing::Int = 0
     sing::Vector{SingType} = SingType[]
@@ -185,7 +176,6 @@ A mutable struct holding internal state variables for stability calculations.
     psilow::Float64 = 0.0   # lower integration bound; raised above the axis by sing_min! when qlow > qmin (RDCON gal)
     qlim::Float64 = 0.0
     q1lim::Float64 = 0.0
-    locstab::FastInterpolations.CubicSeriesInterpolant = cubic_interp(collect(0.0:0.25:1.0), Series(zeros(5, 5)); bc=ZeroCurvBC())
     debug_settings::DebugSettings = DebugSettings()
     wall_settings::Vacuum.WallShapeSettings = Vacuum.WallShapeSettings()
     """
@@ -225,7 +215,8 @@ end
 """
     ForceFreeStatesControl
 
-A mutable struct containing control parameters for stability analysis, set by the user in gpec.toml.
+An immutable struct containing 'ForceFreeStates' parameters set by the user in
+gpec.toml.
 
 ## Fields
 
@@ -268,7 +259,7 @@ A mutable struct containing control parameters for stability analysis, set by th
   - `populate_dense_xi::Bool` - When `use_parallel = true`, append a serial Euler-Lagrange pass at the end of the propagator BVP and let it replace the `odet` returned to the main pipeline.  This populates `u_store` / `du_store` / `xi_s_store` densely in the axis (EL) basis — the only convention the PerturbedEquilibrium / FieldReconstruction downstream code consumes correctly.  Without it the parallel path stores only chunk-endpoint Riccati S matrices with diagnostic derivatives (see Riccati.jl docstring caveats), and HDF5 `integration/xi_psi`/`dxi_psi`/`xi_s` are unusable.  Δ' (`singular/delta_prime_matrix`) is computed from the parallel BVP and is bit-identical between `populate_dense_xi=true` and `false`.  Energies (`vacuum/ep`/`ev`/`et`) are computed by `free_run!` from `odet`, so with `populate_dense_xi=true` they match what a pure serial run (`use_parallel=false`) would produce; with `populate_dense_xi=false` they use the parallel-pass Riccati `odet.u` instead (differs by the ~0.12 % Riccati-vs-axis algorithmic gap on DIIID-class cases).  **Default `false`** to avoid paying the dense-pass cost on Δ'/vacuum/ideal-stability-only runs; **PerturbedEquilibrium-using configs must set `populate_dense_xi = true` explicitly** when `use_parallel = true` (otherwise PE silently reads Riccati-basis garbage).  Auto-disabled when `force_termination = true` regardless of the user setting, since the dense pass has no downstream consumer in that case.  Approximate cost when enabled: one extra serial EL integration (~1× the parallel BVP wall-clock for typical N).
   - `extended_precision_bvp::Bool` - When `true` (default), promote the Δ' BVP linear system to `Complex{Double64}` (~31 digits) for the LU solve and PEST3 combination. Guards against catastrophic cancellation in the PEST3 four-term combination (dp_raw entries can be 10⁴–10⁵× larger than the result; the imaginary part of off-diagonal Δ' is particularly sensitive). Disabling (`false`) saves ~1.5–2× the BVP solve time but on DIIID-class equilibria the imaginary Δ' components can drift by factors of 2–5×; only disable for performance experiments on cases where Float64 has been validated against Double64.
 """
-@kwdef mutable struct ForceFreeStatesControl
+@kwdef struct ForceFreeStatesControl
     verbose::Bool = true
     local_stability_flag::Bool = false
     vac_flag::Bool = false

src/ForceFreeStates/Galerkin/GalerkinSolve.jl

@@ -79,19 +79,18 @@ function galerkin_solve(ctrl::ForceFreeStatesControl, equil, ffit::FourFitVars,
     # right = sig=+1, left = sig=-1, no √det normalization. The Mercier exponent α is a
     # property of the surface, so the left series reuses the right's α (alpha_override). The order is
     # raised to gal_sing_order + ceil(2·Re(α)) for high-Mercier-index surfaces (Fortran sing1_vmat).
-    ctrl_gal = deepcopy(ctrl)
     asymps = GalSingAsymp[]
     for s in sings
-        ctrl_gal.sing_order = ctrl.gal_sing_order
-        ar = compute_sing_asymptotics(s, ctrl_gal, equil, ffit, intr; sig=1.0)
+        sing_order = ctrl.gal_sing_order
+        ar = compute_sing_asymptotics(s, ctrl, equil, ffit, intr; sig=1.0, sing_order=sing_order)
         if ctrl.gal_sing_order_ceiling
             order = ctrl.gal_sing_order + ceil(Int, 2 * real(ar.alpha[1]))
             if order > ctrl.gal_sing_order
-                ctrl_gal.sing_order = order
-                ar = compute_sing_asymptotics(s, ctrl_gal, equil, ffit, intr; sig=1.0)
+                sing_order = order
+                ar = compute_sing_asymptotics(s, ctrl, equil, ffit, intr; sig=1.0, sing_order=sing_order)
             end
         end
-        al = compute_sing_asymptotics(s, ctrl_gal, equil, ffit, intr; sig=-1.0, alpha_override=ar.alpha)
+        al = compute_sing_asymptotics(s, ctrl, equil, ffit, intr; sig=-1.0, alpha_override=ar.alpha, sing_order=sing_order)
         push!(asymps, GalSingAsymp(ar, al))
     end
 

src/ForceFreeStates/Riccati.jl

@@ -1301,19 +1301,19 @@ end
 """
     riccati_eulerlagrange_integration(ctrl, equil, ffit, intr) -> OdeState
 
-Main driver for integrating the dual Riccati ODE across the plasma.
-Functionally identical to `eulerlagrange_integration` except:
-
-1. Uses `riccati_integrate_chunk!`: drives `sing_der!` with `riccati_integrator_callback!`
-   which applies `renormalize_riccati_inplace!` (instead of Gaussian reduction) when
-   column norms exceed ucrit
-2. Uses `riccati_cross_ideal_singular_surf!` instead of `cross_ideal_singular_surf!`:
-   skips Gaussian reduction (avoids near-zero pivot issues when S is small near axis)
-   and renormalizes to (S_new, I) in one step
-3. Skips `transform_u!` — S is already the true solution, no Gaussian-reduction undo needed
-
-Enable via `use_riccati = true` in `[ForceFreeStates]` section of gpec.toml, or by
-setting `ctrl.use_riccati = true` programmatically.
+Integrate the dual Riccati ODE S = U₁·U₂⁻¹ across the plasma (Glasser 2018 Phys. Plasmas 25,
+032507). Reduces stiffness relative to [`serial_eulerlagrange_integration`](@ref), which it
+otherwise mirrors, differing in three places:
+
+1. `riccati_integrate_chunk!` drives `sing_der!` with `riccati_integrator_callback!`, which
+   applies `renormalize_riccati_inplace!` rather than Gaussian reduction when column norms
+   exceed `ctrl.ucrit`
+2. `riccati_cross_ideal_singular_surf!` replaces `cross_ideal_singular_surf!`: it skips
+   Gaussian reduction (avoiding near-zero pivots where S is small near the axis) and
+   renormalizes to (S_new, I) in one step
+3. `transform_u!` is skipped — S is already the true solution, so there is no reduction to undo
+
+Enable via `use_riccati = true` in the `[ForceFreeStates]` section of gpec.toml.
 """
 function riccati_eulerlagrange_integration(
     ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium,
@@ -1605,12 +1605,14 @@ function apply_propagator_inverse!(odet::OdeState, prop::ChunkPropagator)
 end
 
 """
-    parallel_eulerlagrange_integration(ctrl, equil, ffit, intr) -> OdeState
+    parallel_eulerlagrange_integration(ctrl, equil, ffit, intr) -> (odet, propagators, chunks, S_left)
 
-Parallel fundamental matrix (propagator) driver for the EL integration.
+Parallel fundamental matrix (propagator) driver for the EL integration. The trailing three
+return values feed the Δ' BVP in `compute_delta_prime_matrix!`; this is the only branch that
+produces them.
 
-Functionally equivalent to `eulerlagrange_integration`, integrating all bulk chunks
-concurrently using `Threads.@threads`, then re-integrating the outer plasma serially:
+Equivalent to [`serial_eulerlagrange_integration`](@ref), but integrates all bulk chunks
+concurrently using `Threads.@threads`, then re-integrates the outer plasma serially:
 
 1. **Chunk generation**: calls `chunk_el_integration_bounds`, then `balance_integration_chunks`
    to sub-divide chunks for load-balanced parallel execution.
@@ -1626,22 +1628,18 @@ concurrently using `Threads.@threads`, then re-integrating the outer plasma seri
    without renormalization); Riccati integration keeps matrices bounded and provides dense
    checkpoints for `findmax_dW_edge!`.
 
-Enable via `use_parallel = true` in `[ForceFreeStates]` of gpec.toml, or by setting
-`ctrl.use_parallel = true` programmatically. Requires `singfac_min != 0`.
+Enable via `use_parallel = true` in `[ForceFreeStates]` of gpec.toml. Requires `singfac_min != 0`.
 
-**Key differences from standard integration:**
+**Key differences from serial integration:**
 - No Gaussian reduction in the propagator BVP phase (crossings use the
   Riccati-style algorithm, parallel `odet.ifix` stays 0)
 - `transform_u!` is called on the parallel odet but is a no-op (ifix=0)
 - Outer plasma uses serial Riccati integration for numerical stability
-- A serial Euler-Lagrange **dense pass** is appended at the end and
-  replaces the parallel `odet` so that `u_store` / `du_store` are dense and
-  in axis basis — the only convention the PerturbedEquilibrium downstream
-  code consumes correctly.  Δ' (`singular/delta_prime_matrix`) is computed
-  from the parallel BVP and is bit-identical with vs. without this pass.
-  Toggle off with `ctrl.populate_dense_xi = false` if only Δ' / vacuum /
-  energies are needed and the extra serial-EL cost is unwanted (HDF5
-  `integration/xi_*` will then be sparse / zero).
+- When `ctrl.populate_dense_xi` is set, a serial EL dense pass is appended and replaces the
+  parallel `odet`, so `u_store` / `du_store` / `xi_s_store` come back in the axis basis that
+  PerturbedEquilibrium requires. Δ' is computed from the parallel BVP either way and is
+  bit-identical between the two. See the `populate_dense_xi` entry in the
+  [`ForceFreeStatesControl`](@ref) docstring for the cost trade-off.
 
 **Bidirectional integration for large-N accuracy:**
 The crossing chunk (nearest to each rational surface singL[j]) is integrated *backward*
@@ -1927,7 +1925,7 @@ and does NOT populate `delta_prime`; we keep the parallel pass's values
 which `compute_delta_prime_matrix!` uses).
 
 Called from `parallel_eulerlagrange_integration` when
-`ctrl.populate_dense_xi = true` (default).  Approximate cost: one serial
+`ctrl.populate_dense_xi = true`.  Approximate cost: one serial
 EL integration on top of the parallel BVP phase.  Required to make
 `use_parallel = true` produce DCON eigenfunctions usable by the
 PerturbedEquilibrium downstream pipeline.
@@ -1955,27 +1953,12 @@ function _populate_dense_xi_via_serial_el!(
         ) for s in 1:msing],
     )
 
-    # Temporarily switch dispatch flags so `eulerlagrange_integration`
-    # follows the serial EL branch (axis-basis u_store) for this call.
-    saved_use_parallel = ctrl.use_parallel
-    saved_use_riccati  = ctrl.use_riccati
-    saved_verbose      = ctrl.verbose
-    ctrl.use_parallel = false
-    ctrl.use_riccati  = false
-    ctrl.verbose      = false  # suppress duplicate per-chunk logging
-
-    if saved_verbose
+    if ctrl.verbose
         @info "   S → ξ: serial EL dense pass for HDF5 integration/xi_*"
     end
 
-    local fresh_odet::OdeState
-    try
-        fresh_odet, _, _, _ = eulerlagrange_integration(ctrl, equil, ffit, intr)
-    finally
-        ctrl.use_parallel = saved_use_parallel
-        ctrl.use_riccati  = saved_use_riccati
-        ctrl.verbose      = saved_verbose
-    end
+    # Run the serial branch but suppress logging
+    fresh_odet, _, _, _ = serial_eulerlagrange_integration(ctrl, equil, ffit, intr; verbose=false)
 
     # Restore BVP-result fields on `intr`.
     intr.psilim = saved.psilim