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

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 <noreply@anthropic.com>

Author: d-burg

.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

test/runtests.jl

@@ -30,6 +30,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_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