ForceFreeStates - REFACTOR - ⏸️ Single-source the auto chunk-count target

docs/src/stability.md

@@ -108,9 +108,9 @@ nchunks    = 0         # 0 = auto: derived from the singular-surface count alone
 #### Chunking and thread independence
 
 `balance_integration_chunks` splits the base chunks until the count reaches a target derived
-from the number of singular surfaces, `max(2 m_s + 3, 8(m_s + 1) + m_s)`, where ``m_s`` is
-`msing`.  Setting `nchunks` overrides that target; a value below the ``2 m_s + 3`` floor is
-clamped up with a warning.  The target never consults `Threads.nthreads()`, and every chunk
+from the number of singular surfaces alone ([`auto_chunk_target`](@ref)).  Setting `nchunks`
+overrides that target; a value below the [`min_crossing_chunks`](@ref) floor is clamped up with a
+warning.  The target never consults `Threads.nthreads()`, and every chunk
 integrates independently from identity initial conditions, so the results do not depend on how
 many threads `julia -t` provides — threads change wall-clock only.
 

src/ForceFreeStates/EulerLagrange.jl

@@ -37,9 +37,9 @@ end
 # at the interval endpoints. Coefficients are ported from STRIDE's ode_itime cost model
 # (Fortran reference) and unchanged here. Tune only after re-fitting against a per-chunk
 # step-count sweep; touching these affects parallel-chunk load balancing.
-const ODE_COST_AXIS  = (a = 39695.0, b = 212830.0)
-const ODE_COST_RAT   = (a = 17147.0, b = 470710.0)
-const ODE_COST_EDGE  = (a =  1646.0, b =   4683.0)
+const ODE_COST_AXIS = (a=39695.0, b=212830.0)
+const ODE_COST_RAT = (a=17147.0, b=470710.0)
+const ODE_COST_EDGE = (a=1646.0, b=4683.0)
 
 """
     ode_itime_cost(psi1, psi2, intr) -> Float64
@@ -61,6 +61,24 @@ function ode_itime_cost(psi1::Float64, psi2::Float64, intr::ForceFreeStatesInter
     return cost
 end
 
+"""
+    min_crossing_chunks(msing) -> Int
+
+Floor on the chunk count set by the `msing` rational-surface crossings; explicit `nchunks`
+requests below it are clamped up (with a warning) by `balance_integration_chunks`.
+"""
+min_crossing_chunks(msing::Integer) = 2 * msing + 3
+
+"""
+    auto_chunk_target(msing) -> Int
+
+The chunk-count target `balance_integration_chunks` resolves `nchunks = 0` to: 8 non-crossing
+sub-chunks per segment (axis→surf₁, surfᵢ→surfᵢ₊₁, surfₙ→edge) plus the crossing chunks, so that
+`assemble_fm_matrix(condition=true)` can keep the accumulated propagator products well-conditioned
+(STRIDE uses 33 intervals for comparable problems). Always above [`min_crossing_chunks`](@ref).
+"""
+auto_chunk_target(msing::Integer) = 8 * (msing + 1) + msing
+
 """
     balance_integration_chunks(chunks, ctrl, intr) -> Vector{IntegrationChunk}
 
@@ -72,27 +90,21 @@ never from `Threads.nthreads()` — so the chunk list, and hence every Riccati o
 identical whatever thread count `julia -t` provides.
 
 Each split finds the equal-cost midpoint ψ_mid via bisection:
-  ode_itime_cost(psi_start, psi_mid) ≈ ode_itime_cost(psi_start, psi_end) / 2
+ode_itime_cost(psi_start, psi_mid) ≈ ode_itime_cost(psi_start, psi_end) / 2
 
 Sub-chunks inherit `needs_crossing=false` and `ising=0`. Only the LAST sub-chunk of
 each original chunk retains `needs_crossing=true` and the original `ising`, so the
 rational surface crossing still fires at the correct ψ in the serial assembly phase.
 """
 function balance_integration_chunks(chunks::Vector{IntegrationChunk}, ctrl::ForceFreeStatesControl, intr::ForceFreeStatesInternal)
-    min_chunks = 2 * intr.msing + 3
-    # Ensure enough sub-chunks for BVP propagator conditioning: at least 5 non-crossing
-    # sub-chunks per segment (axis→surf₁, surfᵢ→surfᵢ₊₁, surfₙ→edge), plus crossing
-    # chunks. STRIDE uses 33 intervals for comparable problems. Without enough sub-chunks,
-    # assemble_fm_matrix(condition=true) can't keep accumulated products well-conditioned
-    # because single long-span propagators may already have cond ~ 10²⁴.
-    min_bvp_intervals = 8 * (intr.msing + 1) + intr.msing
+    min_chunks = min_crossing_chunks(intr.msing)
     if ctrl.nchunks > 0
         if ctrl.nchunks < min_chunks
             @warn "nchunks = $(ctrl.nchunks) is below the $min_chunks chunks required by $(intr.msing) singular surfaces; clamping up."
         end
         target_n = max(ctrl.nchunks, min_chunks)
     else
-        target_n = max(min_chunks, min_bvp_intervals)
+        target_n = auto_chunk_target(intr.msing)
     end
 
     result = collect(chunks)
@@ -131,10 +143,10 @@ function balance_integration_chunks(chunks::Vector{IntegrationChunk}, ctrl::Forc
         psi_mid = (lo + hi) / 2.0
 
         left = IntegrationChunk(; psi_start=chunk.psi_start, psi_end=psi_mid,
-                                  needs_crossing=false, ising=0, direction=1)
+            needs_crossing=false, ising=0, direction=1)
         right = IntegrationChunk(; psi_start=psi_mid, psi_end=chunk.psi_end,
-                                   needs_crossing=chunk.needs_crossing, ising=chunk.ising,
-                                   direction=chunk.direction)
+            needs_crossing=chunk.needs_crossing, ising=chunk.ising,
+            direction=chunk.direction)
         splice!(result, best_idx, [left, right])
     end
 
@@ -284,20 +296,20 @@ is identified by dominant |U₁| component, giving the physically correct consta
 Frobenius solution and avoiding the spurious logarithmic irregularity.
 """
 function compute_axis_init(ffit::FourFitVars, profiles::Equilibrium.ProfileSplines,
-        intr::ForceFreeStatesInternal, psi_low::Float64)
-    N    = intr.numpert_total
+    intr::ForceFreeStatesInternal, psi_low::Float64)
+    N = intr.numpert_total
     hint = Ref(1)
 
     # Evaluate stability matrices at psi_low
     F_lower = zeros(ComplexF64, N, N)
-    kmat    = zeros(ComplexF64, N, N)
-    gmat    = zeros(ComplexF64, N, N)
+    kmat = zeros(ComplexF64, N, N)
+    gmat = zeros(ComplexF64, N, N)
     ffit.fmats_lower(vec(F_lower), psi_low; hint=hint)
-    ffit.kmats(vec(kmat),          psi_low; hint=hint)
-    ffit.gmats(vec(gmat),          psi_low; hint=hint)
+    ffit.kmats(vec(kmat), psi_low; hint=hint)
+    ffit.gmats(vec(gmat), psi_low; hint=hint)
 
     # singfac[j] = 1 / (m_j − n_j · q) for each mode j
-    q0      = profiles.q_spline(psi_low; hint=hint)
+    q0 = profiles.q_spline(psi_low; hint=hint)
     singfac = vec(1.0 ./ ((intr.mlow:intr.mhigh) .- q0 .* (intr.nlow:intr.nhigh)'))
 
     # F̄⁻¹ = (F_lower · F_lower')⁻¹ via the Cholesky factor
@@ -311,15 +323,15 @@ function compute_axis_init(ffit::FourFitVars, profiles::Equilibrium.ProfileSplin
     for j in 1:N
         sf = singfac[j]
         fi = Finv[j, j]
-        k  = kmat[j, j]
+        k = kmat[j, j]
         kd = conj(k)          # K̄†[j,j]
-        g  = gmat[j, j]
+        g = gmat[j, j]
 
         # 2×2 ODE matrix block for mode j [Glasser 2016 Eq. 22-24, diagonal approximation]
         m11 = -sf * fi * k
-        m12 =  sf^2 * fi
-        m21 =  g - kd * fi * k
-        m22 =  sf * kd * fi
+        m12 = sf^2 * fi
+        m21 = g - kd * fi * k
+        m22 = sf * kd * fi
 
         # Frobenius matrix A₀_j = ψ_low · M_j [Glasser 2016 Eq. 51]
         #! format: off
@@ -346,7 +358,7 @@ function compute_axis_init(ffit::FourFitVars, profiles::Equilibrium.ProfileSplin
         if abs(v2) > Base.sqrt(Base.eps(Float64)) * abs(v1)
             U1_init[j, j] = v1 / v2
         else
-            U1_init[j, j] =  one(ComplexF64)
+            U1_init[j, j] = one(ComplexF64)
             U2_init[j, j] = zero(ComplexF64)
         end
     end
@@ -365,7 +377,7 @@ Formerly `ode_axis_init!`. This now only initializes `psifac`, `ising_start`, an
 Move ising_start logic to chunk_el_integration_bounds?
 """
 function initialize_el_at_axis!(odet::OdeState, ctrl::ForceFreeStatesControl, ffit::FourFitVars,
-        profiles::Equilibrium.ProfileSplines, intr::ForceFreeStatesInternal)
+    profiles::Equilibrium.ProfileSplines, intr::ForceFreeStatesInternal)
 
     # Default psifac to minimum equilibrium psi value
     odet.psifac = profiles.xs[1]
@@ -542,7 +554,7 @@ function chunk_el_integration_bounds(odet::OdeState, ctrl::ForceFreeStatesContro
                 psi_end=psi_end,
                 needs_crossing=true,
                 ising=ising_current,
-                direction = bidirectional ? -1 : 1
+                direction=bidirectional ? -1 : 1
             ))
 
             # After crossing, we jump to the other side of the singular surface
@@ -975,7 +987,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`
@@ -992,7 +1004,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

src/ForceFreeStates/ForceFreeStatesStructs.jl

@@ -262,7 +262,7 @@ gpec.toml.
   - `save_interval::Int` - Save every Nth ODE step (1=all, 10=every 10th). Always saves near rational surfaces. (Same as `euler_step` in the Fortran)
   - `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.
+  - `nchunks::Int` - Target number of Riccati integration chunks. `0` (the default) derives the count from problem structure alone ([`auto_chunk_target`](@ref)), enough sub-chunks per segment to keep the accumulated propagator products well-conditioned. An explicit value below [`min_crossing_chunks`](@ref) 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
@@ -493,71 +493,46 @@ 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)`.
-
   - `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.
-
   - `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
-
   - `unorm::Vector{Float64}` - Current norms of the solution vectors (length `numpert_total`).
-
   - `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)
@@ -566,11 +541,8 @@ 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`).
-
   - `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)`.
-
   - `fixstep::Vector{Int64}` - Step indices (psi step positions) at which normalization/fixups were performed (length `numunorms_init`).
 """
 @kwdef mutable struct OdeState

test/runtests_decomposition_invariance.jl

@@ -49,11 +49,10 @@ end
     @test auto.delta_prime !== nothing
     msing = size(auto.delta_prime.matrix, 1)
 
-    # The nchunks=0 target, mirroring balance_integration_chunks' internal formula (as
-    # runtests_parallel_integration.jl does). Invariance is asserted ABOVE this floor only:
-    # fewer chunks than the msing-derived minimum is not a different decomposition but a
-    # structurally deficient one (the floor gives the rational-surface crossings room).
-    auto_target = max(2 * msing + 3, 8 * (msing + 1) + msing)
+    # The nchunks=0 target, from its single source. Invariance is asserted ABOVE this floor
+    # only: fewer chunks than the msing-derived minimum is not a different decomposition but
+    # a structurally deficient one (the floor gives the rational-surface crossings room).
+    auto_target = GP_TI.ForceFreeStates.auto_chunk_target(msing)
     finer = _solve_at_nchunks(dir, auto_target + 11)
     @test finer.delta_prime !== nothing