From bc9a4085158befd2c5bf5076d55d7b1d6678fc16 Mon Sep 17 00:00:00 2001 From: StuartBenjamin Date: Sat, 31 Oct 2026 17:39:10 -0400 Subject: [PATCH 1/2] ForceFreeStates - BUGFIX - Riccati: count the zero crossings inside each chunk, not only at chunk endpoints MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The Riccati integrator evaluated the fixed-boundary stability criterion only on the stored chunk-endpoint states in the core, so crossings inside a chunk were missed or cancelled (DIII-D 174956 at 1715 ms: 0 at n=1 and n=2 where RDCON finds 2 and 1). Each chunk now counts the crossings of its principal solution at every step of its propagator solve (step_zero_crossings: inertia of U₂ᴴU₁ by Bunch–Kaufman, crit checked only where it changes), and the assembly adds the crossings of the actual solution beyond those from the chunk's end data (crossings_beyond_principal, a Morse index split exact for F̄ > 0). The re-integrated outer plasma keeps its dense check. Same method as the Fortran STRIDE fix. Co-Authored-By: Claude Opus 5.5 --- src/ForceFreeStates/FixedBoundaryStability.jl | 82 ++++++++++++++++--- src/ForceFreeStates/Riccati/Driver.jl | 20 ++++- src/ForceFreeStates/Riccati/Propagators.jl | 15 +++- src/ForceFreeStates/Riccati/Types.jl | 8 +- test/runtests_riccati.jl | 40 +++++++++ 5 files changed, 146 insertions(+), 19 deletions(-) diff --git a/src/ForceFreeStates/FixedBoundaryStability.jl b/src/ForceFreeStates/FixedBoundaryStability.jl index 29896ac3db..7f0d9074bb 100644 --- a/src/ForceFreeStates/FixedBoundaryStability.jl +++ b/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,19 @@ 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 interpolated midpoint confirms it is a genuine crossing +rather than noise. ### 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,7 +51,7 @@ 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 + if istep > first_step && 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 @@ -106,3 +107,64 @@ 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`. +""" +@with_pool pool function step_zero_crossings(u::AbstractArray{ComplexF64,3}, u_prev::AbstractArray{ComplexF64,3}, nneg_prev::Int, direction::Int) + h = acquire!(pool, ComplexF64, size(u, 1), size(u, 2)) + mul!(h, view(u, :, :, 2)', view(u, :, :, 1)) + nneg = count_negative_eigenvalues!(Hermitian(h)) + 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_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 diff --git a/src/ForceFreeStates/Riccati/Driver.jl b/src/ForceFreeStates/Riccati/Driver.jl index ef85b7bfa4..526cd3c617 100644 --- a/src/ForceFreeStates/Riccati/Driver.jl +++ b/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 diff --git a/src/ForceFreeStates/Riccati/Propagators.jl b/src/ForceFreeStates/Riccati/Propagators.jl index 81bd3d831c..2991ec4c95 100644 --- a/src/ForceFreeStates/Riccati/Propagators.jl +++ b/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 diff --git a/src/ForceFreeStates/Riccati/Types.jl b/src/ForceFreeStates/Riccati/Types.jl index 04fd25355a..46d1cddcfd 100644 --- a/src/ForceFreeStates/Riccati/Types.jl +++ b/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) diff --git a/test/runtests_riccati.jl b/test/runtests_riccati.jl index b3904b9702..bf16b6e944 100644 --- a/test/runtests_riccati.jl +++ b/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, From ef25a06b7db01730128ad8c4399bf7a135261fe1 Mon Sep 17 00:00:00 2001 From: StuartBenjamin Date: Sun, 1 Nov 2026 04:35:05 -0500 Subject: [PATCH 2/2] =?UTF-8?q?ForceFreeStates=20-=20BUGFIX=20-=20zero=20c?= =?UTF-8?q?rossings:=20confirm=20a=20crit=20sign=20change=20by=20the=20ine?= =?UTF-8?q?rtia=20of=20W=E2=81=BB=C2=B9,=20not=20an=20interpolated=20midpo?= =?UTF-8?q?int?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The midpoint test rejects a real crossing when crit changed eigenvalue at the step before. A crossing lowers the number of negative eigenvalues of W⁻¹ (F̄ > 0); a switch of the smallest eigenvalue or a passage through ∞ does not. Co-Authored-By: Claude Opus 5.5 --- src/ForceFreeStates/FixedBoundaryStability.jl | 35 ++++++++++--------- test/runtests_eulerlagrange.jl | 19 ++++++++++ 2 files changed, 38 insertions(+), 16 deletions(-) diff --git a/src/ForceFreeStates/FixedBoundaryStability.jl b/src/ForceFreeStates/FixedBoundaryStability.jl index 7f0d9074bb..aacbaa4924 100644 --- a/src/ForceFreeStates/FixedBoundaryStability.jl +++ b/src/ForceFreeStates/FixedBoundaryStability.jl @@ -33,8 +33,9 @@ end 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 (never for steps up to -`first_step`). A sign change counts only if the interpolated midpoint confirms it is a genuine crossing -rather than noise. +`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 @@ -52,17 +53,10 @@ rather than noise. zero_cross = false if istep > first_step && 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 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 @@ -116,10 +110,8 @@ F̄ > 0 every crossing moves an eigenvalue of W⁻¹ upward in ψ through 0, so 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`. """ -@with_pool pool function step_zero_crossings(u::AbstractArray{ComplexF64,3}, u_prev::AbstractArray{ComplexF64,3}, nneg_prev::Int, direction::Int) - h = acquire!(pool, ComplexF64, size(u, 1), size(u, 2)) - mul!(h, view(u, :, :, 2)', view(u, :, :, 1)) - nneg = count_negative_eigenvalues!(Hermitian(h)) +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) @@ -147,6 +139,17 @@ function crossings_beyond_principal(P12::AbstractMatrix{ComplexF64}, u_s::Abstra 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 diff --git a/test/runtests_eulerlagrange.jl b/test/runtests_eulerlagrange.jl index 2fe69f4417..53c71ddd2c 100644 --- a/test/runtests_eulerlagrange.jl +++ b/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