From 9671c3f56ea2c214731aa466e6a342726e297a2f Mon Sep 17 00:00:00 2001 From: d-burg Date: Mon, 28 Sep 2026 11:21:55 -0400 Subject: [PATCH 1/6] ForceFreeStates - BUGFIX - Remove the resonant solution at ideal crossings regardless of reduction timing At an ideal rational-surface crossing the Gaussian reduction ordered columns by growth since the previous reduction and the crossing zeroed the first one. When a ucrit reduction had fired a few steps earlier every column had grown alike, the order was arbitrary, the lead column's pivot fell on a non-resonant harmonic, and the crossing removed the wrong solution. A crossing that landed on the step right after a reduction skipped its own reduction and zeroed a column chosen from the previous one. The ideal crossing now always reduces, leads with the column carrying the most of each resonant harmonic, and pivots it on the resonant row, so that column alone holds the resonant component and is the one zeroed. Kinetic crossings and ucrit reductions are unchanged. Co-Authored-By: Claude Opus 5.5 --- src/ForceFreeStates/EulerLagrange.jl | 63 +++++++++++++++++----------- test/runtests_eulerlagrange.jl | 37 ++++++++++++++++ 2 files changed, 75 insertions(+), 25 deletions(-) diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl index 6245d1565..d8fab3b39 100644 --- a/src/ForceFreeStates/EulerLagrange.jl +++ b/src/ForceFreeStates/EulerLagrange.jl @@ -766,11 +766,14 @@ function cross_ideal_singular_surf!( ising::Int ) - # Fixup solution at singular surface - compute_solution_norms!(odet.u, odet, ctrl, intr, true) + singp = intr.sing[ising] + ipert_res = 1 .+ singp.m .- intr.mlow .+ (singp.n .- intr.nlow) .* intr.mpert + + # Fixup solution at singular surface; the ideal jump reduces on the resonant rows so the + # solution it removes below is the resonant one + compute_solution_norms!(odet.u, odet, ctrl, intr, true; resonant_rows=(ctrl.kinetic_factor == 0 ? ipert_res : Int[])) # Compute direction-specific asymptotic power series for this singular surface - singp = intr.sing[ising] sing_asymp_right = compute_sing_asymptotics(singp, ctrl, equil, mats, intr; sig=1.0) sing_asymp_left = compute_sing_asymptotics(singp, ctrl, equil, mats, intr; sig=-1.0, alpha_override=sing_asymp_right.alpha) dpsi = singp.psifac - odet.psifac # ψ_res - ψ (positive) @@ -779,18 +782,14 @@ function cross_ideal_singular_surf!( ua = sing_get_ua(sing_asymp_left, dpsi) odet.ca_l[:, :, :, ising] .= sing_get_ca(odet.u, ua, intr) - # Single n: remove largest solution and sub in asymptotics on the other side - # Multi-n: if we remove the N largest modes in arbitrary order, we can mess up the - # diagonal structure of the matrix and later calculations. zeroed_idx let's us make sure - # the solution vector we're zeroing corresponds to the same block as the resonant mode we - # introduce. It is also needed when transforming u back to the full solution after integration. - ipert_res = 1 .+ singp.m .- intr.mlow .+ (singp.n .- intr.nlow) .* intr.mpert + # Remove the resonant solution for each resonance and sub in asymptotics on the other side. + # The reduction above placed the column pivoted on ipert_res[i] at index position i, which also + # keeps each zeroed column in the same n-block as the resonant mode introduced below; zeroed_idx + # records the positions for transforming u back to the full solution after integration. if ctrl.kinetic_factor == 0 - # Eliminate the solution with the largest norm (in the same block) for each resonance - odet.zeroed_idx[odet.ifix] = Int[] - for i in eachindex(sing_asymp_right.r1) - push!(odet.zeroed_idx[odet.ifix], findfirst(j -> (ipert_res[i] - 1) ÷ intr.mpert == (odet.index[j, odet.ifix] - 1) ÷ intr.mpert, 1:intr.numpert_total)) - odet.u[:, odet.index[odet.zeroed_idx[odet.ifix][i], odet.ifix], :] .= 0 + odet.zeroed_idx[odet.ifix] = collect(eachindex(ipert_res)) + for i in eachindex(ipert_res) + odet.u[:, odet.index[i, odet.ifix], :] .= 0 end end @@ -960,7 +959,7 @@ function integrate_el_region!( end """ - compute_solution_norms!(u::Array{ComplexF64,3}, odet::OdeState, ctrl::ForceFreeStatesControl, intr::ForceFreeStatesInternal, sing_flag::Bool) + compute_solution_norms!(u::Array{ComplexF64,3}, odet::OdeState, ctrl::ForceFreeStatesControl, intr::ForceFreeStatesInternal, sing_flag::Bool; resonant_rows=Int[]) Computes norms of the solution vectors of the array `u` and normalizes them if this is not the first call after a fixup. Formerly `ode_unorm!`. @@ -979,12 +978,16 @@ operate on `u` directly without extra copies. - u: Current solution vector array, updated in-place if fixfac is called - sing_flag: Indicates if normalization is occuring at a singular surface or not + - resonant_rows: Rows of the resonant harmonics at an ideal crossing. When given, the reduction + always runs, even right after a `ucrit` reduction, and pivots on these rows (see + `apply_gaussian_reduction!`) ### TODOs Add resizing logic for unorm arrays when ifix exceeds allocated size """ -function compute_solution_norms!(u::Array{ComplexF64,3}, odet::OdeState, ctrl::ForceFreeStatesControl, intr::ForceFreeStatesInternal, sing_flag::Bool) +function compute_solution_norms!(u::Array{ComplexF64,3}, odet::OdeState, ctrl::ForceFreeStatesControl, intr::ForceFreeStatesInternal, sing_flag::Bool; + resonant_rows=Int[]) # Compute norms of first solution vectors, abort if any are zero odet.unorm .= norm.(eachcol(u[:, :, 1])) @@ -994,11 +997,11 @@ function compute_solution_norms!(u::Array{ComplexF64,3}, odet::OdeState, ctrl::F end # Normalize unorm and perform Gaussian reduction if required - if odet.new + if odet.new && isempty(resonant_rows) odet.new = false odet.unorm0 .= odet.unorm else - odet.unorm ./= odet.unorm0 + odet.new || (odet.unorm ./= odet.unorm0) uratio = maximum(odet.unorm) / minimum(odet.unorm) if uratio > ctrl.ucrit || sing_flag # TODO: add resizing logic here as well @@ -1008,14 +1011,14 @@ function compute_solution_norms!(u::Array{ComplexF64,3}, odet::OdeState, ctrl::F @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 - apply_gaussian_reduction!(u, odet, intr, sing_flag) + apply_gaussian_reduction!(u, odet, intr, sing_flag; resonant_rows) odet.new = true end end end """ - apply_gaussian_reduction!(u::Array{ComplexF64,3}, odet::OdeState, intr::ForceFreeStatesInternal, sing_flag::Bool) + apply_gaussian_reduction!(u::Array{ComplexF64,3}, odet::OdeState, intr::ForceFreeStatesInternal, sing_flag::Bool; resonant_rows=Int[]) Applies Gaussian reduction to orthogonalize solution vectors in `u`. Formerly `ode_fixup!`. Performs the same function as `ode_fixup` in the Fortran code, @@ -1024,8 +1027,14 @@ Used when the spread in norms exceeds a threshold or when a rational surface is This will update both `u` and relevant fields in `odet` in-place. See the description of `compute_solution_norms!` for more details on the benefits of in-place `u` updates. + +Columns are triangularized in order of their growth since the last reduction. At an ideal +crossing, `resonant_rows` puts first, for each resonant row, the column with the largest +component there and pivots it on that row, so that column alone carries the resonant harmonic. +Growth order cannot identify the resonant solution when a `ucrit` reduction fired only a few +steps before the crossing, because every column has then grown alike. """ -function apply_gaussian_reduction!(u::Array{ComplexF64,3}, odet::OdeState, intr::ForceFreeStatesInternal, sing_flag::Bool) +function apply_gaussian_reduction!(u::Array{ComplexF64,3}, odet::OdeState, intr::ForceFreeStatesInternal, sing_flag::Bool; resonant_rows=Int[]) # Store data for the current fixup ifix = odet.ifix @@ -1038,16 +1047,20 @@ function apply_gaussian_reduction!(u::Array{ComplexF64,3}, odet::OdeState, intr: for isol in 1:intr.numpert_total odet.fixfac[isol, isol, ifix] = 1 end - # Sort unorm in descending order (since we triangularize from largest to smallest) - odet.index[:, ifix] = sortperm(odet.unorm; rev=true) + # Resonant columns first, then the rest in descending unorm (we triangularize from largest to smallest) + lead = Int[] + for r in resonant_rows + push!(lead, argmax(j -> j in lead ? -Inf : abs(u[r, j, 1]), 1:intr.numpert_total)) + end + odet.index[:, ifix] = vcat(lead, filter(j -> !(j in lead), sortperm(odet.unorm; rev=true))) # Triangularize primary solutions mask = trues(2, intr.numpert_total) for isol in 1:intr.numpert_total ksol = odet.index[isol, ifix] mask[2, ksol] = false - # Set pivot row based on max location - @views kpert = argmax(abs.(u[:, ksol, 1]) .* mask[1, :]) + # Set pivot row: the resonant row for a leading resonant column, else the max location + @views kpert = isol <= length(lead) ? resonant_rows[isol] : argmax(abs.(u[:, ksol, 1]) .* mask[1, :]) mask[1, kpert] = false # Eliminate other solution vectors below the pivot for jsol in 1:intr.numpert_total diff --git a/test/runtests_eulerlagrange.jl b/test/runtests_eulerlagrange.jl index 497e1cff0..64ab74d5f 100644 --- a/test/runtests_eulerlagrange.jl +++ b/test/runtests_eulerlagrange.jl @@ -294,6 +294,43 @@ end @test odet.new == true # fixup triggered end + @testset "ideal-crossing reduction isolates the resonant harmonic" begin + FFS = GeneralizedPerturbedEquilibrium.ForceFreeStates + # Column 1 has the largest norm but no resonant (row 3) component; column 2 carries the resonance. + # Growth order alone would lead with column 1 and pivot it on row 1, leaving row 3 in columns 2 and 3. + u0 = zeros(ComplexF64, 3, 3, 2) + u0[:, 1, 1] .= [6.0, 0.5, 0.0] + u0[:, 2, 1] .= [0.4, 0.3, 5.0] + u0[:, 3, 1] .= [0.2, 2.0, 1.5] + u0[:, :, 2] .= 1.0 + intr = FFS.ForceFreeStatesInternal(; numpert_total=3) + ctrl = FFS.ForceFreeStatesControl(; ucrit=1e3) + fresh(u) = (odet = FFS.OdeState(3, 10, 10, 1); odet.u .= u; odet) + + odet = fresh(u0) + odet.new = false + odet.unorm0 .= norm.(eachcol(u0[:, :, 1])) # every column grown alike since the last reduction + FFS.compute_solution_norms!(odet.u, odet, ctrl, intr, true; resonant_rows=[3]) + @test odet.index[1, 1] == 2 + @test odet.u[3, 2, 1] == u0[3, 2, 1] + @test odet.u[3, 1, 1] == 0 && odet.u[3, 3, 1] == 0 # only the lead column keeps the resonant harmonic + @test odet.sing_flag[1] + + # A crossing right after a ucrit reduction (new = true) still reduces on the resonant row. + odet = fresh(u0) + FFS.compute_solution_norms!(odet.u, odet, ctrl, intr, true; resonant_rows=[3]) + @test odet.ifix == 1 && odet.index[1, 1] == 2 && odet.new + @test odet.u[3, 1, 1] == 0 && odet.u[3, 3, 1] == 0 + + # Without resonant rows the growth-ordered reduction is unchanged. + odet = fresh(u0) + odet.new = false + odet.unorm0 .= 1.0 + FFS.compute_solution_norms!(odet.u, odet, ctrl, intr, true) + @test odet.index[1, 1] == 1 + @test odet.u[3, 2, 1] != 0 + end + @testset "OdeState construction" begin # Test basic OdeState initialization numpert_total = 5 From 7de3a61dfebe930d730d43c39eb064a7ae9d7207 Mon Sep 17 00:00:00 2001 From: d-burg Date: Mon, 28 Sep 2026 11:30:13 -0400 Subject: [PATCH 2/6] Regression - TEST - Track Solovev n=1 with reductions just before the rational-surface crossings ucrit = 2 at eulerlagrange_tolerance 1e-10 makes a ucrit reduction fire two steps before the q = 2 crossing, which on develop leaves the crossing removing a non-resonant solution. Co-Authored-By: Claude Opus 5.5 --- .../cases/solovev_n1_frequent_reduction.toml | 48 +++++++++++++++++++ 1 file changed, 48 insertions(+) create mode 100644 regression-harness/cases/solovev_n1_frequent_reduction.toml diff --git a/regression-harness/cases/solovev_n1_frequent_reduction.toml b/regression-harness/cases/solovev_n1_frequent_reduction.toml new file mode 100644 index 000000000..212190f5b --- /dev/null +++ b/regression-harness/cases/solovev_n1_frequent_reduction.toml @@ -0,0 +1,48 @@ +# Regression case: Solovev n=1 ideal stability with Gaussian reductions every few steps. +# Reuses the Solovev_ideal_example deck via [overrides] with ucrit = 2, so a ucrit reduction +# regularly fires a step or two before an ideal rational-surface crossing. The crossing must still +# remove the resonant solution; if it removes another one, et[1] lands far from solovev_n1's value. +# Each [quantities.*] block names an HDF5 path in the run output, how to extract it, +# and the noise floor below which a difference is treated as zero. +[case] +name = "solovev_n1_frequent_reduction" +description = "Solovev analytical equilibrium, n=1, ideal stability, ucrit = 2 so reductions land just before the rational-surface crossings" +example_dir = "examples/Solovev_ideal_example" + +# Reduce whenever the solution-norm spread exceeds 2, at a tight Euler-Lagrange tolerance. +[overrides] +"ForceFreeStates.ucrit" = 2.0 +"ForceFreeStates.eulerlagrange_tolerance" = 1e-10 + +# Energies — leading generalized (W,N) pencil eigenvalues at the final truncation (psilim). +[quantities.et_real] +h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_energies" +type = "complex_vector" +extract = "real_first" +label = "total energy Re(et[1])" +noise_threshold = 1e-10 +order = 10 + +[quantities.et_imag] +h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_energies" +type = "complex_vector" +extract = "imag_first" +label = "total energy Im(et[1])" +noise_threshold = 1e-10 +order = 11 + +[quantities.ep_real] +h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_plasma_energies" +type = "complex_vector" +extract = "real_first" +label = "plasma energy Re(ep[1])" +noise_threshold = 1e-10 +order = 12 + +[quantities.ev_real] +h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_vacuum_energies" +type = "complex_vector" +extract = "real_first" +label = "vacuum energy Re(ev[1])" +noise_threshold = 1e-10 +order = 13 From d3de61ca5f17e4f3faf63f6b8640fbb9a7249cae Mon Sep 17 00:00:00 2001 From: d-burg Date: Sun, 4 Oct 2026 02:21:27 -0400 Subject: [PATCH 3/6] Regression - MINOR - Declare the tolerance class of each tracked quantity Co-Authored-By: Claude Opus 5.5 --- regression-harness/cases/solovev_n1_frequent_reduction.toml | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/regression-harness/cases/solovev_n1_frequent_reduction.toml b/regression-harness/cases/solovev_n1_frequent_reduction.toml index 212190f5b..5daecc6d2 100644 --- a/regression-harness/cases/solovev_n1_frequent_reduction.toml +++ b/regression-harness/cases/solovev_n1_frequent_reduction.toml @@ -22,6 +22,7 @@ extract = "real_first" label = "total energy Re(et[1])" noise_threshold = 1e-10 order = 10 +class = "physics_converged" [quantities.et_imag] h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_energies" @@ -30,6 +31,7 @@ extract = "imag_first" label = "total energy Im(et[1])" noise_threshold = 1e-10 order = 11 +class = "physics_converged" [quantities.ep_real] h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_plasma_energies" @@ -38,6 +40,7 @@ extract = "real_first" label = "plasma energy Re(ep[1])" noise_threshold = 1e-10 order = 12 +class = "physics_converged" [quantities.ev_real] h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_vacuum_energies" @@ -46,3 +49,4 @@ extract = "real_first" label = "vacuum energy Re(ev[1])" noise_threshold = 1e-10 order = 13 +class = "physics_converged" From d309eebf7dad81e84154d9a97103efb7f605862e Mon Sep 17 00:00:00 2001 From: d-burg Date: Fri, 9 Oct 2026 13:01:25 -0400 Subject: [PATCH 4/6] ForceFreeStates - REFACTOR - Drop the unreachable kinetic guards at ideal crossings and flatten the reduction control flow The ideal crossing is only reached with kinetic_factor = 0, so its guards never took the other branch. Review suggestions; no change in behavior intended. Co-Authored-By: Claude Opus 5.5 --- src/ForceFreeStates/EulerLagrange.jl | 74 ++++++++++++++-------------- 1 file changed, 37 insertions(+), 37 deletions(-) diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl index 792d1b435..9596ae262 100644 --- a/src/ForceFreeStates/EulerLagrange.jl +++ b/src/ForceFreeStates/EulerLagrange.jl @@ -133,7 +133,7 @@ and a small set of temporary matrices and factors used to compute singular-layer # Initialization parameters - - `zeroed_idx::Vector{Vector{Int}}` - For each ideal rational surface jump, a vector of indices of solutions that were zeroed. # Data for integrator + - `zeroed_idx::Vector{Vector{Int}}` - For each ideal crossing, the positions in `index` of the zeroed resonant solutions (`1:nres`). # Data for integrator - `fixfac::Array{ComplexF64,3}` - Fix-up factors for Gaussian reduction with shape `(numpert_total, numpert_total, numunorms_init)`. # Data for integrator @@ -765,10 +765,11 @@ function cross_ideal_singular_surf!( singp = intr.sing[ising] ipert_res = 1 .+ singp.m .- intr.mlow .+ (singp.n .- intr.nlow) .* intr.mpert + nres = length(ipert_res) # Fixup solution at singular surface; the ideal jump reduces on the resonant rows so the # solution it removes below is the resonant one - compute_solution_norms!(odet.u, odet, ctrl, intr, true; resonant_rows=(ctrl.kinetic_factor == 0 ? ipert_res : Int[])) + compute_solution_norms!(odet.u, odet, ctrl, intr, true; resonant_rows=ipert_res) # Compute direction-specific asymptotic power series for this singular surface sing_asymp_right = compute_sing_asymptotics(singp, ctrl, equil, mats, intr; sig=1.0) @@ -783,11 +784,9 @@ function cross_ideal_singular_surf!( # The reduction above placed the column pivoted on ipert_res[i] at index position i, which also # keeps each zeroed column in the same n-block as the resonant mode introduced below; zeroed_idx # records the positions for transforming u back to the full solution after integration. - if ctrl.kinetic_factor == 0 - odet.zeroed_idx[odet.ifix] = collect(eachindex(ipert_res)) - for i in eachindex(ipert_res) - odet.u[:, odet.index[i, odet.ifix], :] .= 0 - end + odet.zeroed_idx[odet.ifix] = collect(1:nres) + for i in 1:nres + odet.u[:, odet.index[i, odet.ifix], :] .= 0 end # Re-initialize on opposite side of rational surface by approximating solution @@ -801,13 +800,11 @@ function cross_ideal_singular_surf!( # Apply asymptotic solution on other side of singular surface (right side) ua = sing_get_ua(sing_asymp_right, dpsi) - if ctrl.kinetic_factor == 0 - for i in eachindex(sing_asymp_right.r1) - # Zero out the resonant components - odet.u[ipert_res[i], :, :] .= 0 - # Introduce the small asymptotic resonant solution on the other side of the singular surface - odet.u[:, odet.index[odet.zeroed_idx[odet.ifix][i], odet.ifix], :] .= ua[:, ipert_res[i]+intr.numpert_total, :] - end + for i in 1:nres + # Zero out the resonant components + odet.u[ipert_res[i], :, :] .= 0 + # Introduce the small asymptotic resonant solution on the other side of the singular surface + odet.u[:, odet.index[i, odet.ifix], :] .= ua[:, ipert_res[i]+intr.numpert_total, :] end # Get asymptotic coefficients after crossing rational surface odet.ca_r[:, :, :, ising] .= sing_get_ca(odet.u, ua, intr) @@ -1007,9 +1004,7 @@ operate on `u` directly without extra copies. - u: Current solution vector array, updated in-place if fixfac is called - sing_flag: Indicates if normalization is occuring at a singular surface or not - - resonant_rows: Rows of the resonant harmonics at an ideal crossing. When given, the reduction - always runs, even right after a `ucrit` reduction, and pivots on these rows (see - `apply_gaussian_reduction!`) + - resonant_rows: Resonant rows at an ideal crossing; when given, the reduction always runs and pivots on them ### TODOs @@ -1025,24 +1020,31 @@ 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 - # Normalize unorm and perform Gaussian reduction if required - if odet.new && isempty(resonant_rows) + # 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 odet.new = false odet.unorm0 .= odet.unorm - else - odet.new || (odet.unorm ./= odet.unorm0) - uratio = maximum(odet.unorm) / minimum(odet.unorm) - if uratio > ctrl.ucrit || sing_flag - # TODO: add resizing logic here as well - if odet.ifix < ctrl.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 - apply_gaussian_reduction!(u, odet, intr, sing_flag; resonant_rows) - odet.new = true + return + end + + # Growth since the last reduction (raw norms if a crossing comes right after one) + if !odet.new + odet.unorm ./= odet.unorm0 + end + uratio = maximum(odet.unorm) / minimum(odet.unorm) + + # Perform Gaussian reduction if the ucrit ratio is reached or at a singular surface + if uratio > ctrl.ucrit || sing_flag + # TODO: add resizing logic here as well + if odet.ifix < ctrl.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 + apply_gaussian_reduction!(u, odet, intr, sing_flag; resonant_rows) + odet.new = true end end @@ -1057,11 +1059,9 @@ This will update both `u` and relevant fields in `odet` in-place. See the description of `compute_solution_norms!` for more details on the benefits of in-place `u` updates. -Columns are triangularized in order of their growth since the last reduction. At an ideal -crossing, `resonant_rows` puts first, for each resonant row, the column with the largest -component there and pivots it on that row, so that column alone carries the resonant harmonic. -Growth order cannot identify the resonant solution when a `ucrit` reduction fired only a few -steps before the crossing, because every column has then grown alike. +At an ideal crossing, the first pivots are the `resonant_rows`, each on the remaining column largest +there. Growth order cannot pick the resonant solution if a `ucrit` reduction fired just before the +crossing, since every column has then grown alike. """ function apply_gaussian_reduction!(u::Array{ComplexF64,3}, odet::OdeState, intr::ForceFreeStatesInternal, sing_flag::Bool; resonant_rows=Int[]) From c9cb55d91103825f3ab7e97eebfc3ce39dc57274 Mon Sep 17 00:00:00 2001 From: d-burg Date: Fri, 9 Oct 2026 13:32:51 -0400 Subject: [PATCH 5/6] ForceFreeStates - REFACTOR - Select the resonant pivots inside the reduction loop One loop picks each resonant row's column from the columns still free at that step, then continues in growth order. Identical for a single resonant row; with several, later resonant columns are chosen after the earlier eliminations. Co-Authored-By: Claude Opus 5.5 --- src/ForceFreeStates/EulerLagrange.jl | 29 ++++++++++++++-------------- 1 file changed, 15 insertions(+), 14 deletions(-) diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl index 9596ae262..9a6cb836a 100644 --- a/src/ForceFreeStates/EulerLagrange.jl +++ b/src/ForceFreeStates/EulerLagrange.jl @@ -1076,24 +1076,25 @@ function apply_gaussian_reduction!(u::Array{ComplexF64,3}, odet::OdeState, intr: for isol in 1:intr.numpert_total odet.fixfac[isol, isol, ifix] = 1 end - # Resonant columns first, then the rest in descending unorm (we triangularize from largest to smallest) - lead = Int[] - for r in resonant_rows - push!(lead, argmax(j -> j in lead ? -Inf : abs(u[r, j, 1]), 1:intr.numpert_total)) - end - odet.index[:, ifix] = vcat(lead, filter(j -> !(j in lead), sortperm(odet.unorm; rev=true))) - # Triangularize primary solutions - mask = trues(2, intr.numpert_total) + # Triangularize primary solutions, resonant rows first, then columns from largest to smallest growth + growth_order = sortperm(odet.unorm; rev=true) + row_free = trues(intr.numpert_total) + col_free = trues(intr.numpert_total) for isol in 1:intr.numpert_total - ksol = odet.index[isol, ifix] - mask[2, ksol] = false - # Set pivot row: the resonant row for a leading resonant column, else the max location - @views kpert = isol <= length(lead) ? resonant_rows[isol] : argmax(abs.(u[:, ksol, 1]) .* mask[1, :]) - mask[1, kpert] = false + if isol <= length(resonant_rows) + kpert = resonant_rows[isol] + ksol = argmax(j -> col_free[j] ? abs(u[kpert, j, 1]) : -Inf, 1:intr.numpert_total) + else + ksol = growth_order[findfirst(j -> col_free[j], growth_order)] + @views kpert = argmax(abs.(u[:, ksol, 1]) .* row_free) + end + odet.index[isol, ifix] = ksol + col_free[ksol] = false + row_free[kpert] = false # Eliminate other solution vectors below the pivot for jsol in 1:intr.numpert_total - if mask[2, jsol] + if col_free[jsol] odet.fixfac[ksol, jsol, ifix] = -u[kpert, jsol, 1] / u[kpert, ksol, 1] @. @views u[:, jsol, :] .= u[:, jsol, :] .+ u[:, ksol, :] .* odet.fixfac[ksol, jsol, ifix] u[kpert, jsol, 1] = 0 From 50130bafbe8b949cb87a880a4dbbad7c76009c21 Mon Sep 17 00:00:00 2001 From: d-burg Date: Fri, 9 Oct 2026 13:32:51 -0400 Subject: [PATCH 6/6] ForceFreeStates - TEST - Check that et[1] is independent of the reduction frequency in the full-run tests Runs the coarse Solovev deck at ucrit = 2, 3, 5 and 10 against its default and replaces the solovev_n1_frequent_reduction harness case. Co-Authored-By: Claude Opus 5.5 --- .../cases/solovev_n1_frequent_reduction.toml | 52 ------------------- test/runtests_fullruns.jl | 20 +++++++ 2 files changed, 20 insertions(+), 52 deletions(-) delete mode 100644 regression-harness/cases/solovev_n1_frequent_reduction.toml diff --git a/regression-harness/cases/solovev_n1_frequent_reduction.toml b/regression-harness/cases/solovev_n1_frequent_reduction.toml deleted file mode 100644 index 5daecc6d2..000000000 --- a/regression-harness/cases/solovev_n1_frequent_reduction.toml +++ /dev/null @@ -1,52 +0,0 @@ -# Regression case: Solovev n=1 ideal stability with Gaussian reductions every few steps. -# Reuses the Solovev_ideal_example deck via [overrides] with ucrit = 2, so a ucrit reduction -# regularly fires a step or two before an ideal rational-surface crossing. The crossing must still -# remove the resonant solution; if it removes another one, et[1] lands far from solovev_n1's value. -# Each [quantities.*] block names an HDF5 path in the run output, how to extract it, -# and the noise floor below which a difference is treated as zero. -[case] -name = "solovev_n1_frequent_reduction" -description = "Solovev analytical equilibrium, n=1, ideal stability, ucrit = 2 so reductions land just before the rational-surface crossings" -example_dir = "examples/Solovev_ideal_example" - -# Reduce whenever the solution-norm spread exceeds 2, at a tight Euler-Lagrange tolerance. -[overrides] -"ForceFreeStates.ucrit" = 2.0 -"ForceFreeStates.eulerlagrange_tolerance" = 1e-10 - -# Energies — leading generalized (W,N) pencil eigenvalues at the final truncation (psilim). -[quantities.et_real] -h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_energies" -type = "complex_vector" -extract = "real_first" -label = "total energy Re(et[1])" -noise_threshold = 1e-10 -order = 10 -class = "physics_converged" - -[quantities.et_imag] -h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_energies" -type = "complex_vector" -extract = "imag_first" -label = "total energy Im(et[1])" -noise_threshold = 1e-10 -order = 11 -class = "physics_converged" - -[quantities.ep_real] -h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_plasma_energies" -type = "complex_vector" -extract = "real_first" -label = "plasma energy Re(ep[1])" -noise_threshold = 1e-10 -order = 12 -class = "physics_converged" - -[quantities.ev_real] -h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_vacuum_energies" -type = "complex_vector" -extract = "real_first" -label = "vacuum energy Re(ev[1])" -noise_threshold = 1e-10 -order = 13 -class = "physics_converged" diff --git a/test/runtests_fullruns.jl b/test/runtests_fullruns.jl index d652bb35b..305a23a64 100644 --- a/test/runtests_fullruns.jl +++ b/test/runtests_fullruns.jl @@ -1,4 +1,5 @@ using HDF5 +using TOML # Run GeneralizedPerturbedEquilibrium.main on the provided example directories and assert it completes without throwing. @testset "Full ForceFreeStates runs" begin @@ -9,6 +10,25 @@ using HDF5 true end + # A ucrit reduction landing just before an ideal crossing must not change which solution the crossing removes. + @info "Running Solovev ideal example with frequent Gaussian reductions" + @testset "et[1] is independent of the reduction frequency" begin + deck = TOML.parsefile(joinpath(ex1, "gpec.toml")) + deck["ForceFreeStates"]["write_outputs_to_HDF5"] = true + function total_energy(ucrit) + deck["ForceFreeStates"]["ucrit"] = ucrit + return mktempdir() do workdir + open(io -> TOML.print(io, deck), joinpath(workdir, "gpec.toml"), "w") + GeneralizedPerturbedEquilibrium.main([workdir]) + return h5open(h5 -> read(h5["ForceFreeStates/FreeBoundaryStability/eigenmode_energies"])[1], joinpath(workdir, "gpec.h5"), "r") + end + end + et_ref = total_energy(1e3) + for ucrit in (2.0, 3.0, 5.0, 10.0) + @test total_energy(ucrit) ≈ et_ref rtol = 1e-6 + end + end + ex2 = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example_multi_n") @info "Running Solovev ideal multi-n example" @test begin