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

src/ForceFreeStates/CoreTypes.jl

@@ -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

src/ForceFreeStates/EulerLagrange.jl

@@ -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/Riccati/Driver.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
 

test/runtests.jl

@@ -29,6 +29,7 @@ else
     include("./runtests_coordinate_invariant.jl")
     include("./runtests_eulerlagrange.jl")
     include("./runtests_riccati.jl")
+    include("./runtests_fixed_axis.jl")
     include("./runtests_parallel_integration.jl")
     include("./runtests_result_struct.jl")
     include("./runtests_solve_api.jl")

test/runtests_fixed_axis.jl

@@ -0,0 +1,59 @@
+# Regression test for the axis initial condition of the Euler-Lagrange integration.
+#
+# Free-boundary energies from the Frobenius free-axis start (`fixed_axis = false`, the
+# default before this test existed) were wrong on flat-core diverted equilibria: one
+# eigenvalue of the plasma response matrix W_p sat ~10x off its Fortran DCON value and,
+# where it crossed zero, `et[1]` reported a spurious -1e2 ... -1e5 "instability".  The
+# fixture is a 129x129 TokaMaker geqdsk of a synthetic DIII-D-like ramp-up slice
+# (q0 = 3.42, q95 = 5.9, diverted, no wall, n = 1) on which Fortran DCON (v1.5.5, same
+# mlow/mhigh/psilim/dmlim, mpsi 256) gives a stiff W_p eigenvalue of +1.609e4, no
+# negative W_p eigenvalue, and et[1] = +1.42; the old default gave et[1] = -8.4e4.
+#
+# The `@test_broken` lines pin the DEFECT of the Frobenius start: if `compute_axis_init`
+# is ever repaired to deliver its documented psilow^(|m|/2) limit, they report an
+# "unexpected pass" and should then be promoted to plain `@test`s.
+using LinearAlgebra
+using GeneralizedPerturbedEquilibrium.ForceFreeStates: ForceFreeStatesControl
+
+const FIXED_AXIS_FIXTURE = joinpath(@__DIR__, "test_data", "regression_rampup_fixed_axis")
+# Fortran DCON, GPEC v1.5.5-323, on the 257x257 parent of the fixture at mpsi 256 / mtheta 512.
+const WP_STIFF_DCON = 1.6092e4
+
+function run_fixed_axis_case(fixed_axis::Bool)
+    dir = mktempdir()
+    for f in readdir(FIXED_AXIS_FIXTURE)
+        cp(joinpath(FIXED_AXIS_FIXTURE, f), joinpath(dir, f))
+    end
+    toml = joinpath(dir, "gpec.toml")
+    s = read(toml, String)
+    s = replace(s, r"^fixed_axis.*\n"m => "")
+    s = replace(s, "[ForceFreeStates]\n" => "[ForceFreeStates]\nfixed_axis = $(fixed_axis)\n")
+    write(toml, s)
+    r = GeneralizedPerturbedEquilibrium.main([dir])
+    fb = r.ffs.free_boundary
+    wp = Matrix(fb.wp)
+    wpe = eigvals(Hermitian((wp + wp') / 2))
+    dp = real.(diag(Matrix(r.ffs.delta_prime.matrix)))
+    return (et1 = minimum(real.(fb.et)), wp_min = minimum(wpe), wp_max = maximum(wpe),
+            n_neg = count(<(0), wpe), dprime = dp)
+end
+
+@testset "Fixed-axis default (DCON axis condition)" begin
+    # The default itself is the fix: pin it.
+    @test ForceFreeStatesControl().fixed_axis == true
+
+    fixed = run_fixed_axis_case(true)
+    @test fixed.et1 > 0                      # physical fundamental, not a pole
+    @test fixed.n_neg == 0                   # W_p has no negative eigenvalue (DCON: none)
+    @test isapprox(fixed.wp_max, WP_STIFF_DCON; rtol=0.02)  # stiff eigenvalue on DCON's value (this fixture gives 1.6090e4; 2 % covers 129^2 / mpsi 128 vs DCON's 257^2 / mpsi 256)
+
+    # The Frobenius free-axis start, kept as an opt-in: record its defect.
+    free = run_fixed_axis_case(false)
+    @test_broken free.et1 > 0
+    @test_broken free.n_neg == 0
+    @test_broken isapprox(free.wp_max, WP_STIFF_DCON; rtol=0.02)
+
+    # The Delta' BVP does not depend on the axis condition.
+    @test length(fixed.dprime) == length(free.dprime)
+    @test isapprox(fixed.dprime, free.dprime; rtol=1e-3)
+end

test/test_data/README_test_data.md

@@ -1,3 +1,14 @@
 Put relevant data for various tests here, most likely outputs from the Fortran code used to validate Julia outputs
 
 # TODO: store as hdf5 files instead?
+
+## regression_rampup_fixed_axis/
+
+Fixture for `test/runtests_fixed_axis.jl`: a synthetic, diverted DIII-D-like ramp-up
+slice (TokaMaker free-boundary solve, 129×129 geqdsk, q0 = 3.42, q95 = 5.9, no wall,
+n = 1) on which the Frobenius free-axis start (`fixed_axis = false`) produces a spurious
+`et[1]` ≈ −1e5 while the DCON fixed-axis start (`fixed_axis = true`, the default) gives
+`et[1]` = +1.50 and a plasma-matrix stiff eigenvalue of +1.609e4, matching Fortran DCON
+(v1.5.5-323) on the 257×257 parent equilibrium. `gpec.toml` carries the production settings
+(mpsi 128, mtheta 256, mthvac 480, `set_psilim_via_dmlim = true`); the test injects
+`fixed_axis` itself. Analytic profiles, not experimental data.