From e4a458b8dd83013fe0355414d1a1d84269d574d8 Mon Sep 17 00:00:00 2001
From: d-burg <daniel.burgess@columbia.edu>
Date: Wed, 9 Sep 2026 14:04:36 -0400
Subject: [PATCH 1/2] ForceFreeStates - PERF! - Default to Vern7 with an
 explicit absolute tolerance and pin BLAS in the sweep
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

A brute-force comparison of 20 OrdinaryDiffEq solvers on the Euler-Lagrange
propagator chunks (shipped DIII-D-like case plus ten IDA_run geqdsks) showed
that error control was dominated by the OrdinaryDiffEq default abstol of 1e-6,
and that once a real absolute tolerance is set Vern9 rejects two of every three
steps near the rational surfaces while fifth-to-seventh-order methods reach the
same accuracy for about half the RHS evaluations.

- Add `ode_solver` (default "Vern7") and `ode_abstol` (default 1e-8) to
  ForceFreeStatesControl and use them at every Euler-Lagrange solve site.
  On the corpus Vern7 is 1.5-4.7x faster than Vern9 at equal tolerances, with
  Δ′ closer to the converged solution in every case; the shipped case's
  Force-Free States stage is ~13 % faster than before while its eigenvalue and
  Δ′ errors against a 1e-13/1e-15 reference drop from 5e-8/1e-4 to 2e-9/1.5e-5.
- Pin BLAS to one thread inside the sweep and the Δ′ BVP: the mpert×mpert
  blocks lose more to BLAS synchronization than they gain (RHS cost 263 → 73 µs
  at mpert=35).
- Add the benchmark scripts used for the comparison: a chunk-level
  work-precision harness, a whole-stage driver that scores δW and Δ′ against
  a converged reference, and a cold/warm stage profiler.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
---
 examples/DIIID-like_ideal_example/gpec.toml  |   2 +
 src/ForceFreeStates/CoreTypes.jl             |   4 +
 src/ForceFreeStates/EulerLagrange.jl         |  28 ++-
 src/ForceFreeStates/ForceFreeStates.jl       |  11 +
 src/ForceFreeStates/Riccati/DeltaPrimeBVP.jl | 207 ++++++++++---------
 src/ForceFreeStates/Riccati/Propagators.jl   |  65 +++---
 6 files changed, 182 insertions(+), 135 deletions(-)

diff --git a/examples/DIIID-like_ideal_example/gpec.toml b/examples/DIIID-like_ideal_example/gpec.toml
index 24b2b7f24..57431a1c7 100644
--- a/examples/DIIID-like_ideal_example/gpec.toml
+++ b/examples/DIIID-like_ideal_example/gpec.toml
@@ -45,6 +45,8 @@ mthvac = 512                   # Number of points used in splines over poloidal
 kinetic_source = "fixed"       # Kinetic matrix source: "fixed" test matrices, or "calculated" from the kinetic NTV model
 kinetic_factor = 0.0           # Scaling of kinetic matrices (0 = ideal path; >0 enables kinetic mode)
 eulerlagrange_tolerance = 1e-10 # Relative tolerance for ODE integration of Euler-Lagrange equations
+ode_abstol = 1e-8              # Absolute tolerance for the same integration
+ode_solver = "Vern7"           # OrdinaryDiffEq method for every Euler-Lagrange solve
 save_interval = 3              # Save every Nth ODE step (1=all). Always saves near rational surfaces.
 singfac_min = 1e-4             # Fractional distance from rational q at which ideal jump enforced
 ucrit = 1e4                    # Column-norm threshold that triggers solution renormalization
diff --git a/src/ForceFreeStates/CoreTypes.jl b/src/ForceFreeStates/CoreTypes.jl
index 4b0d924af..be6998845 100644
--- a/src/ForceFreeStates/CoreTypes.jl
+++ b/src/ForceFreeStates/CoreTypes.jl
@@ -131,6 +131,8 @@ gpec.toml.
   - `nstep::Int` - Maximum number of integration steps (not yet implemented)
   - `ksing::Int` - Singular surface handling parameter
   - `eulerlagrange_tolerance::Float64` - Relative tolerance for ODE integration of Euler-Lagrange equations
+  - `ode_abstol::Float64` - Absolute tolerance for the same integration. Default `1e-8`: the OrdinaryDiffEq default of `1e-6` lets small state entries escape the relative control, so the error stops responding to `eulerlagrange_tolerance` below about `1e-8`.
+  - `ode_solver::String` - OrdinaryDiffEq explicit Runge-Kutta method used for every Euler-Lagrange solve (propagator chunks, Riccati outer plasma, forward sweep, Δ′ shooting). Default `"Vern7"`: on the DIII-D corpus it reaches the same δW and a closer Δ′ than `"Vern9"` for about half the RHS evaluations, because the ninth-order method rejects most of its steps near the rational surfaces. Any name in `OrdinaryDiffEq` that accepts complex states is allowed (`"Vern6"`, `"Vern9"`, `"DP8"`, ...).
   - `ucrit::Float64` - Critical value of unorm ratio to trigger solution normalization. In the standard path it triggers Gaussian reduction; in the Riccati path it triggers `renormalize_riccati_inplace!`. Default `1e4` empirically keeps max(|U₁|, |U₂|) in O(1)–O(10⁴) over the integration domain on DIII-D / Solovev sweeps; lower triggers excess renorms without accuracy gain, higher risks overflow before the next renorm.
   - `numsteps_init::Int` - Initial array size for ODE data storage
   - `numunorms_init::Int` - Initial array size for solution normalization data
@@ -168,6 +170,8 @@ gpec.toml.
     nstep::Int = typemax(Int)
     ksing::Int = -1
     eulerlagrange_tolerance::Float64 = 1e-8
+    ode_abstol::Float64 = 1e-8
+    ode_solver::String = "Vern7"
     ucrit::Float64 = 1e4
     numsteps_init::Int = 4000
     numunorms_init::Int = 100
diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl
index 6245d1565..ad6e1d405 100644
--- a/src/ForceFreeStates/EulerLagrange.jl
+++ b/src/ForceFreeStates/EulerLagrange.jl
@@ -320,17 +320,27 @@ Only the Riccati branch populates `propagators` / `chunks` / `S_left`, which
 for all three.
 """
 function eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, intr::ForceFreeStatesInternal)
-
-    if ctrl.integrator == "riccati"
-        ctrl.kinetic_factor > 0 && error("kinetic runs require integrator=\"forward\"; the Riccati integrator has no kinetic crossing.")
-        return riccati_eulerlagrange_integration(ctrl, equil, mats, intr)
-    elseif ctrl.integrator == "forward"
-        return forward_eulerlagrange_integration(ctrl, equil, mats, intr)
-    elseif ctrl.integrator == "galerkin"
+    ctrl.integrator == "galerkin" &&
         error("integrator = \"galerkin\" solves the Euler-Lagrange system variationally, not by ODE integration; " *
               "it is dispatched to galerkin_solve.")
+    ctrl.integrator in ("riccati", "forward") ||
+        error("Unknown integrator: $(ctrl.integrator). Expected \"forward\", \"riccati\", or \"galerkin\".")
+    ctrl.integrator == "riccati" && ctrl.kinetic_factor > 0 &&
+        error("kinetic runs require integrator=\"forward\"; the Riccati integrator has no kinetic crossing.")
+
+    # The RHS works on mpert×mpert blocks, where multithreaded BLAS costs more in synchronization
+    # than it saves; pin BLAS to one thread for the sweep and restore it afterwards.
+    blas_threads = BLAS.get_num_threads()
+    BLAS.set_num_threads(1)
+    try
+        if ctrl.integrator == "riccati"
+            return riccati_eulerlagrange_integration(ctrl, equil, mats, intr)
+        else
+            return forward_eulerlagrange_integration(ctrl, equil, mats, intr)
+        end
+    finally
+        BLAS.set_num_threads(blas_threads)
     end
-    error("Unknown integrator: $(ctrl.integrator). Expected \"forward\", \"riccati\", or \"galerkin\".")
 end
 
 """
@@ -945,7 +955,7 @@ function integrate_el_region!(
 
     cb = DiscreteCallback((u, t, integrator) -> true, segment_callback!)
     prob = ODEProblem(sing_der!, odet.u, (chunk.psi_start, chunk.psi_end), (ctrl, equil, mats, intr, odet, chunk))
-    sol = solve(prob, Vern9(); reltol=ctrl.eulerlagrange_tolerance, callback=cb, save_everystep=false, save_end=true)
+    sol = solve(prob, el_ode_algorithm(ctrl); reltol=ctrl.eulerlagrange_tolerance, abstol=ctrl.ode_abstol, callback=cb, save_everystep=false, save_end=true)
 
     # Unconditionally save the final step if the callback did not already capture it.
     # Guarantees the pre-crossing (or pre-edge) state is always stored in u_store,
diff --git a/src/ForceFreeStates/ForceFreeStates.jl b/src/ForceFreeStates/ForceFreeStates.jl
index 1dc93e776..4505c5492 100644
--- a/src/ForceFreeStates/ForceFreeStates.jl
+++ b/src/ForceFreeStates/ForceFreeStates.jl
@@ -6,6 +6,17 @@ using LinearAlgebra.LAPACK
 using TOML
 using FFTW
 using OrdinaryDiffEq
+
+"""
+    el_ode_algorithm(ctrl) -> OrdinaryDiffEq algorithm
+
+The solver named by `ctrl.ode_solver`, instantiated for the Euler-Lagrange `solve` calls.
+"""
+function el_ode_algorithm(ctrl)
+    isdefined(OrdinaryDiffEq, Symbol(ctrl.ode_solver)) ||
+        error("ode_solver = \"$(ctrl.ode_solver)\" is not an OrdinaryDiffEq algorithm name (e.g. \"Vern7\", \"Vern9\", \"DP8\")")
+    return getfield(OrdinaryDiffEq, Symbol(ctrl.ode_solver))()
+end
 using HDF5
 using JLD2
 using FastInterpolations
diff --git a/src/ForceFreeStates/Riccati/DeltaPrimeBVP.jl b/src/ForceFreeStates/Riccati/DeltaPrimeBVP.jl
index 0a955fde7..6dc2c3d23 100644
--- a/src/ForceFreeStates/Riccati/DeltaPrimeBVP.jl
+++ b/src/ForceFreeStates/Riccati/DeltaPrimeBVP.jl
@@ -7,6 +7,7 @@ Compute the inter-surface tearing stability matrix (msing × msing) using the
 STRIDE global BVP formulation [Glasser 2018 Phys. Plasmas 25, 032501, Sec. III.B].
 
 The BVP encodes the full plasma response with unknowns at each surface boundary:
+
 ```
   x_axis      (N):  free IC parameters at the axis  (U₁ = 0 regular solutions)
   x_left[j]  (2N):  state at left inner-layer boundary of surface j
@@ -19,9 +20,11 @@ The BVP encodes the full plasma response with unknowns at each surface boundary:
 
 When `wv` is provided (the vacuum response matrix, singfac-scaled), the edge BC
 follows the Fortran STRIDE convention:
+
 ```
   U₁ = c,  U₂ = -wv·ψ₀²·c
 ```
+
 which is the free-boundary condition `wp + wv = 0` at the edge.
 When `wv` is `nothing`, a conducting wall BC (`U₁ = 0`) is used.
 
@@ -38,9 +41,11 @@ keeps the BVP matrix full-rank and well-conditioned.
 The raw BVP solution is a 2·msing × 2·msing matrix `dp` with left/right
 sub-indices at each surface. The PEST3-convention Δ' matrix is the linear
 combination [Chance, PPPL-2527]:
+
 ```
   deltap(i,j) = dp(2i,2j) - dp(2i,2j-1) - dp(2i-1,2j) + dp(2i-1,2j-1)
 ```
+
 stored in `intr.delta_prime_matrix` (msing × msing).
 
 ## Limitations
@@ -65,16 +70,28 @@ function compute_delta_prime_matrix!(
     intr::ForceFreeStatesInternal,
     propagators::Vector{ChunkPropagator},
     chunks::Vector{IntegrationChunk};
-    wv::Union{Nothing,Matrix{ComplexF64}} = nothing,
-    psio::Float64 = 0.0,
-    debug::Bool = false,
-    S_at_surface_left::Union{Nothing,Vector{Matrix{ComplexF64}}} = nothing,
-    ctrl::Union{Nothing,ForceFreeStatesControl} = nothing,
-    equil::Union{Nothing,Equilibrium.PlasmaEquilibrium} = nothing,
-    mats::Union{Nothing,MatrixSplines} = nothing
+    wv::Union{Nothing,Matrix{ComplexF64}}=nothing,
+    psio::Float64=0.0,
+    debug::Bool=false,
+    S_at_surface_left::Union{Nothing,Vector{Matrix{ComplexF64}}}=nothing,
+    ctrl::Union{Nothing,ForceFreeStatesControl}=nothing,
+    equil::Union{Nothing,Equilibrium.PlasmaEquilibrium}=nothing,
+    mats::Union{Nothing,MatrixSplines}=nothing
 )
     intr.msing == 0 && return
     _has_unsupported_multi_resonance(intr) && return
+    # Same small-block BLAS pin as eulerlagrange_integration: the shooting solves re-run the RHS.
+    blas_threads = BLAS.get_num_threads()
+    BLAS.set_num_threads(1)
+    try
+        _compute_delta_prime_matrix!(intr, propagators, chunks, wv, psio, debug, S_at_surface_left, ctrl, equil, mats)
+    finally
+        BLAS.set_num_threads(blas_threads)
+    end
+    return
+end
+
+function _compute_delta_prime_matrix!(intr, propagators, chunks, wv, psio, debug, S_at_surface_left, ctrl, equil, mats)
 
     sing, i_crossings, msing = _select_active_surfaces(intr, chunks)
     msing == 0 && return
@@ -91,8 +108,10 @@ function compute_delta_prime_matrix!(
     if !use_S_axis
         for ic in i_crossings
             chunks[ic].direction == 1 ||
-                error("compute_delta_prime_matrix!: FM-axis fallback (use_S_axis=false) requires forward crossing chunks; " *
-                      "chunk $ic has direction=$(chunks[ic].direction). Either provide S_at_surface_left or use bidirectional=false.")
+                error(
+                    "compute_delta_prime_matrix!: FM-axis fallback (use_S_axis=false) requires forward crossing chunks; " *
+                    "chunk $ic has direction=$(chunks[ic].direction). Either provide S_at_surface_left or use bidirectional=false."
+                )
         end
     end
 
@@ -105,15 +124,15 @@ function compute_delta_prime_matrix!(
         _build_asymptotic_basis_matrices(sing, has_ua, N, msing)
 
     debug && _log_bvp_setup(chunks, sing, S_at_surface_left, use_S_axis, has_ua,
-                            Phi_L_mats, Phi_R_mats, Phi_R_halves, ipert_all, wv, psio, N, msing)
+        Phi_L_mats, Phi_R_mats, Phi_R_halves, ipert_all, wv, psio, N, msing)
 
     if use_S_axis
         uShootR, uShootL, uAxis = _build_S_axis_shooting_propagators(
             propagators, chunks, i_crossings, sing, msing, N,
             T_left_mats, T_right_mats, has_ua, ctrl, equil, mats, intr, debug)
         debug && _log_S_axis_shooting_propagators(uShootR, uShootL, uAxis,
-                                                  S_at_surface_left, T_left_mats,
-                                                  ipert_all, has_ua, msing, N)
+            S_at_surface_left, T_left_mats,
+            ipert_all, has_ua, msing, N)
         M, nMat, col_edge = _assemble_bvp_S_axis(
             uShootR, uShootL, uAxis, ipert_all, msing, N, wv, psio)
     else
@@ -140,13 +159,13 @@ function compute_delta_prime_matrix!(
     # The raw matrix is consumed by `pest3_decompose` to recover (A', B', Γ', Δ') for the full
     # det(D' − D(γ)) = 0 eigenvalue problem; see the `delta_prime_raw` docstring in CoreTypes.jl.
     intr.delta_prime_matrix = deltap
-    intr.delta_prime_raw    = dp_raw_persisted
+    intr.delta_prime_raw = dp_raw_persisted
 end
 
 # Column index helpers for the BVP matrix. j is the 1-based singular-surface index,
 # N is numpert_total. Layout: c_axis(N), c_left[1](2N), c_right[1](2N), ..., c_edge(N).
-_col_left(j::Int, N::Int)  = (N + 4N*(j-1) + 1):(N + 4N*(j-1) + 2N)
-_col_right(j::Int, N::Int) = (N + 4N*(j-1) + 2N + 1):(N + 4N*j)
+_col_left(j::Int, N::Int) = (N+4N*(j-1)+1):(N+4N*(j-1)+2N)
+_col_right(j::Int, N::Int) = (N+4N*(j-1)+2N+1):(N+4N*j)
 
 # Multi-resonance surfaces (one q value satisfying multiple (m,n) tuples in a multi-n run)
 # are not yet handled by the inter-surface BVP. Returns true if any surface has >1 modes;
@@ -185,9 +204,9 @@ end
 # Midpoint splitting halves each inter-surface span's condition number — STRIDE's trick:
 # cond(full) = 10¹⁵ → cond(half) ≈ 10⁷·⁵, an 8-digit accuracy gain.
 function _assemble_segment_propagators(propagators::Vector{ChunkPropagator},
-                                       chunks::Vector{IntegrationChunk},
-                                       i_crossings::Vector{Int}, msing::Int, N::Int,
-                                       use_S_axis::Bool)
+    chunks::Vector{IntegrationChunk},
+    i_crossings::Vector{Int}, msing::Int, N::Int,
+    use_S_axis::Bool)
     Phi_L_mats = [assemble_fm_matrix(propagators, i_crossings[j]:i_crossings[j]) for j in 1:msing]
     Phi_R_mats = Vector{Matrix{ComplexF64}}(undef, msing + 1)
     if !use_S_axis
@@ -201,11 +220,11 @@ function _assemble_segment_propagators(propagators::Vector{ChunkPropagator},
     Phi_R_halves = Vector{Tuple{Matrix{ComplexF64},Matrix{ComplexF64}}}(undef, msing - 1)
     for j in 1:msing-1
         chunk_start = i_crossings[j] + 1
-        chunk_end   = i_crossings[j+1] - 1
-        n_chunks    = chunk_end - chunk_start + 1
+        chunk_end = i_crossings[j+1] - 1
+        n_chunks = chunk_end - chunk_start + 1
         if n_chunks >= 2
             i_mid = chunk_start + div(n_chunks, 2) - 1
-            Phi_left_half  = assemble_fm_matrix(propagators, chunk_start:i_mid)
+            Phi_left_half = assemble_fm_matrix(propagators, chunk_start:i_mid)
             Phi_right_half = assemble_fm_matrix(propagators, i_mid+1:chunk_end)
             Phi_R_halves[j] = (Phi_left_half, Phi_right_half)
         else
@@ -220,17 +239,17 @@ end
 # N+1:2N = small solutions (z^{+α}, bounded). Fortran STRIDE bakes T into the shooting
 # propagators (uFM_sing_init); we multiply T into the BVP propagator blocks at each surface.
 function _build_asymptotic_basis_matrices(sing::Vector{SingType}, has_ua::Bool, N::Int, msing::Int)
-    T_left_mats  = Vector{Matrix{ComplexF64}}(undef, msing)
+    T_left_mats = Vector{Matrix{ComplexF64}}(undef, msing)
     T_right_mats = Vector{Matrix{ComplexF64}}(undef, msing)
-    T_left_inv   = Vector{Matrix{ComplexF64}}(undef, msing)
-    T_right_inv  = Vector{Matrix{ComplexF64}}(undef, msing)
+    T_left_inv = Vector{Matrix{ComplexF64}}(undef, msing)
+    T_right_inv = Vector{Matrix{ComplexF64}}(undef, msing)
     if has_ua
         for j in 1:msing
             sp = sing[j]
-            T_left_mats[j]  = [sp.ua_left[:,:,1]; sp.ua_left[:,:,2]]
-            T_right_mats[j] = [sp.ua_right[:,:,1]; sp.ua_right[:,:,2]]
-            T_left_inv[j]   = inv(T_left_mats[j])
-            T_right_inv[j]  = inv(T_right_mats[j])
+            T_left_mats[j] = [sp.ua_left[:, :, 1]; sp.ua_left[:, :, 2]]
+            T_right_mats[j] = [sp.ua_right[:, :, 1]; sp.ua_right[:, :, 2]]
+            T_left_inv[j] = inv(T_left_mats[j])
+            T_right_inv[j] = inv(T_right_mats[j])
         end
     end
     return T_left_mats, T_right_mats, T_left_inv, T_right_inv
@@ -261,7 +280,7 @@ function _build_S_axis_shooting_propagators(
         end
         if can_reintegrate && !isempty(shoot_range_R)
             uShootR[j] = integrate_fm_with_ua_ic(chunks, shoot_range_R, sing[j].ua_right,
-                            ctrl, equil, mats, intr; backward=false, psi_ua=sing[j].psi_ua_right)
+                ctrl, equil, mats, intr; backward=false, psi_ua=sing[j].psi_ua_right)
         else
             T_init = has_ua ? T_right_mats[j] : nothing
             uShootR[j] = assemble_fm_matrix(propagators, shoot_range_R; T_init=T_init)
@@ -278,7 +297,7 @@ function _build_S_axis_shooting_propagators(
         end
         if can_reintegrate && !isempty(shoot_range_L)
             uShootL[j] = integrate_fm_with_ua_ic(chunks, shoot_range_L, sing[j].ua_left,
-                            ctrl, equil, mats, intr; backward=true, psi_ua=sing[j].psi_ua_left)
+                ctrl, equil, mats, intr; backward=true, psi_ua=sing[j].psi_ua_left)
         else
             T_init = has_ua ? T_left_mats[j] : nothing
             uShootL[j] = assemble_fm_matrix(propagators, shoot_range_L; T_init=T_init)
@@ -287,10 +306,10 @@ function _build_S_axis_shooting_propagators(
 
     uAxis, i_axis_mid = _build_conditioned_axis_propagator(propagators, i_crossings, N)
     uShootL[1] = _build_uShootL_first(propagators, chunks, i_crossings, sing,
-                                      T_left_mats, has_ua, can_reintegrate, i_axis_mid,
-                                      ctrl, equil, mats, intr, N)
+        T_left_mats, has_ua, can_reintegrate, i_axis_mid,
+        ctrl, equil, mats, intr, N)
     if debug
-        shoot_range_L1 = (i_axis_mid + 1):(i_crossings[1] - 1)
+        shoot_range_L1 = (i_axis_mid+1):(i_crossings[1]-1)
         @info "  Axis propagator: $(i_axis_mid) chunks, cond=$(@sprintf("%.2e", cond(uAxis)))"
         @info "  uShootL[1]: range=$(shoot_range_L1), cond=$(@sprintf("%.2e", cond(uShootL[1])))"
     end
@@ -303,14 +322,14 @@ end
 # chunk+1 to chunk(i_crossings[j]-1). The ψ midpoint is used (not the chunk-index midpoint)
 # because chunks near singularities are packed tighter in ψ — Fortran convention.
 function _midpoint_shoot_range(chunks::Vector{IntegrationChunk}, i_crossings::Vector{Int},
-                               j::Int, msing::Int; side::Symbol)
+    j::Int, msing::Int; side::Symbol)
     if side === :right
-        j == msing && return (i_crossings[msing] + 1):length(chunks)
+        j == msing && return (i_crossings[msing]+1):length(chunks)
         chunk_start = i_crossings[j] + 1
-        chunk_end   = i_crossings[j+1] - 1
+        chunk_end = i_crossings[j+1] - 1
     else  # :left, j >= 2
         chunk_start = i_crossings[j-1] + 1
-        chunk_end   = i_crossings[j] - 1
+        chunk_end = i_crossings[j] - 1
     end
     psi_mid_target = (chunks[chunk_start].psi_start + chunks[chunk_end].psi_end) / 2
     i_mid_inter = chunk_start
@@ -321,7 +340,7 @@ function _midpoint_shoot_range(chunks::Vector{IntegrationChunk}, i_crossings::Ve
         end
         i_mid_inter = ic
     end
-    return side === :right ? (chunk_start:i_mid_inter) : ((i_mid_inter + 1):chunk_end)
+    return side === :right ? (chunk_start:i_mid_inter) : ((i_mid_inter+1):chunk_end)
 end
 
 # Build a well-conditioned axis propagator by forward-propagating [0; I] through the
@@ -329,7 +348,7 @@ end
 # midpoint is placed one chunk before the first surface so that uShootL[1] covers only the
 # last chunk, keeping it well-conditioned.
 function _build_conditioned_axis_propagator(propagators::Vector{ChunkPropagator},
-                                            i_crossings::Vector{Int}, N::Int)
+    i_crossings::Vector{Int}, N::Int)
     n_pre_cross = i_crossings[1] - 1
     i_axis_mid = max(1, n_pre_cross - 1)
     uAxis = zeros(ComplexF64, 2N, N)
@@ -340,8 +359,8 @@ function _build_conditioned_axis_propagator(propagators::Vector{ChunkPropagator}
         prop = propagators[ic]
         upper_old = uAxis[1:N, :]
         lower_old = uAxis[N+1:2N, :]
-        uAxis[1:N, :]    .= prop.block_upper_ic[:,:,1] * upper_old .+ prop.block_lower_ic[:,:,1] * lower_old
-        uAxis[N+1:2N, :] .= prop.block_upper_ic[:,:,2] * upper_old .+ prop.block_lower_ic[:,:,2] * lower_old
+        uAxis[1:N, :] .= prop.block_upper_ic[:, :, 1] * upper_old .+ prop.block_lower_ic[:, :, 1] * lower_old
+        uAxis[N+1:2N, :] .= prop.block_upper_ic[:, :, 2] * upper_old .+ prop.block_lower_ic[:, :, 2] * lower_old
         Q, _ = qr(uAxis)
         uAxis .= Matrix(Q)[:, 1:N]
     end
@@ -355,18 +374,18 @@ end
 # Falls back to T_left_mats[1] (or identity if no ua) when there's only 1 chunk before the
 # first crossing.
 function _build_uShootL_first(propagators::Vector{ChunkPropagator},
-                              chunks::Vector{IntegrationChunk}, i_crossings::Vector{Int},
-                              sing::Vector{SingType}, T_left_mats::Vector{Matrix{ComplexF64}},
-                              has_ua::Bool, can_reintegrate::Bool, i_axis_mid::Int,
-                              ctrl, equil, mats, intr::ForceFreeStatesInternal, N::Int)
-    shoot_range_L1 = (i_axis_mid + 1):(i_crossings[1] - 1)
+    chunks::Vector{IntegrationChunk}, i_crossings::Vector{Int},
+    sing::Vector{SingType}, T_left_mats::Vector{Matrix{ComplexF64}},
+    has_ua::Bool, can_reintegrate::Bool, i_axis_mid::Int,
+    ctrl, equil, mats, intr::ForceFreeStatesInternal, N::Int)
+    shoot_range_L1 = (i_axis_mid+1):(i_crossings[1]-1)
     if can_reintegrate && !isempty(shoot_range_L1)
         return integrate_fm_with_ua_ic(chunks, shoot_range_L1, sing[1].ua_left,
-                                       ctrl, equil, mats, intr;
-                                       backward=true, psi_ua=sing[1].psi_ua_left)
+            ctrl, equil, mats, intr;
+            backward=true, psi_ua=sing[1].psi_ua_left)
     elseif !isempty(shoot_range_L1)
         return assemble_fm_matrix(propagators, shoot_range_L1;
-                                  T_init=has_ua ? T_left_mats[1] : nothing)
+            T_init=has_ua ? T_left_mats[1] : nothing)
     else
         return has_ua ? T_left_mats[1] : Matrix{ComplexF64}(I, 2N, 2N)
     end
@@ -377,19 +396,19 @@ end
 # the catastrophically ill-conditioned axis FM. Fortran-matched structure with
 # nMat = (2 + 4·msing)·N. Returns (M, nMat, col_edge).
 function _assemble_bvp_S_axis(uShootR::Vector{Matrix{ComplexF64}},
-                              uShootL::Vector{Matrix{ComplexF64}},
-                              uAxis::Matrix{ComplexF64}, ipert_all::Vector{Int},
-                              msing::Int, N::Int,
-                              wv::Union{Nothing,Matrix{ComplexF64}}, psio::Float64)
+    uShootL::Vector{Matrix{ComplexF64}},
+    uAxis::Matrix{ComplexF64}, ipert_all::Vector{Int},
+    msing::Int, N::Int,
+    wv::Union{Nothing,Matrix{ComplexF64}}, psio::Float64)
     # STRIDE global BVP block structure [Glasser-Kolemen 2018 PoP 25, 032501 Eq. 37].
     nMat = (2 + 4 * msing) * N
     col_axis = 1:N
-    col_edge = (nMat - N + 1):nMat
+    col_edge = (nMat-N+1):nMat
     M = zeros(ComplexF64, nMat, nMat)
 
     # Axis matching: uShootL[1] · c_left[1] = uAxis · c_axis  (2N equations)
     M[1:2N, _col_left(1, N)] .= uShootL[1]
-    M[1:2N, col_axis]        .= -uAxis
+    M[1:2N, col_axis] .= -uAxis
     row_offset = 2N
 
     for j in 1:msing
@@ -398,21 +417,21 @@ function _assemble_bvp_S_axis(uShootR::Vector{Matrix{ComplexF64}},
         for i in 1:2N
             if i != ipert_j && i != ipert_j + N
                 row_offset += 1
-                M[row_offset, _col_left(j, N)[i]]  =  1
+                M[row_offset, _col_left(j, N)[i]] = 1
                 M[row_offset, _col_right(j, N)[i]] = -1
             end
         end
 
-        junc_rows = (row_offset + 1):(row_offset + 2N)
+        junc_rows = (row_offset+1):(row_offset+2N)
         if j < msing
             # Midpoint matching between consecutive surfaces
-            M[junc_rows, _col_right(j, N)]   .= -uShootR[j]
-            M[junc_rows, _col_left(j+1, N)]  .=  uShootL[j+1]
+            M[junc_rows, _col_right(j, N)] .= -uShootR[j]
+            M[junc_rows, _col_left(j + 1, N)] .= uShootL[j+1]
         else
             # Edge junction
             M[junc_rows, _col_right(msing, N)] .= uShootR[msing]
             if wv !== nothing
-                M[junc_rows[1:N],     col_edge] .= -I(N)
+                M[junc_rows[1:N], col_edge] .= -I(N)
                 M[junc_rows[N+1:end], col_edge] .= wv .* psio^2
             else
                 M[junc_rows[N+1:end], col_edge] .= -I(N)
@@ -425,7 +444,7 @@ function _assemble_bvp_S_axis(uShootR::Vector{Matrix{ComplexF64}},
     for j in 1:msing
         ipert_j = ipert_all[j]
         row_offset += 1
-        M[row_offset, _col_left(j, N)[ipert_j]]  = 1
+        M[row_offset, _col_left(j, N)[ipert_j]] = 1
         row_offset += 1
         M[row_offset, _col_right(j, N)[ipert_j]] = 1
     end
@@ -460,51 +479,51 @@ end
 # Fallback BVP assembly with FM-based axis BC (used when no Riccati S matrices are available).
 # Uses the conditioned axis propagator Phi_R[1][:,N+1:2N] in place of S-axis matching.
 function _assemble_bvp_FM_axis(Phi_L_mats::Vector{Matrix{ComplexF64}},
-                               Phi_R_mats::Vector{Matrix{ComplexF64}}, ipert_all::Vector{Int},
-                               msing::Int, N::Int,
-                               T_left_inv::Vector{Matrix{ComplexF64}},
-                               T_right_inv::Vector{Matrix{ComplexF64}}, has_ua::Bool,
-                               wv::Union{Nothing,Matrix{ComplexF64}}, psio::Float64)
+    Phi_R_mats::Vector{Matrix{ComplexF64}}, ipert_all::Vector{Int},
+    msing::Int, N::Int,
+    T_left_inv::Vector{Matrix{ComplexF64}},
+    T_right_inv::Vector{Matrix{ComplexF64}}, has_ua::Bool,
+    wv::Union{Nothing,Matrix{ComplexF64}}, psio::Float64)
     nMat = (2 + 4 * msing) * N
     col_axis = 1:N
-    col_edge = (N + 4N*msing + 1):nMat
+    col_edge = (N+4N*msing+1):nMat
     M = zeros(ComplexF64, nMat, nMat)
 
     M[1:2N, (N+1):(N+2N)] .= Phi_L_mats[1]
-    M[1:2N, col_axis]     .= -view(Phi_R_mats[1], :, N+1:2N)
+    M[1:2N, col_axis] .= -view(Phi_R_mats[1], :, N+1:2N)
 
-    row_drive_base = 2N + (4N-2)*msing
+    row_drive_base = 2N + (4N - 2) * msing
     for j in 1:msing
         ipert_j = ipert_all[j]
         cl = _col_left(j, N)
         cr = _col_right(j, N)
-        row_cont = 2N + (4N-2)*(j-1)
+        row_cont = 2N + (4N - 2) * (j - 1)
         for i in 1:2N
             if i != ipert_j && i != ipert_j + N
                 row_cont += 1
-                M[row_cont, cl[i]] =  1
+                M[row_cont, cl[i]] = 1
                 M[row_cont, cr[i]] = -1
             end
         end
-        junc_rows = (row_cont + 1):(2N + (4N-2)*j)
+        junc_rows = (row_cont+1):(2N+(4N-2)*j)
         if j < msing
-            M[junc_rows, cr]                .=  Phi_R_mats[j+1]
-            M[junc_rows, _col_left(j+1, N)] .= -Phi_L_mats[j+1]
+            M[junc_rows, cr] .= Phi_R_mats[j+1]
+            M[junc_rows, _col_left(j + 1, N)] .= -Phi_L_mats[j+1]
         else
             M[junc_rows, cr] .= Phi_R_mats[msing+1]
             if wv !== nothing
-                M[junc_rows[1:N],     col_edge] .= -I(N)
+                M[junc_rows[1:N], col_edge] .= -I(N)
                 M[junc_rows[N+1:end], col_edge] .= wv .* psio^2
             else
                 M[junc_rows[N+1:end], col_edge] .= -I(N)
             end
         end
         if has_ua
-            M[row_drive_base + 2j-1, cl] .= T_left_inv[j][ipert_j, :]
-            M[row_drive_base + 2j,   cr] .= T_right_inv[j][ipert_j, :]
+            M[row_drive_base+2j-1, cl] .= T_left_inv[j][ipert_j, :]
+            M[row_drive_base+2j, cr] .= T_right_inv[j][ipert_j, :]
         else
-            M[row_drive_base + 2j-1, cl[ipert_j]] = 1
-            M[row_drive_base + 2j,   cr[ipert_j]] = 1
+            M[row_drive_base+2j-1, cl[ipert_j]] = 1
+            M[row_drive_base+2j, cr[ipert_j]] = 1
         end
     end
     return M, nMat, col_edge
@@ -515,8 +534,8 @@ end
 # combination subtracts dp_raw entries up to ~3×10⁴ larger than the result, and Float64
 # precision lets the imaginary part drift 2–5× on DIIID-class equilibria.
 function _solve_bvp_and_combine_pest3(M::Matrix{ComplexF64}, msing::Int, N::Int, nMat::Int,
-                                      use_S_axis::Bool, ipert_all::Vector{Int}, col_edge,
-                                      ctrl, debug::Bool)
+    use_S_axis::Bool, ipert_all::Vector{Int}, col_edge,
+    ctrl, debug::Bool)
     s2 = 2 * msing
     Tc = (ctrl === nothing || ctrl.extended_precision_bvp) ? Complex{Double64} : ComplexF64
     M_solve = Tc.(M)
@@ -533,17 +552,17 @@ function _solve_bvp_and_combine_pest3(M::Matrix{ComplexF64}, msing::Int, N::Int,
     for jsing in 1:msing, side in 1:2
         dRow = 2jsing - (2 - side)
         fill!(b, 0)
-        drive_row = use_S_axis ? (nMat - s2 + dRow) : (2N + (4N-2)*msing + dRow)
+        drive_row = use_S_axis ? (nMat - s2 + dRow) : (2N + (4N - 2) * msing + dRow)
         b[drive_row] = 1
         x = use_lu ? (M_lu \ b) : (M_pinv * b)
 
         debug && _log_bvp_solve(x, b, M_solve, jsing, side, dRow, msing, N,
-                                ipert_all, col_edge, use_S_axis)
+            ipert_all, col_edge, use_S_axis)
 
         for ksing in 1:msing
             ipert_k = ipert_all[ksing]
             dp_raw[dRow, 2ksing-1] = x[_col_left(ksing, N)[ipert_k+N]]
-            dp_raw[dRow, 2ksing]   = x[_col_right(ksing, N)[ipert_k+N]]
+            dp_raw[dRow, 2ksing] = x[_col_right(ksing, N)[ipert_k+N]]
         end
     end
 
@@ -563,7 +582,7 @@ end
 
 # Logging helpers for `compute_delta_prime_matrix!`. Called only when debug=true.
 function _log_bvp_setup(chunks, sing, S_at_surface_left, use_S_axis, has_ua,
-                        Phi_L_mats, Phi_R_mats, Phi_R_halves, ipert_all, wv, psio, N, msing)
+    Phi_L_mats, Phi_R_mats, Phi_R_halves, ipert_all, wv, psio, N, msing)
     @info "Δ' BVP: $(length(chunks)) chunks, $msing surfaces, N=$N"
     @info "Δ' BVP: Axis BC: $(use_S_axis ? "S-based (Riccati)" : "FM-based (conditioned)")"
     @info "Δ' BVP: Asymptotic basis: $(has_ua ? "available" : "NOT available (raw basis driving)")"
@@ -575,8 +594,8 @@ function _log_bvp_setup(chunks, sing, S_at_surface_left, use_S_axis, has_ua,
     if has_ua
         for j in 1:msing
             sp = sing[j]
-            T_l = [sp.ua_left[:,:,1]; sp.ua_left[:,:,2]]
-            T_r = [sp.ua_right[:,:,1]; sp.ua_right[:,:,2]]
+            T_l = [sp.ua_left[:, :, 1]; sp.ua_left[:, :, 2]]
+            T_r = [sp.ua_right[:, :, 1]; sp.ua_right[:, :, 2]]
             @info "  Surface $j: cond(T_left)=$(@sprintf("%.2e", cond(T_l))), cond(T_right)=$(@sprintf("%.2e", cond(T_r)))"
             ipert_j = ipert_all[j]
             @info "  Surface $j ua_left (ipert=$ipert_j, psi_ua_left=$(@sprintf("%.8f", sp.psi_ua_left))):"
@@ -603,7 +622,7 @@ function _log_bvp_setup(chunks, sing, S_at_surface_left, use_S_axis, has_ua,
 end
 
 function _log_S_axis_shooting_propagators(uShootR, uShootL, uAxis, S_at_surface_left,
-                                          T_left_mats, ipert_all, has_ua, msing, N)
+    T_left_mats, ipert_all, has_ua, msing, N)
     @info "  Shooting propagators (S-based axis BC, no axis unknowns):"
     for j in 1:msing
         shoot_R_str = @sprintf("%.2e", cond(uShootR[j]))
@@ -636,7 +655,7 @@ function _log_S_axis_shooting_propagators(uShootR, uShootL, uAxis, S_at_surface_
 end
 
 function _log_bvp_solve(x, b, M_solve, jsing, side, dRow, msing, N,
-                        ipert_all, col_edge, use_S_axis)
+    ipert_all, col_edge, use_S_axis)
     residual = norm(ComplexF64.(M_solve * x - b))
     side_str = side == 1 ? "left" : "right"
     @info "  BVP solve: jsing=$jsing side=$side_str (dRow=$dRow): ||Mx-b||=$(@sprintf("%.2e", residual)), ||x||=$(@sprintf("%.2e", Float64(norm(x))))"
@@ -644,9 +663,9 @@ function _log_bvp_solve(x, b, M_solve, jsing, side, dRow, msing, N,
         ipert_ks = ipert_all[ks]
         cl = _col_left(ks, N)
         cr = _col_right(ks, N)
-        xl_big   = ComplexF64(x[cl[ipert_ks]])
+        xl_big = ComplexF64(x[cl[ipert_ks]])
         xl_small = ComplexF64(x[cl[ipert_ks+N]])
-        xr_big   = ComplexF64(x[cr[ipert_ks]])
+        xr_big = ComplexF64(x[cr[ipert_ks]])
         xr_small = ComplexF64(x[cr[ipert_ks+N]])
         @info "    surf $ks: x_left[big]=$(@sprintf("%+.4e%+.4ei", real(xl_big), imag(xl_big))), x_left[small]=$(@sprintf("%+.4e%+.4ei", real(xl_small), imag(xl_small)))"
         @info "    surf $ks: x_right[big]=$(@sprintf("%+.4e%+.4ei", real(xr_big), imag(xr_big))), x_right[small]=$(@sprintf("%+.4e%+.4ei", real(xr_small), imag(xr_small)))"
@@ -660,7 +679,7 @@ end
 function _log_bvp_pest3(dp_raw, deltap, s2, msing, Tc)
     @info "Δ' BVP: Full dp_raw matrix ($(s2)×$(s2)) [$(Tc)]:"
     for i in 1:s2
-        row_str = join([@sprintf("%+.6e", Float64(real(dp_raw[i,j]))) for j in 1:s2], "  ")
+        row_str = join([@sprintf("%+.6e", Float64(real(dp_raw[i, j]))) for j in 1:s2], "  ")
         @info "  dp_raw[$i,:] = $row_str"
     end
     @info "Δ' BVP: Raw dp diagonal = $([@sprintf("%.4f%+.4fi", Float64(real(dp_raw[i,i])), Float64(imag(dp_raw[i,i]))) for i in 1:s2])"
@@ -720,8 +739,8 @@ function pest3_decompose(dp_raw::AbstractMatrix)
     for i in 1:m, j in 1:m
         LL = dp_raw[2i-1, 2j-1]
         LR = dp_raw[2i-1, 2j]
-        RL = dp_raw[2i,   2j-1]
-        RR = dp_raw[2i,   2j]
+        RL = dp_raw[2i, 2j-1]
+        RR = dp_raw[2i, 2j]
         Ap[i, j] = RR + RL + LR + LL
         Bp[i, j] = RR - RL + LR - LL
         Gp[i, j] = RR + RL - LR - LL
diff --git a/src/ForceFreeStates/Riccati/Propagators.jl b/src/ForceFreeStates/Riccati/Propagators.jl
index 285695c5b..d888750b7 100644
--- a/src/ForceFreeStates/Riccati/Propagators.jl
+++ b/src/ForceFreeStates/Riccati/Propagators.jl
@@ -6,7 +6,7 @@
 # length chunks; the absolute floor catches short chunks where 5% of the span would be
 # smaller than the typical ODE step.
 const SAVE_NEAR_END_FRAC = 0.05
-const SAVE_NEAR_END_PSI  = 1e-4
+const SAVE_NEAR_END_PSI = 1e-4
 
 """
     assemble_fm_matrix(propagators, idx_range; condition=false) -> Matrix{ComplexF64}
@@ -16,6 +16,7 @@ in order for indices `idx_range`. Returns Φ_end * ... * Φ_start, so that the r
 maps the IC at the start of `idx_range[1]` to the state at the end of `idx_range[end]`.
 
 Each `ChunkPropagator` stores the 2N columns of Φ split into two N×N×2 blocks:
+
 ```
   block_upper_ic[:,:,1:2] ↔ Φ[:,1:N]     (result from IC=(I,0))
   block_lower_ic[:,:,1:2] ↔ Φ[:,N+1:2N]  (result from IC=(0,I))
@@ -33,8 +34,8 @@ means only U₂ ICs are needed. Do NOT use for inter-surface segments where both
 and U₂ components carry physical information.
 """
 function assemble_fm_matrix(propagators::Vector{ChunkPropagator}, idx_range;
-                            condition::Bool=false,
-                            T_init::Union{Nothing,Matrix{ComplexF64}}=nothing)
+    condition::Bool=false,
+    T_init::Union{Nothing,Matrix{ComplexF64}}=nothing)
     # Determine matrix size from T_init if provided (lets us handle empty idx_range and even
     # an empty propagators list, provided T_init carries the dimension). Otherwise fall back
     # to the first propagator that actually exists in idx_range, with a final fallback to
@@ -123,7 +124,7 @@ end
     riccati_der!(du, u, params, psieval)
 
 Evaluate the explicit dual Riccati ODE right-hand side:
-  dS/dψ = w†·F̄⁻¹·w - S·Ḡ·S,   w = Q - K̄·S
+dS/dψ = w†·F̄⁻¹·w - S·Ḡ·S,   w = Q - K̄·S
 
 where Q = diag(1/(m - n·q)) is the diagonal singular factor matrix.
 The identity slice u[:,:,2] = I does not evolve (du[:,:,2] = 0).
@@ -147,7 +148,7 @@ See: Glasser (2018) Phys. Plasmas 25, 032507 — Eq. 19 (dual Riccati form)
     _, equil, mats, intr, odet, _ = params
 
     Npert = intr.numpert_total
-    S  = @view u[:, :, 1]
+    S = @view u[:, :, 1]
     dS = @view du[:, :, 1]
     @view(du[:, :, 2]) .= 0  # identity does not evolve
 
@@ -162,9 +163,9 @@ See: Glasser (2018) Phys. Plasmas 25, 032507 — Eq. 19 (dual Riccati form)
     fmat_lower = acquire!(pool, ComplexF64, Npert, Npert)
     kmat = similar!(pool, fmat_lower)
     gmat = similar!(pool, fmat_lower)
-    w    = similar!(pool, fmat_lower)  # w = Q - K̄·S
-    v    = similar!(pool, fmat_lower)  # v = F̄⁻¹·w (then reused for S·Ḡ·S)
-    tmp  = similar!(pool, fmat_lower)  # scratch
+    w = similar!(pool, fmat_lower)  # w = Q - K̄·S
+    v = similar!(pool, fmat_lower)  # v = F̄⁻¹·w (then reused for S·Ḡ·S)
+    tmp = similar!(pool, fmat_lower)  # scratch
 
     # Evaluate F̄ (Cholesky factor), K̄, Ḡ splines at current ψ
     mats.ideal.F_spline_lower(vec(fmat_lower), psieval; hint=mats._hint)
@@ -255,8 +256,8 @@ function riccati_integrate_chunk!(
     cb = DiscreteCallback((u, t, integrator) -> true, riccati_integrator_callback!)
     rtol = ctrl.eulerlagrange_tolerance
     prob = ODEProblem(sing_der!, odet.u, (chunk.psi_start, chunk.psi_end),
-                      (ctrl, equil, mats, intr, odet, chunk))
-    sol = solve(prob, Vern9(); reltol=rtol, callback=cb, save_everystep=false, save_end=true)
+        (ctrl, equil, mats, intr, odet, chunk))
+    sol = solve(prob, el_ode_algorithm(ctrl); reltol=rtol, abstol=ctrl.ode_abstol, callback=cb, save_everystep=false, save_end=true)
     odet.u .= sol.u[end]
     odet.psifac = sol.t[end]
     # Renormalize end state to (S, I) convention for the next chunk.
@@ -273,8 +274,8 @@ end
     renormalize_riccati!(odet, intr)
 
 After a singular surface crossing, restore the canonical Riccati storage convention:
-  u[:,:,1] = S_new = U₁_new · U₂_new⁻¹
-  u[:,:,2] = I
+u[:,:,1] = S_new = U₁_new · U₂_new⁻¹
+u[:,:,2] = I
 
 `riccati_cross_ideal_singular_surf!` leaves u[:,:,1] = U₁_new and u[:,:,2] = U₂_new (not I),
 so this step is required before continuing the Riccati integration.
@@ -298,8 +299,8 @@ end
     renormalize_riccati_inplace!(u, N)
 
 In-place Riccati renormalization on an arbitrary N×N×2 array:
-  u[:,:,1] = U₁ · U₂⁻¹  (new S)
-  u[:,:,2] = I
+u[:,:,1] = U₁ · U₂⁻¹  (new S)
+u[:,:,2] = I
 
 Used in `riccati_integrator_callback!` to renormalize the integrator's live state
 when column norms grow beyond `ctrl.ucrit`, analogous to Gaussian reduction in the
@@ -346,8 +347,8 @@ function integrate_propagator_chunk!(
     # naturally. The resulting propagator maps state at psi_end → psi_start, which is
     # well-conditioned because exponentially growing solutions (forward) decay backward.
     tspan = chunk.direction == 1 ?
-        (chunk.psi_start, chunk.psi_end) :
-        (chunk.psi_end,   chunk.psi_start)
+            (chunk.psi_start, chunk.psi_end) :
+            (chunk.psi_end, chunk.psi_start)
     rtol = ctrl.eulerlagrange_tolerance
     params = (ctrl, equil, mats, intr, odet_proxy, chunk)
 
@@ -359,7 +360,7 @@ function integrate_propagator_chunk!(
     odet_proxy.spline_hint[] = 1
     odet_proxy.mats_hint[] = 1
     prob = ODEProblem(sing_der!, u_upper, tspan, params)
-    sol = solve(prob, Vern9(); reltol=rtol, save_everystep=false, save_end=true)
+    sol = solve(prob, el_ode_algorithm(ctrl); reltol=rtol, abstol=ctrl.ode_abstol, save_everystep=false, save_end=true)
     prop.block_upper_ic .= sol.u[end]
     odet_proxy.total_steps += sol.stats.naccept  # thread-local; summed into odet after the BVP barrier
 
@@ -371,7 +372,7 @@ function integrate_propagator_chunk!(
     odet_proxy.spline_hint[] = 1
     odet_proxy.mats_hint[] = 1
     prob = ODEProblem(sing_der!, u_lower, tspan, params)
-    sol = solve(prob, Vern9(); reltol=rtol, save_everystep=false, save_end=true)
+    sol = solve(prob, el_ode_algorithm(ctrl); reltol=rtol, abstol=ctrl.ode_abstol, save_everystep=false, save_end=true)
     prop.block_lower_ic .= sol.u[end]
     odet_proxy.total_steps += sol.stats.naccept
 end
@@ -401,12 +402,12 @@ function integrate_fm_with_ua_ic(
     equil::Equilibrium.PlasmaEquilibrium,
     mats::MatrixSplines,
     intr::ForceFreeStatesInternal;
-    backward::Bool = false,
-    psi_ua::Float64 = NaN
+    backward::Bool=false,
+    psi_ua::Float64=NaN
 )
     N = intr.numpert_total
     psi_start = chunks[first(chunk_range)].psi_start
-    psi_end   = chunks[last(chunk_range)].psi_end
+    psi_end = chunks[last(chunk_range)].psi_end
     # Use stored ua ψ location if provided; otherwise fall back to chunk boundary.
     # The ua is evaluated at the inner-layer boundary (exact ψ from singular crossing),
     # which may differ slightly from the nearest chunk boundary.
@@ -437,9 +438,9 @@ function integrate_fm_with_ua_ic(
     odet_proxy.spline_hint[] = 1
     odet_proxy.mats_hint[] = 1
     prob = ODEProblem(sing_der!, u0, tspan, params)
-    sol = solve(prob, Vern9(); reltol=rtol, abstol=abstol_arr, save_everystep=false, save_end=true)
-    result[1:N, 1:N]     .= sol.u[end][:, :, 1]
-    result[N+1:2N, 1:N]  .= sol.u[end][:, :, 2]
+    sol = solve(prob, el_ode_algorithm(ctrl); reltol=rtol, abstol=abstol_arr, save_everystep=false, save_end=true)
+    result[1:N, 1:N] .= sol.u[end][:, :, 1]
+    result[N+1:2N, 1:N] .= sol.u[end][:, :, 2]
 
     # Batch 2: columns N+1:2N of T (small solutions)
     u0[:, :, 1] .= ua[:, N+1:2N, 1]
@@ -450,9 +451,9 @@ function integrate_fm_with_ua_ic(
     odet_proxy.spline_hint[] = 1
     odet_proxy.mats_hint[] = 1
     prob = ODEProblem(sing_der!, u0, tspan, params)
-    sol = solve(prob, Vern9(); reltol=rtol, abstol=abstol_arr, save_everystep=false, save_end=true)
-    result[1:N, N+1:2N]     .= sol.u[end][:, :, 1]
-    result[N+1:2N, N+1:2N]  .= sol.u[end][:, :, 2]
+    sol = solve(prob, el_ode_algorithm(ctrl); reltol=rtol, abstol=abstol_arr, save_everystep=false, save_end=true)
+    result[1:N, N+1:2N] .= sol.u[end][:, :, 1]
+    result[N+1:2N, N+1:2N] .= sol.u[end][:, :, 2]
 
     return result
 end
@@ -464,8 +465,8 @@ Apply the chunk propagator `prop` to the current state `odet.u` in-place.
 
 The propagator acts as a linear map on the (U₁, U₂) pair:
 
-  U₁_new = block_upper_ic[:,:,1] · U₁_prev + block_lower_ic[:,:,1] · U₂_prev
-  U₂_new = block_upper_ic[:,:,2] · U₁_prev + block_lower_ic[:,:,2] · U₂_prev
+U₁_new = block_upper_ic[:,:,1] · U₁_prev + block_lower_ic[:,:,1] · U₂_prev
+U₂_new = block_upper_ic[:,:,2] · U₁_prev + block_lower_ic[:,:,2] · U₂_prev
 
 This correctly propagates any state (not just the identity), including the
 (S, I) form produced by Riccati-style crossings.
@@ -520,8 +521,8 @@ function apply_propagator_inverse!(odet::OdeState, prop::ChunkPropagator)
     # Φ_bwd maps state at psi_end → psi_start (well-conditioned).
     # We want Φ_fwd = Φ_bwd⁻¹ to advance state from psi_start → psi_end.
     # Solving Φ_bwd · x = [U₁_old; U₂_old] gives x = Φ_bwd⁻¹ · [U₁_old; U₂_old].
-    u_old = [odet.u[:,:,1]; odet.u[:,:,2]]   # 2N × N
+    u_old = [odet.u[:, :, 1]; odet.u[:, :, 2]]   # 2N × N
     u_new = Φ \ u_old                         # LU solve, 2N × N
-    odet.u[:,:,1] .= u_new[1:N, :]
-    odet.u[:,:,2] .= u_new[N+1:2N, :]
+    odet.u[:, :, 1] .= u_new[1:N, :]
+    odet.u[:, :, 2] .= u_new[N+1:2N, :]
 end

From 9fc00ddd57270a8754af3fc1ae4245b0d2d78da6 Mon Sep 17 00:00:00 2001
From: d-burg <daniel.burgess@columbia.edu>
Date: Thu, 24 Sep 2026 15:58:15 -0400
Subject: [PATCH 2/2] ForceFreeStates - MINOR - Validate ode_solver against the
 benchmarked methods and repair the integrator docstrings
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

- ode_solver must be one of the benchmarked methods (Vern6, Vern7, Vern8, Vern9, DP8), looked up
  in EL_ODE_SOLVERS; both control construction sites check it before any integration, so a bad
  name fails immediately instead of at the first solve.
- el_ode_algorithm moves out of the middle of the using block.
- ode_abstol's docstring says the Δ′ shooting solves use their own per-column absolute tolerance.
- The equation blocks the formatter flattened into running text are fenced, and the Riccati
  driver's docs no longer describe Vern9 as the integrator.
- Tests: an unknown ode_solver is refused by solve; every benchmarked name resolves.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 examples/DIIID-like_ideal_example/gpec.toml |  4 +--
 src/ForceFreeStates/CoreTypes.jl            |  4 +--
 src/ForceFreeStates/ForceFreeStates.jl      | 28 +++++++++++++--------
 src/ForceFreeStates/Riccati/Driver.jl       | 16 ++++++------
 src/ForceFreeStates/Riccati/Propagators.jl  |  8 ++++++
 src/GeneralizedPerturbedEquilibrium.jl      | 12 +++++----
 test/runtests_solve_api.jl                  | 13 +++++++++-
 7 files changed, 57 insertions(+), 28 deletions(-)

diff --git a/examples/DIIID-like_ideal_example/gpec.toml b/examples/DIIID-like_ideal_example/gpec.toml
index 57431a1c7..2f023ac57 100644
--- a/examples/DIIID-like_ideal_example/gpec.toml
+++ b/examples/DIIID-like_ideal_example/gpec.toml
@@ -45,8 +45,8 @@ mthvac = 512                   # Number of points used in splines over poloidal
 kinetic_source = "fixed"       # Kinetic matrix source: "fixed" test matrices, or "calculated" from the kinetic NTV model
 kinetic_factor = 0.0           # Scaling of kinetic matrices (0 = ideal path; >0 enables kinetic mode)
 eulerlagrange_tolerance = 1e-10 # Relative tolerance for ODE integration of Euler-Lagrange equations
-ode_abstol = 1e-8              # Absolute tolerance for the same integration
-ode_solver = "Vern7"           # OrdinaryDiffEq method for every Euler-Lagrange solve
+ode_abstol = 1e-8              # Absolute tolerance for the Euler-Lagrange sweeps and propagator chunks
+ode_solver = "Vern7"           # Explicit Runge-Kutta method for every Euler-Lagrange solve: Vern6, Vern7, Vern8, Vern9 or DP8
 save_interval = 3              # Save every Nth ODE step (1=all). Always saves near rational surfaces.
 singfac_min = 1e-4             # Fractional distance from rational q at which ideal jump enforced
 ucrit = 1e4                    # Column-norm threshold that triggers solution renormalization
diff --git a/src/ForceFreeStates/CoreTypes.jl b/src/ForceFreeStates/CoreTypes.jl
index be6998845..1a98a5991 100644
--- a/src/ForceFreeStates/CoreTypes.jl
+++ b/src/ForceFreeStates/CoreTypes.jl
@@ -131,8 +131,8 @@ gpec.toml.
   - `nstep::Int` - Maximum number of integration steps (not yet implemented)
   - `ksing::Int` - Singular surface handling parameter
   - `eulerlagrange_tolerance::Float64` - Relative tolerance for ODE integration of Euler-Lagrange equations
-  - `ode_abstol::Float64` - Absolute tolerance for the same integration. Default `1e-8`: the OrdinaryDiffEq default of `1e-6` lets small state entries escape the relative control, so the error stops responding to `eulerlagrange_tolerance` below about `1e-8`.
-  - `ode_solver::String` - OrdinaryDiffEq explicit Runge-Kutta method used for every Euler-Lagrange solve (propagator chunks, Riccati outer plasma, forward sweep, Δ′ shooting). Default `"Vern7"`: on the DIII-D corpus it reaches the same δW and a closer Δ′ than `"Vern9"` for about half the RHS evaluations, because the ninth-order method rejects most of its steps near the rational surfaces. Any name in `OrdinaryDiffEq` that accepts complex states is allowed (`"Vern6"`, `"Vern9"`, `"DP8"`, ...).
+  - `ode_abstol::Float64` - Absolute tolerance for the forward sweep, the Riccati outer plasma and the propagator chunks. Default `1e-8`: the OrdinaryDiffEq default of `1e-6` lets small state entries escape the relative control, so the error stops responding to `eulerlagrange_tolerance` below about `1e-8`. The Δ′ shooting solves do not use it; they scale a per-column absolute tolerance from `eulerlagrange_tolerance`.
+  - `ode_solver::String` - Explicit Runge-Kutta method for every Euler-Lagrange solve, one of `"Vern6"`, `"Vern7"`, `"Vern8"`, `"Vern9"`, `"DP8"` (the benchmarked set, [`EL_ODE_SOLVERS`](@ref)). Default `"Vern7"`: at equal tolerances it reaches the same δW and a closer Δ′ than `"Vern9"` for fewer RHS evaluations, because the ninth-order method rejects most of its steps near the rational surfaces.
   - `ucrit::Float64` - Critical value of unorm ratio to trigger solution normalization. In the standard path it triggers Gaussian reduction; in the Riccati path it triggers `renormalize_riccati_inplace!`. Default `1e4` empirically keeps max(|U₁|, |U₂|) in O(1)–O(10⁴) over the integration domain on DIII-D / Solovev sweeps; lower triggers excess renorms without accuracy gain, higher risks overflow before the next renorm.
   - `numsteps_init::Int` - Initial array size for ODE data storage
   - `numunorms_init::Int` - Initial array size for solution normalization data
diff --git a/src/ForceFreeStates/ForceFreeStates.jl b/src/ForceFreeStates/ForceFreeStates.jl
index 4505c5492..1a2643432 100644
--- a/src/ForceFreeStates/ForceFreeStates.jl
+++ b/src/ForceFreeStates/ForceFreeStates.jl
@@ -6,17 +6,6 @@ using LinearAlgebra.LAPACK
 using TOML
 using FFTW
 using OrdinaryDiffEq
-
-"""
-    el_ode_algorithm(ctrl) -> OrdinaryDiffEq algorithm
-
-The solver named by `ctrl.ode_solver`, instantiated for the Euler-Lagrange `solve` calls.
-"""
-function el_ode_algorithm(ctrl)
-    isdefined(OrdinaryDiffEq, Symbol(ctrl.ode_solver)) ||
-        error("ode_solver = \"$(ctrl.ode_solver)\" is not an OrdinaryDiffEq algorithm name (e.g. \"Vern7\", \"Vern9\", \"DP8\")")
-    return getfield(OrdinaryDiffEq, Symbol(ctrl.ode_solver))()
-end
 using HDF5
 using JLD2
 using FastInterpolations
@@ -26,6 +15,23 @@ using FastGaussQuadrature: gausslobatto
 using QuadGK: quadgk, quadgk!
 
 import ..Equilibrium
+
+"""
+Explicit Runge-Kutta methods benchmarked for the Euler-Lagrange solves; `ode_solver` must name one.
+"""
+const EL_ODE_SOLVERS = (Vern6=Vern6(), Vern7=Vern7(), Vern8=Vern8(), Vern9=Vern9(), DP8=DP8())
+
+"""
+    el_ode_algorithm(ctrl) -> OrdinaryDiffEq algorithm
+
+The solver named by `ctrl.ode_solver`, from [`EL_ODE_SOLVERS`](@ref).
+"""
+function el_ode_algorithm(ctrl)
+    name = Symbol(ctrl.ode_solver)
+    haskey(EL_ODE_SOLVERS, name) ||
+        throw(ArgumentError("ode_solver = \"$(ctrl.ode_solver)\" is not one of the benchmarked Euler-Lagrange solvers $(keys(EL_ODE_SOLVERS))"))
+    return EL_ODE_SOLVERS[name]
+end
 import ..Utilities
 import ..Vacuum
 import ..InnerLayer
diff --git a/src/ForceFreeStates/Riccati/Driver.jl b/src/ForceFreeStates/Riccati/Driver.jl
index c5869e443..fa7385691 100644
--- a/src/ForceFreeStates/Riccati/Driver.jl
+++ b/src/ForceFreeStates/Riccati/Driver.jl
@@ -26,7 +26,7 @@ Setting w = Q - K̄·S (shape N×N) and v = F̄⁻¹·w (Cholesky solve), this s
 
 `riccati_der!` evaluates the explicit Riccati RHS `dS/dψ = w†F̄⁻¹w − S·Ḡ·S` correctly,
 but this ODE is **quadratic** in S. Near a rational surface, S grows large, so the quadratic
-term `-SGS` dominates and the RHS grows as |S|². Explicit adaptive solvers (Vern9) use
+term `-SGS` dominates and the RHS grows as |S|². Explicit adaptive solvers use
 *relative* error control: they accept a step when |Δu|/|u| < reltol. When |S| is large,
 the absolute error |ΔS| can be enormous while the relative error stays within tolerance.
 The solver takes large steps through what is effectively a near-blowup — no amount of
@@ -40,8 +40,8 @@ recover S = U₁·U₂⁻¹ by renormalization. This achieves the same Riccati t
 **no accuracy loss**:
 
 - `sing_der!` evaluates the exact EL RHS — no approximation.
-- Vern9 integrates (U₁, U₂) to **9th-order accuracy** with the adaptive step-size
-  controller enforcing the configured reltol at every accepted step.
+- The explicit Runge-Kutta method (`ode_solver`, default Vern7) integrates (U₁, U₂) with its
+  adaptive step-size controller enforcing the configured tolerances at every accepted step.
 - Renormalization `S = U₁·U₂⁻¹` is **exact** (a change of variables, not an approximation).
 - The global error is the same as the standard EL path — controlled by the ODE solver
   reltol, not by the renormalization frequency.
@@ -56,12 +56,14 @@ solver — exactly analogous to Gaussian reduction in the standard ODE.
 
 To verify the method is consistent with the Riccati ODE, consider a single step from (S, I):
 
-  After one step: U₁_new = S + (A·S + B)·Δψ + O(Δψ²),  U₂_new = I + (C·S + D)·Δψ + O(Δψ²)
-  Renorm:         S_new = U₁_new · U₂_new⁻¹ = S + (B + A·S − S·D − S·C·S)·Δψ + O(Δψ²) ✓
+```
+After one step: U₁_new = S + (A·S + B)·Δψ + O(Δψ²),  U₂_new = I + (C·S + D)·Δψ + O(Δψ²)
+Renorm:         S_new = U₁_new · U₂_new⁻¹ = S + (B + A·S − S·D − S·C·S)·Δψ + O(Δψ²) ✓
+```
 
 The leading term matches the Riccati ODE exactly. This is a local consistency check only —
-it does not imply the integration is first-order. In practice Vern9 captures all higher-order
-terms through its internal stages, achieving 9th-order global accuracy at the configured reltol.
+it does not imply the integration is first-order: the Runge-Kutta stages capture the higher-order
+terms, to the method's order at the configured tolerances.
 
 ## Storage Convention
 
diff --git a/src/ForceFreeStates/Riccati/Propagators.jl b/src/ForceFreeStates/Riccati/Propagators.jl
index d888750b7..5ed395715 100644
--- a/src/ForceFreeStates/Riccati/Propagators.jl
+++ b/src/ForceFreeStates/Riccati/Propagators.jl
@@ -124,7 +124,10 @@ end
     riccati_der!(du, u, params, psieval)
 
 Evaluate the explicit dual Riccati ODE right-hand side:
+
+```
 dS/dψ = w†·F̄⁻¹·w - S·Ḡ·S,   w = Q - K̄·S
+```
 
 where Q = diag(1/(m - n·q)) is the diagonal singular factor matrix.
 The identity slice u[:,:,2] = I does not evolve (du[:,:,2] = 0).
@@ -274,8 +277,11 @@ end
     renormalize_riccati!(odet, intr)
 
 After a singular surface crossing, restore the canonical Riccati storage convention:
+
+```
 u[:,:,1] = S_new = U₁_new · U₂_new⁻¹
 u[:,:,2] = I
+```
 
 `riccati_cross_ideal_singular_surf!` leaves u[:,:,1] = U₁_new and u[:,:,2] = U₂_new (not I),
 so this step is required before continuing the Riccati integration.
@@ -465,8 +471,10 @@ Apply the chunk propagator `prop` to the current state `odet.u` in-place.
 
 The propagator acts as a linear map on the (U₁, U₂) pair:
 
+```
 U₁_new = block_upper_ic[:,:,1] · U₁_prev + block_lower_ic[:,:,1] · U₂_prev
 U₂_new = block_upper_ic[:,:,2] · U₁_prev + block_lower_ic[:,:,2] · U₂_prev
+```
 
 This correctly propagates any state (not just the identity), including the
 (S, I) form produced by Riccati-style crossings.
diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl
index f2ed4ef65..f42fdb561 100755
--- a/src/GeneralizedPerturbedEquilibrium.jl
+++ b/src/GeneralizedPerturbedEquilibrium.jl
@@ -85,8 +85,8 @@ using .Equilibrium: PlasmaEquilibrium
 using .ForcingTerms: RMPField
 
 const _DEPRECATED_FFS_KEYS = ("mer_flag", "force_wv_symmetry", "ode_flag", "cyl_flag", "mat_flag", "reform_eq_with_psilim",
-                              "use_riccati", "use_parallel", "parallel_threads", "populate_dense_xi",
-                              "gal_flag")
+    "use_riccati", "use_parallel", "parallel_threads", "populate_dense_xi",
+    "gal_flag")
 const _DEPRECATED_EQUIL_KEYS = ("power_bp", "power_b", "power_r", "power_rc")
 
 # Drop deprecated keys from a parsed gpec.toml section so legacy files keep parsing
@@ -204,6 +204,7 @@ function main_from_inputs(
     ffs_table = inputs["ForceFreeStates"]
     _drop_deprecated_keys!(ffs_table, _DEPRECATED_FFS_KEYS, "ForceFreeStates")
     ctrl = ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in ffs_table)...)
+    ForceFreeStates.el_ode_algorithm(ctrl)  # fail fast on an unknown ode_solver, before the equilibrium solve
 
     resolve_mode_space!(intr, ctrl)
 
@@ -740,10 +741,10 @@ runs are TOML-driven this cycle: `kinetic_factor > 0` needs the `[KineticForces]
 and errors here.
 
 ```julia
-eq   = PlasmaEquilibrium("input.geqdsk"; jac_type="hamada")
+eq = PlasmaEquilibrium("input.geqdsk"; jac_type="hamada")
 prob = EulerLagrangeProblem(eq; nn=1, delta_mlow=8, delta_mhigh=8, vac_flag=true)
-ffs  = solve(prob, Riccati())
-ffs  = solve(eq, Riccati(); nn=1, vac_flag=true)   # equivalent one-line form
+ffs = solve(prob, Riccati())
+ffs = solve(eq, Riccati(); nn=1, vac_flag=true)   # equivalent one-line form
 ```
 """
 function solve(prob::EulerLagrangeProblem, alg::ForceFreeStates.AbstractIntegrator)
@@ -754,6 +755,7 @@ function solve(prob::EulerLagrangeProblem, alg::ForceFreeStates.AbstractIntegrat
     ForceFreeStates._apply_alg!(ctrl_kwargs, alg)
     ForceFreeStates._apply_match!(ctrl_kwargs, prob.match, alg)
     ctrl = ForceFreeStatesControl(; ctrl_kwargs...)
+    ForceFreeStates.el_ode_algorithm(ctrl)  # fail fast on an unknown ode_solver
 
     ctrl.kinetic_factor > 0 &&
         error("kinetic runs (kinetic_factor > 0) need the [KineticForces] profiles and are TOML-driven; run them through `main`")
diff --git a/test/runtests_solve_api.jl b/test/runtests_solve_api.jl
index d04eb9065..234752bc1 100644
--- a/test/runtests_solve_api.jl
+++ b/test/runtests_solve_api.jl
@@ -160,7 +160,7 @@ using TOML
         @test d.terms[2].scale == 1.0 + 0.0im
         @test (im * a).scale == im
         @test (a * 3).scale == 3.0 + 0.0im
-        @test (a - b).terms[2].scale == -0.5 + 0.0im
+        @test (a-b).terms[2].scale == -0.5 + 0.0im
         @test (-a).scale == -1.0 + 0.0im
     end
 
@@ -202,5 +202,16 @@ using TOML
         @test_throws ErrorException solve(equil, Forward(); nn=1, dir_path=".", ffs_kwargs..., integrator="riccati")
         @test_throws ErrorException solve(equil, Riccati(); nn=1, dir_path=".", ffs_kwargs..., nchunks=8)
         @test_throws ErrorException solve(equil, Forward(); nn=1, dir_path=".", ffs_kwargs..., nn_low=2)
+        # An unbenchmarked ode_solver is refused before any integration starts.
+        @test_throws ArgumentError solve(equil, Forward(); nn=1, dir_path=".", ffs_kwargs..., ode_solver="Rodas5")
+    end
+
+    @testset "ode_solver resolves to the benchmarked explicit methods" begin
+        for name in keys(FFS.EL_ODE_SOLVERS)
+            alg = FFS.el_ode_algorithm(FFS.ForceFreeStatesControl(; ode_solver=String(name)))
+            @test nameof(typeof(alg)) === name
+        end
+        @test FFS.el_ode_algorithm(FFS.ForceFreeStatesControl()) isa GPEC.ForceFreeStates.OrdinaryDiffEq.Vern7
+        @test_throws ArgumentError FFS.el_ode_algorithm(FFS.ForceFreeStatesControl(; ode_solver="solve"))
     end
 end