diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl index 0c259db9eb..d407000592 100644 --- a/src/ForceFreeStates/EulerLagrange.jl +++ b/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 diff --git a/src/ForceFreeStates/Surfaces/Finding.jl b/src/ForceFreeStates/Surfaces/Finding.jl index 389e7c5431..0ea2b436cb 100644 --- a/src/ForceFreeStates/Surfaces/Finding.jl +++ b/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 diff --git a/test/runtests_eulerlagrange.jl b/test/runtests_eulerlagrange.jl index d818fc9b36..2fe69f4417 100644 --- a/test/runtests_eulerlagrange.jl +++ b/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 diff --git a/test/runtests_layer_overlap_fitzpatrick.jl b/test/runtests_layer_overlap_fitzpatrick.jl index bb1f86a49d..eb8822a42d 100644 --- a/test/runtests_layer_overlap_fitzpatrick.jl +++ b/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