ForceFreeStates - BUGFIX - Count every fixed-boundary zero crossing: inside Riccati chunks, and after a switch of the smallest eigenvalue

src/ForceFreeStates/FixedBoundaryStability.jl

@@ -1,18 +1,18 @@
 """
-    evaluate_stability_criterion!(odet, profiles) -> nzero
+    evaluate_stability_criterion!(odet, profiles; first_step=1) -> nzero
 
 Fill `odet.crit_store` over all saved steps and return the number of zero crossings of the
-critical eigenvalue (each signals a fixed-boundary instability). Warns with a count if W⁻¹ was
-non-Hermitian beyond tolerance at any step.
+critical eigenvalue (each signals a fixed-boundary instability) after step `first_step`. Warns with
+a count if W⁻¹ was non-Hermitian beyond tolerance at any step.
 """
-function evaluate_stability_criterion!(odet::OdeState, profiles::Equilibrium.ProfileSplines)
+function evaluate_stability_criterion!(odet::OdeState, profiles::Equilibrium.ProfileSplines; first_step::Int=1)
 
     resize!(odet.crit_store, odet.step)
     nzero = 0
     nonherm_count = 0
 
     for istep in 1:odet.step
-        zero_cross, nonherm = check_for_zero_crossings!(odet, profiles, istep)
+        zero_cross, nonherm = check_for_zero_crossings!(odet, profiles, istep; first_step)
         if zero_cross
             nzero += 1
         end
@@ -29,18 +29,20 @@ function evaluate_stability_criterion!(odet::OdeState, profiles::Equilibrium.Pro
 end
 
 """
-    check_for_zero_crossings!(odet, profiles, istep) -> zero_cross, nonherm
+    check_for_zero_crossings!(odet, profiles, istep; first_step=1) -> zero_cross, nonherm
 
 Compute `crit` (smallest-magnitude eigenvalue of W⁻¹, scaled by (dV/dψ)²) at step `istep` into
-`odet.crit_store`, and report whether it changed sign since the previous step. A sign change counts
-only if the interpolated midpoint confirms it is a genuine crossing rather than noise.
+`odet.crit_store`, and report whether it changed sign since the previous step (never for steps up to
+`first_step`). A sign change counts only if the number of negative eigenvalues of W⁻¹ also dropped: with
+F̄ > 0 a crossing raises an eigenvalue through 0, while a switch of the smallest eigenvalue or a passage
+through ∞ does not.
 
 ### Returns
 
   - `zero_cross::Bool`: True if a physical zero crossing was detected
   - `nonherm::Bool`: True if W⁻¹ was non-Hermitian beyond tolerance
 """
-@with_pool pool function check_for_zero_crossings!(odet::OdeState, profiles::Equilibrium.ProfileSplines, istep::Int)
+@with_pool pool function check_for_zero_crossings!(odet::OdeState, profiles::Equilibrium.ProfileSplines, istep::Int; first_step::Int=1)
 
     u = acquire!(pool, eltype(odet.u_store), size(odet.u_store)[1:3])
     psi = odet.psi_store[istep]
@@ -50,18 +52,11 @@ only if the interpolated midpoint confirms it is a genuine crossing rather than
     odet.crit_store[istep] = crit_val * dVdpsi^2
 
     zero_cross = false
-    if istep > 1 && odet.crit_store[istep] * odet.crit_store[istep-1] < 0
-        crit = odet.crit_store[istep]
-        crit_prev = odet.crit_store[istep-1]
-        # Reject noise: crit at the interpolated crossing must lie between and well below both neighbors
-        fac = crit / (crit - crit_prev)
-        psi_mid = psi - fac * (psi - odet.psi_store[istep-1])
-        u_mid = u .- fac .* (u .- @view(odet.u_store[:, :, :, istep-1]))
-        dVdpsi = profiles.dVdpsi_spline(psi_mid; hint=odet.spline_hint)
-        crit_mid_val, _ = compute_smallest_eigenvalue(u_mid)
-        crit_mid = crit_mid_val * dVdpsi^2
-        if (crit_mid - crit) * (crit_mid - crit_prev) < 0 && abs(crit_mid) < 0.5 * min(abs(crit), abs(crit_prev))
+    if istep > first_step && odet.crit_store[istep] * odet.crit_store[istep-1] < 0
+        if count_negative_wp_inverse(u) < count_negative_wp_inverse(@view odet.u_store[:, :, :, istep-1])
             zero_cross = true
+            crit, crit_prev = odet.crit_store[istep], odet.crit_store[istep-1]
+            psi_mid = psi - crit / (crit - crit_prev) * (psi - odet.psi_store[istep-1])
             @info "Zero crossing detected at ψ = $(@sprintf("%.3f", psi_mid))"
         end
     end
@@ -106,3 +101,73 @@ smallest magnitude. Also flags whether W⁻¹ was non-Hermitian beyond tolerance
     crit = evals[findmin(abs, evals)[2]]
     return crit, nonherm
 end
+
+"""
+    step_zero_crossings(u, u_prev, nneg_prev, direction) -> ncross, nneg
+
+Zero crossings of a solution over one ODE step from `u_prev` to `u`, taken in ψ direction `direction` (±1). With
+F̄ > 0 every crossing moves an eigenvalue of W⁻¹ upward in ψ through 0, so `nneg`, the number of negative eigenvalues
+of the congruent U₂ᴴU₁, drops by one and `crit` changes sign. A passage through ∞ (det U₂ = 0) also changes `nneg`,
+but not the sign of `crit`, which is evaluated only where `nneg` differs from the previous step's `nneg_prev`.
+"""
+function step_zero_crossings(u::AbstractArray{ComplexF64,3}, u_prev::AbstractArray{ComplexF64,3}, nneg_prev::Int, direction::Int)
+    nneg = count_negative_wp_inverse(u)
+    ncross = 0
+    if nneg != nneg_prev
+        crit, _ = compute_smallest_eigenvalue(u)
+        crit_prev, _ = compute_smallest_eigenvalue(u_prev)
+        if crit * crit_prev < 0
+            ncross = max(0, direction * (nneg_prev - nneg))
+        end
+    end
+    return ncross, nneg
+end
+
+"""
+    crossings_beyond_principal(P12, u_s, u_e, direction) -> Int
+
+Zero crossings of a solution in a chunk beyond those of the chunk's principal solution (U₁ = 0 at its start): the
+negative eigenvalues of `direction`·herm(X_sᴴ P₁₂⁻¹ X_e), a Morse index split exact for F̄ > 0. `u_s` and `u_e` hold
+the solution at the chunk's integration start and end in one column basis, X are their U₁ blocks, and `P12` is the
+U₁ block of the principal solution at the end. Each column is first divided by the norm of its `u_s` column, as their
+spread near the axis would otherwise turn the small eigenvalues into sign noise.
+"""
+function crossings_beyond_principal(P12::AbstractMatrix{ComplexF64}, u_s::AbstractArray{ComplexF64,3}, u_e::AbstractArray{ComplexF64,3}, direction::Int)
+    colnorms = [norm(view(u_s, :, j, :)) for j in axes(u_s, 2)]'
+    X_s = u_s[:, :, 1] ./ colnorms
+    X_e = u_e[:, :, 1] ./ colnorms
+    return count_negative_eigenvalues!(hermitianpart!(direction .* (X_s' * (P12 \ X_e))))
+end
+
+"""
+    count_negative_wp_inverse(u) -> Int
+
+Number of negative eigenvalues of W⁻¹ = U₁U₂⁻¹, from the congruent U₂ᴴU₁.
+"""
+@with_pool pool function count_negative_wp_inverse(u::AbstractArray{ComplexF64,3})
+    h = acquire!(pool, ComplexF64, size(u, 1), size(u, 2))
+    mul!(h, view(u, :, :, 2)', view(u, :, :, 1))
+    return count_negative_eigenvalues!(Hermitian(h))
+end
+
+"""
+    count_negative_eigenvalues!(H) -> Int
+
+Number of negative eigenvalues of the Hermitian `H`, from its Bunch–Kaufman LDLᴴ factorization (overwrites `H`): a
+1×1 pivot counts by its sign, a 2×2 pivot block holds exactly one. Signs are taken as computed, so unlike `inertia`'s
+relative zero tolerance this keeps the small eigenvalues of a graded `H`.
+"""
+function count_negative_eigenvalues!(H::Hermitian{ComplexF64})
+    F = bunchkaufman!(H; check=false)
+    nneg, k = 0, 1
+    while k <= size(H, 1)
+        if F.ipiv[k] < 0  # 2×2 block in rows k, k+1
+            nneg += 1
+            k += 2
+        else
+            nneg += real(F.LD[k, k]) < 0
+            k += 1
+        end
+    end
+    return nneg
+end

src/ForceFreeStates/Riccati/Driver.jl

@@ -37,6 +37,8 @@ Select with `Riccati()`, or `integrator = "Riccati"` in a deck's `[ForceFreeStat
 `singfac_min != 0`. Results do not depend on the thread count; `ctrl.algorithm.nchunks` is the only
 tunable. `odet.u_store` holds chunk-endpoint Riccati states, not the EL basis
 (`odet.u_store_el_basis == false`), so PerturbedEquilibrium and the HDF5 ξ datasets require `Forward()`.
+`odet.nzero` counts the zero crossings of each chunk before the outer plasma from its propagator (principal
+solution plus `crossings_beyond_principal`), and those of the re-integrated outer plasma from its stored steps.
 """
 function riccati_eulerlagrange_integration(
     ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium,
@@ -51,7 +53,7 @@ function riccati_eulerlagrange_integration(
     # Harvest thread-local solver-step counts; the outer re-integration counts via its own callback.
     odet.total_steps += sum(p.total_steps for p in odet_proxies)
 
-    S_at_surface_left, last_crossing_step =
+    S_at_surface_left, last_crossing_step, chunk_nzero =
         _assemble_propagators_serially!(odet, propagators, chunks, ctrl, equil, mats, intr)
 
     _reintegrate_outer_plasma!(odet, last_crossing_step, ctrl, equil, mats, intr)
@@ -63,7 +65,9 @@ function riccati_eulerlagrange_integration(
     if intr.verbose
         @info "Evaluating fixed-boundary stability criterion"
     end
-    odet.nzero = evaluate_stability_criterion!(odet, equil.profiles)
+    # Chunks before the outer re-integration (or a truncated edge) count their own crossings; crit checks the stored steps after.
+    ncore = min(last_crossing_step, odet.step)
+    odet.nzero = sum(@view chunk_nzero[1:ncore]) + evaluate_stability_criterion!(odet, equil.profiles; first_step=ncore)
     transform_u!(odet, intr)  # no-op when ifix=0 (no Gaussian reduction)
 
     return odet, propagators, chunks, S_at_surface_left
@@ -123,7 +127,9 @@ end
 # Apply chunk propagators serially, renormalizing to (S, I) after each (STRIDE's ode_fixup role:
 # keeps conditioning at O(cond_per_chunk) instead of its K-th power). Backward chunks are applied
 # via apply_propagator_inverse!. S_at_surface_left is the Riccati S at each surface's left boundary,
-# the Δ' BVP axis BC. Returns (S_at_surface_left, last_crossing_step).
+# the Δ' BVP axis BC. Each chunk stores one step, so last_crossing_step also indexes the last crossing
+# chunk. chunk_nzero[i] counts the zero crossings of the solution in chunk i.
+# Returns (S_at_surface_left, last_crossing_step, chunk_nzero).
 function _assemble_propagators_serially!(odet::OdeState, propagators::Vector{ChunkPropagator},
     chunks::Vector{IntegrationChunk},
     ctrl::ForceFreeStatesControl,
@@ -132,12 +138,18 @@ function _assemble_propagators_serially!(odet::OdeState, propagators::Vector{Chu
     N = intr.numpert_total
     S_at_surface_left = Matrix{ComplexF64}[]
     last_crossing_step = 1
+    chunk_nzero = zeros(Int, length(chunks))
     for (i, chunk) in enumerate(chunks)
+        # The solution at the chunk's integration start (u_s) and end (u_e), in one column basis until the renorm
+        u_left = copy(odet.u)
         if chunk.direction == -1
             apply_propagator_inverse!(odet, propagators[i])
+            u_s, u_e = odet.u, u_left
         else
             apply_propagator!(odet, propagators[i])
+            u_s, u_e = u_left, odet.u
         end
+        chunk_nzero[i] = propagators[i].principal_nzero + crossings_beyond_principal(view(propagators[i].block_lower_ic, :, :, 1), u_s, u_e, chunk.direction)
         renormalize_riccati_inplace!(odet.u, N)
         odet.psifac = chunk.psi_end
         odet.q = equil.profiles.q_spline(odet.psifac)
@@ -164,7 +176,7 @@ function _assemble_propagators_serially!(odet::OdeState, propagators::Vector{Chu
             odet.step += 1
         end
     end
-    return S_at_surface_left, last_crossing_step
+    return S_at_surface_left, last_crossing_step, chunk_nzero
 end
 
 # Re-integrate last rational surface → psilim with Riccati: FM propagation loses precision there at

src/ForceFreeStates/Riccati/Propagators.jl

@@ -269,6 +269,8 @@ renormalize_riccati!(odet::OdeState, intr::ForceFreeStatesInternal) = renormaliz
 
 Fill `prop` with the chunk's fundamental matrix by integrating the EL ODE from IC (I, 0) into
 `prop.block_upper_ic` and from (0, I) into `prop.block_lower_ic`, without renormalization or storage.
+The (0, I) solve is the chunk's principal solution; its callback also counts that solution's zero crossings
+into `prop.principal_nzero` (`step_zero_crossings`).
 `odet_proxy` holds the `sing_der!` side effects; concurrent calls need distinct proxies.
 """
 function integrate_propagator_chunk!(
@@ -300,15 +302,24 @@ function integrate_propagator_chunk!(
     prop.block_upper_ic .= sol.u[end]
     odet_proxy.total_steps += sol.stats.naccept  # thread-local; summed into odet after the BVP barrier
 
-    # Lower block IC: U₁ = 0, U₂ = I
+    # Lower block IC: U₁ = 0, U₂ = I (the principal solution), counting its zero crossings at every step
     u_lower = zeros(ComplexF64, N, N, 2)
     for i in 1:N
         u_lower[i, i, 2] = 1
     end
+    u_prev, nneg_prev = copy(u_lower), Ref(0)
+    prop.principal_nzero = 0
+    function count_principal_crossings!(integrator)
+        refresh_column_abstol!(integrator)
+        ncross, nneg_prev[] = step_zero_crossings(integrator.u, u_prev, nneg_prev[], chunk.direction)
+        prop.principal_nzero += ncross
+        u_prev .= integrator.u
+    end
+    cb = DiscreteCallback((u, t, integrator) -> true, count_principal_crossings!; save_positions=(false, false))
     odet_proxy.spline_hint[] = 1
     odet_proxy.mats_hint[] = 1
     prob = ODEProblem(sing_der!, u_lower, tspan, params)
-    sol = solve_el(prob, rtol)
+    sol = solve_el(prob, rtol; callback=cb)
     prop.block_lower_ic .= sol.u[end]
     odet_proxy.total_steps += sol.stats.naccept
 end

src/ForceFreeStates/Riccati/Types.jl

@@ -30,7 +30,8 @@ Represents the propagator Φ(ψ₂,ψ₁) computed by integrating the EL ODE fro
 identity-block initial conditions:
 
   - `block_upper_ic`: result of integrating with IC = (I_N, 0_N)  (U₁ = I, U₂ = 0)
-  - `block_lower_ic`: result of integrating with IC = (0_N, I_N)  (U₁ = 0, U₂ = I)
+  - `block_lower_ic`: result of integrating with IC = (0_N, I_N)  (U₁ = 0, U₂ = I), the chunk's principal solution
+  - `principal_nzero`: zero crossings (det U₁ = 0) of the principal solution inside the chunk
 
 Applying the propagator to the current state `u_prev`:
 
@@ -39,8 +40,9 @@ u₂_new = block_upper_ic[:,:,2] · u₁_prev + block_lower_ic[:,:,2] · u₂_pr
 
 Each chunk starts from an identity IC rather than the accumulated state, so chunks integrate in parallel.
 """
-struct ChunkPropagator
+mutable struct ChunkPropagator
     block_upper_ic::Array{ComplexF64,3}   # shape (N, N, 2) — result from IC = (I, 0)
     block_lower_ic::Array{ComplexF64,3}   # shape (N, N, 2) — result from IC = (0, I)
+    principal_nzero::Int
 end
-ChunkPropagator(N::Int) = ChunkPropagator(zeros(ComplexF64, N, N, 2), zeros(ComplexF64, N, N, 2))
+ChunkPropagator(N::Int) = ChunkPropagator(zeros(ComplexF64, N, N, 2), zeros(ComplexF64, N, N, 2), 0)

test/runtests_eulerlagrange.jl

@@ -423,6 +423,25 @@ end
         @test odet.u[3, 2, 1] != 0
     end
 
+    @testset "a zero crossing right after the smallest eigenvalue switches" begin
+        F = GeneralizedPerturbedEquilibrium.ForceFreeStates
+        ex = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example")
+        inputs = TOML.parsefile(joinpath(ex, "gpec.toml"))
+        equil = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(
+            GeneralizedPerturbedEquilibrium.equilibrium_config(inputs, ex),
+            GeneralizedPerturbedEquilibrium.analytic_config(inputs))
+        # W⁻¹ = diag(a, b) per step: the smallest switches from a to b (crit +0.3 → -0.25, no crossing),
+        # then b crosses 0 while crit becomes a (crit -0.25 → +0.6), which an interpolated midpoint rejects.
+        odet = F.OdeState(2, 4, 1, 1)
+        odet.step = 3
+        odet.psi_store[1:3] = [0.2, 0.3, 0.4]
+        for (i, w) in enumerate(([0.3, -0.5], [0.6, -0.25], [0.6, 5.0]))
+            odet.u_store[:, :, 1, i] = diagm(ComplexF64.(w))
+            odet.u_store[:, :, 2, i] = I(2)
+        end
+        @test F.evaluate_stability_criterion!(odet, equil.profiles) == 1
+    end
+
     @testset "OdeState construction" begin
         # Test basic OdeState initialization
         numpert_total = 5

test/runtests_riccati.jl

@@ -80,6 +80,41 @@ end
         @test odet.u[:, :, 1] ≈ S_expected rtol = 1e-12
     end
 
+    @testset "zero crossings of an oscillator: per step and Morse split" begin
+        # U₁' = ω U₂, U₂' = -ω U₁ per mode (B = ω > 0): U₁ = sin(ωψ + φ) vanishes as W⁻¹ = tan(ωψ + φ) rises through 0.
+        # The principal solution (ω = 1, φ = 0) on (0, 7] crosses at π and 2π; its poles at π/2 and 3π/2 do not count.
+        u, u_prev = zeros(ComplexF64, 1, 1, 2), zeros(ComplexF64, 1, 1, 2)
+        u_prev[1, 1, 2] = 1
+        ncross, nneg = 0, 0
+        for ψ in 0.1:0.1:7.0
+            u[1, 1, 1], u[1, 1, 2] = sin(ψ), cos(ψ)
+            k, nneg = FFS.step_zero_crossings(u, u_prev, nneg, 1)
+            ncross += k
+            u_prev .= u
+        end
+        @test ncross == 2
+
+        # Two modes on [0, L]: mode 1 crosses once more than its principal solution, mode 2 as often.
+        ω, φ, L = [1.0, 2.0], [1.5, 0.3], 2.0
+        function state(ψ)
+            y = zeros(ComplexF64, 2, 2, 2)
+            for j in 1:2
+                y[j, j, 1], y[j, j, 2] = sin(ω[j] * ψ + φ[j]), cos(ω[j] * ψ + φ[j])
+            end
+            return y
+        end
+        P12 = diagm(ComplexF64.(sin.(ω .* L)))
+        @test FFS.crossings_beyond_principal(P12, state(0.0), state(L), 1) == 1
+        # Independent of the column basis and its scales
+        T = ComplexF64[2 1; 1 3] * Diagonal([1e-9, 1e9])
+        rebase(y) = cat(y[:, :, 1] * T, y[:, :, 2] * T; dims=3)
+        @test FFS.crossings_beyond_principal(P12, rebase(state(0.0)), rebase(state(L)), 1) == 1
+        # Integrated backward from L, where the backward propagator's P₁₂ is -sin(ωL)
+        @test FFS.crossings_beyond_principal(-P12, state(L), state(0.0), -1) == 1
+        # A 2×2 Bunch–Kaufman pivot block holds one negative eigenvalue
+        @test FFS.count_negative_eigenvalues!(Hermitian(ComplexF64[0 1; 1 0])) == 1
+    end
+
     # ── Shared Solovev setup ──────────────────────────────────────────────────
     #
     # equil (Grad-Shafranov solve) and mats (metric matrices) are expensive and
@@ -122,6 +157,11 @@ end
         # Riccati stores chunk endpoints, the forward path stores every saved ODE step,
         # so the Riccati store is always the sparser of the two.
         @test odet_ric.step <= odet_fwd.step
+
+        # Both find this example fixed-boundary unstable, the Riccati core from its chunk propagators. Only the
+        # verdict is compared: the example is Mercier unstable at q = 2, where the crossings accumulate
+        # log-periodically into the surface, so their number depends on resolution.
+        @test odet_ric.nzero > 0 && odet_fwd.nzero > 0
     end
 
     # Note: a Solovev per-surface Δ' regression testset previously lived here,