ForceFreeStates - BUGFIX - Move a rational qlow start off its surface and grow the reduction store on demand

src/ForceFreeStates/CoreTypes.jl

@@ -54,7 +54,7 @@ Solve-time scratch state for a force-balance solve: the resolved mode space, sin
     kinsing_scan_cond::Vector{Float64} = Float64[]
     kinsing_scan_threshold::Float64 = 0.0
     psilim::Float64 = 0.0
-    psilow::Float64 = 0.0   # lower integration bound; the outermost crossing of q = qlow (sing_min!), else the first grid node
+    psilow::Float64 = 0.0   # lower integration bound; the outermost crossing of q = qlow (sing_min!), moved off a rational it lands on, else the first grid node
     qlim::Float64 = 0.0
     q1lim::Float64 = 0.0
     wall_settings::Union{Nothing,Vacuum.WallSettings} = nothing

src/ForceFreeStates/EulerLagrange.jl

@@ -42,7 +42,7 @@ records needed to undo it, and singular-surface asymptotic coefficients.
 ## Fields
 
   - `numpert_total::Int` - Total number of (m, n) modes.
-  - `numunorms_init::Int` - Allocated number of Gaussian-reduction records.
+  - `numunorms_init::Int` - Capacity of the Gaussian-reduction records; doubled on demand by `grow_reduction_store!`.
   - `msing::Int` - Number of singular surfaces (sizes `ca_l`/`ca_r`).
   - `numsteps_init::Int` - Initial allocation of stored steps.
   - `step::Int` - Index of the next stored step (1-based).
@@ -71,10 +71,10 @@ records needed to undo it, and singular-surface asymptotic coefficients.
   - `unorm::Vector{Float64}` - Current solution-column norms.
   - `unorm0::Vector{Float64}` - Reference norms since the last reduction.
   - `ifix::Int` - Number of Gaussian reductions performed.
-  - `index::Array{Int,2}` - Pivot column order for each reduction, shape `(numpert_total, numunorms_init)`.
+  - `index::Array{Int,2}` - Pivot column order for each reduction, shape `(numpert_total, capacity)`.
   - `sing_flag::Vector{Bool}` - Whether each reduction was at a singular surface.
   - `zeroed_idx::Vector{Vector{Int}}` - For each ideal crossing, the positions in `index` of the zeroed resonant solutions.
-  - `fixfac::Array{ComplexF64,3}` - Gaussian-reduction factors, shape `(numpert_total, numpert_total, numunorms_init)`.
+  - `fixfac::Array{ComplexF64,3}` - Gaussian-reduction factors, shape `(numpert_total, numpert_total, capacity)`.
   - `fixstep::Vector{Int64}` - Last stored step before each reduction.
   - `spline_hint::Base.RefValue{Int}` - Shared interval hint for profile splines evaluated at the same ψ.
   - `mats_hint::Base.RefValue{Int}` - Interval hint for the matrix splines; per `OdeState`, so per-thread copies do not race.
@@ -130,6 +130,31 @@ end
 OdeState(numpert_total::Int, numsteps_init::Int, numunorms_init::Int, msing::Int) =
     OdeState(; numpert_total, numsteps_init, numunorms_init, msing)
 
+"""
+    grow_reduction_store!(odet::OdeState) -> OdeState
+
+Double the capacity of the Gaussian-reduction records (`index`, `sing_flag`, `zeroed_idx`, `fixfac`,
+`fixstep`), keeping the records already written. A reduction fires at every rational surface and at every
+`renorm_threshold` trigger, so high-n runs outgrow `numunorms_init`; without the records `transform_u!`
+cannot rebuild the solution profiles.
+"""
+function grow_reduction_store!(odet::OdeState)
+    cap = length(odet.fixstep)
+    added = max(cap, 1)
+    N = odet.numpert_total
+    index = zeros(Int, N, cap + added)
+    index[:, 1:cap] .= odet.index
+    odet.index = index
+    fixfac = zeros(ComplexF64, N, N, cap + added)
+    fixfac[:, :, 1:cap] .= odet.fixfac
+    odet.fixfac = fixfac
+    append!(odet.sing_flag, falses(added))
+    append!(odet.zeroed_idx, [Int[] for _ in 1:added])
+    append!(odet.fixstep, zeros(Int64, added))
+    odet.numunorms_init = cap + added
+    return odet
+end
+
 """
     compute_delta_prime_from_ca!(odet, intr, equil)
 
@@ -542,7 +567,7 @@ function chunk_el_integration_bounds(odet::OdeState, ctrl::ForceFreeStatesContro
                 error(
                     "Singular surface q=$(sg.q) at ψ=$(sg.psifac) (q'=$(sg.q1)): its crossing standoff " *
                     "singfac_min/|n q'| = $(sg.psifac - psi_end) reaches back past the current chunk start ψ=$psi_current. " *
-                    "Move the domain start (psilow, or a qlow that lands on this surface) off it, or lower singfac_min if a neighboring crossing is too close."
+                    "Move the domain start (psilow) off it, or lower singfac_min if a neighboring crossing is too close."
                 )
             end
             @assert isempty(chunks) || psi_current >= chunks[end].psi_end "Overlapping chunks detected"
@@ -813,12 +838,8 @@ function compute_solution_norms!(u::Array{ComplexF64,3}, odet::OdeState, ctrl::F
     end
 
     if fix_required
-        if odet.ifix < ctrl.algorithm.numunorms_init
-            odet.ifix += 1
-        else
-            @warn "unorm storage reached, no longer saving fixfac data. Stability outputs and unorming will be correct, but cannot reconstruct `u`. \n
-            Increase `numunorms_init` if needed. Automatic resizing will be added in a future version."
-        end
+        odet.ifix == length(odet.fixstep) && grow_reduction_store!(odet)
+        odet.ifix += 1
         apply_gaussian_reduction!(u, odet, intr, sing_flag; resonant_rows)
         odet.new = true
     end

src/ForceFreeStates/Integrators.jl

@@ -16,7 +16,7 @@ Serial forward integration with Gaussian reduction and dense solution profiles.
   - `eulerlagrange_tolerance::Float64` - Relative tolerance for ODE integration.
   - `renorm_threshold::Float64` - Column-norm growth threshold for Gaussian reduction or Riccati renormalization.
   - `numsteps_init::Int` - Initial capacity for stored ODE solution steps.
-  - `numunorms_init::Int` - Initial capacity for solution-normalization records.
+  - `numunorms_init::Int` - Initial capacity for solution-normalization records; grown on demand.
   - `save_interval::Int` - Save every Nth ODE step; always retain steps near rational surfaces.
 """
 @kwdef struct Forward <: AbstractIntegrator
@@ -37,7 +37,7 @@ Chunked Riccati integration and the singular-surface boundary-value solve.
   - `eulerlagrange_tolerance::Float64` - Relative tolerance for ODE integration.
   - `renorm_threshold::Float64` - Column-norm growth threshold for Gaussian reduction or Riccati renormalization.
   - `numsteps_init::Int` - Initial capacity for stored ODE solution steps.
-  - `numunorms_init::Int` - Initial capacity for solution-normalization records.
+  - `numunorms_init::Int` - Initial capacity for solution-normalization records; grown on demand.
   - `save_interval::Int` - Save every Nth ODE step; always retain steps near rational surfaces.
   - `nchunks::Int` - Target chunk count; zero derives the count from the singular surfaces.
   - `extended_precision_bvp::Bool` - Use Complex{Double64} for the Delta-prime BVP and PEST3 combination to reduce cancellation error.

src/ForceFreeStates/Surfaces/Finding.jl

@@ -217,7 +217,8 @@ end
     sing_min!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium)
 
 Set the lower integration bound `intr.psilow` to `qlow_start_psi`, the one bound every
-integrator starts from. Excluding the q < qlow region matters most for the outer-region Galerkin
+integrator starts from, moved off any rational surface it lands on by
+`move_start_off_rational!`. Excluding the q < qlow region matters most for the outer-region Galerkin
 solve, whose Hermite FEM would otherwise integrate through unhandled q ≤ 1 ideal singularities and
 contaminate Δ′ at the innermost kept surface.
 """
@@ -227,6 +228,38 @@ function sing_min!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl,
     intr.psilow > profiles.xs[1] &&
         @info "sing_min: qlow=$(@sprintf("%.3f", ctrl.qlow)); raising psilow from $(@sprintf("%.5f", profiles.xs[1])) " *
               "to $(@sprintf("%.5f", intr.psilow)) (excludes the q < qlow region)"
+    move_start_off_rational!(intr, ctrl)
+    return intr.psilow
+end
+
+"""
+    move_start_off_rational!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl)
+
+Move `intr.psilow` just past any singular surface it lands on. A start on a rational, or inside its
+crossing standoff `singfac_min/|n q'|`, sits where the resonant Euler-Lagrange coefficients diverge
+and the first chunk would have to run backwards to the standoff; the start is placed one standoff
+beyond that surface instead, which then counts as inside the domain start and is not crossed. A
+rational `qlow` lands here for every n that makes `n·qlow` an integer (qlow = 1.5 at n = 2, qlow = 3
+at n = 3, 6, 9, ...).
+"""
+function move_start_off_rational!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl)
+    ctrl.singfac_min > 0 || return intr.psilow   # no standoff, no crossings: nothing to clear
+    for s in intr.sing
+        standoff = ctrl.singfac_min / (minimum(s.n) * abs(s.q1))
+        isfinite(standoff) || continue   # q' = 0 at the surface: not a crossing the chunker stands off from
+        if abs(intr.psilow - s.psifac) <= standoff
+            moved = s.psifac + standoff
+            moved < intr.psilim ||
+                error(
+                    "Integration start ψ_N=$(intr.psilow) lands on the q=$(s.q) rational surface at ψ_N=$(s.psifac), and the standoff " *
+                    "singfac_min/|n q'| = $standoff past it reaches psilim=$(intr.psilim); move psilow or qlow, or lower singfac_min."
+                )
+            @warn "Integration start ψ_N=$(@sprintf("%.7f", intr.psilow)) lands on the q=$(@sprintf("%.4f", s.q)) rational surface " *
+                  "(m/n = $(join(s.m, ","))/$(join(s.n, ","))) at ψ_N=$(@sprintf("%.7f", s.psifac)); " *
+                  "moving the start to ψ_N=$(@sprintf("%.7f", moved)), one crossing standoff singfac_min/|n q'| past it."
+            intr.psilow = moved
+        end
+    end
     return intr.psilow
 end
 

test/runtests_eulerlagrange.jl

@@ -326,6 +326,64 @@ end
         @test c.algorithm.save_interval == 1
     end
 
+    @testset "grow_reduction_store! keeps the written records and doubles the capacity" begin
+        FFS = GeneralizedPerturbedEquilibrium.ForceFreeStates
+        mpert = 2
+        odet = FFS.OdeState(mpert, 10, 2, 1)
+        odet.ifix = 2
+        odet.fixstep[1:2] = [3, 6]
+        odet.sing_flag[1:2] = [false, true]
+        odet.zeroed_idx[2] = [1]
+        for ifix in 1:2
+            odet.fixfac[:, :, ifix] = ComplexF64[1 0.5im*ifix; 0 1]
+            odet.index[:, ifix] = [2, 1]
+        end
+        kept = (copy(odet.fixfac), copy(odet.index), copy(odet.fixstep), copy(odet.sing_flag), deepcopy(odet.zeroed_idx))
+
+        FFS.grow_reduction_store!(odet)
+        @test length(odet.fixstep) == 4 && size(odet.fixfac) == (mpert, mpert, 4) && size(odet.index) == (mpert, 4)
+        @test length(odet.sing_flag) == 4 && length(odet.zeroed_idx) == 4 && odet.numunorms_init == 4
+        @test length(FFS.grow_reduction_store!(FFS.OdeState(mpert, 10, 0, 1)).fixstep) == 1   # a zero capacity still grows
+        @test odet.fixfac[:, :, 1:2] == kept[1] && odet.index[:, 1:2] == kept[2] && odet.fixstep[1:2] == kept[3]
+        @test odet.sing_flag[1:2] == kept[4] && odet.zeroed_idx[1:2] == kept[5] && isempty(odet.zeroed_idx[3])
+
+        # A reduction past the initial capacity grows the store instead of dropping the record
+        intr = FFS.ForceFreeStatesInternal(; mpert=mpert, numpert_total=mpert)
+        ctrl = FFS.ForceFreeStatesControl(; algorithm=Forward(; renorm_threshold=10.0, numunorms_init=1))
+        odet = FFS.OdeState(mpert, 10, 1, 1)
+        odet.u[:, 1, 1] .= [1, 0]
+        odet.u[:, 2, 1] .= [0, 1]
+        for _ in 1:3
+            FFS.compute_solution_norms!(odet.u, odet, ctrl, intr, true)
+        end
+        @test odet.ifix == 3 && length(odet.fixstep) == 4 && all(odet.sing_flag[1:3])
+    end
+
+    @testset "move_start_off_rational! steps the start one standoff past a rational it lands on" begin
+        FFS = GeneralizedPerturbedEquilibrium.ForceFreeStates
+        intr = FFS.ForceFreeStatesInternal(; mpert=2, numpert_total=2)
+        intr.sing = [FFS.SingType(; psifac=p, q=q, q1=2.0, m=[m], n=[2]) for (p, q, m) in ((0.25, 1.5, 3), (0.50, 2.0, 4))]
+        intr.msing = 2
+        intr.psilim = 1.0
+        ctrl = FFS.ForceFreeStatesControl(; algorithm=Forward(), singfac_min=1e-4)
+        standoff = 1e-4 / (2 * 2.0)
+
+        intr.psilow = 0.25
+        @test (@test_logs (:warn, r"lands on the q=1.5000 rational surface") FFS.move_start_off_rational!(intr, ctrl)) ≈ 0.25 + standoff
+        intr.psilow = 0.50 - 0.5 * standoff   # inside the standoff, just short of the surface
+        @test (@test_logs (:warn, r"m/n = 4/2") FFS.move_start_off_rational!(intr, ctrl)) ≈ 0.50 + standoff
+        intr.psilow = 0.30                     # off every rational: untouched, no warning
+        @test (@test_logs FFS.move_start_off_rational!(intr, ctrl)) == 0.30
+        intr.psilow = 0.25                     # no standoff when singfac_min = 0: the start is left where it is
+        @test (@test_logs FFS.move_start_off_rational!(intr, FFS.ForceFreeStatesControl(; algorithm=Forward(), singfac_min=0.0))) == 0.25
+        intr.sing[2].q1 = 0.0                  # q' = 0 at a surface: an infinite standoff is skipped, not applied
+        intr.psilow = 0.50
+        @test (@test_logs FFS.move_start_off_rational!(intr, ctrl)) == 0.50
+        intr.sing[2].q1 = 2.0
+        intr.psilim = 0.50 + 0.5 * standoff    # the moved start would pass psilim: an error, not a silent move
+        @test_throws ErrorException FFS.move_start_off_rational!(intr, ctrl)
+    end
+
     @testset "compute_solution_norms!" begin
         mpert = 2
         odet = GeneralizedPerturbedEquilibrium.ForceFreeStates.OdeState(mpert, 10, 10, 10)