From 2296514686f65124297593f5154dc5f3d2242e5d Mon Sep 17 00:00:00 2001 From: d-burg Date: Wed, 19 Aug 2026 13:53:28 -0400 Subject: [PATCH] FFS - BUG FIX - Per-column absolute tolerance in the Riccati shooting integration MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The shooting solves batched all columns under one uniform absolute tolerance (the ODE library default, 1e-6), so the resonant big-solution column set the error floor of the small-solution column and the extracted Δ' inherited an absolute error that scaled with the matching radius (dpsi^~2α). Scale abstol per column to the column's own magnitude, mirroring the reference implementation's per-column atol policy. Δ' becomes matching-radius independent: on the TJ circular benchmark the 3/1 diagonal moves 0.909 → 1.054 at the default singfac_min=1e-4, in agreement with the Galerkin path (1.057), Fortran STRIDE (1.05), and TJ (1.067); the 2/1 and all shaped-case surfaces move ≤0.1%. Co-Authored-By: Claude Opus 5 --- src/ForceFreeStates/Riccati.jl | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/src/ForceFreeStates/Riccati.jl b/src/ForceFreeStates/Riccati.jl index c1710cc12..b1e29f0a7 100644 --- a/src/ForceFreeStates/Riccati.jl +++ b/src/ForceFreeStates/Riccati.jl @@ -1405,20 +1405,29 @@ function integrate_fm_with_ua_ic( u0 = zeros(ComplexF64, N, N, 2) u0[:, :, 1] .= ua[:, 1:N, 1] u0[:, :, 2] .= ua[:, 1:N, 2] + # Per-column absolute tolerance so a batch's largest column cannot set the error floor of its smallest: + # otherwise the resonant small-solution column inherits an absolute error set by the big solution's magnitude. + abstol_arr = similar(u0, Float64) + for j in 1:N + abstol_arr[:, j, :] .= max(maximum(abs, @view u0[:, j, :]), 1e-30) * rtol + end odet_proxy.spline_hint[] = 1 odet_proxy.ffit_hint[] = 1 prob = ODEProblem(sing_der!, u0, tspan, params) - sol = solve(prob, Vern9(); reltol=rtol, save_everystep=false, save_end=true) + sol = solve(prob, Vern9(); reltol=rtol, abstol=abstol_arr, save_everystep=false, save_end=true) result[1:N, 1:N] .= sol.u[end][:, :, 1] result[N+1:2N, 1:N] .= sol.u[end][:, :, 2] # Batch 2: columns N+1:2N of T (small solutions) u0[:, :, 1] .= ua[:, N+1:2N, 1] u0[:, :, 2] .= ua[:, N+1:2N, 2] + for j in 1:N + abstol_arr[:, j, :] .= max(maximum(abs, @view u0[:, j, :]), 1e-30) * rtol + end odet_proxy.spline_hint[] = 1 odet_proxy.ffit_hint[] = 1 prob = ODEProblem(sing_der!, u0, tspan, params) - sol = solve(prob, Vern9(); reltol=rtol, save_everystep=false, save_end=true) + sol = solve(prob, Vern9(); reltol=rtol, abstol=abstol_arr, save_everystep=false, save_end=true) result[1:N, N+1:2N] .= sol.u[end][:, :, 1] result[N+1:2N, N+1:2N] .= sol.u[end][:, :, 2]