ForceFreeStates - BUGFIX! - Force the reduction at kinetic singular crossings in the previous pivot order

src/ForceFreeStates/EulerLagrange.jl

@@ -779,15 +779,15 @@ end
     compute_solution_norms!(u::Array{ComplexF64,3}, odet::OdeState, ctrl::ForceFreeStatesControl, intr::ForceFreeStatesInternal, sing_flag::Bool; resonant_rows=Int[])
 
 Compute the column norms of u₁ and apply `apply_gaussian_reduction!` to `u` in place when their
-growth since the last reduction exceeds `ctrl.algorithm.renorm_threshold`, or when `sing_flag` is set.
-The first call after a reduction only records reference norms. Errors if any norm is zero.
-Port of Fortran `ode_unorm`.
+growth since the last reduction exceeds `ctrl.algorithm.renorm_threshold`, and always when `sing_flag` is set.
+The first call after a reduction records reference norms; if `sing_flag` forces a reduction there, it keeps
+the previous reduction's pivot order. Errors if any norm is zero. Port of Fortran `ode_unorm`.
 
 ### Arguments
 
   - `u` - Solution array, updated in place (the integrator state, to avoid copies)
   - `sing_flag` - Whether this call is at a singular surface
-  - `resonant_rows` - Resonant rows at an ideal crossing; when given, the reduction always runs and pivots on them
+  - `resonant_rows` - Resonant rows at an ideal crossing, pivoted on first
 """
 function compute_solution_norms!(u::Array{ComplexF64,3}, odet::OdeState, ctrl::ForceFreeStatesControl, intr::ForceFreeStatesInternal, sing_flag::Bool;
     resonant_rows=Int[])
@@ -798,21 +798,21 @@ function compute_solution_norms!(u::Array{ComplexF64,3}, odet::OdeState, ctrl::F
         error("One of the first solution vector norms unorm(1,$jmax) = 0")
     end
 
-    # The first call after a reduction records the reference norms; an ideal crossing reduces regardless
-    at_ideal_crossing = !isempty(resonant_rows)
-    if odet.new && !at_ideal_crossing
+    # A singular surface always reduces; otherwise reduce when the growth since the last reduction is too uneven
+    fix_required = sing_flag
+    if odet.new
         odet.new = false
         odet.unorm0 .= odet.unorm
-        return
-    end
-
-    # Growth since the last reduction (raw norms if a crossing comes right after one)
-    if !odet.new
+        # A forced reduction right after another keeps that reduction's pivot order (none to keep before the first)
+        if sing_flag && odet.ifix > 0
+            @views odet.unorm[odet.index[:, odet.ifix]] .= (intr.numpert_total-1):-1:0
+        end
+    else
         odet.unorm ./= odet.unorm0
+        fix_required |= maximum(odet.unorm) / minimum(odet.unorm) > ctrl.algorithm.renorm_threshold
     end
-    uratio = maximum(odet.unorm) / minimum(odet.unorm)
 
-    if uratio > ctrl.algorithm.renorm_threshold || sing_flag
+    if fix_required
         if odet.ifix < ctrl.algorithm.numunorms_init
             odet.ifix += 1
         else

src/ForceFreeStates/Surfaces/Finding.jl

@@ -152,6 +152,14 @@ function sing_lim!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl,
         end
     end
 
+    # A target below q_min has no psilim: say why instead of failing in the root find, as Fortran's "Can't find psilim"
+    intr.qlim >= equil.params.qmin ||
+        error(
+            "sing_lim!: the integration limit qlim = $(@sprintf("%.3f", intr.qlim)) lies below q_min = $(@sprintf("%.3f", equil.params.qmin)), " *
+            "so no flux surface reaches it. With truncate_at_rational_offset no rational surface of n = $(intr.nlow) is left to truncate at; " *
+            "widen the domain (qhigh, the layer-overlap cap) or set truncate_at_rational_offset = false."
+        )
+
     # qlim below qmax: Newton-solve q(psilim) = qlim
     if intr.qlim < equil.params.qmax
         # Seed from the edge inward at the outermost node with q ≤ qlim, so a reversed-shear profile, which crosses

test/runtests_eulerlagrange.jl

@@ -369,6 +369,21 @@ end
         odet.u[:, 2, 1] .= [1, 0]
         GeneralizedPerturbedEquilibrium.ForceFreeStates.compute_solution_norms!(odet.u, odet, ctrl, intr, true)
         @test odet.new == true  # fixup triggered
+
+        # Case 6: a singular call right after a reduction still reduces, keeping that reduction's pivot order
+        # rather than the raw-norm order [1, 2, 3]
+        FFS = GeneralizedPerturbedEquilibrium.ForceFreeStates
+        odet = FFS.OdeState(3, 10, 10, 1)
+        odet.u[:, :, 1] .= [6.0 0.4 0.2; 0.5 0.3 2.0; 0.1 5.0 1.5]
+        odet.u[:, :, 2] .= 1.0
+        odet.ifix = 1
+        odet.index[:, 1] .= [3, 1, 2]
+        odet.new = true
+        FFS.compute_solution_norms!(odet.u, odet, ctrl, FFS.ForceFreeStatesInternal(; numpert_total=3), true)
+        @test odet.ifix == 2
+        @test odet.new
+        @test odet.sing_flag[2]
+        @test odet.index[:, 2] == [3, 1, 2]
     end
 
     @testset "ideal-crossing reduction isolates the resonant harmonic" begin

test/runtests_layer_overlap_fitzpatrick.jl

@@ -196,6 +196,15 @@ end
         @test i.qlim <= qh
     end
 
+    # No rational surface far enough below the bound leaves nothing to truncate at: a clear error, not a failed root find.
+    # The offset puts the last rational's truncation point above qhigh and the one before it below q_min.
+    qmin = equil.params.qmin
+    c_none = ForceFreeStatesControl(; truncate_at_rational_offset=true, rational_offset_fraction=mod(qmin, 1) + 0.04,
+        qhigh=qmin + 0.02)
+    @test_throws r"lies below q_min" sing_lim!(_fresh(), c_none, equil)
+    # Without truncation, a qhigh below q_min fails the same way.
+    @test_throws r"lies below q_min" sing_lim!(_fresh(), ForceFreeStatesControl(; truncate_at_rational_offset=false, qhigh=qmin - 0.1), equil)
+
     # A cap placed exactly on a rational must exclude it on either side of the float knife-edge.
     q_of = equil.profiles.q_spline
     lo, hi = equil.profiles.xs[1], equil.params.psihigh_resolved