|
| 1 | +using Test |
| 2 | +using TOML |
| 3 | + |
| 4 | +# Thread-invariance of the parallel FM/BVP path. |
| 5 | +# |
| 6 | +# `ForceFreeStatesControl` documents that the parallel path produces bit-identical Δ′ across |
| 7 | +# thread counts, and the example decks disagree about whether to rely on it: the SLAYER deck |
| 8 | +# pins `parallel_threads = 1` for reproducibility while the DIII-D-like ideal deck — whose Δ′ |
| 9 | +# the regression harness tracks — runs `parallel_threads = 2`. These tests hold the source and |
| 10 | +# the deck fixed and vary only the BVP thread cap, so a golden Δ′ is a property of the physics |
| 11 | +# rather than of the machine it was measured on. |
| 12 | +# |
| 13 | +# Two axes are reachable through `parallel_threads`: |
| 14 | +# |
| 15 | +# - scheduling: `Threads.@threads` over chunks vs a serial loop (any cap ≥ 2) |
| 16 | +# - decomposition: `balance_integration_chunks` targets |
| 17 | +# `max(2·msing+3, 4·effective_threads, 8·(msing+1)+msing)` sub-chunks, so once |
| 18 | +# `4·effective_threads` exceeds the `min_bvp_intervals` floor the chunk boundaries |
| 19 | +# themselves move and the propagator products reassociate |
| 20 | +# |
| 21 | +# The decomposition axis needs `effective_threads > (9·msing+8)/4` — about 14 threads for the |
| 22 | +# DIII-D-like deck's msing=5 — so it is out of reach of a typical CI runner and is exercised by |
| 23 | +# the nightly harness instead. Measured on this deck at parallel_threads = 1, 2, 4 and 16 |
| 24 | +# (53 vs 64 chunks, confirmed different boundaries): Δ′ and et[1] were bit-identical throughout. |
| 25 | +# |
| 26 | +# `effective_threads = min(Threads.nthreads(), parallel_threads)` collapses every cap to 1 in a |
| 27 | +# single-threaded session, which would make these comparisons trivially true, so they are skipped |
| 28 | +# there rather than passing vacuously. |
| 29 | + |
| 30 | +const GP_TI = GeneralizedPerturbedEquilibrium |
| 31 | + |
| 32 | +""" |
| 33 | +Run the ideal stability pipeline on `dir` at a given BVP thread cap. |
| 34 | +
|
| 35 | +Mirrors the standalone setup used by the parallel-integration tests: build the equilibrium |
| 36 | +(applying the two-pass auto grid when the deck asks for it), integrate the Euler-Lagrange |
| 37 | +system, then assemble the STRIDE BVP Δ′ matrix. Returns the Δ′ matrix, the leading energy |
| 38 | +eigenvalue, and the chunk boundaries, so a caller can tell scheduling changes from |
| 39 | +decomposition changes. |
| 40 | +""" |
| 41 | +function _run_at_thread_cap(dir::String, parallel_threads::Int) |
| 42 | + inputs = TOML.parsefile(joinpath(dir, "gpec.toml")) |
| 43 | + inputs["ForceFreeStates"]["verbose"] = false |
| 44 | + inputs["ForceFreeStates"]["use_parallel"] = true |
| 45 | + inputs["ForceFreeStates"]["write_outputs_to_HDF5"] = false |
| 46 | + inputs["ForceFreeStates"]["parallel_threads"] = parallel_threads |
| 47 | + |
| 48 | + intr = GP_TI.ForceFreeStates.ForceFreeStatesInternal(; dir_path=dir) |
| 49 | + ctrl = GP_TI.ForceFreeStates.ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...) |
| 50 | + eq_config = GP_TI.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], dir) |
| 51 | + sol_config = haskey(inputs, "SOL_INPUT") ? GP_TI.Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) : nothing |
| 52 | + equil = GP_TI.Equilibrium.setup_equilibrium(eq_config, sol_config) |
| 53 | + if GP_TI.Equilibrium.wants_two_pass(eq_config) |
| 54 | + mand = GP_TI.ForceFreeStates.rational_psi_nodes(equil; nlow=ctrl.nn_low, nhigh=ctrl.nn_high) |
| 55 | + psi_nodes = GP_TI.Equilibrium.refined_psi_grid(equil; tau=eq_config.psi_accuracy, mandatory=mand) |
| 56 | + rerun_input = GP_TI.Equilibrium.build_direct_from_ingest(eq_config, equil.ingest) |
| 57 | + equil = GP_TI.Equilibrium.setup_equilibrium(eq_config, rerun_input; override_psi_nodes=psi_nodes) |
| 58 | + end |
| 59 | + intr.wall_settings = GP_TI.Vacuum.WallShapeSettings(; (Symbol(k) => v for (k, v) in inputs["Wall"])...) |
| 60 | + GP_TI.ForceFreeStates.sing_lim!(intr, ctrl, equil) |
| 61 | + intr.nlow = ctrl.nn_low |
| 62 | + intr.nhigh = ctrl.nn_high |
| 63 | + intr.npert = 1 |
| 64 | + GP_TI.ForceFreeStates.sing_find!(intr, equil) |
| 65 | + intr.mlow = min(intr.nlow * equil.params.qmin, 0) - 4 - ctrl.delta_mlow |
| 66 | + intr.mhigh = trunc(Int, intr.nhigh * equil.params.qmax) + ctrl.delta_mhigh |
| 67 | + intr.mpert = intr.mhigh - intr.mlow + 1 |
| 68 | + intr.numpert_total = intr.mpert * intr.npert |
| 69 | + metric = GP_TI.ForceFreeStates.make_metric(equil, intr.mpert) |
| 70 | + ffit = GP_TI.ForceFreeStates.make_matrix(equil, intr, metric) |
| 71 | + odet, fm_propagators, fm_chunks, fm_S_left = GP_TI.ForceFreeStates.eulerlagrange_integration(ctrl, equil, ffit, intr) |
| 72 | + vac = GP_TI.ForceFreeStates.free_run!(odet, ctrl, equil, ffit, intr) |
| 73 | + GP_TI.ForceFreeStates.compute_delta_prime_matrix!(intr, fm_propagators, fm_chunks; |
| 74 | + wv=vac.wv, psio=equil.psio, S_at_surface_left=fm_S_left, ctrl=ctrl, equil=equil, ffit=ffit) |
| 75 | + return (dpm=copy(intr.delta_prime_matrix), et1=vac.et[1], msing=intr.msing, |
| 76 | + bounds=[(c.psi_start, c.psi_end) for c in fm_chunks]) |
| 77 | +end |
| 78 | + |
| 79 | +@testset "Thread invariance of the parallel BVP path" begin |
| 80 | + if Threads.nthreads() < 2 |
| 81 | + @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)." |
| 82 | + @test true |
| 83 | + else |
| 84 | + @testset "Solovev — leading eigenvalue" begin |
| 85 | + dir = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example") |
| 86 | + serial = _run_at_thread_cap(dir, 1) |
| 87 | + threaded = _run_at_thread_cap(dir, 2) |
| 88 | + # Bit-identical, not approximate: the parallel path claims exactness, and a |
| 89 | + # tolerance here would hide precisely the reassociation this test exists to catch. |
| 90 | + @test threaded.et1 === serial.et1 |
| 91 | + end |
| 92 | + |
| 93 | + @testset "DIII-D-like — Δ′ diagonal and leading eigenvalue" begin |
| 94 | + # The deck whose Δ′ the regression harness pins, and where the BVP Δ′ is |
| 95 | + # well-conditioned (Solovev sits near marginal stability and its BVP Δ′ is not). |
| 96 | + dir = joinpath(@__DIR__, "..", "examples", "DIIID-like_ideal_example") |
| 97 | + serial = _run_at_thread_cap(dir, 1) |
| 98 | + threaded = _run_at_thread_cap(dir, 2) |
| 99 | + |
| 100 | + @test threaded.msing == serial.msing |
| 101 | + @test size(threaded.dpm) == size(serial.dpm) |
| 102 | + @test threaded.et1 === serial.et1 |
| 103 | + for j in 1:serial.msing |
| 104 | + @test threaded.dpm[j, j] === serial.dpm[j, j] |
| 105 | + end |
| 106 | + @test threaded.dpm == serial.dpm |
| 107 | + |
| 108 | + # At these caps the min_bvp_intervals floor fixes the chunk count, so the |
| 109 | + # boundaries should be untouched and only scheduling differs. If this fails the |
| 110 | + # decomposition moved and the Δ′ comparison above became a stronger claim than |
| 111 | + # the one this testset intends to make. |
| 112 | + @test threaded.bounds == serial.bounds |
| 113 | + end |
| 114 | + end |
| 115 | +end |
0 commit comments