ForceFreeStates - BUGFIX! - Use trace-aware Frobenius exponents at resonant surfaces

src/ForceFreeStates/CoreTypes.jl

@@ -153,7 +153,6 @@ gpec.toml.
   - `force_termination::Bool` - Terminate after force-free states (skip perturbed equilibrium calculations)
   - `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.
 """
 @kwdef struct ForceFreeStatesControl
     verbose::Bool = true
@@ -191,7 +190,6 @@ gpec.toml.
     force_termination::Bool = false
     integrator::String = "riccati"   # Default: unlocks SingularSurfaces/Delta_prime_matrix (STRIDE BVP Δ′ matrix) used by SLAYER/GGJ downstream. Use "forward" for dense ξ (PerturbedEquilibrium) or kinetic runs.
     nchunks::Int = 0                 # Riccati chunk-count target; 0 = auto (derived from msing alone, never from Threads.nthreads()).
-    extended_precision_bvp::Bool = true   # Promote Δ' BVP to Complex{Double64}; default on (Float64 drifts the imaginary Δ' by 2–5× on DIIID-class cases).
 
     # --- RDCON outer-region Galerkin Δ′ solver (gal_solve port); selected by integrator = "galerkin" ---
     gal_solver::String = "auto"     # "auto" (banded Cholesky if the matrix is positive definite, else LU), "LU" (zgbtrf/zgbtrs) or "cholesky" (zpbtrf/zpbtrs, error if not positive definite)

src/ForceFreeStates/EulerLagrange.jl

@@ -769,7 +769,7 @@ function cross_ideal_singular_surf!(
     # Compute direction-specific asymptotic power series for this singular surface
     singp = intr.sing[ising]
     sing_asymp_right = compute_sing_asymptotics(singp, ctrl, equil, mats, intr; sig=1.0)
-    sing_asymp_left = compute_sing_asymptotics(singp, ctrl, equil, mats, intr; sig=-1.0, alpha_override=sing_asymp_right.alpha)
+    sing_asymp_left = compute_sing_asymptotics(singp, ctrl, equil, mats, intr; sig=-1.0)
     dpsi = singp.psifac - odet.psifac # ψ_res - ψ (positive)
 
     # Get asymptotic coefficients before crossing (left side)

src/ForceFreeStates/Galerkin/GalerkinSolve.jl

@@ -39,7 +39,9 @@ function gal_make_arrays!(ws::GalWorkspace, ctrl::ForceFreeStatesControl, equil,
     return ws
 end
 
-"""Empty `GalerkinResult` for a domain with no resonant surfaces."""
+"""
+Empty `GalerkinResult` for a domain with no resonant surfaces.
+"""
 function empty_galerkin_result()
     return GalerkinResult(0, Float64[], Float64[], Int[], Int[], Float64[], ComplexF64[], nothing, nothing)
 end
@@ -75,22 +77,20 @@ function galerkin_solve(ctrl::ForceFreeStatesControl, equil, mats::MatrixSplines
 
     ctrl.verbose && @info "Starting outer-region Galerkin Δ′ solve (msing=$msing, solver=$(ctrl.gal_solver))"
 
-    # Per-surface two-sided asymptotic series matching Fortran sing.f vmatr/vmatl:
-    # right = sig=+1, left = sig=-1, no √det normalization. The Mercier exponent α is a
-    # property of the surface, so the left series reuses the right's α (alpha_override). The order is
-    # raised to gal_sing_order + ceil(2·Re(α)) for high-Mercier-index surfaces (Fortran sing1_vmat).
+    # Per-surface two-sided asymptotic series (sig = ±1). Each side uses its own resonant block.
     asymps = GalSingAsymp[]
     for s in sings
         sing_order = ctrl.gal_sing_order
         ar = compute_sing_asymptotics(s, ctrl, equil, mats, intr; sig=1.0, sing_order=sing_order)
         if ctrl.gal_sing_order_ceiling
+            # The order is raised for high-Mercier-index surfaces
             order = ctrl.gal_sing_order + ceil(Int, 2 * real(ar.alpha[1]))
             if order > ctrl.gal_sing_order
                 sing_order = order
                 ar = compute_sing_asymptotics(s, ctrl, equil, mats, intr; sig=1.0, sing_order=sing_order)
             end
         end
-        al = compute_sing_asymptotics(s, ctrl, equil, mats, intr; sig=-1.0, alpha_override=ar.alpha, sing_order=sing_order)
+        al = compute_sing_asymptotics(s, ctrl, equil, mats, intr; sig=-1.0, sing_order=sing_order)
         push!(asymps, GalSingAsymp(ar, al))
     end
 
@@ -158,8 +158,7 @@ function galerkin_solve(ctrl::ForceFreeStatesControl, equil, mats::MatrixSplines
 
     Ap, Bp, Gammap, Deltap = gal_pest3_blocks(delta, msing)
 
-    # Mercier index D_I = -Re(α²); α is the small/large solution exponent (a surface property, taken
-    # from the right series — the left series reuses the same α via alpha_override above).
+    # Mercier index D_I = -Re(α²). α is the half-separation of the resonant exponents.
     di = [real(-asymps[i].right.alpha[1]^2) for i in 1:msing]
     alpha = [asymps[i].right.alpha[1] for i in 1:msing]
 
@@ -206,7 +205,7 @@ function gal_scaled_lu_solve!(mat::Matrix{ComplexF64}, sol::AbstractVecOrMat{Com
     n = size(mat, 2)
     off = kl + ku + 1  # band row holding the diagonal
     d = [iszero(mat[off, j]) ? 1.0 : 1 / sqrt(abs(mat[off, j])) for j in 1:n]
-    for j in 1:n, i in max(1, j - ku):min(n, j + kl)
+    for j in 1:n, i in max(1, j-ku):min(n, j+kl)
         mat[off+i-j, j] *= d[i] * d[j]
     end
     sol .*= d
@@ -311,11 +310,19 @@ const GALERKIN_H5_ANNOTATIONS = [
     "ForceFreeStates/Solutions/GalerkinIntegration/q" =>
         (; long_name="safety factor q on the Galerkin solution grid", dims=("psi",), attach=(1 => "ForceFreeStates/Solutions/GalerkinIntegration/psi",)),
     "ForceFreeStates/Solutions/GalerkinIntegration/xi_psi" =>
-        (; long_name="closed axis-to-edge ξ^ψ profiles (ideal or inner-layer closure; identity-at-edge basis)", dims=("mode", "solution", "psi"), attach=(3 => "ForceFreeStates/Solutions/GalerkinIntegration/psi",)),
+        (;
+            long_name="closed axis-to-edge ξ^ψ profiles (ideal or inner-layer closure; identity-at-edge basis)",
+            dims=("mode", "solution", "psi"),
+            attach=(3 => "ForceFreeStates/Solutions/GalerkinIntegration/psi",)
+        ),
     "ForceFreeStates/Solutions/GalerkinIntegration/dxi_psidpsi" =>
         (; long_name="analytic ψ_N derivative of the closed ξ^ψ profiles", dims=("mode", "solution", "psi"), attach=(3 => "ForceFreeStates/Solutions/GalerkinIntegration/psi",)),
     "ForceFreeStates/Solutions/GalerkinIntegration/xi_s" =>
-        (; long_name="surface displacement ξ_s from the outer ideal-MHD relation", dims=("mode", "solution", "psi"), attach=(3 => "ForceFreeStates/Solutions/GalerkinIntegration/psi",)),
+        (;
+            long_name="surface displacement ξ_s from the outer ideal-MHD relation",
+            dims=("mode", "solution", "psi"),
+            attach=(3 => "ForceFreeStates/Solutions/GalerkinIntegration/psi",)
+        ),
     "ForceFreeStates/Solutions/GalerkinIntegration/Basis/psi" => (; long_name="normalized poloidal flux ψ_N grid of the raw Galerkin basis", scale="psi_gal_basis"),
     "ForceFreeStates/Solutions/GalerkinIntegration/Basis/is_rational" =>
         (; long_name="flag: grid node lies on a rational surface", dims=("psi",), attach=(1 => "ForceFreeStates/Solutions/GalerkinIntegration/Basis/psi",)),

src/ForceFreeStates/Riccati/Crossings.jl

@@ -6,14 +6,14 @@
 Cross a singular surface for the Riccati formulation. Replaces `cross_ideal_singular_surf!`
 for the Riccati integration path with two key differences:
 
-1. **No Gaussian reduction**: `cross_ideal_singular_surf!` calls `compute_solution_norms!`
-   which applies Gaussian reduction to (S, I). This divides by pivot elements of S, which
-   can be near-zero (S = 0 at axis and grows slowly), producing NaN/Inf in U₂. For Riccati,
-   S is bounded so Gaussian reduction is unnecessary.
+ 1. **No Gaussian reduction**: `cross_ideal_singular_surf!` calls `compute_solution_norms!`
+    which applies Gaussian reduction to (S, I). This divides by pivot elements of S, which
+    can be near-zero (S = 0 at axis and grows slowly), producing NaN/Inf in U₂. For Riccati,
+    S is bounded so Gaussian reduction is unnecessary.
 
-2. **Direct column zeroing**: Instead of using the GR-sorted `odet.index` to identify the
-   column to zero, we use `ipert_res` directly (the resonant mode index). This is valid since
-   without GR there is no permutation applied to the columns of S.
+ 2. **Direct column zeroing**: Instead of using the GR-sorted `odet.index` to identify the
+    column to zero, we use `ipert_res` directly (the resonant mode index). This is valid since
+    without GR there is no permutation applied to the columns of S.
 
 **Δ' normalization**: This function expects `odet.u` in the bounded (U₁, U₂) form produced by
 `riccati_integrate_chunk!` with `needs_crossing=true` (final renorm skipped). ca_l is computed
@@ -54,16 +54,13 @@ end
 """
     _two_sided_singular_asymptotics(singp, ctrl, equil, mats, intr) -> (left, right)
 
-Compute left- (`sig=-1`) and right- (`sig=+1`) side singular asymptotics matching
-Fortran STRIDE's separate vmatl/vmatr (sing_vmat). Alpha is taken from the right
-side and shared with the left.
+Left (`sig=-1`) and right (`sig=+1`) singular asymptotics, each from its own resonant block.
 """
 function _two_sided_singular_asymptotics(singp::SingType, ctrl::ForceFreeStatesControl,
-                                         equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines,
-                                         intr::ForceFreeStatesInternal)
+    equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines,
+    intr::ForceFreeStatesInternal)
     sing_asymp_right = compute_sing_asymptotics(singp, ctrl, equil, mats, intr; sig=1.0)
-    sing_asymp_left  = compute_sing_asymptotics(singp, ctrl, equil, mats, intr; sig=-1.0,
-                                                alpha_override=sing_asymp_right.alpha)
+    sing_asymp_left = compute_sing_asymptotics(singp, ctrl, equil, mats, intr; sig=-1.0)
     return sing_asymp_left, sing_asymp_right
 end
 
@@ -83,7 +80,7 @@ end
 
 # Capture left-side asymptotic data into odet.ca_l and singp.ua_left/psi_ua_left.
 function _capture_left_crossing_data!(odet::OdeState, singp::SingType, sing_asymp_left,
-                                      dpsi::Float64, intr::ForceFreeStatesInternal, ising::Int)
+    dpsi::Float64, intr::ForceFreeStatesInternal, ising::Int)
     ua = sing_get_ua(sing_asymp_left, dpsi)
     singp.ua_left = copy(ua)
     singp.psi_ua_left = odet.psifac
@@ -95,9 +92,9 @@ end
 # odet.psifac to the right side. The zeroed columns stay zero through the predictor
 # since du[:, ipert_res, :] = 0 when u[:, ipert_res, :] = 0.
 function _predict_across_singular_surface!(odet::OdeState, ctrl::ForceFreeStatesControl,
-                                           equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines,
-                                           intr::ForceFreeStatesInternal, ising::Int,
-                                           ipert_res, dpsi::Float64, sing_asymp_right)
+    equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines,
+    intr::ForceFreeStatesInternal, ising::Int,
+    ipert_res, dpsi::Float64, sing_asymp_right)
     if ctrl.kinetic_factor == 0
         for i in eachindex(sing_asymp_right.r1)
             odet.u[:, ipert_res[i], :] .= 0
@@ -117,8 +114,8 @@ end
 # Column ipert_res of [U₁_new; U₂_new] = ua[:, ipert_res+N, :] (the introduced small asymptotic),
 # so ca_r[ipert_res, ipert_res, 2] = 1 regardless of other columns' normalization.
 function _capture_right_crossing_data!(odet::OdeState, singp::SingType, sing_asymp_right,
-                                       dpsi::Float64, intr::ForceFreeStatesInternal, ising::Int,
-                                       ipert_res, ctrl::ForceFreeStatesControl)
+    dpsi::Float64, intr::ForceFreeStatesInternal, ising::Int,
+    ipert_res, ctrl::ForceFreeStatesControl)
     ua = sing_get_ua(sing_asymp_right, dpsi)
     singp.ua_right = copy(ua)
     singp.psi_ua_right = odet.psifac
@@ -135,9 +132,9 @@ end
 # STUB: per-surface ca-based Δ' (not physically valid; see SingType.delta_prime docstring).
 # The canonical Δ' is intr.delta_prime_matrix from compute_delta_prime_matrix!.
 function _stash_per_surface_delta_prime_stub!(odet::OdeState, intr::ForceFreeStatesInternal,
-                                              ising::Int, ipert_res, sing_asymp_right,
-                                              equil::Equilibrium.PlasmaEquilibrium,
-                                              ctrl::ForceFreeStatesControl)
+    ising::Int, ipert_res, sing_asymp_right,
+    equil::Equilibrium.PlasmaEquilibrium,
+    ctrl::ForceFreeStatesControl)
     ctrl.kinetic_factor == 0 || return
     denom = (2π)^2 * equil.psio
     n_res = length(sing_asymp_right.r1)