diff --git a/src/ForceFreeStates/CoreTypes.jl b/src/ForceFreeStates/CoreTypes.jl index 1cb794f65..062b3490c 100644 --- a/src/ForceFreeStates/CoreTypes.jl +++ b/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 diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl index d40700059..beed25f0c 100644 --- a/src/ForceFreeStates/EulerLagrange.jl +++ b/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 diff --git a/src/ForceFreeStates/Integrators.jl b/src/ForceFreeStates/Integrators.jl index 8210f4f8c..470c0030f 100644 --- a/src/ForceFreeStates/Integrators.jl +++ b/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. diff --git a/src/ForceFreeStates/Surfaces/Finding.jl b/src/ForceFreeStates/Surfaces/Finding.jl index 8bf9f6d3e..aac509cce 100644 --- a/src/ForceFreeStates/Surfaces/Finding.jl +++ b/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 diff --git a/test/runtests_eulerlagrange.jl b/test/runtests_eulerlagrange.jl index 2fe69f441..b2627b704 100644 --- a/test/runtests_eulerlagrange.jl +++ b/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)