ForceFreeStates: Dual Riccati + Parallel FM integration + Δ' output

[logan-nc]

Summary

Implements two new integration paths for the Euler-Lagrange ODE and adds full Δ' (tearing stability parameter) computation to the Riccati and parallel FM paths.

• Dual Riccati integration (use_riccati = true): Sequential reformulation as a matrix Riccati ODE S = U₁·U₂⁻¹ with periodic renormalization to maintain bounded (U₁, U₂). Validated on Solovev (0.006% energy error) and DIIID n=1 (0.12% error). Faster than standard for large N.

• Parallel FM integration (use_parallel = true): Each integration chunk is solved independently from identity ICs using Threads.@threads, then assembled serially. Uses bidirectional integration — crossing chunks (near rational surfaces) are integrated backward to keep FM propagators well-conditioned. The well-conditioned backward FM inverse is applied via LU solve during serial assembly. Fixes a ~10% energy error for DIIID N=26 that existed with the all-forward approach.

• Per-surface Δ': Both paths compute intr.sing[s].delta_prime (scalar) and delta_prime_col (full N-vector of off-diagonal coupling). Written to singular/delta_prime and singular/delta_prime_col in HDF5.

• Inter-surface Δ' matrix (parallel FM only): intr.delta_prime_matrix (2·msing × 2·msing) via the STRIDE global BVP [Glasser 2018 Phys. Plasmas 25, 032501]. Written to singular/delta_prime_matrix in HDF5.

• Edge-dW scan decoupled from integration truncation (new, commit 5d5b8ee): The edge-dW scan over ψ ∈ [psiedge, psilim] is now diagnostic-only by default. It no longer silently moves psilim/qlim/u to the dW peak, which was corrupting Δ' of the outermost rational by tens of percent and making the answer depend on a subtle, hidden heuristic. Legacy behaviour is preserved behind an opt-in truncate_at_dW_peak = true flag (marked in the docstring as "Δ'/δW unreliable in this mode"). psihigh / qhigh / dmlim are now the only controls over the integration domain. Validated against Fortran STRIDE on 42+56-point TJ β and ε scans: Δ'(2/1), Δ'(3/1), δWp, δWv, δWt all match within ~1% (see LAR_beta_scan and LAR_epsilon_scan comparison plots on the companion perf/riccati-full-geqdsk-scans branch).

---

Accuracy

Example	Path	et[1]	Error vs std
Solovev N=8	standard	16.56600	—
Solovev N=8	riccati	16.56691	0.006%
Solovev N=8	parallel	16.56691	0.006%
DIIID N=26	standard	1.21237	—
DIIID N=26	riccati	1.21386	0.123%
DIIID N=26	parallel	1.21386	0.123%

The ~0.12% gap on DIIID is algorithmic (different crossing convention: Riccati-style vs Gaussian Reduction), not ODE tolerance. Both paths are physically correct; 0.12% is well within physical uncertainty.

---

Performance (4 threads)

Example	Path	Time	Speedup
Solovev N=8	standard	0.25s	1.00×
Solovev N=8	riccati	0.09s	2.83×
Solovev N=8	parallel	0.14s	1.80×
DIIID N=26	standard	1.17s	1.00×
DIIID N=26	riccati	0.91s	1.29×
DIIID N=26	parallel	0.90s	1.30×

---

Tests

Suite	Pass	Total
runtests_riccati.jl	23	23
runtests_parallel_integration.jl	56	56
runtests_fullruns.jl	9	9

Riccati tests include unit tests for: renormalize_riccati_inplace!, renormalize_riccati!, Riccati end state U₂≈I, Δ' regression values, delta_prime_col shape/diagonal consistency, riccati_der! formula (Glasser 2018 Eq. 19), and compute_delta_prime_from_ca! bit-identical consistency.

Parallel tests include: apply_propagator! identity and linearity, apply_propagator_inverse! round-trip, chunk balance target count, direction field propagation through bidirectional=true, Solovev energy accuracy, DIIID energy accuracy (key bidirectional fix regression), ode_itime_cost additivity, and Δ' matrix regression for both Solovev and DIIID.

runtests_fullruns.jl updated: regression_solovev_kinetic_multi_n/et[1] expected from −0.01248 → −0.19359 with inline comment. The old value reflected the truncated-integration behaviour from the edge-dW heuristic; the new value reflects the full-domain answer controlled solely by psihigh/qhigh.

---

Breaking behaviour change

Edge-dW scan is no longer a silent psilim truncation. Users who relied on psiedge < psilim to shrink the integration domain must now either:

1. Set psihigh or qhigh (or dmlim) explicitly to the desired outer boundary (recommended — the whole point of this change is that the edge is user-controlled and visible), or
2. Add truncate_at_dW_peak = true to their [ForceFreeStates] TOML section to restore the legacy heuristic (flagged as unreliable for Δ' and δW; preserved for experimental work on improved edge-mode filters).

The edge_scan/ HDF5 group is still written with psi, q, plasma_energy, vacuum_energy, total_energy, vacuum_eigenvalue for whatever post-processing previously used it as a diagnostic.

---

Regression harness snapshot (develop @ e0ee9df vs perf/riccati @ 685a92a)

Run via regress --cases diiid_n1,solovev_n1,solovev_multi_n,tj_epsilon_pole --refs develop,local --force.

• solovev_n1: 16 CHANGED, 5 unchanged
  • et[1]: −0.2646 → +0.0999 (sign flip — driven by the edge-dW decoupling plus pre-existing Riccati/parallel-FM work on this branch)
  • vacuum matrix min eigenvalue: 1.905 → 2.284 (new HDF5 quantity)
  • ODE steps (saved): 607 → 273 (55 % fewer due to the psihigh refactor on develop between 04-09 and 04-12, unrelated to this PR's edits)
• diiid_n1: 21 CHANGED, 6 unchanged
  • Leading-eigenvalue shifts track the energy-principle refactor
  • Large shifts in delta prime, island half-widths, Chirikov parameter, ||resonant flux|| are pre-existing on develop (||resonant flux|| went 452 → 3.26 × 10⁶ between develop @ 04-09 and develop @ 04-12 already — see ticket below) and not caused by this PR. The in-PR perturbed-equilibrium changes are expected to re-baseline these; a separate investigation should confirm the DIIID singular-coupling numbers are physically sane.
• solovev_multi_n: 11 CHANGED, 4 unchanged — qualitatively the same pattern as solovev_n1
• tj_epsilon_pole: new case on this branch, not present on develop

---

Files changed

File	Change
src/ForceFreeStates/ForceFreeStatesStructs.jl	IntegrationChunk.direction, SingType.delta_prime_col, ForceFreeStatesInternal.delta_prime_matrix; new truncate_at_dW_peak::Bool = false control with docstring warning
src/ForceFreeStates/EulerLagrange.jl	chunk_el_integration_bounds(bidirectional=false); balance_integration_chunks direction propagation; edge-dW scan split into diagnostic (default) vs legacy-truncation paths
src/ForceFreeStates/Riccati.jl	Full Riccati path; integrate_propagator_chunk! bidirectional tspan; apply_propagator_inverse!; bidirectional serial assembly; compute_delta_prime_matrix! with Phi_L/Phi_R BVP; assemble_fm_matrix empty-range guard; single-resonance assertion; edge-dW decoupling at both the sequential Riccati and parallel FM exits
src/GeneralizedPerturbedEquilibrium.jl	HDF5 output for singular/delta_prime, delta_prime_col, delta_prime_matrix; edge_scan/ diagnostic group preserved
test/runtests_riccati.jl	23 tests, shared setup, riccati_der! formula test, compute_delta_prime_from_ca! test
test/runtests_parallel_integration.jl	56 tests including apply_propagator_inverse!, direction field, DIIID Δ' matrix
test/runtests_fullruns.jl	Solovev kinetic multi-n et[1] expected value updated to reflect full-domain integration
examples/LAR_{beta,epsilon}_scan/gpec.toml	Drop stale kin_flag/con_flag keys (removed from the control struct during the develop merge)
examples/TJ_epsilon_pole_example/gpec.toml	Same stale-flag cleanup
benchmarks/benchmark_threads.jl	Thread-scaling benchmark
benchmarks/benchmark_riccati_der.jl	riccati_der! sanity check
benchmarks/benchmark_delta_prime_methods.jl	compute_delta_prime_from_ca! sanity check

---

Code review responses

Fixed: assemble_fm_matrix empty-range guard · single-resonance assertion in compute_delta_prime_matrix! · stale comment in parallel tests · dead parallel_threads field removed · stale kin_flag/con_flag in LAR/TJ example TOMLs (would crash at runtime) · edge-dW silent psilim truncation (was corrupting Δ' by >50 % at some ε-scan points)

Deferred: apply_propagator! per-call allocations · test boilerplate in runtests_parallel_integration.jl · DIIID perturbed-equilibrium coupling-metric magnitudes (pre-existing on develop, separate investigation)

---

🤖 Generated with Claude Code

[logan-nc]

@claude review this

[claude]

Claude finished @logan-nc's task in 3m 55s —— View job

---

Code Review

This is a well-designed PR implementing three complementary integration paths. The theory is clearly traced to Glasser (2016/2018), the docstrings are excellent, and the bidirectional fix for large-N conditioning is elegant. Below are issues found, ordered by severity.

---

Bug: assemble_fm_matrix panics on empty range

src/ForceFreeStates/Riccati.jl:77-87

function assemble_fm_matrix(propagators::Vector{ChunkPropagator}, idx_range)
    N = size(propagators[1].block_upper_ic, 1)   # ← accesses propagators[1] unconditionally
    Phi = Matrix{ComplexF64}(I, 2N, 2N)
    for i in idx_range
        Phi = Phi_i * Phi
    end
    return Phi
end

Called in compute_delta_prime_matrix! as:

Phi_R_mats[1] = assemble_fm_matrix(propagators, 1:i_crossings[1]-1)

If the first chunk is already a crossing chunk (i_crossings[1] == 1), the range is 1:0 (empty). propagators[1] still exists but the indexing to determine N is a fragile pattern — and more importantly, the result should be the identity matrix, which requires knowing N from somewhere else when the range is empty. The function should guard against this:

function assemble_fm_matrix(propagators::Vector{ChunkPropagator}, idx_range)
    isempty(idx_range) && return Matrix{ComplexF64}(I, 2*size(propagators[1].block_upper_ic,1), 2*size(propagators[1].block_upper_ic,1))
    N = size(propagators[first(idx_range)].block_upper_ic, 1)
    ...
end

Or better: pass N as an explicit argument. Fix this →

---

Minor: apply_propagator! allocates on every call in the serial assembly loop

src/ForceFreeStates/Riccati.jl:739-758

function apply_propagator!(odet::OdeState, prop::ChunkPropagator)
    u1_prev = copy(@view odet.u[:, :, 1])   # allocation
    u2_prev = copy(@view odet.u[:, :, 2])   # allocation
    tmp = similar(u1_prev)                   # allocation
    ...
end

This is called once per chunk in the serial assembly phase, so for a typical DIIID run (say 50 chunks) the overhead is ~150 allocations of N×N complex matrices. For N=26 this is negligible, but for the design intent of supporting large N it is worth pre-allocating a work buffer in parallel_eulerlagrange_integration and threading it through. Same for apply_propagator_inverse!. Fix this →

---

Minor: Stale comment in runtests_parallel_integration.jl

test/runtests_parallel_integration.jl:125-133

# Note: this test uses the Solovev example (N=8 modes) where FM propagators
# are well-conditioned. For large-N problems (N ≳ 20, e.g. DIIID with N=26),
# FM propagator ill-conditioning leads to ~10% energy error with no speedup
# over the serial Riccati path. See parallel_eulerlagrange_integration docstring
# for details and deferred fix approaches (bidirectional integration / continuous QR).

This comment predates the bidirectional fix that is the main contribution of this PR. The large-N issue is now resolved — the DIIID accuracy is within 0.24% (see test at line 190). This comment should be updated to reflect the current state. Fix this →

---

Minor: Dead code — riccati_der! and compute_delta_prime_from_ca!

Both functions have docstrings acknowledging they are unused ("retained for reference"):

• riccati_der! (Riccati.jl:281): uses the @with_pool pool macro, tying it to an execution context. If this function is never called, it should either be tested (e.g., as a sanity check that the explicit form matches sing_der!), or removed. The docstring says it could be used with an implicit solver — that's a legitimate reason to keep it, but at minimum add a comment linking to the issue/TODO tracking that future work.

• compute_delta_prime_from_ca! (EulerLagrange.jl:22): the docstring explicitly says "no longer called from any integration driver". If it's retained for testing, add a test. Otherwise remove.

---

Minor: compute_delta_prime_matrix! single-resonance assumption not asserted early

Riccati.jl:167-172

ipert_all = [begin
    sp = intr.sing[j]
    idx = 1 + sp.m[1] - intr.mlow + (sp.n[1] - intr.nlow) * intr.mpert
    @assert 1 <= idx <= N "Resonant mode index out of range"
    idx
end for j in 1:msing]

The function uses sp.m[1] / sp.n[1], silently assuming exactly one resonant mode per surface. The limitation is documented in the docstring, but multi-resonance surfaces will silently give wrong results (only the first resonant mode is coupled). A stronger guard at the function entry would be safer:

@assert all(j -> length(intr.sing[j].m) == 1, 1:msing) "compute_delta_prime_matrix! only supports single-resonance surfaces; got multi-resonance surface"

Fix this →

---

Minor: Test boilerplate duplication

The 15-line equilibrium setup block (parse TOML → setup_equilibrium → sing_find! → set modes) is repeated verbatim across 7+ test cases in runtests_riccati.jl and runtests_parallel_integration.jl. Extract a helper function:

function setup_solovev_intr_ctrl(; use_riccati=false, use_parallel=false)
    ex = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example")
    inputs = TOML.parsefile(joinpath(ex, "jpec.toml"))
    inputs["ForceFreeStates"]["verbose"] = false
    inputs["ForceFreeStates"]["use_riccati"] = use_riccati
    inputs["ForceFreeStates"]["use_parallel"] = use_parallel
    ...
    return intr, ctrl, equil, ffit
end

This would cut the test file length by ~40% and make failures easier to localize. Fix this →

---

Minor: parallel_threads field now superseded

ForceFreeStatesControl.parallel_threads (already present, default 1) was previously documented as "not yet implemented". Now that use_parallel routes through Threads.nthreads(), parallel_threads is effectively dead. Either: (a) wire parallel_threads into parallel_eulerlagrange_integration via @threads :static with thread count control, or (b) document it as deprecated in favour of JULIA_NUM_THREADS env variable.

---

Positive observations

• The bidirectional integration rationale is crystal clear. The docstring on parallel_eulerlagrange_integration (and the STRIDE BVP comments in compute_delta_prime_matrix!) trace every design decision back to the specific equations in Glasser (2018). This level of traceability is excellent.

• The outer plasma re-integration (Riccati rather than FM for the post-last-crossing region) is a thoughtful decision: it avoids FM ill-conditioning in the outer plasma and provides dense checkpoints for findmax_dW_edge!. The comment explaining why is clear.

• renormalize_riccati_inplace! after every chunk in the serial assembly loop is the Julia analogue of ode_fixup in STRIDE — this connection is explicitly called out in the code comment, which is appreciated.

• integrate_propagator_chunk! with no callback is correctly justified: identity ICs keep the solution bounded within each chunk (the chunk is designed to stop before the rational surface), so normalization is unnecessary.

• Δ' normalization convention is carefully documented in riccati_cross_ideal_singular_surf!: the skip-renorm-before-crossing approach ensuring ca_r[ipert_res, ipert_res, 2] = 1 is a subtle but correct choice that makes the left/right Δ' comparison gauge-consistent.

• Test coverage is solid. Unit tests for renormalize_riccati_inplace! (idempotent case, random case), apply_propagator! (identity, linearity), balance_integration_chunks (chunk count, contiguity, crossing count), and end-to-end regression tests for both Solovev and DIIID.

---

[logan-nc]

Post-merge update — after develop (GPEC rename) merge

Merged origin/develop (GPEC name-change PR #168) into this branch. One conflict resolved in write_outputs_to_HDF5: kept the new Δ' HDF5 outputs from this branch, adopted develop's updated vacuum field names (vac_data, plasma_pts/wall_pts, 3D Cartesian). Updated all test and benchmark files to reference GeneralizedPerturbedEquilibrium and gpec.toml.

Test results (post-merge)

Suite	Pass	Total	Time
runtests_riccati.jl	23	23	1m27s
runtests_parallel_integration.jl	56	56	1m31s

Benchmark results (post-merge, 1 thread)

Example	Path	et[1]	Time	Speedup vs std	Error vs std
Solovev N=8	standard	16.56600	0.09s	1.00×	—
Solovev N=8	riccati	16.56691	0.17s	0.51×	0.0055%
Solovev N=8	parallel	16.56691	0.13s	0.66×	0.0055%
DIIID N=26	standard	1.21237	1.26s	1.00×	—
DIIID N=26	riccati	1.21386	0.63s	1.98×	0.1230%
DIIID N=26	parallel	1.21386	1.19s	1.06×	0.1231%

Benchmark results (post-merge, 4 threads)

Example	Path	et[1]	Time	Speedup vs std	Error vs std
Solovev N=8	standard	16.56600	0.25s	1.00×	—
Solovev N=8	riccati	16.56691	0.09s	2.83×	0.0055%
Solovev N=8	parallel	16.56691	0.14s	1.80×	0.0055%
DIIID N=26	standard	1.21237	1.17s	1.00×	—
DIIID N=26	riccati	1.21386	0.91s	1.29×	0.1230%
DIIID N=26	parallel	1.21386	0.90s	1.30×	0.1231%

Note: 1-thread Solovev timings have high variance (sub-100ms runs) and are not meaningful for speedup comparisons. DIIID (1.2s+) is the reliable benchmark. At 4 threads, Riccati achieves 1.29× and parallel achieves 1.30× on DIIID. Energy agreement is unchanged: <0.13% for both paths vs standard.

[logan-nc]

Pre-merge cleanup (commits 1515591..725a527)

Three rounds of fixes applied after a code review pass:

CI fixes (1515591)

• Random stdlib missing from Project.toml — runtests_riccati.jl uses using Random but the package was not declared in [deps], causing CI to fail with Package Random not found. Added to both [deps] and [compat].
• Broken markdown links in docstrings — adding @autodocs in stability.md surfaced docstring patterns like x_left[j] (2N) that CommonMark parses as [link text](url) even with spaces between ] and (. Fixed by wrapping affected notation blocks in code fences (assemble_fm_matrix, compute_delta_prime_matrix!, VacuumData field docs).

Documentation (142a79c)

• New docs/src/stability.md page covering the ForceFreeStates module: Newcomb/DCON stability criterion with paper citations (Glasser 2016/2018a/2018b), all three integration drivers, Δ' computation (per-surface and STRIDE BVP matrix), configuration reference, @autodocs API block, and example usage.
• Added to docs/make.jl navigation; cross-linked from equilibrium.md.

Pre-merge code review fixes (725a527)

Three issues flagged as must-fix or high-priority:

1. Threads.@threads :static — since Julia 1.7, the default :dynamic scheduler can migrate tasks between OS threads, making Threads.threadid() unreliable for indexing odet_proxies. :static guarantees a 1:1 task-to-thread mapping.

2. outer_chunk psi_end guard — the outer-plasma re-integration now uses psilim * (1 - eps) to match the guard applied by chunk_el_integration_bounds, avoiding a potential boundary evaluation at exactly psilim.

3. Consistent logging — all println calls in Riccati.jl replaced with @info macros, matching the pattern in EulerLagrange.jl and enabling log-level filtering/suppression in tests.

All 23 Riccati tests and 56 parallel FM tests pass locally after these changes.

[logan-nc]

@jhalpern30 will clean this up by moving U solutions to be 2D instead of 3D after this merges (#89)

[d-burg]

Final checks from 04/22 meeting:

• Check how truncate_at_dW_peak is changing Δ', specifically order of finishing integration vs Δ' calculation
• Double δW use_parallel=False
• Add all TJ equil parameters to gpec.toml
• Check energy U_store and Ud_store in use_parallel=True, these are DCON ξ functions. Essential for perturbed equilibrium. Compare against DCON and add to benchmark

[matt-pharr]

@d-burg status on this?

[d-burg]

perf/riccati: parallel ξ + Δ' from one run, TJ-analytic rename, multi-resonance skip

@logan-nc @matt-pharr what do you think? I believe we're just about ready to merge, pending your recommendations.

Summary

Brings the parallel FM-propagator BVP path to feature-parity with the
serial DCON path so a single use_parallel = true run produces both
the inter-surface Δ' matrix and dense DCON ξ functions in axis basis.
Renames GPEC's analytic-equilibrium adaptation of R. Fitzpatrick's TJ code
to tj_analytic (with citation everywhere it's defined or used), turns the
LAR scan examples into self-describing single-file gpec.toml, and
gracefully handles multi-n runs where one q value is rational for
multiple (m, n) tuples.

Commits since c6c845f

Hash	Title
54d12fe2	ForceFreeStates - NEW FEATURE - Port set_psilim_via_dmlim + default tightenings
4845ec80	ForceFreeStates - IMPROVEMENT - Default use_parallel=true, singfac_min=1e-4
db7c490a	ForceFreeStates - BUG FIX - Wire ctrl.parallel_threads into BVP path; default 1 (serial) eliminates DIII-D 147131 thread-race
7ac87c8e	ForceFreeStates - PERFORMANCE - parallel_threads default 1 → 2 (≈20 % BVP speedup; bit-identical Δ' in 15-trial DIII-D 147131 sweep)
3c8130da	ForceFreeStates - BUG FIX + EXAMPLES - truncate_at_dW_peak self-consistent Δ' + LAR TJ TOML refactor + parallel ξ benchmark
5acf1478	ForceFreeStates - NEW FEATURE - Dense ξ in parallel BVP path + bit-identical regression + pinned Δ'
6d07c07d	EQUIL - REFACTOR - Rename TJ → TJ-like with Fitzpatrick citation everywhere
085de133	EXAMPLES - CLEANUP - LAR scans: single-file gpec.toml, per-line annotation, TJ → TJ-like
5a4c2c29	EXAMPLES - CLEANUP - Remove TJ_epsilon_pole_example and its regression case
cc2affeb	EQUIL - REFACTOR - Rename tj_like → tj_analytic (cleaner, less hedge-y)
8073c126	ForceFreeStates - BUG FIX - Preserve Riccati-gauge ca_l/ca_r across dense ξ pass

Total: 20 files changed (+1355, −718) over the 5 new commits, on top of
the 5 already-pushed perf/riccati commits.

"What's in the box" (a proper Claude-ism)

ForceFreeStates — dense ξ in the parallel path (5acf1478)

Three coupled changes so the parallel FM-propagator BVP delivers both Δ'
and DCON eigenfunctions usable by downstream PerturbedEquilibrium:

1. Dense ξ pass. parallel_eulerlagrange_integration now appends a
   serial Euler-Lagrange dense pass (helper
   _populate_dense_xi_via_serial_el!) that replaces the propagator-BVP
   odet with a fresh serial-EL odet whose u_store / ud_store are
   dense and in axis basis — the only convention PerturbedEquilibrium
   / FieldReconstruction consumes correctly. All BVP-relevant fields
   (intr.psilim, intr.qlim, intr.sing[*].delta_prime,
   delta_prime_col, ua_left, psi_ua_left) are saved/restored
   across the pass. Gated by the new
   ctrl.populate_dense_xi::Bool = true (default on).

2. Multi-resonance skip in compute_delta_prime_matrix!. Replace
   the hard @assert (which crashed multi-n runs whose q value was
   rational for two distinct (m, n) tuples) with an early return +
   warning. Per-surface scalar Δ' from
   riccati_cross_ideal_singular_surf! and HDF5
   singular/delta_prime / singular/delta_prime_col remain
   populated; only the inter-surface BVP singular/delta_prime_matrix
   is omitted in that regime.

3. Tests + benchmark.

   • New @testset "ξ functions bit-identical between use_parallel modes (populate_dense_xi)" proves psi_store / q_store / u_store / ud_store / crit_store / step / nzero from
     use_parallel=true; populate_dense_xi=true are byte-for-byte
     identical to a use_parallel=false run on both Solovev (small N)
     and DIIID-like (large N). A populate_dense_xi=false control
     assertion proves the result isn't trivially passing through
     accidental sparse storage on both sides.
   • Pinned per-surface intr.sing[s].delta_prime values added to both
     "Parallel FM integration matches standard ODE" testsets
     (rtol = 5 %, matching the existing et_par ≈ 1.29 style).
   • Pinned diagonal delta_prime_matrix values added to both STRIDE
     BVP testsets (rtol = 5 %).
   • benchmarks/benchmark_xi_parallel_vs_serial.jl rewritten:
     accepts any example dir (defaults to Solovev + DIIID-like),
     overlays all resonant modes on log-y, adds a right-column residual
     panel. Confirms exact (max|Δ| == 0) match on both equilibria.

   Net: runtests_parallel_integration.jl grew 113 → 127 tests
   (≈13 s extra per CI matrix entry); runtests_fullruns.jl went from
   8/9 (pre-existing multi-n crash) to 9/9 pass.

EQUIL — TJ → TJ-analytic rename + Fitzpatrick citation (6d07c07d, cc2affeb)

GPEC's adaptation of R. Fitzpatrick's TJ code analytic-profile family
(https://github.com/rfitzp/TJ) is now consistently named "TJ-analytic" in
identifiers and prose. The upstream TJ code is cited at every
definition and use site. Why two commits: the first pass introduced
tj_like as the new identifier, but on review it read as a weak hedge
("kinda like the real one"); tj_analytic says exactly what this is —
the analytic-profile model from Fitzpatrick's TJ, fed through GPEC's
pipeline rather than the full TJ code.

Identifier renames:

• TJConfig → TJAnalyticConfig (both file-path and dict constructors)
• tj_run / tj_run_direct → tj_analytic_run / tj_analytic_run_direct
• tj_f1 / tj_f1p → tj_analytic_f1 / tj_analytic_f1p
• tj_shape_rhs! / _initial / _solve → tj_analytic_shape_*
• tj_find_nu → tj_analytic_find_nu
• TJShapeParams → TJAnalyticShapeParams
• Local parameter tj::TJAnalyticConfig (same short tj as before;
  the type signature carries the disambiguation).

Config / user-facing renames:

• eq_type values: "tj" → "tj_analytic", "tj_direct" → "tj_analytic_direct"
• Embedded TOML section: [TJ_INPUT] → [TJ_ANALYTIC_INPUT]
• EquilibriumConfig now makes eq_filename optional when an embedded
  [TJ_ANALYTIC_INPUT] / [SOL_INPUT] / [LAR_INPUT] section is present.
• Dropped a stale sigma_type="tj" reference on
  LargeAspectRatioConfig.qa.

Tests: test/runtests_tj_analytic.jl (16/16 pass). Coincidental
matches in Vacuum/Field.jl (fintjj) and
InnerLayer/GGJ/{Shooting,InnerAsymptotics}.jl (_build_tjmat,
inps_tjmat, loop-local tj) are intentionally left alone — they
have nothing to do with Fitzpatrick's TJ code.

Prose convention: the -analytic suffix is reserved for our model
class. Bare "TJ" refers to the upstream code (e.g. "TJ-analytic
GetPSIvac (cf. Fitzpatrick's TJ)" in comments referring to functions
we ported the logic of).

EXAMPLES — LAR scans (085de133) + TJ_epsilon_pole removal (5a4c2c29)

LAR_beta_scan and LAR_epsilon_scan each now consist of a single
self-describing gpec.toml plus run_scan.jl. No more tj.toml
side-cars: all TJ-analytic equilibrium parameters live in an embedded
[TJ_ANALYTIC_INPUT] section. Every field across [Equilibrium],
[TJ_ANALYTIC_INPUT], [Wall], and [ForceFreeStates] has a
one-liner comment describing what it is (not just a label) — e.g.
"Number of radial spline nodes used to discretize ψ" rather than
"Radial grid points". Each header cites
https://github.com/rfitzp/TJ.

run_scan.jl scripts updated:

• import TJAnalyticConfig (was TJConfig)
• override config["TJ_ANALYTIC_INPUT"][...] (was config["TJ_INPUT"][...])
• LAR_epsilon_scan flips eq_type → "tj_analytic_direct" per point
• banners say "TJ-analytic β scan" / "TJ-analytic ε scan"

diagnose_profiles.jl docstring clarified that its "TJ" geqdsk
comparison data are produced by Fitzpatrick's external TJ code, not
GPEC's internal tj_analytic model.

examples/TJ_epsilon_pole_example/ and the matching
regression-harness/cases/tj_epsilon_pole.toml removed. The
ε ≈ 0.66 near-pole physics it exercised is already covered by the
ε-scan in examples/LAR_epsilon_scan/ (which sweeps ε up to 0.660
along the same kink-pole approach) and by the
"tj_analytic_run_direct (Option B) — pole-approach physics at ε = 0.60"
testset in test/runtests_tj_analytic.jl.

Default behaviour after merge

With the shipped example tomls (ode_flag=true, vac_flag=true,
use_parallel=true, populate_dense_xi=true), one main() call
produces:

Quantity	HDF5 path	Required flags
DCON ξ functions, dense in axis basis	integration/xi_psi, dxi_psi, xi_s, psi, q	ode_flag, populate_dense_xi
Per-surface Δ'	singular/delta_prime, singular/delta_prime_col	ode_flag, use_parallel
Inter-surface PEST3 Δ' matrix	singular/delta_prime_matrix	ode_flag, use_parallel, vac_flag, single-resonance surfaces only
Free-boundary energies	vacuum/ep, vacuum/ev, vacuum/et	ode_flag, vac_flag

Tests

File	Result
runtests_utilities.jl	11 / 11 ✓
runtests_vacuum.jl	252 / 252 ✓
runtests_equil.jl	240 / 240 ✓
runtests_eulerlagrange.jl	92 / 92 ✓
runtests_riccati.jl	24 / 24 ✓
runtests_parallel_integration.jl	127 / 127 ✓ (was 113 before this PR; +14 from pinned Δ' + bit-identical tests)
runtests_sing.jl	75 / 75 ✓
runtests_tj_analytic.jl	16 / 16 ✓
runtests_fullruns.jl	9 / 9 ✓ (was 8/9 errored on perf/riccati tip due to pre-existing multi-n assertion crash; now passes thanks to the skip in 5acf1478)
Total	846 / 846 ✓

End-to-end --test runs of both LAR scans complete with Δ' values
bit-identical to pre-rename outputs:

Scan	Δ' values
LAR_beta_scan --test	dp21 = {+10.0150, +15.7659, +292.6038} for pc ∈ {0.001, 0.10, 0.17}
LAR_epsilon_scan --test	dp21 = {+9.2087, +5.5595, +2.4427} for ε ∈ {0.2495, 0.4072, 0.5510}

Regression-harness sweep against develop

Ran regress --cases diiid_n1,solovev_n1,solovev_multi_n --refs develop,local
(develop @ e0ee9dff, 2026-04-12). All differences are accounted for:

Expected shifts (introduced earlier in perf/riccati, intentional)

These show up across all three cases and reflect documented behavior changes:

Class	Cause	Approx scale
ODE step counts (integration/nstep_*) ↓ ~45-55 %	singfac_min = 1e-4 (was 0) tightens inner-layer cutoff (4845ec80); dense ξ pass uses serial-EL adaptive stepping (5acf1478)	DIIID 1526 → 756
Energy eigenvalues (et[1], ep[1], ev[1]) shift O(1)-O(100)	free_run! now runs on the dense-pass serial-EL odet (axis basis) rather than the parallel Riccati odet (renormalized basis)	DIIID et[1] 1.51 → −7.80
mpert change (DIIID 26 → 35)	delta_mhigh doubling in main pipeline + tightenings (54d12fe2)	—
Δ' diagonal magnitudes	Now computed by riccati_cross_ideal_singular_surf! in (S, I) Riccati gauge (4845ec80); old default path used a different normalization	within-branch bit-identical (pinned in tests)
q profile (checksum), pressure profile (checksum)	set_psilim_via_dmlim moved the integration limit (54d12fe2)	—

||resonant flux|| blow-up — diagnosed and fixed in this PR (8073c126)

Pre-fix the dense ξ pass produced max|resonant_flux| = 2.85e+23 on
diiid_n1 (vs 5.81e-03 on the pre-dense-pass perf/riccati tip
3c8130da). Cascading downstream: delta_prime, island_half_widths,
Chirikov parameter, resonant_current, penetrated_field all blew up.

Root cause: _populate_dense_xi_via_serial_el! (in 5acf1478)
replaced the parallel BVP's odet with a fresh serial-EL odet but
preserved only intr.sing[*] fields — silently letting the fresh EL
pass's axis-basis ca_l / ca_r overwrite the parallel BVP's
(S, I) Riccati-gauge values. PE's SingularCoupling.jl reads
ca_l[resnum, resnum, 2, s] / ca_r[resnum, resnum, 2, s] and is
written against the Riccati convention; in EL axis basis these
matrices carry the exponentially-growing U₁ at the inner-layer
boundary, inflating downstream resonant flux by ~25 orders of magnitude.

Fix (8073c126): also save odet.ca_l / odet.ca_r in the
helper's saved tuple and copy them onto fresh_odet after the dense
pass returns. The fresh EL odet's own ca_l/ca_r are discarded —
they're not needed because compute_delta_prime_matrix! uses
propagators/chunks (not ca) and the bit-identical-ξ test in
runtests_parallel_integration.jl only compares
u_store/ud_store/psi_store/q_store/crit_store (not ca).

Verification: after the fix, max|resonant_flux| = 0.005813333867190486
on diiid_n1 — bit-identical to 3c8130da (the perf/riccati
pre-dense-pass tip). Full test suite still 846/846.

solovev_n1 Ref 1 replay failure

The cached develop @ e0ee9dff entry for solovev_n1 fails to replay
inside the harness (runner.jl:384). Not a code issue — the harness
cache for that one case has gone stale. Local ran cleanly with
ODE steps = 274, et[1] = +9.985e-02, qmax = 3.147 (consistent
with the LAR scan results).

Known issues / open questions

• Multi-resonance surface support in compute_delta_prime_matrix!
  (this PR makes the path graceful but doesn't fully extend it). When a
  multi-n run produces a singular surface satisfying two distinct
  (m, n) tuples (e.g. q = 2 with both (m=2, n=1) and (m=4, n=2)),
  the inter-surface Δ' BVP is now skipped with a warning instead of
  crashing. Per-surface Δ' (singular/delta_prime, delta_prime_col)
  is still produced. Tracked as a follow-up: matrix shape becomes
  (n_res_total, n_res_total) with a (surface, mode, side) ↔ BVP-row
  map; the BVP matrix size stays (2 + 4·msing)·N so the algorithmic
  change is bookkeeping, not new math (or at least I don't think so?).

• Reverse-shear q-profile handling in sing_find! (pre-existing,
  separate code path from the multi-resonance issue above; not fixed by
  this PR). Question
  for review: can DCON or STRIDE type codes be run on a reverse shear equilibrium with multiple e.g. q = 2 crossings? If so, such an equilibrium where q
  crosses a rational value at two distinct ψ locations,
  the existing-surface check in sing_find! (Sing.jl:32) keys only
  on s.q, so the second crossing's (m, n) is pushed into the
  existing SingType's m/n vectors — silently losing the second ψ
  location. Two distinct singular surfaces collapse into one with
  multiplicity-2 m/n vectors.

Pre-merge checklist

• Full test suite green (846/846) on macOS-aarch64 Julia 1.11
• Per-line annotation pass on examples/LAR_*_scan/gpec.toml
• Fitzpatrick TJ citation at every TJ-analytic definition + use site
• Bit-identical ξ regression test on Solovev + DIIID-like
• Pinned Δ' values regression test on Solovev + DIIID-like
• TJ_epsilon_pole_example + its regression-harness case removed
• No trailing whitespace, LF line endings, TOML/YAML syntax valid
• No .h5 / figures/ / outputs/ artifacts committed
• Regression harness sweep against develop (diiid_n1, solovev_n1,
  solovev_multi_n) — all diffs accounted for; see "Regression-harness
  sweep" section above. The ||resonant flux|| blow-up surfaced by
  the sweep has been fixed in 8073c126.
• Multi-resonance Δ' BVP extension — separate follow-up PR
• Reverse-shear sing_find! dedup fix — possible separate PR?

[logan-nc]

@claude review this PR. Be succinct. Do not summarize. Just point out any bugs or potential improvements needed prior to merging.

[logan-nc]

@d-burg thanks for getting this to the finish line! I am out of the office today but will give it a read through as soon as I can. I really appreciate the great effort here. This is a 🏆 feature.

[claude]

Code Review — PR #178: Dual Riccati + Parallel FM integration + Δ' output

This PR implements two new integration paths and full Δ' computation. The physics is sound and the validation against Fortran STRIDE is well-documented. Issues below range from blocking to advisory.

---

Breaking default changes (most critical)

use_parallel = true is the new default. This is the single biggest concern. The parallel path:

• Runs the Δ' STRIDE BVP matrix on every run via compute_delta_prime_matrix! — extra cost not reflected in the benchmark table.
• With populate_dense_xi = true (also default), appends a full serial EL pass after the parallel BVP. The PR description benchmarks the Riccati path at 2.83× speedup, but the actual default runtime with use_parallel + populate_dense_xi is closer to 2× the serial path (parallel BVP ~1× + serial EL ~1×), which is a net regression, not a speedup. This should be measured and documented honestly before making it default.
• Is thread-nondeterministic. The test comment in runtests_fullruns.jl explicitly says: "rtol is loose because this result is thread-count sensitive (drifts ~15% between single- and multi-threaded invocations)." A 15% energy variation based on OS scheduler decisions is a correctness concern, not just a performance note.

Other changed defaults that will silently alter existing user results without opt-in:

• eulerlagrange_tolerance: 1e-7 → 1e-8 (tighter, but silent)
• singfac_min: 0.0 → 1e-4 (required nonzero for parallel path, but changes standard path behavior too)
• sing_order: 2 → 6 (more expensive asymptotic expansion everywhere)
• parallel_threads: 1 → 2

Recommendation: Default use_parallel = false until the thread-nondeterminism is resolved and the performance accounting includes the populate_dense_xi overhead.

---

use_double64_bvp missing from docstring

use_double64_bvp::Bool = true is added to ForceFreeStatesControl but not listed in the struct's docstring. Every other field has a bullet; this one is silently absent.

---

Fullruns test tolerance is not a meaningful regression test

In runtests_fullruns.jl, the test uses rtol = 0.2 (±20%) on an expected value of -0.18. The comment says the result drifts ±15% by thread count, so the test can pass and fail on the same code depending on how many Julia threads are available. This is a compile check, not a regression guard. The thread-nondeterminism needs to be resolved before this test can have a meaningful tolerance.

The Δ' regression tests in runtests_parallel_integration.jl use rtol=0.05. For the Solovev case, some Δ' values have imaginary parts 40–500× larger than the real parts (e.g., -7.278138e+00 + 4.172681e+03im). The tests should either target the physically meaningful component or explain why the complex value is the right quantity to pin.

---

benchmark_integration_paths.jl mutates production example files

The benchmark writes a modified gpec.toml in place and restores it in a finally block. A Julia crash (OOM, segfault, SIGKILL) will leave the example's gpec.toml modified. The benchmark should copy the example to a mktempdir and run there, as benchmark_xi_parallel_vs_serial.jl already does.

---

diagnose_profiles.jl hard-codes a private git branch

The file calls git show perf/riccati-full-geqdsk-scans:examples/LAR_epsilon_scan/equilibria/$fname. This branch is not part of the main repo and will fail for all contributors who don't have it locally. Document the prerequisite clearly or make the file self-contained.

---

Docs files not wired into Documenter

docs/delta_prime_numerical_analysis.md and docs/stride_delta_prime_validation.md are placed in docs/ (not docs/src/) and are not added to docs/make.jl. They won't appear in the public documentation site. Move them to docs/src/ and add them to make.jl if they are permanent reference material.

---

Minor issues

• balance_integration_chunks test: @test length(balanced) >= min(target_n, length(base_chunks) * 50) — in the common case this asserts >= target_n. The function produces at most target_n chunks; the test should use == or justify >=.

• Edge-scan restore logic: After the diagnostic-only scan, the code saves/restores odet.psifac and odet.u, but findmax_dW_edge! also writes to odet.edge_scan. The comment should document all mutated state so future maintainers don't inadvertently add more.

• Regression harness: The PR description compares develop @ e0ee9dff vs local, which conflates develop-side changes with PR changes. Per CLAUDE.md, the harness should compare develop vs the PR branch directly.

• benchmark_xi_parallel_vs_serial.jl calls Pkg.activate inside the script, overriding the project environment when the script is already expected to run with --project=... Remove the Pkg.activate call.

---

What works well

The bidirectional integration approach for crossing chunks is well-motivated: the direction field cleanly tracks backward chunks, balance_integration_chunks preserves it correctly when splitting, and apply_propagator_inverse! provides the matching inverse. The edge-dW decoupling from integration truncation is a genuine correctness fix for Δ' at the outermost surface, and the opt-in truncate_at_dW_peak flag preserves backward compatibility without silent behavior changes. The docs/src/stability.md documentation is thorough and useful.

[logan-nc]

@d-burg I am ok with changing defaults, as we are still in the chaos of active dev.
But I think the rest of the claude review points are worth addressing. Especially the danger of the gpec.toml restore approach. Can you knock 'em out?

[d-burg]

@logan-nc on it right now! These were helpful catches

[d-burg]

Proposal: top-level [Run] section with three physics flags

Currently what GPEC computes is gated by a scattered set of low-level flags plus implicit conditions: two same-named force_termination flags in different structs, use_parallel implicitly unlocking the inter-surface Δ' matrix, and [PerturbedEquilibrium] section-presence as the PE on/off switch. Before final cleanup I'd like to add a brief user-facing layer:

[Run]
compute_ideal_stability       = true   # δW + per-surface Δ'
compute_deltaprime_matrix     = false  # Inter-surface Δ' matrix (STRIDE BVP, Glasser 2018)
compute_perturbed_equilibrium = false  # GPEC response, island widths, Chirikov, resonant flux

The driver auto-derives ode_flag, vac_flag, use_parallel, populate_dense_xi, and the PE entry point from these three. Explicit low-level flags still work as overrides.

Mode coverage:

Mode	ideal_stability	deltaprime_matrix	perturbed_equilibrium
Equilibrium only	F	F	F
Ideal stability + δW	T	F	F
Δ' matrix (tearing study, no PE)	T	T	F
GPEC response (no Δ' matrix)	T	F	T
Full pipeline	T	T	T

When deltaprime_matrix=true the parallel FM path is used, which leaves odet.u_store in (S, I) Riccati basis and ud_store ≈ 0 – unusable by PerturbedEquilibrium.FieldReconstruction, which expects axis-basis (EL basis) dense ξ. The dense-ξ pass (_populate_dense_xi_via_serial_el!) re-runs a serial Euler–Lagrange integration after the parallel BVP to repopulate ξ in axis basis. So whenever deltaprime_matrix=true AND perturbed_equilibrium=true, the driver must force populate_dense_xi=true – otherwise PE silently reads Riccati-basis garbage. The proposed auto-derivation does exactly this. With only one of the two set, populate_dense_xi stays off (no wasted serial pass for tearing-only studies; PE-only runs use the standard serial path and don't need it either). Kinetic damping stays orthogonal (continuous kinetic_factor).

What do you think @logan-nc @matt-pharr ?

[logan-nc]

@d-burg did you address the review issues? Please push so I can do a final pass

[d-burg]

Two new commits address the latest pre-merge audit and bring the branch up to date with develop.

a89bd5d2

• Project.toml: Random moved to [extras]/[targets].test; DoubleFloats opt-in via ctrl.extended_precision_bvp (default true — Float64 drifts imaginary Δ' by 2–5× on DIIID).
• Stubs: SingType.delta_prime / delta_prime_col and compute_delta_prime_from_ca! explicitly annotated as stubs; canonical Δ' is ForceFreeStatesInternal.delta_prime_matrix (STRIDE BVP).
• Decompositions: compute_delta_prime_matrix! 540 → 63 LOC (+11 helpers), parallel_eulerlagrange_integration 281 → 36 LOC (+7), riccati_cross_ideal_singular_surf! 122 → 20 LOC (+6). All bit-identical.
• Logging: Per-crossing vmat/asymptotic dumps demoted from @info to @debug.
• Correctness: Assertion guards FM-axis BVP fallback against direction=-1 crossing chunks.
• Inline citations: Glasser 2018 Eqs. 19 / 29 / 31 / 33 / 37; STRIDE sing_vmat reference in Sing.jl.
• Naming: SAVE_NEAR_END_FRAC, SAVE_NEAR_END_PSI, ODE_COST_AXIS/RAT/EDGE constants; ucrit=1e4 rationale documented.
• Performance: Auto-skip the populate_dense_xi serial-EL pass when force_termination=true (eliminates ~1× wasted runtime in non-PE runs).
• Tests: Tightened Riccati Solovev rtol 1e-2 → 1e-4 (matches PR's 0.006% claim). Δ' matrix tolerances split by entry magnitude / |Im|/|Re| ratio (V4). Platform-tolerance brackets on FP-sensitive et[1] values. Fixed sing_lim! NaN when nn_low ≤ 0.

ace5e21c — Merge origin/develop

7 conflict files resolved:

• Project.toml — kept julia=1.11 (develop) + [extras] (ours).
• EquilibriumTypes.jl, GeneralizedPerturbedEquilibrium.jl, EulerLagrange.jl — preserved both sides.
• SingularCoupling.jl, Analysis/PerturbedEquilibrium.jl, PerturbedEquilibrium/Utils.jl — took develop's PE refactor in full; clarified that ForceFreeStates' BVP delta_prime_matrix (equilibrium-intrinsic) and PE's state.delta_prime (forcing-driven) are distinct physical quantities sharing only a symbol — both now written independently to HDF5.

Three latent runtime crashes from auto-merge gaps also patched: missing use_galgrid struct field; missing bt_sign arg in DirectRunInput calls (AnalyticEquilibrium.jl, ReadEquilibrium.jl IMAS reader); stale state.coupling_coefficient / state.resonant_amplitude HDF5 writes (fields removed by develop's refactor).

Verification

• Pkg.test() passes the full suite on the merged tree (Apple aarch64 / Julia 1.11).
• Local regression-harness runs cleanly on all three standard cases (diiid_n1, solovev_n1, solovev_multi_n); develop-side comparison currently blocked by a harness/Pkg.instantiate FastInterpolations-BC issue independent of this PR.

[d-burg]

Pre-merge audit cleanup — 4fbd882
Five small fixes addressing the most recent audit, all green on 19 testsets.

B1 — Per-thread ffit_hint in sing_der!. 21 hot-path callsites converted from shared ffit._hint::Ref{Int} to odet.ffit_hint (per-thread via odet_proxies). Honest framing: I built a reproducer on Apple M1 Max — 19 runs across julia_nthreads ∈ {1,4,8} × parallel_threads ∈ {2,8} were all bit-identical to -0.193593591803846. The race exists in source but FastInterpolations.RefHint is stale-tolerant so it doesn't produce drift in practice. Fix is defensive (completes the per-thread isolation pattern, removes a future-flake source if FastInterpolations' hint semantics change), not a drift fix.

B3 — assemble_fm_matrix size inference. Empty-range guard ordering fixed; size now inferred from T_init or propagators[first(idx_range)] instead of unconditionally dereferencing propagators[1].

H1 — parallel_threads honored in balance_integration_chunks. Uses min(nthreads, ctrl.parallel_threads); a julia -t 16 user setting parallel_threads=2 no longer pays for 8× the requested sub-chunks.

H2 / H3 — Comment cleanup. Drops re-introduced [Fortran sing.f lines XXXX] brackets in Sing.jl (Logan 2015 App. C eq. refs already carry the provenance); compresses the 7-line historical-narration block in the HDF5 writer to one line per CLAUDE.md.

populate_dense_xi default flip (true → false) — clarity, not perf rescue. The audit's "75 % regression" framing was extrapolated from a small-N benchmark and doesn't generalize. Empirical on examples/Solovev_ideal_example (delta_m=8): parallel+populate ≈ 97 s, serial only ≈ 494 s — parallel wins by 5× even with populate=true. The flip stands as a UI clarification (PE-using TOMLs opt in explicitly; non-PE TOMLs save ~10-30 % parallel-BVP wall-clock). 2 PE examples set populate_dense_xi = true explicitly with a strengthened comment; 4 non-PE examples flip to false; 4 test fixtures untouched to preserve pinned et[1] values.

Tightened kinetic multi-n rtol. runtests_fullruns.jl: decade-wide bracket → isapprox(et, -0.193593591803846; rtol=1e-3). The "Apple silicon drifts ~20 %" warning in the old comment doesn't reproduce on the current code path.