diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl index 88c1d9e4a..de8d8a743 100644 --- a/src/ForceFreeStates/EulerLagrange.jl +++ b/src/ForceFreeStates/EulerLagrange.jl @@ -61,6 +61,29 @@ function ode_itime_cost(psi1::Float64, psi2::Float64, intr::ForceFreeStatesInter return cost end +""" + min_crossing_chunks(msing) -> Int + +Floor on the chunk count set by the `msing` rational-surface crossings; explicit `nchunks` +requests below it are clamped up (with a warning) by `balance_integration_chunks`. +""" +min_crossing_chunks(msing::Integer) = 2 * msing + 3 + +""" + auto_chunk_target(msing) -> Int + +The chunk-count target `balance_integration_chunks` resolves `nchunks = 0` to. Beyond the +crossing floor, BVP propagator conditioning needs at least 8 non-crossing sub-chunks per +segment (axis→surf₁, surfᵢ→surfᵢ₊₁, surfₙ→edge) plus the crossing chunks — without them, +`assemble_fm_matrix(condition=true)` cannot keep accumulated products well-conditioned +because single long-span propagators may already have cond ~ 10²⁴ (STRIDE uses 33 intervals +for comparable problems). Derived from `msing` alone — never from the thread count — which +is what keeps the chunk list, and hence every Riccati output, identical whatever `julia -t` +provides. Tests that steer decompositions relative to the auto target call this instead of +mirroring the formula. +""" +auto_chunk_target(msing::Integer) = max(min_crossing_chunks(msing), 8 * (msing + 1) + msing) + """ balance_integration_chunks(chunks, ctrl, intr) -> Vector{IntegrationChunk} @@ -79,20 +102,14 @@ each original chunk retains `needs_crossing=true` and the original `ising`, so t rational surface crossing still fires at the correct ψ in the serial assembly phase. """ function balance_integration_chunks(chunks::Vector{IntegrationChunk}, ctrl::ForceFreeStatesControl, intr::ForceFreeStatesInternal) - min_chunks = 2 * intr.msing + 3 - # Ensure enough sub-chunks for BVP propagator conditioning: at least 5 non-crossing - # sub-chunks per segment (axis→surf₁, surfᵢ→surfᵢ₊₁, surfₙ→edge), plus crossing - # chunks. STRIDE uses 33 intervals for comparable problems. Without enough sub-chunks, - # assemble_fm_matrix(condition=true) can't keep accumulated products well-conditioned - # because single long-span propagators may already have cond ~ 10²⁴. - min_bvp_intervals = 8 * (intr.msing + 1) + intr.msing + min_chunks = min_crossing_chunks(intr.msing) if ctrl.nchunks > 0 if ctrl.nchunks < min_chunks @warn "nchunks = $(ctrl.nchunks) is below the $min_chunks chunks required by $(intr.msing) singular surfaces; clamping up." end target_n = max(ctrl.nchunks, min_chunks) else - target_n = max(min_chunks, min_bvp_intervals) + target_n = auto_chunk_target(intr.msing) end result = collect(chunks) diff --git a/test/runtests_decomposition_invariance.jl b/test/runtests_decomposition_invariance.jl index 097af281d..84bcd9052 100644 --- a/test/runtests_decomposition_invariance.jl +++ b/test/runtests_decomposition_invariance.jl @@ -49,11 +49,10 @@ end @test auto.delta_prime !== nothing msing = size(auto.delta_prime.matrix, 1) - # The nchunks=0 target, mirroring balance_integration_chunks' internal formula (as - # runtests_parallel_integration.jl does). Invariance is asserted ABOVE this floor only: - # fewer chunks than the msing-derived minimum is not a different decomposition but a - # structurally deficient one (the floor gives the rational-surface crossings room). - auto_target = max(2 * msing + 3, 8 * (msing + 1) + msing) + # The nchunks=0 target, from its single source. Invariance is asserted ABOVE this floor + # only: fewer chunks than the msing-derived minimum is not a different decomposition but + # a structurally deficient one (the floor gives the rational-surface crossings room). + auto_target = GP_TI.ForceFreeStates.auto_chunk_target(msing) finer = _solve_at_nchunks(dir, auto_target + 11) @test finer.delta_prime !== nothing diff --git a/test/runtests_parallel_integration.jl b/test/runtests_parallel_integration.jl index 7e7737e06..0b3f4fe28 100644 --- a/test/runtests_parallel_integration.jl +++ b/test/runtests_parallel_integration.jl @@ -134,11 +134,11 @@ using TOML base_chunks = GeneralizedPerturbedEquilibrium.ForceFreeStates.chunk_el_integration_bounds(odet, ctrl, intr) balanced = GeneralizedPerturbedEquilibrium.ForceFreeStates.balance_integration_chunks(base_chunks, ctrl, intr) - # Must mirror balance_integration_chunks' internal target_n formula for nchunks = 0 - # (src/ForceFreeStates/EulerLagrange.jl). Keep this in sync. The formula reads only - # intr.msing — no thread count enters it, which is what makes Riccati outputs - # independent of how many threads `julia -t` provides. - target_n = max(2 * intr.msing + 3, 8 * (intr.msing + 1) + intr.msing) + # The nchunks = 0 target, from its single source. It reads only intr.msing — no + # thread count enters it, which is what makes Riccati outputs independent of how + # many threads `julia -t` provides. + target_n = GeneralizedPerturbedEquilibrium.ForceFreeStates.auto_chunk_target(intr.msing) + @test target_n == max(2 * intr.msing + 3, 8 * (intr.msing + 1) + intr.msing) # After balancing, chunk count equals target_n: the while-loop adds exactly one # chunk per iteration (a bisection split) and exits when length(result) >= target_n, @@ -183,7 +183,7 @@ using TOML @test length(balanced_more) == target_n + 7 # An nchunks below the singular-surface floor is clamped up, with a warning. - min_chunks = 2 * intr.msing + 3 + min_chunks = GeneralizedPerturbedEquilibrium.ForceFreeStates.min_crossing_chunks(intr.msing) ctrl_few = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])..., nchunks=1) balanced_few = @test_logs (:warn,) match_mode=:any GeneralizedPerturbedEquilibrium.ForceFreeStates.balance_integration_chunks(base_chunks, ctrl_few, intr)