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.

---

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

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.

---

Files changed

File	Change
src/ForceFreeStates/ForceFreeStatesStructs.jl	IntegrationChunk.direction, SingType.delta_prime_col, ForceFreeStatesInternal.delta_prime_matrix; removed dead parallel_threads field
src/ForceFreeStates/EulerLagrange.jl	chunk_el_integration_bounds(bidirectional=false), balance_integration_chunks direction propagation
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
src/GeneralizedPerturbedEquilibrium.jl	HDF5 output for singular/delta_prime, delta_prime_col, delta_prime_matrix
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
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

Deferred: apply_propagator! per-call allocations · test boilerplate in runtests_parallel_integration.jl

---

🤖 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.