Test - TEST - Decomposition invariance of the unified Riccati Delta-prime path

.github/workflows/test.yaml

@@ -86,7 +86,11 @@ 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 }}
@@ -99,8 +103,19 @@ jobs:
           - '1.x'  # latest (currently 1.12)
         os:
           - ubuntu-latest
+        threads:
+          - '1'
+        include:
+          # 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'
 
     env:
+      JULIA_NUM_THREADS: ${{ matrix.threads }}
       DEPOT_PATHS: |
         ~/.julia/artifacts
         ~/.julia/packages

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_decomposition_invariance.jl  # Riccati Δ' chunk-decomposition invariance
+# - test/runtests_fullruns.jl           # End-to-end tests
 ```
 
 ### Building Documentation

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_decomposition_invariance.jl")
     include("./runtests_sing.jl")
     include("./runtests_innerlayer.jl")
     include("./runtests_tj_analytic.jl")

test/runtests_decomposition_invariance.jl

@@ -0,0 +1,86 @@
+using Test
+using TOML
+
+# 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 is bit-identical when the same
+# integration is cut into a genuinely different set of chunks.
+#
+# 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
+
+"""
+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 _solve_at_nchunks(dir::String, nchunks::Int)
+    inputs = TOML.parsefile(joinpath(dir, "gpec.toml"))
+    ffs_in = inputs["ForceFreeStates"]
+    eq_config = GP_TI.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], dir)
+    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=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
+
+    # 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 = _solve_at_nchunks(dir, 0)
+    @test auto.delta_prime !== nothing
+    msing = size(auto.delta_prime.matrix, 1)
+
+    # 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
+
+    # 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. 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.delta_prime.matrix == auto.delta_prime.matrix
+end