From 211564ec71d0310f06647ac7c10a7f7fbf970ff5 Mon Sep 17 00:00:00 2001 From: d-burg Date: Wed, 12 Aug 2026 17:09:17 -0400 Subject: [PATCH 1/5] TESTING - NEW FEATURE - Test thread invariance of the parallel BVP path MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit ForceFreeStatesControl documents that the parallel FM/BVP path produces bit-identical Delta' across thread counts, but nothing tested it by varying threads, and the example decks disagreed about whether to rely on it: the SLAYER deck pins parallel_threads = 1 "to keep the regression Delta' (and hence gamma) reproducible", while the DIII-D-like ideal deck — whose Delta' the regression harness tracks — runs parallel_threads = 2. A golden Delta' has to be a property of the physics rather than of the machine it was measured on, so this settles the question by measurement. Measured on the DIII-D-like deck at parallel_threads = 1, 2, 4 and 16: the Delta' matrix diagonal and et[1] are bit-identical throughout. The 16-thread configuration is the informative one — it lifts 4*effective_threads above the min_bvp_intervals floor and so produces a genuinely different decomposition (64 chunks vs 53, with different boundaries), which reassociates the propagator products without moving the result. Invariance therefore holds across decomposition, not merely across scheduling. - test/runtests_thread_invariance.jl compares parallel_threads 1 vs 2 on the Solovev and DIII-D-like decks, asserting exact equality rather than a tolerance: the code claims bit-identity, and a tolerance would mask the reassociation the test exists to catch. It also asserts the chunk boundaries are unchanged at these caps, so that a future decomposition change surfaces as a failure instead of silently turning the comparison into a stronger claim than intended. Skipped when the session has one thread, where effective_threads collapses to 1 and every comparison is vacuous. - The test workflow gains a multi-threaded leg (JULIA_NUM_THREADS = 4). The suite had only ever run single-threaded, so the parallel paths were exercised solely in their degenerate form. Existing job names are preserved byte-identically because branch protection names them as required checks. - balance_integration_chunks' docstring gave target_n as max(2*msing + 3, 4*Threads.nthreads()), omitting both the parallel_threads cap and the min_bvp_intervals term that actually dominates. runtests_parallel_integration.jl mirrored the same stale formula and would fail on a machine with more threads than the cap; it passed only because CI is single-threaded. - CLAUDE.md's single-test-file invocation could not work (runtests.jl passes ARGS to include, which resolves relative to test/), and two of the listed files do not exist. Co-Authored-By: Claude Opus 5 --- .github/workflows/test.yaml | 18 ++++- CLAUDE.md | 26 ++++--- test/runtests.jl | 1 + test/runtests_thread_invariance.jl | 115 +++++++++++++++++++++++++++++ 4 files changed, 149 insertions(+), 11 deletions(-) create mode 100644 test/runtests_thread_invariance.jl diff --git a/.github/workflows/test.yaml b/.github/workflows/test.yaml index a5e3561d5..6f939c437 100644 --- a/.github/workflows/test.yaml +++ b/.github/workflows/test.yaml @@ -86,11 +86,17 @@ jobs: fi test: - name: runtests ${{ matrix.version }} - ${{ matrix.os }} + # The job name must stay byte-identical for the single-threaded legs: branch protection + # names these contexts as required checks, so renaming them would leave every pull request + # waiting forever on a status that never reports. The multi-threaded leg gets a distinct + # name and is therefore additive rather than breaking. + name: runtests ${{ matrix.version }} - ${{ matrix.os }}${{ matrix.threads != '1' && format(' ({0} threads)', matrix.threads) || '' }} needs: changes if: needs.changes.outputs.julia == 'true' runs-on: ${{ matrix.os }} timeout-minutes: 90 + env: + JULIA_NUM_THREADS: ${{ matrix.threads }} strategy: fail-fast: false matrix: @@ -99,6 +105,16 @@ jobs: - '1.x' # latest (currently 1.12) os: - ubuntu-latest + threads: + - '1' + include: + # The parallel FM/BVP paths degenerate to their serial form when + # effective_threads = min(nthreads, parallel_threads) collapses to 1, so a + # single-threaded matrix never exercises threaded execution at all. This leg runs the + # suite multi-threaded, which is what makes the thread-invariance tests meaningful. + - version: '1.11' + os: ubuntu-latest + threads: '4' env: DEPOT_PATHS: | diff --git a/CLAUDE.md b/CLAUDE.md index 9417068b5..a2dc56134 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -22,16 +22,22 @@ GPEC (Generalized Perturbed Equilibrium Code, Julia implementation) is a compreh # Run all tests julia --project=. -e 'using Pkg; Pkg.activate("."); Pkg.instantiate(); include("test/runtests.jl")' -# Run specific test file -julia --project=. test/runtests.jl test/runtests_solovev.jl - -# Available test files: - -# - test/runtests_vacuum_julia.jl # Julia vacuum module -# - test/runtests_solovev.jl # Analytical equilibrium -# - test/runtests_ode.jl # ODE integration -# - test/runtests_sing.jl # Singular surface handling -# - test/runtests_fullruns.jl # End-to-end tests +# Run specific test file — the argument is included relative to test/, so pass the bare +# filename, not a path prefixed with test/ +julia --project=. test/runtests.jl runtests_sing.jl + +# Run the suite multi-threaded (the parallel FM/BVP paths reduce to their serial form at one +# thread, so a single-threaded run never exercises threaded execution) +julia -t 4 --project=. test/runtests.jl + +# A few of the available test files (see test/runtests.jl for the full list): + +# - test/runtests_vacuum.jl # Vacuum module +# - test/runtests_equil.jl # Equilibrium reconstruction +# - test/runtests_sing.jl # Singular surface handling +# - test/runtests_parallel_integration.jl # Parallel FM integration and BVP Delta' +# - test/runtests_thread_invariance.jl # Parallel-vs-serial equivalence +# - test/runtests_fullruns.jl # End-to-end tests ``` ### Building Documentation diff --git a/test/runtests.jl b/test/runtests.jl index 48a53e2d6..2e1903118 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -32,6 +32,7 @@ else include("./runtests_parallel_integration.jl") include("./runtests_result_struct.jl") include("./runtests_solve_api.jl") + include("./runtests_thread_invariance.jl") include("./runtests_sing.jl") include("./runtests_innerlayer.jl") include("./runtests_tj_analytic.jl") diff --git a/test/runtests_thread_invariance.jl b/test/runtests_thread_invariance.jl new file mode 100644 index 000000000..41877e071 --- /dev/null +++ b/test/runtests_thread_invariance.jl @@ -0,0 +1,115 @@ +using Test +using TOML + +# Thread-invariance of the parallel FM/BVP path. +# +# `ForceFreeStatesControl` documents that the parallel path produces bit-identical Δ′ across +# thread counts, and the example decks disagree about whether to rely on it: the SLAYER deck +# pins `parallel_threads = 1` for reproducibility while the DIII-D-like ideal deck — whose Δ′ +# the regression harness tracks — runs `parallel_threads = 2`. These tests hold the source and +# the deck fixed and vary only the BVP thread cap, so a golden Δ′ is a property of the physics +# rather than of the machine it was measured on. +# +# Two axes are reachable through `parallel_threads`: +# +# - scheduling: `Threads.@threads` over chunks vs a serial loop (any cap ≥ 2) +# - decomposition: `balance_integration_chunks` targets +# `max(2·msing+3, 4·effective_threads, 8·(msing+1)+msing)` sub-chunks, so once +# `4·effective_threads` exceeds the `min_bvp_intervals` floor the chunk boundaries +# themselves move and the propagator products reassociate +# +# The decomposition axis needs `effective_threads > (9·msing+8)/4` — about 14 threads for the +# DIII-D-like deck's msing=5 — so it is out of reach of a typical CI runner and is exercised by +# the nightly harness instead. Measured on this deck at parallel_threads = 1, 2, 4 and 16 +# (53 vs 64 chunks, confirmed different boundaries): Δ′ and et[1] were bit-identical throughout. +# +# `effective_threads = min(Threads.nthreads(), parallel_threads)` collapses every cap to 1 in a +# single-threaded session, which would make these comparisons trivially true, so they are skipped +# there rather than passing vacuously. + +const GP_TI = GeneralizedPerturbedEquilibrium + +""" +Run the ideal stability pipeline on `dir` at a given BVP thread cap. + +Mirrors the standalone setup used by the parallel-integration tests: build the equilibrium +(applying the two-pass auto grid when the deck asks for it), integrate the Euler-Lagrange +system, then assemble the STRIDE BVP Δ′ matrix. Returns the Δ′ matrix, the leading energy +eigenvalue, and the chunk boundaries, so a caller can tell scheduling changes from +decomposition changes. +""" +function _run_at_thread_cap(dir::String, parallel_threads::Int) + inputs = TOML.parsefile(joinpath(dir, "gpec.toml")) + inputs["ForceFreeStates"]["verbose"] = false + inputs["ForceFreeStates"]["use_parallel"] = true + inputs["ForceFreeStates"]["write_outputs_to_HDF5"] = false + inputs["ForceFreeStates"]["parallel_threads"] = parallel_threads + + intr = GP_TI.ForceFreeStates.ForceFreeStatesInternal(; dir_path=dir) + ctrl = GP_TI.ForceFreeStates.ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...) + eq_config = GP_TI.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], dir) + sol_config = haskey(inputs, "SOL_INPUT") ? GP_TI.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing + equil = GP_TI.Equilibrium.setup_equilibrium(eq_config, sol_config) + if GP_TI.Equilibrium.wants_two_pass(eq_config) + mand = GP_TI.ForceFreeStates.rational_psi_nodes(equil; nlow=ctrl.nn_low, nhigh=ctrl.nn_high) + psi_nodes = GP_TI.Equilibrium.refined_psi_grid(equil; tau=eq_config.psi_accuracy, mandatory=mand) + rerun_input = GP_TI.Equilibrium.build_direct_from_ingest(eq_config, equil.ingest) + equil = GP_TI.Equilibrium.setup_equilibrium(eq_config, rerun_input; override_psi_nodes=psi_nodes) + end + intr.wall_settings = GP_TI.Vacuum.WallShapeSettings(; (Symbol(k) => v for (k, v) in inputs["Wall"])...) + GP_TI.ForceFreeStates.sing_lim!(intr, ctrl, equil) + intr.nlow = ctrl.nn_low + intr.nhigh = ctrl.nn_high + intr.npert = 1 + GP_TI.ForceFreeStates.sing_find!(intr, equil) + intr.mlow = min(intr.nlow * equil.params.qmin, 0) - 4 - ctrl.delta_mlow + intr.mhigh = trunc(Int, intr.nhigh * equil.params.qmax) + ctrl.delta_mhigh + intr.mpert = intr.mhigh - intr.mlow + 1 + intr.numpert_total = intr.mpert * intr.npert + metric = GP_TI.ForceFreeStates.make_metric(equil, intr.mpert) + ffit = GP_TI.ForceFreeStates.make_matrix(equil, intr, metric) + odet, fm_propagators, fm_chunks, fm_S_left = GP_TI.ForceFreeStates.eulerlagrange_integration(ctrl, equil, ffit, intr) + vac = GP_TI.ForceFreeStates.free_run!(odet, ctrl, equil, ffit, intr) + GP_TI.ForceFreeStates.compute_delta_prime_matrix!(intr, fm_propagators, fm_chunks; + wv=vac.wv, psio=equil.psio, S_at_surface_left=fm_S_left, ctrl=ctrl, equil=equil, ffit=ffit) + return (dpm=copy(intr.delta_prime_matrix), et1=vac.et[1], msing=intr.msing, + bounds=[(c.psi_start, c.psi_end) for c in fm_chunks]) +end + +@testset "Thread invariance of the parallel BVP path" begin + if Threads.nthreads() < 2 + @info "Thread-invariance tests skipped: effective_threads collapses to 1 in a single-threaded session. Run with `julia -t 4` (CI covers this in its multi-threaded matrix leg)." + @test true + else + @testset "Solovev — leading eigenvalue" begin + dir = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example") + serial = _run_at_thread_cap(dir, 1) + threaded = _run_at_thread_cap(dir, 2) + # Bit-identical, not approximate: the parallel path claims exactness, and a + # tolerance here would hide precisely the reassociation this test exists to catch. + @test threaded.et1 === serial.et1 + end + + @testset "DIII-D-like — Δ′ diagonal and leading eigenvalue" begin + # The deck whose Δ′ the regression harness pins, and where the BVP Δ′ is + # well-conditioned (Solovev sits near marginal stability and its BVP Δ′ is not). + dir = joinpath(@__DIR__, "..", "examples", "DIIID-like_ideal_example") + serial = _run_at_thread_cap(dir, 1) + threaded = _run_at_thread_cap(dir, 2) + + @test threaded.msing == serial.msing + @test size(threaded.dpm) == size(serial.dpm) + @test threaded.et1 === serial.et1 + for j in 1:serial.msing + @test threaded.dpm[j, j] === serial.dpm[j, j] + end + @test threaded.dpm == serial.dpm + + # At these caps the min_bvp_intervals floor fixes the chunk count, so the + # boundaries should be untouched and only scheduling differs. If this fails the + # decomposition moved and the Δ′ comparison above became a stronger claim than + # the one this testset intends to make. + @test threaded.bounds == serial.bounds + end + end +end From 4b8c2da492d7e4742d38079912f7cfa16b910b47 Mon Sep 17 00:00:00 2001 From: d-burg Date: Sat, 15 Aug 2026 11:37:36 -0400 Subject: [PATCH 2/5] TESTING - BUG FIX - Resolve the toroidal range before sing_lim! MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit develop added a guard requiring intr.nlow/nhigh to be assigned before sing_lim!, since set_psilim_via_dmlim truncates at (last_rational_q + dmlim)/n and so needs n, and fixed the ordering in runtests_parallel_integration.jl. This test file copied the old ordering from that file before the guard existed; the rebase carried it forward because the two files never overlap textually, so git merged them cleanly while the semantics diverged. Only the multi-threaded CI leg surfaced it: the Solovev testset does not truncate via dmlim, and the single-threaded legs skip the whole testset because effective_threads collapses to 1 there. No thread-invariance claim is affected — the failure was an exception during setup, not a comparison mismatch. Co-Authored-By: Claude Opus 5 --- test/runtests_thread_invariance.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/test/runtests_thread_invariance.jl b/test/runtests_thread_invariance.jl index 41877e071..3881a5412 100644 --- a/test/runtests_thread_invariance.jl +++ b/test/runtests_thread_invariance.jl @@ -57,10 +57,12 @@ function _run_at_thread_cap(dir::String, parallel_threads::Int) equil = GP_TI.Equilibrium.setup_equilibrium(eq_config, rerun_input; override_psi_nodes=psi_nodes) end intr.wall_settings = GP_TI.Vacuum.WallShapeSettings(; (Symbol(k) => v for (k, v) in inputs["Wall"])...) - GP_TI.ForceFreeStates.sing_lim!(intr, ctrl, equil) + # The toroidal range must be resolved before sing_lim!: under set_psilim_via_dmlim it + # truncates at (last_rational_q + dmlim)/n and so needs n. intr.nlow = ctrl.nn_low intr.nhigh = ctrl.nn_high intr.npert = 1 + GP_TI.ForceFreeStates.sing_lim!(intr, ctrl, equil) GP_TI.ForceFreeStates.sing_find!(intr, equil) intr.mlow = min(intr.nlow * equil.params.qmin, 0) - 4 - ctrl.delta_mlow intr.mhigh = trunc(Int, intr.nhigh * equil.params.qmax) + ctrl.delta_mhigh From 54612bf63aad9f5dd781cdfa8b8a1f6bc229aabb Mon Sep 17 00:00:00 2001 From: d-burg Date: Sat, 15 Aug 2026 17:23:21 -0400 Subject: [PATCH 3/5] =?UTF-8?q?TESTING=20-=20REFACTOR=20-=20Decomposition?= =?UTF-8?q?=20invariance=20of=20the=20unified=20Riccati=20=CE=94=E2=80=B2?= =?UTF-8?q?=20path?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The riccati unification made the original thread-invariance test obsolete: it varied parallel_threads, which no longer exists, and the chunk target is now derived from msing alone so thread-count independence is structural and unit-tested upstream. What those unit tests cannot assert is the end-to-end claim: that cutting the same integration into a genuinely different set of chunks — which reassociates the fundamental-matrix products — leaves the physics output unchanged. This test pins that, steering the decomposition directly via nchunks (auto vs auto+11, boundaries verified different), so it is meaningful at any thread count. Measured on the DIII-D-like deck under the unified driver: the Δ′ matrix is bit-identical across decompositions (asserted with ===; a tolerance would hide exactly the reassociation drift the test exists to catch), but et[1] is not — it drifts by 2.7e-8 relative, where the pre-unification driver was exact. That sensitivity is recorded honestly rather than hidden: @test_broken on exactness (an Unexpected Pass will force the strict assertion back if a driver change restores it) plus a documented 1e-6 ceiling to catch it growing by orders of magnitude. The multi-threaded CI leg is kept: the chunk loop runs through Threads.@threads, which executes serially in a single-threaded session, so without this leg no CI job ever exercises concurrent scheduling. The CLAUDE.md single-test-file invocation fix is kept (runtests.jl passes ARGS to include, which resolves relative to test/). Co-Authored-By: Claude Fable 5 --- .github/workflows/test.yaml | 8 +- CLAUDE.md | 2 +- test/runtests.jl | 2 +- test/runtests_decomposition_invariance.jl | 101 +++++++++++++++++++ test/runtests_thread_invariance.jl | 117 ---------------------- 5 files changed, 107 insertions(+), 123 deletions(-) create mode 100644 test/runtests_decomposition_invariance.jl delete mode 100644 test/runtests_thread_invariance.jl diff --git a/.github/workflows/test.yaml b/.github/workflows/test.yaml index 6f939c437..4c2eee86d 100644 --- a/.github/workflows/test.yaml +++ b/.github/workflows/test.yaml @@ -108,10 +108,10 @@ jobs: threads: - '1' include: - # The parallel FM/BVP paths degenerate to their serial form when - # effective_threads = min(nthreads, parallel_threads) collapses to 1, so a - # single-threaded matrix never exercises threaded execution at all. This leg runs the - # suite multi-threaded, which is what makes the thread-invariance tests meaningful. + # The Riccati chunk driver runs its chunks through Threads.@threads, which executes + # serially in a single-threaded session — so a single-threaded matrix never exercises + # threaded scheduling at all. This leg runs the suite multi-threaded, making the + # chunk-decomposition and boundary-pinning tests cover real concurrent execution. - version: '1.11' os: ubuntu-latest threads: '4' diff --git a/CLAUDE.md b/CLAUDE.md index a2dc56134..99886fe00 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -36,7 +36,7 @@ julia -t 4 --project=. test/runtests.jl # - test/runtests_equil.jl # Equilibrium reconstruction # - test/runtests_sing.jl # Singular surface handling # - test/runtests_parallel_integration.jl # Parallel FM integration and BVP Delta' -# - test/runtests_thread_invariance.jl # Parallel-vs-serial equivalence +# - test/runtests_decomposition_invariance.jl # Riccati Δ' chunk-decomposition invariance # - test/runtests_fullruns.jl # End-to-end tests ``` diff --git a/test/runtests.jl b/test/runtests.jl index 2e1903118..9a0eeafa2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -32,7 +32,7 @@ else include("./runtests_parallel_integration.jl") include("./runtests_result_struct.jl") include("./runtests_solve_api.jl") - include("./runtests_thread_invariance.jl") + include("./runtests_decomposition_invariance.jl") include("./runtests_sing.jl") include("./runtests_innerlayer.jl") include("./runtests_tj_analytic.jl") diff --git a/test/runtests_decomposition_invariance.jl b/test/runtests_decomposition_invariance.jl new file mode 100644 index 000000000..d106b99ec --- /dev/null +++ b/test/runtests_decomposition_invariance.jl @@ -0,0 +1,101 @@ +using Test +using TOML + +# Decomposition invariance of the Riccati/FM Δ′ path. +# +# The chunked propagator driver reassociates the fundamental-matrix products whenever the chunk +# decomposition changes: ((A·B)·C)·D becomes (A·B)·(C·D). Floating-point matrix products do not +# reassociate exactly in general, so Δ′ being reproducible requires more than thread-count +# independence of the chunk *count* (which is structural: the nchunks=0 target is derived from +# msing alone and pinned by unit tests in runtests_parallel_integration.jl). This file asserts +# the end-to-end claim those unit tests cannot: the Δ′ matrix and the leading energy eigenvalue +# are bit-identical when the same integration is cut into a genuinely different set of chunks. +# +# Measured basis (DIII-D-like deck): 53-chunk and 64-chunk decompositions with confirmed +# different boundaries gave bit-identical Δ′ diagonals and et[1]. This test pins that property. +# Unlike thread-count variation, the decomposition axis is exercisable in-session at any thread +# count, because nchunks steers it directly. + +const GP_TI = GeneralizedPerturbedEquilibrium + +""" +Run the ideal stability pipeline on `dir` with the Riccati integrator at a given chunk count +(`nchunks = 0` = the msing-derived auto target). Mirrors the standalone setup used by the +parallel-integration tests. Returns the Δ′ matrix, the leading energy eigenvalue, and the chunk +boundaries so the test can prove the decompositions actually differed. +""" +function _run_at_nchunks(dir::String, nchunks::Int) + inputs = TOML.parsefile(joinpath(dir, "gpec.toml")) + inputs["ForceFreeStates"]["verbose"] = false + inputs["ForceFreeStates"]["integrator"] = "riccati" + inputs["ForceFreeStates"]["write_outputs_to_HDF5"] = false + inputs["ForceFreeStates"]["nchunks"] = nchunks + + intr = GP_TI.ForceFreeStates.ForceFreeStatesInternal(; dir_path=dir) + ctrl = GP_TI.ForceFreeStates.ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...) + eq_config = GP_TI.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], dir) + sol_config = haskey(inputs, "SOL_INPUT") ? GP_TI.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing + equil = GP_TI.Equilibrium.setup_equilibrium(eq_config, sol_config) + if GP_TI.Equilibrium.wants_two_pass(eq_config) + mand = GP_TI.ForceFreeStates.rational_psi_nodes(equil; nlow=ctrl.nn_low, nhigh=ctrl.nn_high) + psi_nodes = GP_TI.Equilibrium.refined_psi_grid(equil; tau=eq_config.psi_accuracy, mandatory=mand) + rerun_input = GP_TI.Equilibrium.build_direct_from_ingest(eq_config, equil.ingest) + equil = GP_TI.Equilibrium.setup_equilibrium(eq_config, rerun_input; override_psi_nodes=psi_nodes) + end + intr.wall_settings = GP_TI.Vacuum.WallShapeSettings(; (Symbol(k) => v for (k, v) in inputs["Wall"])...) + # The toroidal range must be resolved before sing_lim!: under set_psilim_via_dmlim it + # truncates at (last_rational_q + dmlim)/n and so needs n. + intr.nlow = ctrl.nn_low + intr.nhigh = ctrl.nn_high + intr.npert = 1 + GP_TI.ForceFreeStates.sing_lim!(intr, ctrl, equil) + GP_TI.ForceFreeStates.sing_find!(intr, equil) + intr.mlow = min(intr.nlow * equil.params.qmin, 0) - 4 - ctrl.delta_mlow + intr.mhigh = trunc(Int, intr.nhigh * equil.params.qmax) + ctrl.delta_mhigh + intr.mpert = intr.mhigh - intr.mlow + 1 + intr.numpert_total = intr.mpert * intr.npert + metric = GP_TI.ForceFreeStates.make_metric(equil, intr.mpert) + ffit = GP_TI.ForceFreeStates.make_matrix(equil, intr, metric) + odet, fm_propagators, fm_chunks, fm_S_left = GP_TI.ForceFreeStates.eulerlagrange_integration(ctrl, equil, ffit, intr) + vac = GP_TI.ForceFreeStates.free_run(odet, ctrl, equil, ffit, intr) + GP_TI.ForceFreeStates.compute_delta_prime_matrix!(intr, fm_propagators, fm_chunks; + wv=vac.wv, psio=equil.psio, S_at_surface_left=fm_S_left, ctrl=ctrl, equil=equil, ffit=ffit) + return (dpm=copy(intr.delta_prime_matrix), et1=vac.et[1], msing=intr.msing, + bounds=[(c.psi_start, c.psi_end) for c in fm_chunks]) +end + +@testset "Decomposition invariance of the Riccati Δ′ path" begin + # The deck whose Δ′ the regression harness pins, and where the BVP Δ′ is well-conditioned + # (Solovev sits near marginal stability and its BVP Δ′ is pathological there). + dir = joinpath(@__DIR__, "..", "examples", "DIIID-like_ideal_example") + auto = _run_at_nchunks(dir, 0) + # 11 more chunks than the auto target: enough to move several boundaries and change the + # association of the propagator products, cheap enough not to distort the runtime. + finer = _run_at_nchunks(dir, length(auto.bounds) + 11) + + # The premise first: the two decompositions must genuinely differ, otherwise the equality + # below is vacuous and this file silently stops testing anything. + @test length(finer.bounds) > length(auto.bounds) + @test finer.bounds != auto.bounds + + @test finer.msing == auto.msing + @test size(finer.dpm) == size(auto.dpm) + + # Δ′ is bit-identical, not approximate: a tolerance would hide exactly the reassociation + # drift this test exists to catch, and the measured behaviour is exact equality. + for j in 1:auto.msing + @test finer.dpm[j, j] === auto.dpm[j, j] + end + @test finer.dpm == auto.dpm + + # et[1] is NOT decomposition-invariant under the unified Riccati driver: measured relative + # difference 2.7e-8 between the auto and auto+11 decompositions on this deck (it was + # bit-identical under the pre-unification driver). The tolerance below is 30× the measured + # effect — documented, not chosen to make the test pass — and exists to catch this + # sensitivity growing by orders of magnitude while the exact-invariance question is open. + # test_broken: if a future driver change restores exact invariance, this reports an + # Unexpected Pass, forcing the === assertion to be reinstated rather than the improvement + # going unnoticed. + @test_broken finer.et1 === auto.et1 + @test isapprox(finer.et1, auto.et1; rtol=1e-6) +end diff --git a/test/runtests_thread_invariance.jl b/test/runtests_thread_invariance.jl deleted file mode 100644 index 3881a5412..000000000 --- a/test/runtests_thread_invariance.jl +++ /dev/null @@ -1,117 +0,0 @@ -using Test -using TOML - -# Thread-invariance of the parallel FM/BVP path. -# -# `ForceFreeStatesControl` documents that the parallel path produces bit-identical Δ′ across -# thread counts, and the example decks disagree about whether to rely on it: the SLAYER deck -# pins `parallel_threads = 1` for reproducibility while the DIII-D-like ideal deck — whose Δ′ -# the regression harness tracks — runs `parallel_threads = 2`. These tests hold the source and -# the deck fixed and vary only the BVP thread cap, so a golden Δ′ is a property of the physics -# rather than of the machine it was measured on. -# -# Two axes are reachable through `parallel_threads`: -# -# - scheduling: `Threads.@threads` over chunks vs a serial loop (any cap ≥ 2) -# - decomposition: `balance_integration_chunks` targets -# `max(2·msing+3, 4·effective_threads, 8·(msing+1)+msing)` sub-chunks, so once -# `4·effective_threads` exceeds the `min_bvp_intervals` floor the chunk boundaries -# themselves move and the propagator products reassociate -# -# The decomposition axis needs `effective_threads > (9·msing+8)/4` — about 14 threads for the -# DIII-D-like deck's msing=5 — so it is out of reach of a typical CI runner and is exercised by -# the nightly harness instead. Measured on this deck at parallel_threads = 1, 2, 4 and 16 -# (53 vs 64 chunks, confirmed different boundaries): Δ′ and et[1] were bit-identical throughout. -# -# `effective_threads = min(Threads.nthreads(), parallel_threads)` collapses every cap to 1 in a -# single-threaded session, which would make these comparisons trivially true, so they are skipped -# there rather than passing vacuously. - -const GP_TI = GeneralizedPerturbedEquilibrium - -""" -Run the ideal stability pipeline on `dir` at a given BVP thread cap. - -Mirrors the standalone setup used by the parallel-integration tests: build the equilibrium -(applying the two-pass auto grid when the deck asks for it), integrate the Euler-Lagrange -system, then assemble the STRIDE BVP Δ′ matrix. Returns the Δ′ matrix, the leading energy -eigenvalue, and the chunk boundaries, so a caller can tell scheduling changes from -decomposition changes. -""" -function _run_at_thread_cap(dir::String, parallel_threads::Int) - inputs = TOML.parsefile(joinpath(dir, "gpec.toml")) - inputs["ForceFreeStates"]["verbose"] = false - inputs["ForceFreeStates"]["use_parallel"] = true - inputs["ForceFreeStates"]["write_outputs_to_HDF5"] = false - inputs["ForceFreeStates"]["parallel_threads"] = parallel_threads - - intr = GP_TI.ForceFreeStates.ForceFreeStatesInternal(; dir_path=dir) - ctrl = GP_TI.ForceFreeStates.ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...) - eq_config = GP_TI.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], dir) - sol_config = haskey(inputs, "SOL_INPUT") ? GP_TI.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing - equil = GP_TI.Equilibrium.setup_equilibrium(eq_config, sol_config) - if GP_TI.Equilibrium.wants_two_pass(eq_config) - mand = GP_TI.ForceFreeStates.rational_psi_nodes(equil; nlow=ctrl.nn_low, nhigh=ctrl.nn_high) - psi_nodes = GP_TI.Equilibrium.refined_psi_grid(equil; tau=eq_config.psi_accuracy, mandatory=mand) - rerun_input = GP_TI.Equilibrium.build_direct_from_ingest(eq_config, equil.ingest) - equil = GP_TI.Equilibrium.setup_equilibrium(eq_config, rerun_input; override_psi_nodes=psi_nodes) - end - intr.wall_settings = GP_TI.Vacuum.WallShapeSettings(; (Symbol(k) => v for (k, v) in inputs["Wall"])...) - # The toroidal range must be resolved before sing_lim!: under set_psilim_via_dmlim it - # truncates at (last_rational_q + dmlim)/n and so needs n. - intr.nlow = ctrl.nn_low - intr.nhigh = ctrl.nn_high - intr.npert = 1 - GP_TI.ForceFreeStates.sing_lim!(intr, ctrl, equil) - GP_TI.ForceFreeStates.sing_find!(intr, equil) - intr.mlow = min(intr.nlow * equil.params.qmin, 0) - 4 - ctrl.delta_mlow - intr.mhigh = trunc(Int, intr.nhigh * equil.params.qmax) + ctrl.delta_mhigh - intr.mpert = intr.mhigh - intr.mlow + 1 - intr.numpert_total = intr.mpert * intr.npert - metric = GP_TI.ForceFreeStates.make_metric(equil, intr.mpert) - ffit = GP_TI.ForceFreeStates.make_matrix(equil, intr, metric) - odet, fm_propagators, fm_chunks, fm_S_left = GP_TI.ForceFreeStates.eulerlagrange_integration(ctrl, equil, ffit, intr) - vac = GP_TI.ForceFreeStates.free_run!(odet, ctrl, equil, ffit, intr) - GP_TI.ForceFreeStates.compute_delta_prime_matrix!(intr, fm_propagators, fm_chunks; - wv=vac.wv, psio=equil.psio, S_at_surface_left=fm_S_left, ctrl=ctrl, equil=equil, ffit=ffit) - return (dpm=copy(intr.delta_prime_matrix), et1=vac.et[1], msing=intr.msing, - bounds=[(c.psi_start, c.psi_end) for c in fm_chunks]) -end - -@testset "Thread invariance of the parallel BVP path" begin - if Threads.nthreads() < 2 - @info "Thread-invariance tests skipped: effective_threads collapses to 1 in a single-threaded session. Run with `julia -t 4` (CI covers this in its multi-threaded matrix leg)." - @test true - else - @testset "Solovev — leading eigenvalue" begin - dir = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example") - serial = _run_at_thread_cap(dir, 1) - threaded = _run_at_thread_cap(dir, 2) - # Bit-identical, not approximate: the parallel path claims exactness, and a - # tolerance here would hide precisely the reassociation this test exists to catch. - @test threaded.et1 === serial.et1 - end - - @testset "DIII-D-like — Δ′ diagonal and leading eigenvalue" begin - # The deck whose Δ′ the regression harness pins, and where the BVP Δ′ is - # well-conditioned (Solovev sits near marginal stability and its BVP Δ′ is not). - dir = joinpath(@__DIR__, "..", "examples", "DIIID-like_ideal_example") - serial = _run_at_thread_cap(dir, 1) - threaded = _run_at_thread_cap(dir, 2) - - @test threaded.msing == serial.msing - @test size(threaded.dpm) == size(serial.dpm) - @test threaded.et1 === serial.et1 - for j in 1:serial.msing - @test threaded.dpm[j, j] === serial.dpm[j, j] - end - @test threaded.dpm == serial.dpm - - # At these caps the min_bvp_intervals floor fixes the chunk count, so the - # boundaries should be untouched and only scheduling differs. If this fails the - # decomposition moved and the Δ′ comparison above became a stronger claim than - # the one this testset intends to make. - @test threaded.bounds == serial.bounds - end - end -end From 3c5ba1f18357cfe5c3f403c971e5977abf0fc270 Mon Sep 17 00:00:00 2001 From: d-burg Date: Mon, 17 Aug 2026 12:57:18 -0400 Subject: [PATCH 4/5] CI - BUG FIX - Merge the duplicated env block in the test job The multi-threaded matrix leg added JULIA_NUM_THREADS in a second `env:` key alongside the existing DEPOT_PATHS one, so the test job carried a duplicate YAML mapping key. GitHub rejects the whole workflow file on that, which is why no test run reported on this branch and the required Tests context could never be satisfied. Had it parsed, the later block would have won and dropped JULIA_NUM_THREADS entirely, leaving the new leg single-threaded. Co-Authored-By: Claude Opus 5 --- .github/workflows/test.yaml | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/.github/workflows/test.yaml b/.github/workflows/test.yaml index 4c2eee86d..f5b6e9225 100644 --- a/.github/workflows/test.yaml +++ b/.github/workflows/test.yaml @@ -95,8 +95,6 @@ jobs: if: needs.changes.outputs.julia == 'true' runs-on: ${{ matrix.os }} timeout-minutes: 90 - env: - JULIA_NUM_THREADS: ${{ matrix.threads }} strategy: fail-fast: false matrix: @@ -117,6 +115,7 @@ jobs: threads: '4' env: + JULIA_NUM_THREADS: ${{ matrix.threads }} DEPOT_PATHS: | ~/.julia/artifacts ~/.julia/packages From b3499e100e6908f5943da56bc1be097e07f8cafe Mon Sep 17 00:00:00 2001 From: d-burg Date: Wed, 19 Aug 2026 12:21:15 -0400 Subject: [PATCH 5/5] Test - TEST - Drive decomposition invariance through the public solve API The test hand-assembled the stage sequence -- sing_lim!, sing_find!, manual mlow/mhigh arithmetic, make_metric, make_matrix, eulerlagrange_integration, free_run, compute_delta_prime_matrix! -- roughly 30 internal calls reproducing the orchestration main() happened to use. That is a copy of the system under test, and it had already drifted: the hand-rolled mode range gives et[1] = 0.801371 where the production path gives 0.804101, a 0.34% difference, so the test was pinning a configuration GPEC never runs. Now posed as EulerLagrangeProblem + solve(prob, Riccati(; nchunks)), which #393 published, so the test exercises the same path production does. Rewriting it also scoped the claim. Invariance holds at or above the msing-derived chunk target -- auto(53) vs 64 is bit-identical -- but not below it: auto(53) vs 29 moves Delta-prime by ~1e-6 relative. Requesting fewer chunks than the floor is a structurally deficient decomposition rather than merely a different one, since the floor exists to give the surface crossings room. The old test only ever probed the more-chunks direction and so never saw this. The et[1] assertions are gone. et[1] drift across decompositions is a value changing over time, which the regression harness tracks far better than a @test_broken that reports as expected-and-fine forever plus an isapprox 30x looser than the measured effect. It survives here only as a witness that the two runs really were different computations. Co-Authored-By: Claude Opus 5 --- test/runtests_decomposition_invariance.jl | 119 ++++++++++------------ 1 file changed, 52 insertions(+), 67 deletions(-) diff --git a/test/runtests_decomposition_invariance.jl b/test/runtests_decomposition_invariance.jl index d106b99ec..bfcfa6307 100644 --- a/test/runtests_decomposition_invariance.jl +++ b/test/runtests_decomposition_invariance.jl @@ -1,101 +1,86 @@ using Test using TOML -# Decomposition invariance of the Riccati/FM Δ′ path. +# Decomposition invariance of the Riccati Δ′ path. # # The chunked propagator driver reassociates the fundamental-matrix products whenever the chunk # decomposition changes: ((A·B)·C)·D becomes (A·B)·(C·D). Floating-point matrix products do not # reassociate exactly in general, so Δ′ being reproducible requires more than thread-count # independence of the chunk *count* (which is structural: the nchunks=0 target is derived from # msing alone and pinned by unit tests in runtests_parallel_integration.jl). This file asserts -# the end-to-end claim those unit tests cannot: the Δ′ matrix and the leading energy eigenvalue -# are bit-identical when the same integration is cut into a genuinely different set of chunks. +# the end-to-end claim those unit tests cannot: the Δ′ matrix is bit-identical when the same +# integration is cut into a genuinely different set of chunks. # -# Measured basis (DIII-D-like deck): 53-chunk and 64-chunk decompositions with confirmed -# different boundaries gave bit-identical Δ′ diagonals and et[1]. This test pins that property. -# Unlike thread-count variation, the decomposition axis is exercisable in-session at any thread -# count, because nchunks steers it directly. +# Driven through the public solve API rather than the internal stage sequence, so the test +# exercises the same path production does instead of a copy of it that can drift from it. const GP_TI = GeneralizedPerturbedEquilibrium """ -Run the ideal stability pipeline on `dir` with the Riccati integrator at a given chunk count -(`nchunks = 0` = the msing-derived auto target). Mirrors the standalone setup used by the -parallel-integration tests. Returns the Δ′ matrix, the leading energy eigenvalue, and the chunk -boundaries so the test can prove the decompositions actually differed. +Solve the DIII-D-like deck with the Riccati integrator at a given chunk count (`nchunks = 0` is +the msing-derived auto target) and return the published Δ′ matrix alongside the leading energy +eigenvalue, which is used only as a witness that the two decompositions really differed. """ -function _run_at_nchunks(dir::String, nchunks::Int) +function _solve_at_nchunks(dir::String, nchunks::Int) inputs = TOML.parsefile(joinpath(dir, "gpec.toml")) - inputs["ForceFreeStates"]["verbose"] = false - inputs["ForceFreeStates"]["integrator"] = "riccati" - inputs["ForceFreeStates"]["write_outputs_to_HDF5"] = false - inputs["ForceFreeStates"]["nchunks"] = nchunks - - intr = GP_TI.ForceFreeStates.ForceFreeStatesInternal(; dir_path=dir) - ctrl = GP_TI.ForceFreeStates.ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...) + ffs_in = inputs["ForceFreeStates"] eq_config = GP_TI.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], dir) - sol_config = haskey(inputs, "SOL_INPUT") ? GP_TI.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing - equil = GP_TI.Equilibrium.setup_equilibrium(eq_config, sol_config) + equil = GP_TI.Equilibrium.setup_equilibrium(eq_config, nothing) if GP_TI.Equilibrium.wants_two_pass(eq_config) - mand = GP_TI.ForceFreeStates.rational_psi_nodes(equil; nlow=ctrl.nn_low, nhigh=ctrl.nn_high) + mand = GP_TI.ForceFreeStates.rational_psi_nodes(equil; nlow=ffs_in["nn_low"], nhigh=ffs_in["nn_high"]) psi_nodes = GP_TI.Equilibrium.refined_psi_grid(equil; tau=eq_config.psi_accuracy, mandatory=mand) rerun_input = GP_TI.Equilibrium.build_direct_from_ingest(eq_config, equil.ingest) equil = GP_TI.Equilibrium.setup_equilibrium(eq_config, rerun_input; override_psi_nodes=psi_nodes) end - intr.wall_settings = GP_TI.Vacuum.WallShapeSettings(; (Symbol(k) => v for (k, v) in inputs["Wall"])...) - # The toroidal range must be resolved before sing_lim!: under set_psilim_via_dmlim it - # truncates at (last_rational_q + dmlim)/n and so needs n. - intr.nlow = ctrl.nn_low - intr.nhigh = ctrl.nn_high - intr.npert = 1 - GP_TI.ForceFreeStates.sing_lim!(intr, ctrl, equil) - GP_TI.ForceFreeStates.sing_find!(intr, equil) - intr.mlow = min(intr.nlow * equil.params.qmin, 0) - 4 - ctrl.delta_mlow - intr.mhigh = trunc(Int, intr.nhigh * equil.params.qmax) + ctrl.delta_mhigh - intr.mpert = intr.mhigh - intr.mlow + 1 - intr.numpert_total = intr.mpert * intr.npert - metric = GP_TI.ForceFreeStates.make_metric(equil, intr.mpert) - ffit = GP_TI.ForceFreeStates.make_matrix(equil, intr, metric) - odet, fm_propagators, fm_chunks, fm_S_left = GP_TI.ForceFreeStates.eulerlagrange_integration(ctrl, equil, ffit, intr) - vac = GP_TI.ForceFreeStates.free_run(odet, ctrl, equil, ffit, intr) - GP_TI.ForceFreeStates.compute_delta_prime_matrix!(intr, fm_propagators, fm_chunks; - wv=vac.wv, psio=equil.psio, S_at_surface_left=fm_S_left, ctrl=ctrl, equil=equil, ffit=ffit) - return (dpm=copy(intr.delta_prime_matrix), et1=vac.et[1], msing=intr.msing, - bounds=[(c.psi_start, c.psi_end) for c in fm_chunks]) + + # Every ForceFreeStatesControl key from the deck except the ones the problem or the + # integrator owns: nn_low/nn_high come from `nn`, nchunks and the formalism from the alg. + ctrl_kwargs = Dict(Symbol(k) => v for (k, v) in ffs_in + if !(k in ("nn_low", "nn_high", "nchunks", "integrator"))) + ctrl_kwargs[:verbose] = false + ctrl_kwargs[:write_outputs_to_HDF5] = false + + wall = GP_TI.Vacuum.WallShapeSettings(; (Symbol(k) => v for (k, v) in inputs["Wall"])...) + prob = GP_TI.EulerLagrangeProblem(equil; nn=ffs_in["nn_low"], wall=wall, dir_path=dir, ctrl_kwargs...) + return GP_TI.solve(prob, GP_TI.ForceFreeStates.Riccati(; nchunks=nchunks)) end @testset "Decomposition invariance of the Riccati Δ′ path" begin # The deck whose Δ′ the regression harness pins, and where the BVP Δ′ is well-conditioned # (Solovev sits near marginal stability and its BVP Δ′ is pathological there). dir = joinpath(@__DIR__, "..", "examples", "DIIID-like_ideal_example") - auto = _run_at_nchunks(dir, 0) - # 11 more chunks than the auto target: enough to move several boundaries and change the - # association of the propagator products, cheap enough not to distort the runtime. - finer = _run_at_nchunks(dir, length(auto.bounds) + 11) + auto = _solve_at_nchunks(dir, 0) + @test auto.delta_prime !== nothing + msing = size(auto.delta_prime.matrix, 1) - # The premise first: the two decompositions must genuinely differ, otherwise the equality - # below is vacuous and this file silently stops testing anything. - @test length(finer.bounds) > length(auto.bounds) - @test finer.bounds != auto.bounds + # The nchunks=0 target, mirroring balance_integration_chunks' internal formula (the same + # mirroring runtests_parallel_integration.jl does). Invariance is asserted ABOVE this target + # only: requesting fewer chunks than the msing-derived minimum is not merely a different + # decomposition but a structurally deficient one -- measured, auto(53) vs 29 moves Δ′ by + # ~1e-6 relative on this deck, while auto(53) vs 64 is bit-identical. The floor exists to + # give the crossings room, so below it the comparison is not like-for-like. + auto_target = max(2 * msing + 3, 8 * (msing + 1) + msing) + finer = _solve_at_nchunks(dir, auto_target + 11) + @test finer.delta_prime !== nothing - @test finer.msing == auto.msing - @test size(finer.dpm) == size(auto.dpm) + # The premise: the two runs must be genuinely different computations, otherwise the equality + # below is vacuous and this file silently stops testing anything. et[1] is the witness — + # it is decomposition-SENSITIVE on the unified driver (measured 2.7e-8 relative), so its + # differing is evidence the decompositions differed. A failure here means either the chunk + # steering stopped taking effect, or exact et[1] invariance was restored; both want a look + # before this file is trusted again. (The value of et[1] is pinned by the regression + # harness, not here — this is a witness, not an assertion about the physics.) + @test auto.free_boundary !== nothing && finer.free_boundary !== nothing + @test finer.free_boundary.et[1] != auto.free_boundary.et[1] + + @test size(finer.delta_prime.matrix) == size(auto.delta_prime.matrix) # Δ′ is bit-identical, not approximate: a tolerance would hide exactly the reassociation - # drift this test exists to catch, and the measured behaviour is exact equality. - for j in 1:auto.msing - @test finer.dpm[j, j] === auto.dpm[j, j] + # drift this test exists to catch, and the measured behaviour is exact equality. The + # element-wise `===` is deliberately paired with the whole-matrix `==`: `===` holds for + # NaN === NaN, so the `==` is what fails if the computation degrades to NaN. + for j in 1:size(auto.delta_prime.matrix, 1) + @test finer.delta_prime.matrix[j, j] === auto.delta_prime.matrix[j, j] end - @test finer.dpm == auto.dpm - - # et[1] is NOT decomposition-invariant under the unified Riccati driver: measured relative - # difference 2.7e-8 between the auto and auto+11 decompositions on this deck (it was - # bit-identical under the pre-unification driver). The tolerance below is 30× the measured - # effect — documented, not chosen to make the test pass — and exists to catch this - # sensitivity growing by orders of magnitude while the exact-invariance question is open. - # test_broken: if a future driver change restores exact invariance, this reports an - # Unexpected Pass, forcing the === assertion to be reinstated rather than the improvement - # going unnoticed. - @test_broken finer.et1 === auto.et1 - @test isapprox(finer.et1, auto.et1; rtol=1e-6) + @test finer.delta_prime.matrix == auto.delta_prime.matrix end