TESTING - NEW FEATURE - Test thread invariance of the parallel BVP path

.github/workflows/test.yaml

@@ -87,11 +87,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:
@@ -100,6 +106,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'
 
     steps:
       - name: Checkout repository

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

src/ForceFreeStates/EulerLagrange.jl

@@ -67,7 +67,13 @@ end
 Sub-divide integration chunks to produce a load-balanced set for parallel execution.
 Starts from the output of `chunk_el_integration_bounds` and iteratively splits the
 highest-cost chunk (by `ode_itime_cost`) until the total chunk count reaches
-`max(2*msing + 3, 4 * Threads.nthreads())`.
+`max(2*msing + 3, 4*effective_threads, 8*(msing + 1) + msing)`, where
+`effective_threads = min(Threads.nthreads(), ctrl.parallel_threads)`.
+
+The last term (BVP propagator conditioning) dominates for every realistic thread count, so
+the decomposition — and hence the floating-point association of the propagator products — is
+normally independent of the thread count. It stops dominating only when
+`4*effective_threads` exceeds it, i.e. beyond roughly `(9*msing + 8)/4` threads.
 
 Each split finds the equal-cost midpoint ψ_mid via bisection:
   ode_itime_cost(psi_start, psi_mid) ≈ ode_itime_cost(psi_start, psi_end) / 2

test/runtests.jl

@@ -29,6 +29,7 @@ else
     include("./runtests_eulerlagrange.jl")
     include("./runtests_riccati.jl")
     include("./runtests_parallel_integration.jl")
+    include("./runtests_thread_invariance.jl")
     include("./runtests_sing.jl")
     include("./runtests_innerlayer.jl")
     include("./runtests_tj_analytic.jl")

test/runtests_parallel_integration.jl

@@ -135,8 +135,11 @@ using TOML
         balanced = GeneralizedPerturbedEquilibrium.ForceFreeStates.balance_integration_chunks(base_chunks, ctrl, intr)
 
         # Must mirror balance_integration_chunks' internal target_n formula
-        # (src/ForceFreeStates/EulerLagrange.jl). Keep this in sync.
-        target_n = max(2 * intr.msing + 3, 4 * Threads.nthreads(), 8 * (intr.msing + 1) + intr.msing)
+        # (src/ForceFreeStates/EulerLagrange.jl). Keep this in sync. The parallel width is
+        # capped by ctrl.parallel_threads, not the raw thread count: using Threads.nthreads()
+        # here overestimates the target on a machine with more threads than the cap.
+        effective_threads = min(Threads.nthreads(), max(ctrl.parallel_threads, 1))
+        target_n = max(2 * intr.msing + 3, 4 * effective_threads, 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,

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