matt-pharr · Oct 9, 2026 · Oct 8, 2026
diff --git a/examples/DIIID-like_ideal_multi_n_example/gpec.toml b/examples/DIIID-like_ideal_multi_n_example/gpec.toml
@@ -4,9 +4,6 @@
 # small edge singular factor. The q = 2 and q = 3 surfaces resonate at both n, so SingularCoupling
 # repeats them across rows. psihigh and delta_mhigh keep the case small for the regression harness;
 # every n block spans m = -8..20, which a single-n n = 1 run matches only with delta_mhigh = 8.
-# frobenius_psi_max = 0 uses the fixed axis start: the Frobenius start's m = 0 branch is picked by
-# roundoff, which makes the n = 1 anti-Hermitian part of W (and so the torque) differ between
-# otherwise identical runs.
 
 [Equilibrium]
 eq_filename = "../DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode.geqdsk" # Path to equilibrium file

diff --git a/src/ForceFreeStates/CoreTypes.jl b/src/ForceFreeStates/CoreTypes.jl
@@ -123,7 +123,7 @@ gpec.toml.
   - `mthvac::Int` - Number of vacuum poloidal grid points (corresponds to `mtheta` in VacuumInput)
   - `nzvac::Int` - Number of vacuum toroidal grid points (corresponds to `nzeta` in VacuumInput3D)
   - `sing_start::Int` - Start integration at the `sing_start`-th singular surface
-  - `frobenius_psi_max::Float64` - Largest starting ψ_N at which the Frobenius start is used. The Frobenius start [Glasser 2016 Eq. 51] is a power series about the magnetic axis, so an integration beginning above this value (an interior start from a raised `psilow`, or `qlow` above q₀) instead starts every solution with zero displacement and unit "momentum" (U₁ = 0, U₂ = I), the initialization Fortran DCON always uses, and warns that it did. `0` selects the fixed start everywhere, silently. Default 0.01, which keeps every axis-started case on the Frobenius start.
+  - `frobenius_psi_max::Float64` - Largest starting ψ_N at which the opt-in Frobenius start is used. Default `0` selects the fixed start everywhere and silently: U₁ = 0, U₂ = I, the initialization Fortran DCON uses. A positive value uses the diagonal Frobenius 2×2 when the integration starts at or below it, and warns that this start is currently unreliable. A start above the threshold uses the fixed start with no warning. That 2×2 does not recover a displacement that vanishes as a positive power of the starting ψ, and its m = 0 branch is an imaginary pair whose sign follows roundoff. Leave this at 0 unless the Frobenius start itself is under test.
   - `nn_low::Int` - Lower bound for toroidal modes
   - `nn_high::Int` - Upper bound for toroidal modes
   - `delta_mlow::Int` - Expands lower bound of Fourier harmonics by delta_mlow
@@ -221,5 +221,5 @@ gpec.toml.
     gal_rho::Vector{Float64} = Float64[]      # per-surface mass density ρ [kg/m³] (length msing, core→edge); Fortran rmatch `massden`
     gal_rotation::Vector{Float64} = Float64[] # per-surface rotation frequency f [Hz] (length msing, core→edge); forced eigenvalue γ_s = 2πi·n·f. Fortran rmatch `rotation`
     gal_gamma::Float64 = 5 / 3       # ratio of specific heats Γ for the resistive-layer coefficients (resist_eval G term)
-    frobenius_psi_max::Float64 = 0.01
+    frobenius_psi_max::Float64 = 0.0
 end
diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl
@@ -45,49 +45,81 @@ and a small set of temporary matrices and factors used to compute singular-layer
   - `numpert_total::Int` - Total number of Fourier mode combinations (m × n) used in the calculation.
 
   - `numunorms_init::Int` - Initial allocation size for the number of normalization operations recorded.
+
   - `msing::Int` - Number of singular surfaces in the equilibrium (used to size asymptotic coefficient arrays).
+
   - `numsteps_init::Int` - Initial allocation size for the number of integration steps to store.
+
   - `step::Int` - Current integration step index (1-based, like `istep` in the original Fortran).
+
   - `psi_store::Vector{Float64}` - Stored psi values at each saved integration step (length `numsteps_init`).
+
   - `q_store::Vector{Float64}` - Stored q values at each saved integration step (length `numsteps_init`).
+
   - `u_store::Array{ComplexF64,4}` - Stored solution arrays at each saved step with shape
     `(numpert_total, numpert_total, 2, numsteps_init)` (complex solution state used by the solver).
+
   - `du_store::Array{ComplexF64,3}` - dΞ_ψ/dψ (the u₁ block only) at each saved step, shape
     `(numpert_total, numpert_total, step)`. Empty until `materialize_derivative_stores!` fills it,
     except on the galerkin-matched path which supplies the analytic derivative at construction.
     du₂/dψ is never stored densely — its only consumer evaluates it on demand at bracket nodes.
+
   - `xi_s_store::Array{ComplexF64,3}` - Clebsch displacement Ξ_s at each saved step, eq. 18 of Glasser 2016,
     shape `(numpert_total, numpert_total, step)`. Empty until materialized, same as `du_store`.
+
   - `u_store_el_basis::Bool` - True when `u_store` holds the Euler-Lagrange state `(u₁, u₂)`, so the
     derivative kernel can be re-applied to it. False on the sparse parallel path, whose stored columns
     are chunk-endpoint Riccati matrices; `materialize_derivative_stores!` refuses to run there.
+
   - `du_store_populated::Bool` - True once `du_store`/`xi_s_store` hold valid data in the final
     (post-transform, post-normalization) basis. Set by `materialize_derivative_stores!` or by the
     galerkin-matched constructor; stays false where the stores cannot be materialized, e.g. the
     sparse parallel path whose solution is in the Riccati basis.
+
   - `crit_store::Vector{Float64}` - Stored crit parameter values (smallest eigenvalue of W⁻ꜝ) (length `numsteps_init`).
+
   - `ca_r::Array{ComplexF64,4}` - Asymptotic coefficients just to the right of each singular surface
     with shape `(numpert_total, numpert_total, 2, msing)`.
+
   - `ca_l::Array{ComplexF64,4}` - Asymptotic coefficients just to the left of each singular surface
     with shape `(numpert_total, numpert_total, 2, msing)`.
+
   - `ca_populated::Bool` - True once an ideal singular-surface crossing has filled `ca_l`/`ca_r`; kinetic and
     galerkin-matched runs never populate them and leave this false, and the HDF5 writer then emits zero-extent
     `ca_left`/`ca_right` datasets instead of unpopulated arrays.
+
   - `edge_scan::EdgeScanState` - Edge dW scan state and results. Initialized as a disabled sentinel (N_edge=0) and replaced by `findmax_dW_edge!` when a scan runs.
+
   - `psifac::Float64` - Current normalized flux coordinate for the integrator.
+
   - `q::Float64` - Safety factor value at `psifac` (current q during integration).
+
   - `u::Array{ComplexF64,3}` - Current working solution arrays with shape `(numpert_total, numpert_total, 2)`.
+
   - `ising_start::Int` - Index of the starting singular surface to be crossed during integration.
+
   - `psimax::Float64` - Maximum psi value for which the integrator is allowed to run in next integration region.
+
   - `needs_crossing::Bool` - Flag indicating whether a rational surface needs to be crossed after the current integration region.
+
+    # Initialization parameters
+
   - `nzero::Int` - Count of detected zero crossings (used for diagnostics).
+
   - `new::Bool` - Flag indicating whether a new `unorm0` should be computed after a fixup.
 
     # Initialization parameters
+
+    # Saved data throughout integration
+
   - `unorm::Vector{Float64}` - Current norms of the solution vectors (length `numpert_total`).
+
+# Total ODE solver steps taken (all steps, not just saved ones)
+
   - `unorm0::Vector{Float64}` - Reference/initial norms of the solution vectors (length `numpert_total`).
 
     # Saved data throughout integration
+
   - `ifix::Int` - Number of normalization operations performed (index into normalization arrays).
 
 # Total ODE solver steps taken (all steps, not just saved ones)
@@ -96,8 +128,15 @@ and a small set of temporary matrices and factors used to compute singular-layer
 
   - `sing_flag::Vector{Bool}` - Boolean flags indicating which stored normalizations correspond to singular solutions    # Edge dW scan state and results (disabled sentinel when psiedge >= psilim, i.e. no edge scan)
     (length `numunorms_init`).
+
+    # Edge dW scan state and results (disabled sentinel when psiedge >= psilim, i.e. no edge scan)
+
+    # Initialization parameters
+
   - `zeroed_idx::Vector{Vector{Int}}` - For each ideal rational surface jump, a vector of indices of solutions that were zeroed.    # Data for integrator
-  - `fixfac::Array{ComplexF64,3}` - Fix-up factors for Gaussian reduction with shape `(numpert_total, numpert_total, numunorms_init)`.
+
+  - `fixfac::Array{ComplexF64,3}` - Fix-up factors for Gaussian reduction with shape `(numpert_total, numpert_total, numunorms_init)`.    # Data for integrator
+
   - `fixstep::Vector{Int64}` - Step indices (psi step positions) at which normalization/fixups were performed (length `numunorms_init`).
 """
 @kwdef mutable struct OdeState
@@ -404,20 +443,15 @@ end
 """
     compute_axis_init(mats, profiles, intr, psi_low) -> (U1_init, U2_init)
 
-Compute axis initial conditions for the Euler-Lagrange ODE via the Frobenius
-leading-coefficient eigenvalue problem [Glasser Phys. Plasmas 2016 112506 Eq. 51]:
-
-    lim_{ψ→0} [ψ M(ψ) − a I] v = 0
-
-For each mode j, solves the 2×2 Frobenius eigenvalue problem for the diagonal block of
-A₀ = ψ_low · M(ψ_low). The eigenvector with Re(a) ≥ 0 is the regular (non-singular)
-Frobenius solution. Returns U₁_init and U₂_init normalized so that U₂_init = I (consistent
-with the N independent solutions convention).
+Opt-in diagonal approximation to the axis Frobenius problem [Glasser Phys. Plasmas 2016 112506 Eq. 51],
+evaluated at finite `psi_low` rather than in the limit ψ → 0. Used only when `frobenius_psi_max` is
+positive and the integration starts at or below it.
 
-For m≠0 the regular eigenvector has a negligible U₁ component (~ψ_low^(|m|/2)), recovering
-the Glasser [0, I] limit as ψ_low → 0. For m=0 (degenerate a≈0), the regular eigenvector
-is identified by dominant |U₁| component, giving the physically correct constant-displacement
-Frobenius solution and avoiding the spurious logarithmic irregularity.
+For each mode, forms the 2×2 diagonal block of ψ_low · M from the diagonal entries of F̄⁻¹, K̄, and Ḡ,
+discards mode coupling, and returns that block's eigenvector normalized so that U₂ = I (or U₁ = 1, U₂ = 0
+when the momentum component underflows). The component ratio does not scale as a positive power of
+`psi_low`. A negative Ḡ diagonal, including the m = 0 block, gives an imaginary conjugate pair, so the
+returned U₁ is non-real and which branch is stored follows the eigen-solver's column order.
 """
 function compute_axis_init(mats::MatrixSplines, profiles::Equilibrium.ProfileSplines,
     intr::ForceFreeStatesInternal, psi_low::Float64)
@@ -495,12 +529,9 @@ end
 
 Initialize the OdeState struct for the case of sing_start = 0 (axis initialization).
 Formerly `ode_axis_init!`. This now only initializes `psifac`, `ising_start`, and `u`.
-The Frobenius start is used only when the starting surface lies at or below `ctrl.frobenius_psi_max`;
-an interior start (or `ctrl.frobenius_psi_max = 0`) uses U₁ = 0, U₂ = I.
-
-### TODOs
-
-Move ising_start logic to chunk_el_integration_bounds?
+The default `ctrl.frobenius_psi_max = 0` uses the fixed start, U₁ = 0, U₂ = I, with no warning. The
+Frobenius start is used only when `frobenius_psi_max` is positive and the starting surface lies at or
+below it.
 """
 function initialize_el_at_axis!(odet::OdeState, ctrl::ForceFreeStatesControl, mats::MatrixSplines,
     profiles::Equilibrium.ProfileSplines, intr::ForceFreeStatesInternal)
@@ -533,23 +564,19 @@ function initialize_el_at_axis!(odet::OdeState, ctrl::ForceFreeStatesControl, ma
         odet.ising_start = searchsortedfirst(getfield.(intr.sing, :psifac), odet.psifac) - 1
     end
 
-    # The Frobenius start is a series about the axis; an interior start uses the fixed start Fortran DCON always uses.
+    # Default (frobenius_psi_max = 0) is the fixed start Fortran DCON uses. A positive threshold opts into the diagonal Frobenius start.
     fixed_start = odet.psifac > ctrl.frobenius_psi_max
-    if fixed_start && ctrl.frobenius_psi_max > 0
-        @warn "Integration starts at ψ_N = $(round(odet.psifac; sigdigits=4)), above frobenius_psi_max = $(ctrl.frobenius_psi_max): " *
-              "the axis Frobenius start is only valid near the axis, so the fixed start (U₁ = 0, U₂ = I) is used, as in Fortran DCON. " *
-              "Set frobenius_psi_max = 0 to select it silently, or raise it to keep the Frobenius start."
-    end
     if fixed_start
         # Original Glasser initialization: U₁=0, U₂=I [Glasser 2016 §VI].
         # Constrains the displacement ξ^ψ=0 for all modes at the starting surface (fixed axis or rigid core).
         for ipert in 1:intr.numpert_total
             odet.u[ipert, ipert, 2] = 1
         end
     else
-        # Frobenius initialization [Glasser 2016 §VI Eq. 51]: selects the regular
-        # (non-logarithmic) solution for each mode, including the correct constant
-        # displacement solution for the degenerate m=0 case (free magnetic axis).
+        # Opt-in diagonal Frobenius start. See compute_axis_init.
+        @warn "Frobenius axis start at ψ_N = $(round(odet.psifac; sigdigits=4)). " *
+              "This start is currently unreliable - use only if explicitly requested." *
+              "Set frobenius_psi_max = 0 for the fixed start (U₁ = 0, U₂ = I)."
         U1_init, U2_init = compute_axis_init(mats, profiles, intr, odet.psifac)
         odet.u[:, :, 1] .= U1_init
         odet.u[:, :, 2] .= U2_init
@@ -897,7 +924,6 @@ making it clear what region is being integrated.
   - `mats::MatrixSplines` - Fourier fit variables
   - `intr::ForceFreeStatesInternal` - Internal data
   - `chunk::IntegrationChunk` - Integration chunk containing start and end ψ for integration
-
 """
 function integrate_el_region!(
     odet::OdeState,
@@ -1177,7 +1203,7 @@ function transform_u!(odet::OdeState, intr::ForceFreeStatesInternal)
                     temp[ksol, jsol] = odet.fixfac[ksol, jsol, ifix]
                 end
             end
-            mul!(gauss_buffer, view(gauss, :, :, ifix), temp)
+            mul!(gauss_buffer, view(gauss,:,:,ifix), temp)
             gauss[:, :, ifix] .= gauss_buffer
         end
         # Account for zeroed indices at singular surfaces in `ode_ideal_cross`
@@ -1194,7 +1220,7 @@ function transform_u!(odet::OdeState, intr::ForceFreeStatesInternal)
     # and mfix + 1 is the for the region after the last fixup and before the edge
     transforms[:, :, end] .= identity
     for ifix in odet.ifix:-1:1
-        mul!(view(transforms, :, :, ifix), view(gauss, :, :, ifix), view(transforms, :, :, (ifix + 1)))
+        mul!(view(transforms,:,:,ifix), view(gauss,:,:,ifix), view(transforms,:,:,(ifix+1)))
     end
 
     # Now that we have the transform matrices, we can apply them to the solution vectors

diff --git a/src/ForceFreeStates/Riccati/Driver.jl b/src/ForceFreeStates/Riccati/Driver.jl
@@ -76,10 +76,10 @@ After renormalization (at crossing or when norms exceed ucrit):
 This is compatible with downstream code (which uses U₁/U₂ ratio):
   - Free.jl:     wp = u[:,:,2] / u[:,:,1] = I · S⁻¹ = P  ✓  (post-renorm)
   - FixedBoundaryStability.jl: crit = min_eigval(u[:,:,1] / u[:,:,2]) = min_eigval(S)  ✓
-  - Axis init:   determined by `ctrl.frobenius_psi_max`. When the start lies above it (or it is 0), U₁=0, U₂=I → S(ψ₀)=0 (original
-    Glasser fixed-axis BC). Otherwise (default near the axis), Frobenius eigenvalue init [Glasser 2016 Eq. 51]
-    sets U₂=I and U₁ to the regular Frobenius eigenvector per mode → S(ψ₀) = U₁_Frobenius is
-    nonzero in general. Riccati S-evolution remains well-defined either way.
+  - Axis init:   determined by `ctrl.frobenius_psi_max`. The default 0, and any start above a positive threshold,
+    uses U₁=0, U₂=I → S(ψ₀)=0 (Fortran DCON's fixed-axis BC). A positive threshold at or above the start
+    opts into the diagonal Frobenius init, which sets U₂=I and a per-mode U₁, so S(ψ₀) is nonzero.
+    Riccati S-evolution remains well-defined either way.
 
 ## Key Differences from Standard Integration
 

diff --git a/test/runtests_eulerlagrange.jl b/test/runtests_eulerlagrange.jl
@@ -136,8 +136,10 @@ end
     @testset "truncate_chunks! cuts the Riccati chunks back to the edge" begin
         F = GeneralizedPerturbedEquilibrium.ForceFreeStates
         # Crossing chunks end just short of surfaces 1-3; the edge drops surface 3 (msing = 2).
-        chunks() = [F.IntegrationChunk(; psi_start=a, psi_end=b, needs_crossing=i <= 3, ising=i <= 3 ? i : 0)
-                    for (i, (a, b)) in enumerate(((0.10, 0.29), (0.31, 0.59), (0.61, 0.89), (0.91, 0.99)))]
+        chunks() = [
+            F.IntegrationChunk(; psi_start=a, psi_end=b, needs_crossing=i <= 3, ising=i <= 3 ? i : 0)
+            for (i, (a, b)) in enumerate(((0.10, 0.29), (0.31, 0.59), (0.61, 0.89), (0.91, 0.99)))
+        ]
         cs = chunks()
         @test F.truncate_chunks!(cs, 0.75, 2) == 3  # straddling chunk shortened, so re-integrate it
         @test length(cs) == 3 && cs[3].psi_end == 0.75 && count(c -> c.needs_crossing, cs) == 2
@@ -405,7 +407,7 @@ end
         @test length(odet.unorm0) == numpert_total
     end
 
-    @testset "interior start falls back to the fixed initialization" begin
+    @testset "axis start defaults to the fixed initialization" begin
         FFS = GeneralizedPerturbedEquilibrium.ForceFreeStates
         ex = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example")
         inputs = TOML.parsefile(joinpath(ex, "gpec.toml"))
@@ -430,19 +432,21 @@ end
             odet = FFS.OdeState(intr.numpert_total, ctrl.numsteps_init, ctrl.numunorms_init, intr.msing)
             return odet, ctrl, mats, equil, intr
         end
-        # Near the axis the Frobenius start gives the regular solution: U₂ = I with a nonzero U₁.
+        # The default threshold is 0: near the axis the fixed start is used, with no warning.
         odet, ctrl, mats, equil, intr = axis_state(1e-4)
-        FFS.initialize_el_at_axis!(odet, ctrl, mats, equil.profiles, intr)
-        @test odet.u[:, :, 2] ≈ I
-        @test any(!iszero, odet.u[:, :, 1])
-        # An interior start switches to the fixed start (U₁ = 0, U₂ = I) and says so.
-        odet, ctrl, mats, equil, intr = axis_state(0.3)
-        @test_logs (:warn, r"fixed start") FFS.initialize_el_at_axis!(odet, ctrl, mats, equil.profiles, intr)
+        @test ctrl.frobenius_psi_max == 0
+        @test_logs FFS.initialize_el_at_axis!(odet, ctrl, mats, equil.profiles, intr)
         @test iszero(odet.u[:, :, 1]) && odet.u[:, :, 2] ≈ I
-        # The threshold is a control: raising it keeps the Frobenius start, and zero selects the fixed start silently.
-        odet, ctrl, mats, equil, intr = axis_state(0.3; frobenius_psi_max=0.5)
+        # A start above a positive threshold uses the fixed start, still with no warning.
+        odet, ctrl, mats, equil, intr = axis_state(0.3; frobenius_psi_max=0.01)
         @test_logs FFS.initialize_el_at_axis!(odet, ctrl, mats, equil.profiles, intr)
+        @test iszero(odet.u[:, :, 1]) && odet.u[:, :, 2] ≈ I
+        # Raising the threshold to cover the start opts into the Frobenius start and warns.
+        odet, ctrl, mats, equil, intr = axis_state(0.3; frobenius_psi_max=0.5)
+        @test_logs (:warn, r"Frobenius axis start") FFS.initialize_el_at_axis!(odet, ctrl, mats, equil.profiles, intr)
+        @test odet.u[:, :, 2] ≈ I
         @test any(!iszero, odet.u[:, :, 1])
+        # An explicit zero selects the fixed start silently even near the axis.
         odet, ctrl, mats, equil, intr = axis_state(1e-4; frobenius_psi_max=0.0)
         @test_logs FFS.initialize_el_at_axis!(odet, ctrl, mats, equil.profiles, intr)
         @test iszero(odet.u[:, :, 1]) && odet.u[:, :, 2] ≈ I
@@ -613,7 +617,7 @@ end
         # materialize after the Gaussian fixups and free-boundary normalization rather than
         # transforming stored derivatives alongside u_store.
         npert = intr.numpert_total
-        T = Matrix{ComplexF64}(I, npert, npert) .+ 0.25 .* ComplexF64.(reshape(sin.(1:npert^2), npert, npert))
+        T = Matrix{ComplexF64}(I, npert, npert) .+ 0.25 .* ComplexF64.(reshape(sin.(1:(npert^2)), npert, npert))
         odet_t = deepcopy(odet_pristine)
         for istep in 1:odet_t.step
             odet_t.u_store[:, :, 1, istep] = odet_t.u_store[:, :, 1, istep] * T