ForceFreeStates - BUGFIX! - Default to the DCON fixed-axis initial condition

`fixed_axis` now defaults to `true`: the Euler-Lagrange integration starts from
U1 = 0, U2 = I at psilow (xi^psi = 0 at the magnetic axis), the condition used by
Fortran DCON (dcon/ode.f, ode_axis_init) and Glasser 2016 Sec. VI.

The previous default (`false`) started from the free-axis Frobenius state of
`compute_axis_init`. At practical psilow (1e-2 ... 1e-4) that state has U1 entries
of order 0.1-10 with mixed signs instead of the documented psilow^(|m|/2) limit.
Started from it, the signed critical eigenvalue of W_p^-1 is displaced from the
first stored step onward and never recovers: one eigenvalue of the plasma
response matrix W_p sits ~10x off its DCON value in every equilibrium tested,
and where it crosses zero the free-boundary energy et[1] reports a spurious
-1e2 ... -1e5 "instability". Verified against Fortran DCON on identical
geqdsks (diverted DIII-D-like, n = 1): with the fixed-axis start crit(psi)
agrees to a median relative error of 8e-5, the stiff W_p eigenvalue to ~1e-4,
and et[1] returns to the physical fundamental. The Delta' BVP is unaffected
(<= 1e-4 relative). Integrator (riccati/forward), tolerances, coefficient
matrices and the free-boundary assembly were all exonerated by measurement.

The Frobenius start is retained behind `fixed_axis = false` for comparison; its
docstring now records the measured discrepancy. Two figures documenting the
defect and the fix are added under docs/src/assets/.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>

Author: d-burg

docs/src/assets/fixed_axis_crit_trace_vs_dcon.png

@@ -152,6 +152,7 @@ gpec.toml.
   - `integrator::String` - Which formalism integrates the Euler-Lagrange system. `"forward"` sweeps the plasma serially with Gaussian reduction and returns `u_store` / `du_store` / `xi_s_store` dense in the axis (EL) basis — the only convention PerturbedEquilibrium and FieldReconstruction consume correctly, and the only path that supports `kinetic_factor > 0`. `"riccati"` (default) runs the chunked fundamental-matrix propagator driver (Glasser 2018 Phys. Plasmas 25, 032507): chunks are integrated independently from identity initial conditions and assembled serially with Riccati-style crossings, which is the only way to obtain the singular-surface Δ' matrix for the tearing-mode solvers downstream, but leaves `u_store` as sparse chunk-endpoint Riccati states, so dense ξ profiles are unavailable. `"galerkin"` solves the same Euler-Lagrange system variationally instead of by radial ODE integration — the RDCON outer-region singular Galerkin method (Glasser, Wang & Park 2016 Phys. Plasmas 23, 112506), which discretizes the displacement on packed Hermite-cubic elements and solves one global banded system — producing the resistive Δ′ matrix and, when `gal_match_flag` is set, the RPEC inner-layer-matched ξ; it computes its own vacuum response and returns no free-boundary energies, and does not support `kinetic_factor > 0`. Requires `singfac_min != 0` for `"riccati"`.
   - `nchunks::Int` - Target number of Riccati integration chunks. `0` (the default) derives the count from problem structure alone: `max(2·msing + 3, 8·(msing + 1) + msing)`, enough sub-chunks per segment to keep the accumulated propagator products well-conditioned. An explicit value below `2·msing + 3` is clamped up with a warning. Chunk sizing never consults `Threads.nthreads()`, so Riccati outputs are identical whatever thread count `julia -t` provides; threads only change wall-clock.
   - `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.
+  - `fixed_axis::Bool` - Axis initial condition of the Euler-Lagrange integration. `true` (default) starts from U₁ = 0, U₂ = I at `psilow`, i.e. ξ^ψ = 0 at the magnetic axis — the condition used by Fortran DCON (`dcon/ode.f`, `ode_axis_init`) and by Glasser 2016 §VI. `false` starts from the free-axis Frobenius state returned by [`compute_axis_init`](@ref) (Glasser 2016 Eq. 51). The Frobenius state is meant to reduce to the [0, I] limit as `psilow → 0`, but at the `psilow` values used in practice (1e-2 … 1e-4) it returns U₁ entries of order 0.1–10 with mixed signs. Started from that state, the signed critical eigenvalue of W_p⁻¹ (the DCON `crit`) is displaced from the axis onward and never recovers: one eigenvalue of the plasma response matrix W_p stays ~10× off its DCON value in every equilibrium tested, and where it crosses zero the free-boundary energy `et[1]` reports a spurious −10²…−10⁵ "instability" (diverted DIII-D-like n = 1 scans, verified against Fortran DCON on identical geqdsks). With `true` the `crit(ψ)` trace, W_p spectrum and `et` agree with DCON. The Δ′ BVP is unaffected either way (≤ 1e-4 relative). Set `false` only to reproduce pre-existing free-axis results.
 """
 @kwdef struct ForceFreeStatesControl
     verbose::Bool = true
@@ -220,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)
-    fixed_axis::Bool = false
+    fixed_axis::Bool = true          # DCON axis condition U₁=0, U₂=I (ξ^ψ=0 at the axis). `false` = Frobenius free-axis init (compute_axis_init), opt-in only; see docstring.
 end

docs/src/assets/fixed_axis_scan_et1_wp.png

@@ -477,6 +477,15 @@ For m≠0 the regular eigenvector has a negligible U₁ component (~ψ_low^(|m|/
 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.
+
+!!! warning "Not the default"
+    Measured at `psilow = 0.01` on a diverted DIII-D-like equilibrium (n = 1, m = −14…27) the
+    returned U₁ diagonal is 0.11…0.23 for m = −14…−3 and +6.1, −9.8, −13.7 for m = 3, 4, 5 —
+    not the ~ψ_low^(|m|/2) limit stated above — and the downstream signed critical eigenvalue
+    of W_p⁻¹ then disagrees with Fortran DCON from the first stored step onward (see the
+    `fixed_axis` docstring in `CoreTypes.jl`). Whether the discrepancy is in the diagonal 2×2
+    truncation of A₀, in the regular-branch selection, or in the asymptotics claimed here has
+    not been established. Used only when `fixed_axis = false`.
 """
 function compute_axis_init(mats::MatrixSplines, profiles::Equilibrium.ProfileSplines,
         intr::ForceFreeStatesInternal, psi_low::Float64)
@@ -591,16 +600,20 @@ function initialize_el_at_axis!(odet::OdeState, ctrl::ForceFreeStatesControl, ma
     end
 
     if ctrl.fixed_axis
-        # Original Glasser initialization: U₁=0, U₂=I [Glasser 2016 §VI].
-        # Constrains the axis displacement ξ^ψ=0 for all modes (fixed magnetic axis).
-        # Retained as a reference/comparison option; the default (fixed_axis=false) is Frobenius.
+        # Default. Glasser initialization: U₁=0, U₂=I [Glasser 2016 §VI] — the DCON axis
+        # condition (dcon/ode.f, ode_axis_init): ξ^ψ=0 for all modes (fixed magnetic axis).
+        # This is the condition against which the free-boundary energies are validated;
+        # see the `fixed_axis` docstring in CoreTypes.jl for why the Frobenius state below
+        # is not the default.
         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. Frobenius initialization [Glasser 2016 §VI Eq. 51]: selects the regular
+        # (non-logarithmic) solution for each mode, including the constant-displacement
+        # solution for the degenerate m=0 case (free magnetic axis). At practical psilow
+        # the returned U₁ is O(0.1–10), not the documented ψ_low^(|m|/2) limit, and the
+        # resulting crit(ψ) / W_p disagree with DCON — retained for comparison only.
         U1_init, U2_init = compute_axis_init(mats, profiles, intr, odet.psifac)
         odet.u[:, :, 1] .= U1_init
         odet.u[:, :, 2] .= U2_init

src/ForceFreeStates/CoreTypes.jl

@@ -76,10 +76,11 @@ 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.fixed_axis`. When `true`, U₁=0, U₂=I → S(ψ₀)=0 (original
-    Glasser fixed-axis BC). When `false` (default), 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.fixed_axis`. When `true` (default), U₁=0, U₂=I → S(ψ₀)=0
+    (Glasser fixed-axis BC, identical to Fortran DCON). When `false` (opt-in), 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, but only the fixed-axis start reproduces DCON's crit(ψ) and W_p.
 
 ## Key Differences from Standard Integration