From ab58849823e2065302bff825208eb4de4503d5c1 Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Fri, 14 Aug 2026 22:14:14 -0400 Subject: [PATCH 1/3] FFS - REFACTOR - Unify Riccati integrator behind integrator= key Replace the use_riccati / use_parallel / parallel_threads / populate_dense_xi flag set with a single integrator = "forward" | "riccati" | "galerkin" key plus an nchunks knob, per the ForceFreeStates modularization plan (PR 1 of 5). - Delete the standalone serial-Riccati driver and the populate_dense_xi dense-xi re-run; no integration result is merged with another. - Rename serial_eulerlagrange_integration -> forward_eulerlagrange_integration and give the STRIDE propagator driver the freed riccati_ name. - Chunk sizing no longer reads Threads.nthreads(), so Riccati outputs are identical at any thread count; nchunks is the only tunable and clamps up to the singular-surface floor with a warning. - Kinetic runs now error on the Riccati path rather than silently mis-crossing. - The four retired keys join _DEPRECATED_FFS_KEYS and the toml-no-deprecated-keys hook, so old decks and gpec.h5 replays warn and continue. - New examples/DIIID-like_riccati_deltaprime_example plus the diiid_n1_riccati regression case keep the delta-prime matrix covered once the DIIID-like deck moves to the forward integrator. Verified on origin/develop: full suite 59 testsets / 2044 assertions / 0 failures; docs build clean; harness diiid_n1, solovev_n1 and diiid_slayer_n1 all unchanged; delta-prime bit-identical between -t 1 and -t 4. REFACTOR_PLAN.md is amended in the same commit: #363 and #345 merged before this sequence branched, so the plan now records that it builds on the new HDF5 schema. Co-Authored-By: Claude Opus 5 (1M context) --- .pre-commit-config.yaml | 2 +- REFACTOR_PLAN.md | 39 ++- benchmarks/benchmark_delta_prime_methods.jl | 4 +- benchmarks/benchmark_threads.jl | 32 +-- benchmarks/compare_gal_vs_el.jl | 2 +- benchmarks/compare_jbgradpsi_m2.jl | 14 +- benchmarks/plot_xi_eigenmode.jl | 2 +- benchmarks/scan_resistivity_m2.jl | 10 +- benchmarks/scan_rotation_m2.jl | 8 +- docs/src/citations.md | 2 +- docs/src/stability.md | 103 ++++--- examples/DIIID-like_SLAYER_example/gpec.toml | 6 +- .../gpec.toml | 6 +- .../gpec.toml | 6 +- examples/DIIID-like_ideal_example/gpec.toml | 6 +- .../gpec.toml | 57 ++++ examples/LAR_beta_scan/gpec.toml | 4 +- examples/LAR_epsilon_scan/gpec.toml | 4 +- examples/LAR_ideal_match_test/gpec.toml | 4 +- examples/LAR_resistive_match_test/gpec.toml | 4 +- examples/Solovev_ideal_example/gpec.toml | 6 +- examples/Solovev_ideal_example_3D/gpec.toml | 6 +- .../Solovev_ideal_example_multi_n/gpec.toml | 7 +- .../Solovev_kinetic_NTV_example/gpec.toml | 6 +- .../gpec.toml | 1 + examples/a10_kinetic_example/gpec.toml | 1 + regression-harness/cases/diiid_n1.toml | 22 +- .../cases/diiid_n1_riccati.toml | 163 +++++++++++ src/ForceFreeStates/EulerLagrange.jl | 59 ++-- src/ForceFreeStates/ForceFreeStatesStructs.jl | 16 +- src/ForceFreeStates/Galerkin/GalerkinMatch.jl | 2 +- src/ForceFreeStates/Riccati.jl | 267 +++--------------- src/GeneralizedPerturbedEquilibrium.jl | 7 +- src/PerturbedEquilibrium/SingularCoupling.jl | 6 +- src/Tearing/Dispersion/CoupledFullMatch.jl | 2 +- test/runtests_eulerlagrange.jl | 2 +- test/runtests_parallel_integration.jl | 172 +++++------ test/runtests_rerun_from_h5.jl | 17 ++ test/runtests_riccati.jl | 37 +-- .../gpec.toml | 6 +- .../gpec.toml | 6 +- .../gpec.toml | 1 + .../gpec.toml | 6 +- .../gpec.toml | 6 +- .../gpec.toml | 1 + 45 files changed, 560 insertions(+), 580 deletions(-) create mode 100644 examples/DIIID-like_riccati_deltaprime_example/gpec.toml create mode 100644 regression-harness/cases/diiid_n1_riccati.toml diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 54e179385..241049189 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -69,5 +69,5 @@ repos: - id: toml-no-deprecated-keys name: 'TOML conventions: no deprecated config keys' language: pygrep - entry: '^(mer_flag|force_wv_symmetry|ode_flag|cyl_flag|mat_flag|power_bp|power_b|power_r|power_rc)\s*=' + entry: '^(mer_flag|force_wv_symmetry|ode_flag|cyl_flag|mat_flag|use_riccati|use_parallel|parallel_threads|populate_dense_xi|power_bp|power_b|power_r|power_rc)\s*=' files: ^(examples/.*\.toml|test/test_data/.*\.toml)$ diff --git a/REFACTOR_PLAN.md b/REFACTOR_PLAN.md index ba6fa153b..0f5bdc00b 100644 --- a/REFACTOR_PLAN.md +++ b/REFACTOR_PLAN.md @@ -8,11 +8,17 @@ > changes ride with it) so everyone with open PRs can see what is coming and where it > will touch their work. Key coordination points: > -> - The PR sequence below assumes **nothing else merges into `develop` before it -> finishes**. If your PR must land mid-sequence, talk to Matthew first. -> - The HDF5 schema PRs (#363/#364) are treated as guidance on the final HDF5 shape; -> writers refactored here keep today's dataset paths and will be re-targeted by -> those PRs afterward. #367 (immutable control structs) merges after this sequence. +> - The PR sequence below assumes **nothing else merges into `develop` mid-sequence**. +> If your PR must land before it finishes, talk to Matthew first. +> - **Amendment (post-PR-1):** the module-mirroring HDF5 schema (#363) and the vacuum +> surface-inductance migration (#345) merged into develop before PR 1 branched, so +> the sequence is built on the NEW schema (`ForceFreeStates/…`, `SingularSurfaces/…`, +> `LocalStability/…`, `Equilibrium/…`). Writer refactors in PR 3/4 therefore target +> that schema directly; dataset-path references below have been updated accordingly. +> Only #364 (self-describing metadata) remains as a later re-target, and #367 +> (immutable control structs) still merges after this sequence. Neither merged PR +> changes any decision: #345 left `calc_surface_inductance` in PerturbedEquilibrium +> (only its Vacuum support types moved), so the PE consumer map stands. > - The regression harness will be re-baselined during this work; a fresh > Fortran-agreement comparison is the final validation gate for the whole sequence. > - The "Decision record" and "Verified code facts" sections are settled — please do @@ -53,8 +59,8 @@ pe = perturbed_equilibrium(ffs, rmp) | D7 | Public API via **CommonSolve.jl**: `solve(eq::PlasmaEquilibrium, alg; kwargs...)`. `PlasmaEquilibrium(path; kwargs...)` constructor. Module names unchanged. | | D8 | Galerkin standalone computes its own vacuum `wv` (no ODE state needed — verified); its result has `free_boundary = nothing`. | | D9 | TOML: new `integrator = "forward"|"riccati"|"galerkin"` key. Old keys (`use_riccati`, `use_parallel`, `parallel_threads`, `populate_dense_xi`, later `gal_flag`) go to the `_DEPRECATED_FFS_KEYS` warn-and-ignore list AND the `toml-no-deprecated-keys` pre-commit hook pattern. | -| D10 | No back-compat burden; examples/fixtures updated freely; regression re-baselining accepted (HDF5 schema PRs #363/#364 churn it anyway). Final validation = fresh Fortran comparison after the sequence. Assume nothing else merges first. | -| D11 | New structs immutable from day one (eases the later #367 merge). HDF5 writers become functions on result structs (schema itself unchanged here; #363/#364 re-target them later). | +| D10 | No back-compat burden; examples/fixtures updated freely; regression re-baselining accepted. Final validation = fresh Fortran comparison after the sequence. Nothing else merges mid-sequence without coordination (#363/#345 landed before PR 1 and are absorbed — see header amendment). | +| D11 | New structs immutable from day one (eases the later #367 merge). HDF5 writers become functions on result structs, keeping the merged #363 schema paths unchanged; #364 (metadata) re-targets them later. | | D12 | Analysis module reads HDF5 files, not live structs — untouched except where dataset names would change (they don't in this plan). | ### Verified code facts the workers must not re-derive @@ -227,7 +233,7 @@ All code must be JuliaFormatter-clean per `.JuliaFormatter.toml` before commit. `integrator="riccati"`) so the canonical Δ′-matrix fixture survives the DIIID-like_ideal switch to forward. Add a matching regression case `regression-harness/cases/diiid_n1_riccati.toml` tracking - `singular/delta_prime_matrix`-derived quantities (mirror the Δ′ entries of the + `SingularSurfaces/delta_prime_matrix`-derived quantities (mirror the Δ′ entries of the existing `diiid_n1` case; ξ/PE quantities stay on `diiid_n1`). - `benchmarks/benchmark_threads.jl`, `benchmarks/benchmark_delta_prime_methods.jl`: update flag names (`use_riccati`/`parallel_threads` → `integrator`/`nchunks`); @@ -323,7 +329,7 @@ All code must be JuliaFormatter-clean per `.JuliaFormatter.toml` before commit. Full test suite; targeted `runtests_resist_eval.jl`, `runtests_fullruns.jl` (exercises `local_stability_flag=true` decks); docs build (missing-docs gate); -harness `--cases diiid_n1 --refs develop,local` (locstab datasets must be identical). +harness `--cases diiid_n1 --refs develop,local` (`LocalStability/*` datasets must be identical; note the #363 group name already matches the new module name). --- @@ -370,8 +376,8 @@ end Contract of `solution`/`solution_basis`: - Forward → dense EL-basis odet, `:el_axis`. PE-usable. - Riccati → its odet IS carried (psi/q/crit/edge_scan are valid, u_store is chunk - snapshots), `:riccati`. NOT PE-usable; HDF5 `integration/psi|crit|EdgeScan` still - written from it, `integration/xi_*` written empty. + snapshots), `:riccati`. NOT PE-usable; HDF5 `ForceFreeStates/Solutions/ForwardIntegration/psi|crit` and `EdgeScan/*` still + written from it, `.../ForwardIntegration/xi_*` written empty. - Galerkin with `gal_match_flag` → `gal_matched_odestate(...)`, `:gal_native` (PE-usable; `du_store_populated=true` analytic derivatives). Without match → `nothing`, `:none`. @@ -454,7 +460,7 @@ crosses module boundaries after `build_result`. ### 5.4 Verification Full suite; `runtests_fullruns.jl` (all decks — forward decks produce identical -HDF5 vs pre-PR, riccati decks now emit empty `integration/xi_*` + PE-skip warnings); +HDF5 vs pre-PR, riccati decks now emit empty `ForwardIntegration/xi_*` + PE-skip warnings); harness `--cases diiid_n1,solovev_n1 --refs develop,local` (forward-deck tracked quantities unchanged); docs build. @@ -512,14 +518,14 @@ logic stays a pre-FFS stage but is owned by the FFS-facing function disappear from the HDF5 — accepted per D10; gal datasets identical because `galerkin_solve` inputs are unchanged). `gal_*` sub-knobs stay (they become `Galerkin(...)` fields in PR 5). -- FFS writer: tolerate `solution === nothing` (write empty `integration/*` datasets — +- FFS writer: tolerate `solution === nothing` (write empty `ForwardIntegration/*` datasets — extend the existing empty-fallback pattern). ### 6.3 Tests / verification - Update `runtests_fullruns.jl` gal decks' expectations (gal-only HDF5). - New testset (in `runtests_fullruns.jl` or the gal tests): `integrator="galerkin"` - on `LAR_ideal_match_test` produces `galerkin/delta` identical to the PR-3 additive + on `LAR_ideal_match_test` produces the Galerkin `delta` dataset identical to the PR-3 additive run (same `wv` by construction — assert against a stored reference or a paired riccati+gal_flag run on the pre-PR commit during development). - Full suite; harness (gal cases if present, plus diiid/solovev); docs. @@ -660,7 +666,10 @@ manual smoke: run the 4-line UX from the Context section in a REPL against ## 10. Progress -- [ ] PR 1 — `refactor/riccati-unification` +- [ ] PR 1 — `refactor/riccati-unification` — **implemented, in review.** Two deltas + from the §3 spec, both improvements: the new Δ′ example references the DIIID geqdsk + by relative path instead of copying it, and the TOML sweep covered six regression + fixtures (two more had landed on develop since the plan was written), all `forward`. - [ ] PR 2 — `refactor/local-stability-module` - [ ] PR 3 — `refactor/forcefreestates-result` - [ ] PR 4 — `refactor/staged-main` diff --git a/benchmarks/benchmark_delta_prime_methods.jl b/benchmarks/benchmark_delta_prime_methods.jl index 4bf77f179..917e7c9c5 100644 --- a/benchmarks/benchmark_delta_prime_methods.jl +++ b/benchmarks/benchmark_delta_prime_methods.jl @@ -23,7 +23,7 @@ function setup_and_run_solovev() ex = joinpath(@__DIR__, "..", "test", "test_data", "regression_solovev_ideal_example") inputs = TOML.parsefile(joinpath(ex, "gpec.toml")) inputs["ForceFreeStates"]["verbose"] = false - inputs["ForceFreeStates"]["use_riccati"] = true + inputs["ForceFreeStates"]["integrator"] = "riccati" intr = FFS.ForceFreeStatesInternal(; dir_path=ex) ctrl = FFS.ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...) @@ -42,7 +42,7 @@ function setup_and_run_solovev() intr.numpert_total = intr.mpert * intr.npert metric = FFS.make_metric(equil, intr.mpert) ffit = FFS.make_matrix(equil, intr, metric) - odet = FFS.riccati_eulerlagrange_integration(ctrl, equil, ffit, intr) + odet, _, _, _ = FFS.riccati_eulerlagrange_integration(ctrl, equil, ffit, intr) return ctrl, equil, ffit, intr, odet end diff --git a/benchmarks/benchmark_threads.jl b/benchmarks/benchmark_threads.jl index 048c64c66..96e37236a 100644 --- a/benchmarks/benchmark_threads.jl +++ b/benchmarks/benchmark_threads.jl @@ -1,17 +1,17 @@ -# Thread-scaling benchmark for the bidirectional parallel FM integration. -# Runs the Solovev (N=8) and DIIID-like (N=26) examples with use_parallel=true -# across 1, 2, 4, 8 threads and compares against the serial Riccati path. +# Thread-scaling benchmark for the Riccati chunked propagator integration. +# Runs the Solovev (N=8) and DIIID-like (N=26) examples with integrator="riccati" +# across 1, 2, 4, 8 threads and compares against the forward path. The chunk list is +# thread-independent, so et[1] must not move with the thread count — only wall-clock. # # Usage (from JPEC_main root): # for t in 1 2 4 8; do julia -t $t --project=. benchmarks/benchmark_threads.jl; done using GeneralizedPerturbedEquilibrium, TOML, Printf, Statistics -function run_ffs(ex; use_parallel, use_riccati=false) +function run_ffs(ex; integrator) inputs = TOML.parsefile(joinpath(ex, "gpec.toml")) inputs["ForceFreeStates"]["verbose"] = false - inputs["ForceFreeStates"]["use_parallel"] = use_parallel - inputs["ForceFreeStates"]["use_riccati"] = use_riccati + inputs["ForceFreeStates"]["integrator"] = integrator inputs["ForceFreeStates"]["write_outputs_to_HDF5"] = false intr = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesInternal(; dir_path=ex) ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; @@ -36,17 +36,17 @@ function run_ffs(ex; use_parallel, use_riccati=false) return real(vac.et[1]), intr.numpert_total end -function timed_run(ex; use_parallel, use_riccati=false, nwarm=1, nrep=2) +function timed_run(ex; integrator, nwarm=1, nrep=2) # Warmup for _ in 1:nwarm - run_ffs(ex; use_parallel, use_riccati) + run_ffs(ex; integrator) end # Timed runs times = Float64[] local et1, N for _ in 1:nrep t0 = time() - et1, N = run_ffs(ex; use_parallel, use_riccati) + et1, N = run_ffs(ex; integrator) push!(times, time() - t0) end return mean(times), et1, N @@ -60,18 +60,14 @@ diiid_ex = joinpath(root, "examples", "DIIID-like_ideal_example") println("\n=== Thread-scaling benchmark ($(nthreads) thread(s)) ===\n") for (label, ex) in [("Solovev", sol_ex), ("DIIID-like", diiid_ex)] - t_std, et_std, N = timed_run(ex; use_parallel=false, use_riccati=false) - t_ric, et_ric, _ = timed_run(ex; use_parallel=false, use_riccati=true) - t_par, et_par, _ = timed_run(ex; use_parallel=true, use_riccati=false) + t_fwd, et_fwd, N = timed_run(ex; integrator="forward") + t_ric, et_ric, _ = timed_run(ex; integrator="riccati") - err_ric = abs(et_ric - et_std) / abs(et_std) * 100 - err_par = abs(et_par - et_std) / abs(et_std) * 100 + err_ric = abs(et_ric - et_fwd) / abs(et_fwd) * 100 println("$label (N=$N, nthreads=$nthreads)") - @printf(" standard et[1]=%.5f t=%.2fs speedup=1.00×\n", et_std, t_std) + @printf(" forward et[1]=%.5f t=%.2fs speedup=1.00×\n", et_fwd, t_fwd) @printf(" riccati et[1]=%.5f t=%.2fs speedup=%.2f× err=%.4f%%\n", - et_ric, t_ric, t_std/t_ric, err_ric) - @printf(" parallel et[1]=%.5f t=%.2fs speedup=%.2f× err=%.4f%%\n", - et_par, t_par, t_std/t_par, err_par) + et_ric, t_ric, t_fwd/t_ric, err_ric) println() end diff --git a/benchmarks/compare_gal_vs_el.jl b/benchmarks/compare_gal_vs_el.jl index cff8b26d4..55f747f16 100644 --- a/benchmarks/compare_gal_vs_el.jl +++ b/benchmarks/compare_gal_vs_el.jl @@ -5,7 +5,7 @@ # gal: ξ_gal(ψ) = U_gal(ψ) · w (identity-at-edge ⇒ coefficient is w itself) # w = eigenvector of the total energy operator W = W_plasma + W_vacuum (FreeBoundaryStability). # -# Needs ONE gpec.h5 from a run with populate_dense_xi=true (EL dense u_store), gal_match_flag=true, +# Needs ONE gpec.h5 from a run with integrator="forward" (EL dense u_store), gal_match_flag=true, # gal_ideal_flag=true (gal ideal matched), vac_flag=true (energy operator). # Usage: julia --project=. benchmarks/compare_gal_vs_el.jl [gpec.h5] [out.png] [mode] # mode = "highest" (default, most stable), "lowest" (most unstable), or an integer eigenmode index. diff --git a/benchmarks/compare_jbgradpsi_m2.jl b/benchmarks/compare_jbgradpsi_m2.jl index dd30afbd2..15047e6fd 100644 --- a/benchmarks/compare_jbgradpsi_m2.jl +++ b/benchmarks/compare_jbgradpsi_m2.jl @@ -1,16 +1,16 @@ # Compare the area-normalized b^ψ (PerturbedEquilibrium/Response/psi_area = b^ψ/⟨J·|∇ψ|⟩_θ) for one # poloidal harmonic between two GPEC runs that are identical except for which ξ feeds PerturbedEquilibrium: # (1) IDEAL galerkin matched ξ (gal_match_flag=true, gal_ideal_flag=true) -# (2) SHOOTING ξ (gal_match_flag=false) +# (2) FORWARD ξ (gal_match_flag=false) # # PE writes no ψ grid, so it's reconstructed: gal-ideal → ForceFreeStates/Solutions/GalerkinIntegration/Solution/psi minus issing points; -# shooting → ForceFreeStates/Solutions/ForwardIntegration/psi. +# forward → ForceFreeStates/Solutions/ForwardIntegration/psi. # Usage: julia --project=. benchmarks/compare_jbgradpsi_m2.jl [gal_h5] [shoot_h5] [out.png] [m] using HDF5, Plots, Printf gal_h5 = length(ARGS) >= 1 ? ARGS[1] : "/tmp/gal_ideal_test/gpec.h5" -sh_h5 = length(ARGS) >= 2 ? ARGS[2] : "/tmp/shooting_test/gpec.h5" +sh_h5 = length(ARGS) >= 2 ? ARGS[2] : "/tmp/forward_test/gpec.h5" outpng = length(ARGS) >= 3 ? ARGS[3] : joinpath(@__DIR__, "cmp_jbgradpsi_m2.png") mtarget = length(ARGS) >= 4 ? parse(Int, ARGS[4]) : 2 @@ -23,7 +23,7 @@ pa_g, psi_g, mlow, sing_psi, sing_m = h5open(gal_h5) do f (pa, read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/psi"])[.!iss], read(f["Info/mlow"]), read(f["SingularSurfaces/GalerkinDeltaPrime/sing_psi"]), read(f["SingularSurfaces/GalerkinDeltaPrime/sing_m"])) end -# shooting run: PE grid = ForceFreeStates/Solutions/ForwardIntegration/psi +# forward run: PE grid = ForceFreeStates/Solutions/ForwardIntegration/psi pa_s, psi_s = h5open(sh_h5) do f (to_c(read(f["PerturbedEquilibrium/Response/psi_area"])), read(f["ForceFreeStates/Solutions/ForwardIntegration/psi"])) end @@ -43,15 +43,15 @@ if !isnan(psi_res) ig = argmin(abs.(psi_g .- psi_res)) is = argmin(abs.(psi_s .- psi_res)) @printf("m=%d resonant surface ψ=%.4f\n", mtarget, psi_res) - @printf(" nearest |b^ψ_area|: gal-ideal=%.4e (ψ=%.4f) shooting=%.4e (ψ=%.4f) ratio=%.3f\n", + @printf(" nearest |b^ψ_area|: gal-ideal=%.4e (ψ=%.4f) forward=%.4e (ψ=%.4f) ratio=%.3f\n", abs(g[ig]), psi_g[ig], abs(s[is]), psi_s[is], abs(g[ig]) / abs(s[is])) end plt = plot(; size=(1000, 620), xlabel="ψ_N", ylabel="|b^ψ / ⟨J·|∇ψ|⟩| (area-normalized)", - title="m=$mtarget perturbed normal field — gal-ideal vs shooting PE", legend=:topleft, + title="m=$mtarget perturbed normal field — gal-ideal vs forward PE", legend=:topleft, left_margin=13Plots.mm, bottom_margin=5Plots.mm, right_margin=4Plots.mm) plot!(plt, psi_g, abs.(g); lw=2.2, color=1, label="gal-ideal ξ → PE") -plot!(plt, psi_s, abs.(s); lw=1.8, color=2, ls=:dash, label="shooting ξ → PE") +plot!(plt, psi_s, abs.(s); lw=1.8, color=2, ls=:dash, label="forward ξ → PE") for (j, ps) in enumerate(sing_psi) ism2 = sing_m[j] == mtarget vline!(plt, [ps]; color=(ism2 ? :red : :black), ls=:dot, lw=(ism2 ? 1.8 : 1), diff --git a/benchmarks/plot_xi_eigenmode.jl b/benchmarks/plot_xi_eigenmode.jl index c4be18af7..faea03202 100644 --- a/benchmarks/plot_xi_eigenmode.jl +++ b/benchmarks/plot_xi_eigenmode.jl @@ -5,7 +5,7 @@ # largest Re(eigenvalue) and reconstructs its radial profile by projecting the EL fundamental matrix # (ForceFreeStates/Solutions/ForwardIntegration/xi_psi) onto that edge eigenvector: c = U_edge \ w, ξ(ψ) = U(ψ)·c. # -# Requires a run with populate_dense_xi=true so ForceFreeStates/Solutions/ForwardIntegration/xi_psi is the dense axis-basis fundamental +# Requires a run with integrator="forward" so ForceFreeStates/Solutions/ForwardIntegration/xi_psi is the dense axis-basis fundamental # matrix (not the Riccati S-matrices). # Usage: julia --project=. benchmarks/plot_xi_eigenmode.jl [path/to/gpec.h5] [out.png] diff --git a/benchmarks/scan_resistivity_m2.jl b/benchmarks/scan_resistivity_m2.jl index 6a5695efd..83908201d 100644 --- a/benchmarks/scan_resistivity_m2.jl +++ b/benchmarks/scan_resistivity_m2.jl @@ -1,6 +1,6 @@ # Plot the m=2 area-normalized b^ψ (PerturbedEquilibrium/Response/psi_area) across a resistivity scan # of the RESISTIVE gal matched PE runs (gal_match_flag=true, gal_ideal_flag=false), one curve per η. -# Overlays the shooting (ideal, η→0) reference. The η-scan dirs are produced by the bash loop over +# Overlays the forward (ideal, η→0) reference. The η-scan dirs are produced by the bash loop over # /tmp/etascan_ (each a copy of the 0.993 config with gal_eta scaled). # Usage: julia --project=. benchmarks/scan_resistivity_m2.jl [out.png] [m] @@ -46,10 +46,10 @@ for (i, (d, e)) in enumerate(zip(scandirs, etas)) psi, v = read_m2(joinpath(d, "gpec.h5"); gal=true) plot!(plt, psi, abs.(v); color=cols[i], lw=2, label=@sprintf("η = %g×η_ref", round(e / eta_ref; sigdigits=2))) end -# shooting (ideal, η→0) reference -if isfile("/tmp/shooting_test/gpec.h5") - psis, vs = read_m2("/tmp/shooting_test/gpec.h5"; gal=false) - plot!(plt, psis, abs.(vs); color=:black, ls=:dash, lw=2, label="shooting (ideal, η→0)") +# forward (ideal, η→0) reference +if isfile("/tmp/forward_test/gpec.h5") + psis, vs = read_m2("/tmp/forward_test/gpec.h5"; gal=false) + plot!(plt, psis, abs.(vs); color=:black, ls=:dash, lw=2, label="forward (ideal, η→0)") end isnan(psi_res) || vline!(plt, [psi_res]; color=:red, ls=:dot, lw=1.6, label="q=$mtarget surface") diff --git a/benchmarks/scan_rotation_m2.jl b/benchmarks/scan_rotation_m2.jl index bd01dfe4e..70acebea7 100644 --- a/benchmarks/scan_rotation_m2.jl +++ b/benchmarks/scan_rotation_m2.jl @@ -1,7 +1,7 @@ # Plot the m=2 area-normalized b^ψ (PerturbedEquilibrium/Response/psi_area) across a ROTATION scan # of the resistive gal matched PE runs (gal_match_flag=true, gal_ideal_flag=false), fixed η=8e-8, # rotation f = 1,2,4,8,16 Hz (forced eigenvalue γ_s = 2πi·n·f). One curve per rotation; overlays the -# shooting (ideal) reference. Scan dirs produced by the bash loop over /tmp/rotscan_. +# forward (ideal) reference. Scan dirs produced by the bash loop over /tmp/rotscan_. # Usage: julia --project=. benchmarks/scan_rotation_m2.jl [out.png] [m] using HDF5, Plots, Printf, TOML @@ -41,9 +41,9 @@ for (i, (d, r)) in enumerate(zip(scandirs, rots)) psi, v = read_m2(joinpath(d, "gpec.h5"); gal=true) plot!(plt, psi, abs.(v); color=cols[i], lw=2, label=@sprintf("f = %g Hz", r)) end -if isfile("/tmp/shooting_test/gpec.h5") - psis, vs = read_m2("/tmp/shooting_test/gpec.h5"; gal=false) - plot!(plt, psis, abs.(vs); color=:black, ls=:dash, lw=2, label="shooting (ideal)") +if isfile("/tmp/forward_test/gpec.h5") + psis, vs = read_m2("/tmp/forward_test/gpec.h5"; gal=false) + plot!(plt, psis, abs.(vs); color=:black, ls=:dash, lw=2, label="forward (ideal)") end isnan(psi_res) || vline!(plt, [psi_res]; color=:red, ls=:dot, lw=1.6, label="q=$mtarget surface") diff --git a/docs/src/citations.md b/docs/src/citations.md index 143b9a493..ea5ddaee3 100644 --- a/docs/src/citations.md +++ b/docs/src/citations.md @@ -36,7 +36,7 @@ The primary reference for the `ForceFreeStates` module. Derives the Euler-Lagran > *Physics of Plasmas* **25**, 032507 (2018). > DOI: [10.1063/1.5007042](https://doi.org/10.1063/1.5007042) -Reformulates the DCON eigenvalue problem as a Riccati matrix ODE, enabling parallel integration across singular surfaces and faster computation. Implemented in `src/ForceFreeStates/Riccati.jl` and enabled via `use_riccati = true` in `[ForceFreeStates]`. +Reformulates the DCON eigenvalue problem as a Riccati matrix ODE, enabling parallel integration across singular surfaces and faster computation. Implemented in `src/ForceFreeStates/Riccati.jl` and enabled via `integrator = "riccati"` in `[ForceFreeStates]`. --- diff --git a/docs/src/stability.md b/docs/src/stability.md index 3e0474482..0e8b72880 100644 --- a/docs/src/stability.md +++ b/docs/src/stability.md @@ -45,66 +45,70 @@ finite across every rational surface. ## Integration methods -Three integration drivers are available, all solving the same EL system but with different -numerical strategies. +Two integration drivers are available. Both solve the same EL system, but they differ in +numerical strategy and in what they leave behind for the rest of the pipeline: the forward +driver returns dense displacement profiles, the Riccati driver returns the inter-surface +``\Delta'`` matrix. -### Standard integration +### Forward integration -`eulerlagrange_integration` is the baseline driver. It integrates the EL ODE directly in -``(U_1, U_2)`` using Tsit5 with adaptive step control. Near each rational surface the +`forward_eulerlagrange_integration` is the baseline driver. It integrates the EL ODE directly +in ``(U_1, U_2)`` using Tsit5 with adaptive step control. Near each rational surface the columns of ``U_2`` that correspond to resonant modes are zeroed via Gaussian reduction (GR), -keeping the solution bounded. This is the reference path for correctness comparisons. +keeping the solution bounded; `transform_u!` undoes the reduction at the end, so `u_store` +comes back dense in the axis (Euler-Lagrange) basis. This is the reference path for +correctness comparisons, the only path whose solution the perturbed-equilibrium stage can +consume, and the only path that supports `kinetic_factor > 0`. -Enable with (default): +Enable with: ```toml [ForceFreeStates] -use_riccati = false -use_parallel = false +integrator = "forward" ``` ### Riccati integration -`riccati_eulerlagrange_integration` reformulates the problem in terms of the dual Riccati -matrix ``S = U_1 \cdot U_2^{-1}`` [Glasser 2018a, Eq. 19]: +`riccati_eulerlagrange_integration` (the default) decomposes the radial domain into +independent chunks, integrates each chunk's fundamental-matrix (FM) propagator in parallel +using `Threads.@threads`, then multiplies the propagators in order and applies each +singular-surface crossing serially. It is the only driver that produces the inter-surface +``\Delta'`` matrix. Because chunk endpoints are all it stores, `u_store` is sparse and stays +in the Riccati basis — dense ``\xi`` profiles and the +`ForceFreeStates/Solutions/ForwardIntegration/xi_*` datasets require the forward driver. + +Crossings and the outer-plasma re-integration use the dual Riccati matrix +``S = U_1 \cdot U_2^{-1}`` [Glasser 2018a, Eq. 19]: ```math \frac{dS}{d\psi} = w^\dagger \bar{F}^{-1} w - S\bar{G}S, \qquad w = Q - \bar{K}S. ``` -``S`` remains bounded near rational surfaces (where ``U_1, U_2`` grow exponentially), so the -solver takes fewer steps. Rather than integrating the quadratic Riccati ODE directly (which -blows up when ``|S|`` is large), the code integrates the linear EL system with -`sing_der!` as the RHS and recovers ``S = U_1 U_2^{-1}`` via periodic renormalization — an -approach that is mathematically equivalent to O(Δψ) but uses the ODE solver's full 5th-order -accuracy. - -Renormalization is triggered whenever ``\max(|U_1|)`` or ``\max(|U_2|)`` exceeds the -threshold `ucrit` (default 1e6), and is forced at the end of each chunk. At singular surface -crossings, `riccati_cross_ideal_singular_surf!` applies the small-asymptotic matching -directly in column `ipert_res` — without Gaussian reduction — and renormalizes to ``(S, I)``. +``S`` remains bounded near rational surfaces (where ``U_1, U_2`` grow exponentially). Rather +than integrating the quadratic Riccati ODE directly (which blows up when ``|S|`` is large), +the code integrates the linear EL system with `sing_der!` as the RHS and recovers +``S = U_1 U_2^{-1}`` via periodic renormalization — an approach that is mathematically +equivalent to O(Δψ) but uses the ODE solver's full 5th-order accuracy. Renormalization is +triggered whenever ``\max(|U_1|)`` or ``\max(|U_2|)`` exceeds the threshold `ucrit`, and is +forced at the end of each chunk. At singular surface crossings, +`riccati_cross_ideal_singular_surf!` applies the small-asymptotic matching directly in column +`ipert_res` — without Gaussian reduction — and renormalizes to ``(S, I)``. Enable with: ```toml [ForceFreeStates] -use_riccati = true -use_parallel = false +integrator = "riccati" +nchunks = 0 # 0 = auto: derived from the singular-surface count alone ``` -**Speedup** (benchmarked on reference examples): - -| Example | N modes | Speedup vs standard | -|---------|---------|---------------------| -| Solovev | 8 | ~1.6× (1 thread), ~2.8× (4 threads) | -| DIIID | 26 | ~2.0× (1 thread), ~1.3× (4 threads) | - -### Parallel fundamental-matrix (FM) integration +#### Chunking and thread independence -`parallel_eulerlagrange_integration` decomposes the radial domain into independent chunks and -integrates each chunk in parallel using `Threads.@threads`. Each chunk produces a -fundamental-matrix (FM) propagator. Serial post-processing multiplies the propagators in -order and applies each singular-surface crossing, recovering the same EL trajectory as the -Riccati path. +`balance_integration_chunks` splits the base chunks until the count reaches a target derived +from the number of singular surfaces, `max(2 m_s + 3, 8(m_s + 1) + m_s)`, where ``m_s`` is +`msing`. Setting `nchunks` overrides that target; a value below the ``2 m_s + 3`` floor is +clamped up with a warning. The target never consults `Threads.nthreads()`, and every chunk +integrates independently from identity initial conditions, so the results do not depend on how +many threads `julia -t` provides — threads change wall-clock only. #### Bidirectional integration for large N @@ -126,15 +130,9 @@ The implementation uses a `direction` field on `IntegrationChunk`: crossing chunk. `balance_integration_chunks` preserves this: the sub-chunk closest to the rational surface inherits `direction`, while the earlier sub-chunk always gets `direction=+1`. -Enable with: -```toml -[ForceFreeStates] -use_parallel = true -``` - -**Accuracy** (N=26, DIIID-like example): energy eigenvalue within 2% of standard path. +**Accuracy** (N=26, DIIID-like example): energy eigenvalue within 2% of the forward path. The residual ~2% gap comes from the different crossing convention (Riccati-style direct -zeroing vs GR), not from ODE tolerance; it is present in both 1-thread and 4-thread runs. +zeroing vs GR), not from ODE tolerance; it is present at every thread count. ## Local stability: Mercier and ballooning (s–α) @@ -230,7 +228,7 @@ propagator blocks from bidirectional integration rather than the monolithic forw where ``\Phi_R[j]`` is the forward FM product from ``\psi_{R,j-1}`` to the junction, and ``\Phi_L[j]`` is the backward crossing FM from ``\psi_{L,j}`` to the junction. -The matrix is only populated by the parallel FM path and is written to the HDF5 output +The matrix is only populated by the Riccati path and is written to the HDF5 output under `SingularSurfaces/delta_prime_matrix`. ## Configuration reference @@ -240,8 +238,8 @@ All `ForceFreeStates` options are set in the `[ForceFreeStates]` section of `gpe ```toml [ForceFreeStates] # Integration driver -use_riccati = false # true: Riccati path (faster, same accuracy) -use_parallel = false # true: parallel FM path (multi-thread, large N) +integrator = "riccati" # "forward" for dense xi profiles and kinetic runs +nchunks = 0 # Riccati chunk-count target (0 = auto, from msing alone) # Mode space nn_low = 1 # lowest toroidal mode number @@ -307,8 +305,9 @@ metric = FFS.make_metric(equil, intr.mpert) ffit = FFS.make_matrix(equil, intr, metric) # Choose integration driver. The top-level `eulerlagrange_integration` dispatches -# to the parallel or Riccati path based on ctrl.use_parallel / ctrl.use_riccati, -# and always returns a 4-tuple (odet, propagators, chunks, S_at_surface_left). +# on ctrl.integrator and always returns a 4-tuple +# (odet, propagators, chunks, S_at_surface_left). The trailing three are `nothing` +# on the forward path. odet, _, _, _ = FFS.eulerlagrange_integration(ctrl, equil, ffit, intr) vac = FFS.free_run(odet, ctrl, equil, ffit, intr) @@ -325,10 +324,10 @@ for s in 1:intr.msing end ``` -### Access inter-surface Δ' matrix (parallel FM path) +### Access inter-surface Δ' matrix (Riccati path) ```julia -# intr.delta_prime_matrix is msing × msing after parallel_eulerlagrange_integration. +# intr.delta_prime_matrix is msing × msing after riccati_eulerlagrange_integration. # Internally the solver builds a 2·msing × 2·msing raw matrix; the stored Δ' is # the PEST3 four-term combination that folds the raw block into a per-surface # tearing parameter. diff --git a/examples/DIIID-like_SLAYER_example/gpec.toml b/examples/DIIID-like_SLAYER_example/gpec.toml index a89b382f6..7065439c5 100644 --- a/examples/DIIID-like_SLAYER_example/gpec.toml +++ b/examples/DIIID-like_SLAYER_example/gpec.toml @@ -51,10 +51,8 @@ save_interval = 3 # Save every Nth ODE step (1=all, 10=every 10th). singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced ucrit = 1e4 # Maximum fraction of solutions allowed before re-normalized -# Δ' BVP + parallel integration (see ForceFreeStatesControl docstring for details) -use_parallel = true # Run parallel FM-propagator BVP path (unlocks SingularSurfaces/delta_prime_matrix) -parallel_threads = 1 # serial/bit-deterministic BVP — keeps the regression Δ' (and hence γ) reproducible -populate_dense_xi = false # No PerturbedEquilibrium here; the dense EL pass has no consumer (auto-disabled under force_termination anyway). SLAYER needs only delta_prime_matrix from the parallel BVP. +# Integrator selection and domain truncation (see ForceFreeStatesControl docstring for details) +integrator = "riccati" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic) or "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix) set_psilim_via_dmlim = true # TRUE for diverted geqdsks — q → ∞ at separatrix, so dmlim truncation avoids the δW kink instability at negligible domain cost dmlim = 0.2 # Truncate integration at (last_rational_q + dmlim) / n diff --git a/examples/DIIID-like_gal_resistive_example/gpec.toml b/examples/DIIID-like_gal_resistive_example/gpec.toml index f84d8c52a..cbb702db5 100644 --- a/examples/DIIID-like_gal_resistive_example/gpec.toml +++ b/examples/DIIID-like_gal_resistive_example/gpec.toml @@ -50,10 +50,8 @@ save_interval = 3 # Save every Nth ODE step (1=all, 10=every 10th). singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced ucrit = 1e4 # Maximum fraction of solutions allowed before re-normalized -# Δ' BVP + parallel integration (see ForceFreeStatesControl docstring for details) -use_parallel = true # Run parallel FM-propagator BVP path (unlocks SingularSurfaces/delta_prime_matrix) -parallel_threads = 2 # BVP thread cap (1 = serial/bit-deterministic; 2 ≈ +20% speedup; ≥3 saturates) -populate_dense_xi = false # Dense axis-basis ξ for the FFS HDF5 output. Not needed here: no [PerturbedEquilibrium] section, and the gal-matched path builds its own dense ξ. Set true only for a shooting-fed PE run. +# Integrator selection and domain truncation (see ForceFreeStatesControl docstring for details) +integrator = "riccati" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic) or "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix) set_psilim_via_dmlim = false # Keep psilim at psihigh (do not truncate at last_rational_q + dmlim) dmlim = 0.2 # Truncate integration at (last_rational_q + dmlim) / n (only used when set_psilim_via_dmlim = true) diff --git a/examples/DIIID-like_gal_resistive_pe_example/gpec.toml b/examples/DIIID-like_gal_resistive_pe_example/gpec.toml index e5536c75e..82001c077 100644 --- a/examples/DIIID-like_gal_resistive_pe_example/gpec.toml +++ b/examples/DIIID-like_gal_resistive_pe_example/gpec.toml @@ -49,10 +49,8 @@ save_interval = 3 # Save every Nth ODE step (1=all, 10=every 10th). singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced ucrit = 1e4 # Maximum fraction of solutions allowed before re-normalized -# Δ' BVP + parallel integration (see ForceFreeStatesControl docstring for details) -use_parallel = true # Run parallel FM-propagator BVP path (unlocks SingularSurfaces/delta_prime_matrix) -parallel_threads = 2 # BVP thread cap (1 = serial/bit-deterministic; 2 ≈ +20% speedup; ≥3 saturates) -populate_dense_xi = false # Dense axis-basis ξ for the FFS HDF5 output. Not needed here: no [PerturbedEquilibrium] section, and the gal-matched path builds its own dense ξ. Set true only for a shooting-fed PE run. +# Integrator selection and domain truncation (see ForceFreeStatesControl docstring for details) +integrator = "riccati" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic) or "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix) set_psilim_via_dmlim = false # Keep psilim at psihigh (do not truncate at last_rational_q + dmlim) dmlim = 0.2 # Truncate integration at (last_rational_q + dmlim) / n (only used when set_psilim_via_dmlim = true) diff --git a/examples/DIIID-like_ideal_example/gpec.toml b/examples/DIIID-like_ideal_example/gpec.toml index 9a3bf9bc5..068432501 100644 --- a/examples/DIIID-like_ideal_example/gpec.toml +++ b/examples/DIIID-like_ideal_example/gpec.toml @@ -49,10 +49,8 @@ save_interval = 3 # Save every Nth ODE step (1=all). Always saves n singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced ucrit = 1e4 # Column-norm threshold that triggers solution renormalization -# Δ' BVP + parallel integration (see ForceFreeStatesControl docstring for details) -use_parallel = true # Run parallel FM-propagator BVP path (unlocks SingularSurfaces/delta_prime_matrix) -parallel_threads = 2 # BVP thread cap (1 = serial/bit-deterministic; 2 ≈ +20% speedup; ≥3 saturates) -populate_dense_xi = true # Append serial-EL pass so dense ξ is stored — REQUIRED with a [PerturbedEquilibrium] section +# Integrator selection and domain truncation (see ForceFreeStatesControl docstring for details) +integrator = "forward" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic) or "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix) set_psilim_via_dmlim = true # Truncate at (last_rational_q + dmlim)/n — TRUE for diverted equilibria (q → ∞ at separatrix) dmlim = 0.2 # Truncate integration at (last_rational_q + dmlim)/n (used when set_psilim_via_dmlim = true) diff --git a/examples/DIIID-like_riccati_deltaprime_example/gpec.toml b/examples/DIIID-like_riccati_deltaprime_example/gpec.toml new file mode 100644 index 000000000..17c79ff51 --- /dev/null +++ b/examples/DIIID-like_riccati_deltaprime_example/gpec.toml @@ -0,0 +1,57 @@ +# DIII-D-like H-mode equilibrium — n=1 Δ'-matrix case on the Riccati integrator. +# Reuses the equilibrium of the sibling DIIID-like_ideal_example (geqdsk referenced by +# relative path; not duplicated), and is the canonical SingularSurfaces/delta_prime_matrix +# fixture: the Riccati chunked propagator BVP is the only integrator that produces it. +# No forcing or perturbed equilibrium, so no dense ξ profiles are needed. + +[Equilibrium] +eq_filename = "../DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode.geqdsk" # Path to equilibrium file +eq_type = "efit" # Type of the input 2D equilibrium file +jac_type = "hamada" # Coordinate system (hamada, pest, boozer, equal_arc, park, custom) +grid_type = "auto" # Radial grid packing type ("auto" = two-pass measured-curvature grid) +psilow = 1e-4 # Lower limit of normalized poloidal flux +psihigh = 0.995 # Upper limit of normalized poloidal flux (captures q=6 at this Ip) +mpsi = 0 # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy) +psi_accuracy = 0.001 # Target relative accuracy of splined profile derivatives for the auto grid +mtheta = 256 # Number of poloidal grid points +newq0 = 0 # Override for on-axis safety factor (0 = use input value) +etol = 1e-10 # Error tolerance for equilibrium solver +force_termination = false # Terminate after equilibrium setup (skip stability calculations) + +[Wall] +shape = "nowall" # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath) +a = 0.2415 # Distance from plasma (conformal) or shape parameter +aw = 0.05 # Half-thickness parameter for Dee-shaped walls +bw = 1.5 # Elongation parameter for wall shapes +cw = 0 # Offset of wall center from major radius +dw = 0.5 # Triangularity parameter for wall shapes +tw = 0.05 # Sharpness of wall corners (try 0.05 as initial value) +equal_arc_wall = true # Equal arc length distribution of nodes on wall + +[ForceFreeStates] +local_stability_flag = true # Perform local stability analysis (Mercier and ballooning) across the ψ profile +vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes +force_termination = true # Terminate after force-free states (skip perturbed equilibrium calculations) + +psiedge = 0.99 # Edge dW(ψ) diagnostic scan band [psiedge, psilim]; set ≥ psilim to disable +qlow = 1.02 # Integration initiated at q determined by min(q0, qlow) +qhigh = 1e3 # Integration terminated at q limit determined by min(qa, qhigh) +sing_start = 0 # Start integration at the sing_start'th rational from the axis (psilow) + +nn_low = 1 # Smallest toroidal mode number to include +nn_high = 1 # Largest toroidal mode number to include +delta_mlow = 8 # Expands lower bound of Fourier harmonics +delta_mhigh = 8 # Expands upper bound of Fourier harmonics +mthvac = 512 # Number of points used in splines over poloidal angle at the plasma-vacuum interface + +kinetic_source = "fixed" # Kinetic matrix source: "fixed" test matrices, or "calculated" from the kinetic NTV model +kinetic_factor = 0.0 # Scaling of kinetic matrices (0 = ideal path; >0 enables kinetic mode) +eulerlagrange_tolerance = 1e-10 # Relative tolerance for ODE integration of Euler-Lagrange equations +save_interval = 3 # Save every Nth ODE step (1=all). Always saves near rational surfaces. +singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced +ucrit = 1e4 # Column-norm threshold that triggers solution renormalization + +# Integrator selection and domain truncation (see ForceFreeStatesControl docstring for details) +integrator = "riccati" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic) or "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix) +set_psilim_via_dmlim = true # Truncate at (last_rational_q + dmlim)/n — TRUE for diverted equilibria (q → ∞ at separatrix) +dmlim = 0.2 # Truncate integration at (last_rational_q + dmlim)/n (used when set_psilim_via_dmlim = true) diff --git a/examples/LAR_beta_scan/gpec.toml b/examples/LAR_beta_scan/gpec.toml index 8ac649d06..4ab001876 100644 --- a/examples/LAR_beta_scan/gpec.toml +++ b/examples/LAR_beta_scan/gpec.toml @@ -52,9 +52,7 @@ singfac_min = 1e-4 # Fractional distance from rational q at which id ucrit = 1e4 # Column-norm threshold that triggers solution renormalization sing_order = 6 # Order of the singular-surface (Frobenius) series expansion -use_parallel = true # Run parallel FM-propagator BVP path (unlocks SingularSurfaces/delta_prime_matrix) -parallel_threads = 2 # BVP thread cap (1 = serial/bit-deterministic; 2 ≈ +20% speedup; ≥3 saturates) -populate_dense_xi = false # Append serial-EL pass for dense ξ; not needed without [PerturbedEquilibrium] (default false) +integrator = "riccati" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic) or "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix) set_psilim_via_dmlim = false # FALSE for limited/analytical equilibria — rationals sparse, dmlim would chop too much edge dmlim = 0.2 # Truncate integration at (last_rational_q + dmlim)/n (used when set_psilim_via_dmlim = true) force_termination = true # Stop after force-free states (skip perturbed equilibrium) diff --git a/examples/LAR_epsilon_scan/gpec.toml b/examples/LAR_epsilon_scan/gpec.toml index f0d174ecb..2dc5f069b 100644 --- a/examples/LAR_epsilon_scan/gpec.toml +++ b/examples/LAR_epsilon_scan/gpec.toml @@ -53,9 +53,7 @@ singfac_min = 1e-4 # Fractional distance from rational q at which id ucrit = 1e4 # Column-norm threshold that triggers solution renormalization sing_order = 6 # Order of the singular-surface (Frobenius) series expansion -use_parallel = true # Run parallel FM-propagator BVP path (unlocks SingularSurfaces/delta_prime_matrix) -parallel_threads = 2 # BVP thread cap (1 = serial/bit-deterministic; 2 ≈ +20% speedup; ≥3 saturates) -populate_dense_xi = false # Append serial-EL pass for dense ξ; not needed without [PerturbedEquilibrium] (default false) +integrator = "riccati" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic) or "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix) set_psilim_via_dmlim = false # FALSE for limited/analytical equilibria — rationals sparse, dmlim would chop too much edge dmlim = 0.2 # Truncate integration at (last_rational_q + dmlim)/n (used when set_psilim_via_dmlim = true) force_termination = true # Stop after force-free states (skip perturbed equilibrium) diff --git a/examples/LAR_ideal_match_test/gpec.toml b/examples/LAR_ideal_match_test/gpec.toml index d150cbac2..c53d472eb 100644 --- a/examples/LAR_ideal_match_test/gpec.toml +++ b/examples/LAR_ideal_match_test/gpec.toml @@ -50,9 +50,7 @@ ucrit = 1e4 # Maximum fraction of solutions allowed before r sing_order = 6 # Power-series order for the singular-surface asymptotics save_interval = 3 # Save every Nth ODE step (1=all, 10=every 10th). Always saves near rational surfaces. -use_parallel = true # Run parallel FM-propagator BVP path (unlocks SingularSurfaces/delta_prime_matrix) -parallel_threads = 2 # BVP thread cap (1 = serial/bit-deterministic; 2 is about +20% speedup) -populate_dense_xi = false # Dense axis-basis xi for the FFS HDF5 output. Not needed here: the gal-matched path builds its own dense xi. +integrator = "riccati" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic) or "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix) set_psilim_via_dmlim = false # Keep psilim at psihigh (do not truncate at last_rational_q + dmlim) dmlim = 0.2 # Truncate integration at (last_rational_q + dmlim) / n (only used when set_psilim_via_dmlim = true) force_termination = true # Terminate after the stability stage (no perturbed-equilibrium section here) diff --git a/examples/LAR_resistive_match_test/gpec.toml b/examples/LAR_resistive_match_test/gpec.toml index 9900ad8e7..de8d0e5d1 100644 --- a/examples/LAR_resistive_match_test/gpec.toml +++ b/examples/LAR_resistive_match_test/gpec.toml @@ -51,9 +51,7 @@ ucrit = 1e4 # Maximum fraction of solutions allowed before r sing_order = 6 # Power-series order for the singular-surface asymptotics save_interval = 3 # Save every Nth ODE step (1=all, 10=every 10th). Always saves near rational surfaces. -use_parallel = true # Run parallel FM-propagator BVP path (unlocks SingularSurfaces/delta_prime_matrix) -parallel_threads = 2 # BVP thread cap (1 = serial/bit-deterministic; 2 is about +20% speedup) -populate_dense_xi = false # Dense axis-basis xi for the FFS HDF5 output. Not needed here: the gal-matched path builds its own dense xi. +integrator = "riccati" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic) or "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix) set_psilim_via_dmlim = false # Keep psilim at psihigh (do not truncate at last_rational_q + dmlim) dmlim = 0.2 # Truncate integration at (last_rational_q + dmlim) / n (only used when set_psilim_via_dmlim = true) force_termination = true # Terminate after the stability stage (no perturbed-equilibrium section here) diff --git a/examples/Solovev_ideal_example/gpec.toml b/examples/Solovev_ideal_example/gpec.toml index 66ba0d48b..14c192bd9 100644 --- a/examples/Solovev_ideal_example/gpec.toml +++ b/examples/Solovev_ideal_example/gpec.toml @@ -71,10 +71,8 @@ singfac_min = 1e-4 # Fractional distance from rational q at which ide ucrit = 1e3 # Column-norm threshold that triggers solution renormalization save_interval = 3 # Save every Nth ODE step (1=all). Always saves near rational surfaces. -# Δ' BVP + parallel integration (see ForceFreeStatesControl docstring for details) -use_parallel = true # Run parallel FM-propagator BVP path (unlocks SingularSurfaces/delta_prime_matrix) -parallel_threads = 2 # BVP thread cap (1 = serial/bit-deterministic; 2 ≈ +20% speedup; ≥3 saturates) -populate_dense_xi = true # Append serial-EL pass so dense ξ is stored — REQUIRED with a [PerturbedEquilibrium] section +# Integrator selection and domain truncation (see ForceFreeStatesControl docstring for details) +integrator = "forward" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic) or "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix) set_psilim_via_dmlim = false # FALSE for limited/analytical equilibria — rationals sparse, dmlim would chop too much edge dmlim = 0.2 # Truncate integration at (last_rational_q + dmlim)/n (used when set_psilim_via_dmlim = true) # Solovev analytic equilibrium parameters (eq_type = "sol"); see SolovevConfig in src/Equilibrium. diff --git a/examples/Solovev_ideal_example_3D/gpec.toml b/examples/Solovev_ideal_example_3D/gpec.toml index 70123040b..4d9c89086 100644 --- a/examples/Solovev_ideal_example_3D/gpec.toml +++ b/examples/Solovev_ideal_example_3D/gpec.toml @@ -37,10 +37,8 @@ singfac_min = 1e-4 # Fractional distance from rational q at which ide ucrit = 1e3 # Column-norm threshold that triggers solution renormalization save_interval = 3 # Save every Nth ODE step (1=all). Always saves near rational surfaces. -# Δ' BVP + parallel integration (see ForceFreeStatesControl docstring for details) -use_parallel = true # Run parallel FM-propagator BVP path (unlocks SingularSurfaces/delta_prime_matrix) -parallel_threads = 2 # BVP thread cap (1 = serial/bit-deterministic; 2 ≈ +20% speedup; ≥3 saturates) -populate_dense_xi = false # Append serial-EL pass for dense ξ; not needed without [PerturbedEquilibrium] (default false) +# Integrator selection and domain truncation (see ForceFreeStatesControl docstring for details) +integrator = "riccati" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic) or "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix) set_psilim_via_dmlim = false # FALSE for limited/analytical equilibria — rationals sparse, dmlim would chop too much edge dmlim = 0.2 # Truncate integration at (last_rational_q + dmlim)/n (used when set_psilim_via_dmlim = true) diff --git a/examples/Solovev_ideal_example_multi_n/gpec.toml b/examples/Solovev_ideal_example_multi_n/gpec.toml index 3fb3f7b06..5286b5302 100644 --- a/examples/Solovev_ideal_example_multi_n/gpec.toml +++ b/examples/Solovev_ideal_example_multi_n/gpec.toml @@ -47,10 +47,9 @@ eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced ucrit = 1e3 # Column-norm threshold that triggers solution renormalization -# Δ' BVP + parallel integration (see ForceFreeStatesControl docstring for details) -use_parallel = true # Run parallel FM-propagator BVP path (multi-n Δ' matrix has open issues — sing_lim! warns and skips — but ξ and energies are valid) -parallel_threads = 2 # BVP thread cap (1 = serial/bit-deterministic; 2 ≈ +20% speedup; ≥3 saturates) -populate_dense_xi = false # Append serial-EL pass for dense ξ; not needed without [PerturbedEquilibrium] (default false) +# Integrator selection and domain truncation (see ForceFreeStatesControl docstring for details). +# The multi-n Δ' matrix has open issues and is skipped with a warning; the energies stay valid. +integrator = "riccati" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic) or "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix) set_psilim_via_dmlim = false # FALSE for multi-n — dmlim truncation is ambiguous when n varies (sing_lim! skips anyway) dmlim = 0.2 # Truncate integration at (last_rational_q + dmlim)/n (used when set_psilim_via_dmlim = true) # Solovev analytic equilibrium parameters (eq_type = "sol"); see SolovevConfig in src/Equilibrium. diff --git a/examples/Solovev_kinetic_NTV_example/gpec.toml b/examples/Solovev_kinetic_NTV_example/gpec.toml index 733a92c31..8d3f90d3f 100644 --- a/examples/Solovev_kinetic_NTV_example/gpec.toml +++ b/examples/Solovev_kinetic_NTV_example/gpec.toml @@ -68,10 +68,8 @@ singfac_min = 1e-4 # Fractional distance from rational q at which ide ucrit = 1e3 # Column-norm threshold that triggers solution renormalization save_interval = 3 # Save every Nth ODE step (1=all). Always saves near rational surfaces. -# Δ' BVP + parallel integration (see ForceFreeStatesControl docstring for details) -use_parallel = true # Run parallel FM-propagator BVP path (unlocks SingularSurfaces/delta_prime_matrix) -parallel_threads = 2 # BVP thread cap (1 = serial/bit-deterministic; 2 ≈ +20% speedup; ≥3 saturates) -populate_dense_xi = true # Append serial-EL pass so dense ξ is stored — REQUIRED with a [PerturbedEquilibrium] section +# Integrator selection and domain truncation (see ForceFreeStatesControl docstring for details) +integrator = "forward" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic) or "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix) set_psilim_via_dmlim = false # FALSE for limited/analytical equilibria — rationals sparse, dmlim would chop too much edge dmlim = 0.2 # Truncate integration at (last_rational_q + dmlim)/n (used when set_psilim_via_dmlim = true) # Solovev analytic equilibrium parameters (eq_type = "sol"); see SolovevConfig in src/Equilibrium. diff --git a/examples/Solovev_kinetic_calculated_example/gpec.toml b/examples/Solovev_kinetic_calculated_example/gpec.toml index 149f40af7..760ba8e7d 100644 --- a/examples/Solovev_kinetic_calculated_example/gpec.toml +++ b/examples/Solovev_kinetic_calculated_example/gpec.toml @@ -42,6 +42,7 @@ mthvac = 64 # Number of points used in splines over poloidal a kinetic_source = "calculated" # Kinetic matrix source — exercises KineticForces.compute_calculated_kinetic_matrices callback with real physics kinetic_factor = 1.0 # Full-strength kinetic matrices (the "calculated" path is the real physics; no perturbation scaling) +integrator = "forward" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic) or "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix) eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced ucrit = 1e3 # Maximum fraction of solutions allowed before re-normalized diff --git a/examples/a10_kinetic_example/gpec.toml b/examples/a10_kinetic_example/gpec.toml index d3f61b391..c6a6e9639 100644 --- a/examples/a10_kinetic_example/gpec.toml +++ b/examples/a10_kinetic_example/gpec.toml @@ -39,6 +39,7 @@ mthvac = 512 # Number of points used in splines over poloidal kinetic_source = "calculated" # Kinetic matrix source: "fixed" test matrices, or "calculated" from the kinetic NTV model kinetic_factor = 1.0 # Scaling of kinetic matrices (0 = ideal path; >0 enables kinetic mode; 1.0 = full strength) +integrator = "forward" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic) or "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix) eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced ucrit = 1e4 # Column-norm threshold that triggers solution renormalization diff --git a/regression-harness/cases/diiid_n1.toml b/regression-harness/cases/diiid_n1.toml index 6b1161ab6..da546a115 100644 --- a/regression-harness/cases/diiid_n1.toml +++ b/regression-harness/cases/diiid_n1.toml @@ -1,8 +1,10 @@ # Regression case: DIII-D-like n=1 ideal stability + perturbed equilibrium. -# Tracks energies, singular-surface metrics, Δ' (BVP diagonal), and resonant-flux -# diagnostics for the canonical DIIID-like fixture. Each [quantities.*] block names an -# HDF5 path in the run output, how to extract it, and the noise floor below which a -# difference is treated as zero. Run via the regression-harness (see CLAUDE.md). +# Tracks energies, singular-surface metrics, and resonant-flux diagnostics for the canonical +# DIIID-like fixture, which runs the forward integrator for its dense ξ profiles. The BVP Δ' +# matrix comes from the Riccati integrator only and is tracked by the sibling case +# diiid_n1_riccati. Each [quantities.*] block names an HDF5 path in the run output, how to +# extract it, and the noise floor below which a difference is treated as zero. Run via the +# regression-harness (see CLAUDE.md). [case] name = "diiid_n1" description = "DIII-D-like equilibrium, n=1, ideal + perturbed equilibrium" @@ -270,18 +272,6 @@ label = "npert" noise_threshold = 0 order = 61 -# Tearing stability Δ' — canonical STRIDE BVP matrix diagonal (replaces the -# previous `PerturbedEquilibrium/SingularCoupling/delta_prime` track, which -# was a per-surface stub computed by SingularCoupling from (rbwp1-lbwp1)/(2π·χ'). -# Per-surface Δ' is de-emphasized; SingularCoupling reads this BVP matrix diagonal. -[quantities.delta_prime] -h5path = "SingularSurfaces/delta_prime_matrix" -type = "complex_matrix" -extract = "diagonal_complex" -label = "delta prime (BVP diagonal)" -noise_threshold = 1e-8 -order = 80 - [quantities.island_half_width] h5path = "PerturbedEquilibrium/SingularCoupling/island_half_width" type = "real_vector" diff --git a/regression-harness/cases/diiid_n1_riccati.toml b/regression-harness/cases/diiid_n1_riccati.toml new file mode 100644 index 000000000..a5ae92cc3 --- /dev/null +++ b/regression-harness/cases/diiid_n1_riccati.toml @@ -0,0 +1,163 @@ +# Regression case: DIII-D-like n=1 Δ' matrix from the Riccati chunked propagator BVP. +# The Riccati integrator is the only one that produces SingularSurfaces/delta_prime_matrix, so +# this case is the canonical Δ' track (its sibling diiid_n1 runs the forward integrator and +# carries the ξ/perturbed-equilibrium quantities instead). Also tracks the free-boundary +# energies and singular-surface locations the BVP depends on, so a Δ' move can be attributed. +# Each [quantities.*] block names an HDF5 path in the run output, how to extract it, and the +# noise floor below which a difference is treated as zero. +[case] +name = "diiid_n1_riccati" +description = "DIII-D-like equilibrium, n=1, Riccati integrator Δ' matrix" +example_dir = "examples/DIIID-like_riccati_deltaprime_example" + +# Tearing stability Δ' — the msing×msing BVP matrix diagonal, and the raw side-major +# outer-region matrix it is projected from. +[quantities.delta_prime] +h5path = "SingularSurfaces/delta_prime_matrix" +type = "complex_matrix" +extract = "diagonal_complex" +label = "delta prime (BVP diagonal)" +noise_threshold = 1e-8 +order = 10 + +[quantities.delta_prime_raw] +h5path = "SingularSurfaces/delta_prime_raw" +type = "complex_matrix" +extract = "all_complex" +label = "delta prime (raw side-major)" +noise_threshold = 1e-8 +order = 11 + +[quantities.delta_coil] +h5path = "SingularSurfaces/delta_coil" +type = "complex_matrix" +extract = "all_complex" +label = "edge coil response delta_coil" +noise_threshold = 1e-8 +order = 12 + +# Energies — leading generalized (W,N) pencil eigenvalues at the final truncation (psilim): power-normalized (⟨|ξ|²⟩ = 1 metric) and invariant to the working-coordinate Jacobian. +[quantities.et_real] +h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_energies" +type = "complex_vector" +extract = "real_first" +label = "total energy Re(et[1])" +noise_threshold = 1e-10 +order = 20 + +[quantities.ep_real] +h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_plasma_energies" +type = "complex_vector" +extract = "real_first" +label = "plasma energy Re(ep[1])" +noise_threshold = 1e-10 +order = 21 + +[quantities.ev_real] +h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_vacuum_energies" +type = "complex_vector" +extract = "real_first" +label = "vacuum energy Re(ev[1])" +noise_threshold = 1e-10 +order = 22 + +[quantities.et_all] +h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_energies" +type = "complex_vector" +extract = "all_complex" +label = "total energy (all)" +noise_threshold = 1e-10 +order = 23 + +# Singular surfaces — the Δ' matrix is indexed by these, so a shift here explains a Δ' shift. +[quantities.msing] +h5path = "SingularSurfaces/msing" +type = "int_scalar" +extract = "value" +label = "# singular surfaces" +noise_threshold = 0 +order = 30 + +[quantities.sing_psi] +h5path = "SingularSurfaces/psi" +type = "real_vector" +extract = "all_real" +label = "singular psi locations" +noise_threshold = 1e-8 +order = 31 + +[quantities.sing_q] +h5path = "SingularSurfaces/q" +type = "real_vector" +extract = "all_real" +label = "singular q values" +noise_threshold = 1e-8 +order = 32 + +# Integration — chunk-endpoint step counts move when the chunk sizing changes. +[quantities.nstep] +h5path = "ForceFreeStates/Solutions/ForwardIntegration/nstep" +type = "int_scalar" +extract = "value" +label = "ODE steps (saved)" +noise_threshold = 0 +order = 40 + +[quantities.nstep_total] +h5path = "ForceFreeStates/Solutions/ForwardIntegration/nstep_total" +type = "int_scalar" +extract = "value" +label = "ODE steps (total)" +noise_threshold = 0 +order = 41 + +# Mode numbers +[quantities.mpert] +h5path = "Info/mpert" +type = "int_scalar" +extract = "value" +label = "mpert" +noise_threshold = 0 +order = 50 + +[quantities.npert] +h5path = "Info/npert" +type = "int_scalar" +extract = "value" +label = "npert" +noise_threshold = 0 +order = 51 + +# Equilibrium — shared with diiid_n1; confirms the two cases sit on the same equilibrium. +[quantities.q0] +h5path = "Equilibrium/q0" +type = "real_scalar" +extract = "value" +label = "q0" +noise_threshold = 1e-12 +order = 60 + +[quantities.q95] +h5path = "Equilibrium/q95" +type = "real_scalar" +extract = "value" +label = "q95" +noise_threshold = 1e-12 +order = 61 + +[quantities.betan] +h5path = "Equilibrium/betan" +type = "real_scalar" +extract = "value" +label = "beta_n" +noise_threshold = 1e-12 +order = 62 + +# Runtime +[quantities.runtime] +h5path = "" +type = "runtime" +extract = "value" +label = "Runtime (s)" +noise_threshold = 0.0 +order = 999 diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl index f563dca45..8d0d024f6 100644 --- a/src/ForceFreeStates/EulerLagrange.jl +++ b/src/ForceFreeStates/EulerLagrange.jl @@ -64,10 +64,12 @@ end """ balance_integration_chunks(chunks, ctrl, intr) -> Vector{IntegrationChunk} -Sub-divide integration chunks to produce a load-balanced set for parallel execution. +Sub-divide integration chunks to produce a load-balanced set for the Riccati BVP. Starts from the output of `chunk_el_integration_bounds` and iteratively splits the -highest-cost chunk (by `ode_itime_cost`) until the total chunk count reaches -`max(2*msing + 3, 4 * Threads.nthreads())`. +highest-cost chunk (by `ode_itime_cost`) until the total chunk count reaches the target +set by `ctrl.nchunks` (`0` = auto). The target is derived from problem structure only — +never from `Threads.nthreads()` — so the chunk list, and hence every Riccati output, is +identical whatever thread count `julia -t` provides. Each split finds the equal-cost midpoint ψ_mid via bisection: ode_itime_cost(psi_start, psi_mid) ≈ ode_itime_cost(psi_start, psi_end) / 2 @@ -84,11 +86,14 @@ function balance_integration_chunks(chunks::Vector{IntegrationChunk}, ctrl::Forc # assemble_fm_matrix(condition=true) can't keep accumulated products well-conditioned # because single long-span propagators may already have cond ~ 10²⁴. min_bvp_intervals = 8 * (intr.msing + 1) + intr.msing - # Use the effective parallel width (capped by ctrl.parallel_threads) rather than - # Threads.nthreads() — otherwise a user on `julia -t 16` who sets parallel_threads=2 - # for determinism still pays for 4× the requested sub-chunk count. - effective_threads = min(Threads.nthreads(), max(ctrl.parallel_threads, 1)) - target_n = max(min_chunks, 4 * effective_threads, min_bvp_intervals) + if ctrl.nchunks > 0 + if ctrl.nchunks < min_chunks + @warn "nchunks = $(ctrl.nchunks) is below the $min_chunks chunks required by $(intr.msing) singular surfaces; clamping up." + end + target_n = max(ctrl.nchunks, min_chunks) + else + target_n = max(min_chunks, min_bvp_intervals) + end result = collect(chunks) @@ -140,36 +145,38 @@ end eulerlagrange_integration(ctrl, equil, ffit, intr) -> (odet, propagators, chunks, S_left) Integrate the Euler-Lagrange equations from the axis to `intr.psilim`, crossing each singular -surface on the way (Fortran `ode_run`). Dispatches on `ctrl` to the parallel propagator BVP -(`use_parallel`), the dual Riccati formulation (`use_riccati`), or -[`serial_eulerlagrange_integration`](@ref). +surface on the way (Fortran `ode_run`). Dispatches on `ctrl.integrator` to +[`riccati_eulerlagrange_integration`](@ref) (the chunked propagator BVP) or +[`forward_eulerlagrange_integration`](@ref). -Only the parallel branch populates `propagators` / `chunks` / `S_left`, which -`compute_delta_prime_matrix!` consumes for the Δ' BVP; the other two return `nothing` for all -three. +Only the Riccati branch populates `propagators` / `chunks` / `S_left`, which +`compute_delta_prime_matrix!` consumes for the Δ' BVP; the forward branch returns `nothing` +for all three. """ function eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal) - # Dispatch to parallel or Riccati solver if requested. - # Parallel path returns (odet, propagators, chunks, S_at_surface_left) for deferred Δ' BVP. - if ctrl.use_parallel - return parallel_eulerlagrange_integration(ctrl, equil, ffit, intr) - elseif ctrl.use_riccati - return (riccati_eulerlagrange_integration(ctrl, equil, ffit, intr), nothing, nothing, nothing) + if ctrl.integrator == "riccati" + ctrl.kinetic_factor > 0 && error("kinetic runs require integrator=\"forward\"; the Riccati integrator has no kinetic crossing.") + # Riccati path returns (odet, propagators, chunks, S_at_surface_left) for the deferred Δ' BVP. + return riccati_eulerlagrange_integration(ctrl, equil, ffit, intr) + elseif ctrl.integrator == "forward" + return forward_eulerlagrange_integration(ctrl, equil, ffit, intr) + elseif ctrl.integrator == "galerkin" + error("integrator = \"galerkin\" is not yet a standalone integrator — use gal_flag = true alongside integrator = \"forward\" or \"riccati\".") end - return serial_eulerlagrange_integration(ctrl, equil, ffit, intr) + error("Unknown integrator: $(ctrl.integrator). Expected \"forward\", \"riccati\", or \"galerkin\".") end """ - serial_eulerlagrange_integration(ctrl, equil, ffit, intr; verbose=ctrl.verbose) -> (odet, nothing, nothing, nothing) + forward_eulerlagrange_integration(ctrl, equil, ffit, intr; verbose=ctrl.verbose) -> (odet, nothing, nothing, nothing) -Serial shooting branch of [`eulerlagrange_integration`](@ref): integrates chunk by chunk, +Forward branch of [`eulerlagrange_integration`](@ref): integrates chunk by chunk from the axis, applying Gaussian reduction whenever a solution norm ratio exceeds `ctrl.ucrit` and undoing it via `transform_u!` at the end, so `odet.u_store` comes back dense in the axis basis. Call -directly to force this branch regardless of `ctrl.use_parallel` / `ctrl.use_riccati`; `verbose` -overrides `ctrl.verbose` for progress logging. +directly to force this branch regardless of `ctrl.integrator`; `verbose` overrides +`ctrl.verbose` for progress logging. """ -function serial_eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal; +function forward_eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal; verbose::Bool=ctrl.verbose) # Initialization diff --git a/src/ForceFreeStates/ForceFreeStatesStructs.jl b/src/ForceFreeStates/ForceFreeStatesStructs.jl index 00eaf0ba7..45843bbfa 100644 --- a/src/ForceFreeStates/ForceFreeStatesStructs.jl +++ b/src/ForceFreeStates/ForceFreeStatesStructs.jl @@ -202,7 +202,7 @@ A mutable struct holding internal state variables for stability calculations. via `pest3_decompose(dp_raw)` — needed for the full det(D' − D(γ)) = 0 eigenvalue problem with Glasser stabilization. - Empty unless `ctrl.use_parallel` is true. No ½ prefactor is applied (matches + Empty unless the Riccati integrator was used. No ½ prefactor is applied (matches Fortran rdcon; Pletzer–Dewar paper multiplies by ½). """ delta_prime_raw::Matrix{ComplexF64} = Matrix{ComplexF64}(undef, 0, 0) @@ -249,10 +249,8 @@ gpec.toml. - `HDF5_filename::String` - Name of HDF5 output file - `save_interval::Int` - Save every Nth ODE step (1=all, 10=every 10th). Always saves near rational surfaces. (Same as `euler_step` in the Fortran) - `force_termination::Bool` - Terminate after force-free states (skip perturbed equilibrium calculations) - - `use_riccati::Bool` - Use the dual Riccati reformulation S = U₁·U₂⁻¹ instead of the standard U₁/U₂ ODE. Reduces stiffness for faster integration. See Glasser (2018) Phys. Plasmas 25, 032507. - - `use_parallel::Bool` - Parallel fundamental matrix (propagator) integration using `Threads.@threads`. Each chunk is integrated independently from identity IC and assembled serially. Requires `singfac_min != 0`. Uses the same chunk bounds as the standard path but sub-divides chunks for load balancing. Crossings use the Riccati-style algorithm (no Gaussian reduction). - - `parallel_threads::Int` - Cap on the number of threads the parallel BVP uses. **Default `2`** parallelises the FM chunks across two threads (the BVP has ~10 chunks; 2 threads is enough to amortize them — speedup saturates here, raising to 4 adds scheduling overhead). Set `parallel_threads = 1` to run the FM chunks SERIALLY (no `Threads.@threads`), which is bit-deterministic and immune to the thread-schedule sensitivity that can cause intermittent BVP divergence on numerically delicate equilibria. The parallel path produces bit-identical Δ′ across thread counts; `parallel_threads = 2` is about 20% faster than serial and saturates the speedup. If a parallel run diverges, drop to `parallel_threads = 1` rather than switching `use_parallel = false` — the latter is silently wrong. Capped at `Threads.nthreads()`. - - `populate_dense_xi::Bool` - When `use_parallel = true`, append a serial Euler-Lagrange pass after the propagator BVP so the returned `odet` carries dense axis-basis `u_store`/`du_store`/`xi_s_store` — the only convention PerturbedEquilibrium / FieldReconstruction consume correctly, and what fills HDF5 `ForceFreeStates/Solutions/ForwardIntegration/xi_*`. Δ' (`SingularSurfaces/delta_prime_matrix`) comes from the BVP and is identical either way; free-boundary energies (`ForceFreeStates/FreeBoundaryStability/eigenmode_*`) match a pure serial run when enabled (~0.12 % Riccati-vs-axis gap otherwise). **Default `false`** (skips the extra serial pass, ~1× BVP wall-clock); **PerturbedEquilibrium-using configs must set it `true`** when `use_parallel = true`, else PE reads Riccati-basis garbage. Auto-disabled when `force_termination = true`. + - `integrator::String` - Which formalism integrates the Euler-Lagrange system. `"forward"` sweeps the plasma serially with Gaussian reduction and returns `u_store` / `du_store` / `xi_s_store` dense in the axis (EL) basis — the only convention PerturbedEquilibrium and FieldReconstruction consume correctly, and the only path that supports `kinetic_factor > 0`. `"riccati"` (default) runs the chunked fundamental-matrix propagator driver (Glasser 2018 Phys. Plasmas 25, 032507): chunks are integrated independently from identity initial conditions and assembled serially with Riccati-style crossings, which is the only way to obtain the singular-surface Δ' matrix for the tearing-mode solvers downstream, but leaves `u_store` as sparse chunk-endpoint Riccati states, so dense ξ profiles are unavailable. `"galerkin"` is not yet a standalone integrator and currently errors — use `gal_flag = true` alongside another integrator. Requires `singfac_min != 0` for `"riccati"`. + - `nchunks::Int` - Target number of Riccati integration chunks. `0` (the default) derives the count from problem structure alone: `max(2·msing + 3, 8·(msing + 1) + msing)`, enough sub-chunks per segment to keep the accumulated propagator products well-conditioned. An explicit value below `2·msing + 3` is clamped up with a warning. Chunk sizing never consults `Threads.nthreads()`, so Riccati outputs are identical whatever thread count `julia -t` provides; threads only change wall-clock. - `extended_precision_bvp::Bool` - When `true` (default), promote the Δ' BVP linear system to `Complex{Double64}` (~31 digits) for the LU solve and PEST3 combination. Guards against catastrophic cancellation in the PEST3 four-term combination (dp_raw entries can be 10⁴–10⁵× larger than the result; the imaginary part of off-diagonal Δ' is particularly sensitive). Disabling (`false`) saves ~1.5–2× the BVP solve time but on DIIID-class equilibria the imaginary Δ' components can drift by factors of 2–5×; only disable for performance experiments on cases where Float64 has been validated against Double64. """ @kwdef struct ForceFreeStatesControl @@ -272,7 +270,7 @@ gpec.toml. ucrit::Float64 = 1e4 numsteps_init::Int = 4000 numunorms_init::Int = 100 - singfac_min::Float64 = 1e-4 # Matches Fortran STRIDE; required nonzero for use_parallel path. + singfac_min::Float64 = 1e-4 # Matches Fortran STRIDE; required nonzero for the Riccati path. set_psilim_via_dmlim::Bool = true # Safe default for diverted equilibria (most production use); set false for limited/analytical (LAR, Solovev). Auto-skipped for multi-n. See docstring. dmlim::Float64 = 0.2 sing_order::Int = 6 @@ -283,16 +281,14 @@ gpec.toml. reform_eq_with_psilim::Bool = false psiedge::Float64 = 0.99 truncate_at_dW_peak::Bool = false # Edge-dW peak becomes new physical edge; Δ' BVP made self-consistent. See docstring. - parallel_threads::Int = 2 diagnose::Bool = false diagnose_ca::Bool = false write_outputs_to_HDF5::Bool = true HDF5_filename::String = "gpec.h5" save_interval::Int = 3 force_termination::Bool = false - use_riccati::Bool = false - use_parallel::Bool = true # Default on: unlocks SingularSurfaces/delta_prime_matrix (STRIDE BVP Δ′ matrix) used by SLAYER/GGJ downstream. - populate_dense_xi::Bool = false # When use_parallel=true, set to true ONLY if a PerturbedEquilibrium pipeline will consume dense ξ. Default false avoids the ~1× parallel-BVP serial-EL re-run for non-PE runs (Δ'/vacuum/ideal-stability only). See ForceFreeStatesControl docstring for the full trade-off (et[1] convention differs by ~0.12% on DIIID between populate=true vs false). + integrator::String = "riccati" # Default: unlocks SingularSurfaces/delta_prime_matrix (STRIDE BVP Δ′ matrix) used by SLAYER/GGJ downstream. Use "forward" for dense ξ (PerturbedEquilibrium) or kinetic runs. + nchunks::Int = 0 # Riccati chunk-count target; 0 = auto (derived from msing alone, never from Threads.nthreads()). extended_precision_bvp::Bool = true # Promote Δ' BVP to Complex{Double64}; default on (Float64 drifts the imaginary Δ' by 2–5× on DIIID-class cases). # --- RDCON outer-region Galerkin Δ′ solver (gal_solve port) --- diff --git a/src/ForceFreeStates/Galerkin/GalerkinMatch.jl b/src/ForceFreeStates/Galerkin/GalerkinMatch.jl index cecbb6330..d9c7a4a81 100644 --- a/src/ForceFreeStates/Galerkin/GalerkinMatch.jl +++ b/src/ForceFreeStates/Galerkin/GalerkinMatch.jl @@ -237,7 +237,7 @@ end """ gal_matched_odestate(gal_result, ffit, intr) -> OdeState -Pack the RPEC-matched outer solution into an `OdeState` shaped exactly like the shooting integrator's, +Pack the RPEC-matched outer solution into an `OdeState` shaped exactly like the forward integrator's, so `PerturbedEquilibrium` consumes it unchanged. Mirrors Fortran `idcon_build`'s gal branch (idcon.f) and `globalsol.bin` (the on-surface `issing` points are dropped, match.f): diff --git a/src/ForceFreeStates/Riccati.jl b/src/ForceFreeStates/Riccati.jl index 812bd5c6a..c1710cc12 100644 --- a/src/ForceFreeStates/Riccati.jl +++ b/src/ForceFreeStates/Riccati.jl @@ -1292,110 +1292,6 @@ function _store_crossing_step!(odet::OdeState) store_ode_data!(odet, odet.psifac, odet.u) end -""" - riccati_eulerlagrange_integration(ctrl, equil, ffit, intr) -> OdeState - -Integrate the dual Riccati ODE S = U₁·U₂⁻¹ across the plasma (Glasser 2018 Phys. Plasmas 25, -032507). Reduces stiffness relative to [`serial_eulerlagrange_integration`](@ref), which it -otherwise mirrors, differing in three places: - -1. `riccati_integrate_chunk!` drives `sing_der!` with `riccati_integrator_callback!`, which - applies `renormalize_riccati_inplace!` rather than Gaussian reduction when column norms - exceed `ctrl.ucrit` -2. `riccati_cross_ideal_singular_surf!` replaces `cross_ideal_singular_surf!`: it skips - Gaussian reduction (avoiding near-zero pivots where S is small near the axis) and - renormalizes to (S_new, I) in one step -3. `transform_u!` is skipped — S is already the true solution, so there is no reduction to undo - -Enable via `use_riccati = true` in the `[ForceFreeStates]` section of gpec.toml. -""" -function riccati_eulerlagrange_integration( - ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, - ffit::FourFitVars, intr::ForceFreeStatesInternal -) - # Initialization — same as eulerlagrange_integration - odet = OdeState(intr.numpert_total, ctrl.numsteps_init, ctrl.numunorms_init, intr.msing) - if ctrl.sing_start <= 0 - initialize_el_at_axis!(odet, ctrl, ffit, equil.profiles, intr) - elseif ctrl.sing_start <= intr.msing - error("sing_start > 0 not implemented yet!") - else - error("Invalid value for sing_start: $(ctrl.sing_start) > msing = $(intr.msing)") - end - - chunks = chunk_el_integration_bounds(odet, ctrl, intr) - - # Prime odet.new = false so that compute_solution_norms! (if called elsewhere) - # does not skip Gaussian reduction on first invocation. Also initialize unorm0 - # to safe defaults since the Riccati callback never calls compute_solution_norms!. - odet.new = false - fill!(odet.unorm0, 1.0) - - if ctrl.verbose - @info " ψ = $((@sprintf "%.3f" odet.psifac)), q = $((@sprintf "%.3f" equil.profiles.q_spline(odet.psifac)))" - end - - for chunk in chunks - # Integrate this chunk using the Riccati ODE (Riccati callback skips Gaussian reduction) - riccati_integrate_chunk!(odet, ctrl, equil, ffit, intr, chunk) - if ctrl.verbose - @info " ψ = $((@sprintf "%.3f" odet.psifac)), q= $((@sprintf "%.3f" odet.q)), max(S) = $((@sprintf "%.2e" maximum(abs, odet.u[:,:,1]))), steps = $(odet.step-1)" - end - - # Cross rational surface (Riccati crossing skips GR, uses ipert_res directly) - if chunk.needs_crossing - if ctrl.kinetic_factor > 0 - error("kinetic_factor > 0 not implemented yet in Riccati!") - else - riccati_cross_ideal_singular_surf!(odet, ctrl, equil, ffit, intr, chunk.ising) - # renormalize_riccati! is called inside riccati_cross_ideal_singular_surf! - end - end - end - - # Edge-dW scan over [psiedge, psilim] — populates odet.edge_scan for HDF5 output. - # See EulerLagrange.jl counterpart and ForceFreeStatesControl docstring for the - # diagnostic vs legacy-truncation semantics and reliability caveats on - # truncate_at_dW_peak=true. - odet.step -= 1 - trim_storage!(odet) - if ctrl.psiedge < intr.psilim - saved_psifac, saved_u = odet.psifac, copy(odet.u) - peak_step = findmax_dW_edge!(odet, ctrl, equil, ffit, intr) - if ctrl.truncate_at_dW_peak - # Legacy: truncate integration data to dW peak (corrupts Δ' and δW). - odet.step = peak_step - trim_storage!(odet) - intr.psilim = odet.psi_store[end] - intr.qlim = odet.q_store[end] - odet.u .= odet.u_store[:, :, :, end] - if ctrl.verbose - @info "Truncating integration at peak edge dW (LEGACY — Δ'/δW unreliable): ψ = $((@sprintf "%.2f" odet.psi_store[odet.step])), q = $((@sprintf "%.2f" odet.q_store[odet.step]))" - end - else - odet.psifac = saved_psifac - odet.u .= saved_u - if ctrl.verbose - @info "Edge-dW peak (diagnostic): ψ = $((@sprintf "%.2f" odet.psi_store[peak_step])), q = $((@sprintf "%.2f" odet.q_store[peak_step])); integration domain unchanged" - end - end - end - - # Evaluate fixed-boundary stability criterion - if ctrl.verbose - @info "Evaluating fixed-boundary stability criterion" - end - odet.nzero = evaluate_stability_criterion!(odet, equil.profiles) - - # Note: transform_u! is intentionally skipped. - # S is already the true solution (invariant under Gaussian reduction), - # and u_store entries have u[:,:,1]=S, u[:,:,2]=I throughout integration. - # At crossing steps, u_store has U₁_new/U₂_new which compute_smallest_eigenvalue - # correctly resolves to S_new via rdiv. No transformation is needed. - - return odet -end - """ integrate_propagator_chunk!(prop, chunk, ctrl, equil, ffit, intr, odet_proxy) @@ -1599,41 +1495,42 @@ function apply_propagator_inverse!(odet::OdeState, prop::ChunkPropagator) end """ - parallel_eulerlagrange_integration(ctrl, equil, ffit, intr) -> (odet, propagators, chunks, S_left) + riccati_eulerlagrange_integration(ctrl, equil, ffit, intr) -> (odet, propagators, chunks, S_left) -Parallel fundamental matrix (propagator) driver for the EL integration. The trailing three -return values feed the Δ' BVP in `compute_delta_prime_matrix!`; this is the only branch that -produces them. +The Riccati/STRIDE integrator: a chunked fundamental matrix (propagator) driver for the EL +integration. The trailing three return values feed the Δ' BVP in `compute_delta_prime_matrix!`; +this is the only branch that produces them. -Equivalent to [`serial_eulerlagrange_integration`](@ref), but integrates all bulk chunks -concurrently using `Threads.@threads`, then re-integrates the outer plasma serially: +Solves the same system as [`forward_eulerlagrange_integration`](@ref), but integrates all bulk +chunks concurrently using `Threads.@threads`, then re-integrates the outer plasma serially: 1. **Chunk generation**: calls `chunk_el_integration_bounds`, then `balance_integration_chunks` - to sub-divide chunks for load-balanced parallel execution. -2. **Parallel phase**: `integrate_propagator_chunk!` integrates each chunk independently + to sub-divide chunks for load balancing. The chunk count depends only on `intr.msing` and + `ctrl.nchunks`, never on the thread count, so results are thread-independent. +2. **Propagator phase**: `integrate_propagator_chunk!` integrates each chunk independently from identity initial conditions (no accumulated state, no normalization/callback). Each thread uses a private `OdeState` proxy for `sing_der!` side effects. 3. **Serial assembly**: propagators are applied sequentially with `apply_propagator!`. Rational surface crossings use `riccati_cross_ideal_singular_surf!` (no Gaussian - reduction) matching the Riccati path convention. + reduction). 4. **Outer plasma re-integration**: after the last rational surface crossing, the outer plasma (from last ψ_s to psilim) is re-integrated using `riccati_integrate_chunk!`. FM propagation in this region is prone to precision loss for high N (exponential growth without renormalization); Riccati integration keeps matrices bounded and provides dense checkpoints for `findmax_dW_edge!`. -Enable via `use_parallel = true` in `[ForceFreeStates]` of gpec.toml. Requires `singfac_min != 0`. +Select via `integrator = "riccati"` in `[ForceFreeStates]` of gpec.toml. Requires +`singfac_min != 0`. Uses whatever threads `julia -t` provides; `ctrl.nchunks` is the only +tunable. -**Key differences from serial integration:** +**Key differences from the forward integrator:** - No Gaussian reduction in the propagator BVP phase (crossings use the - Riccati-style algorithm, parallel `odet.ifix` stays 0) -- `transform_u!` is called on the parallel odet but is a no-op (ifix=0) + Riccati-style algorithm, `odet.ifix` stays 0) +- `transform_u!` is called on the odet but is a no-op (ifix=0) - Outer plasma uses serial Riccati integration for numerical stability -- When `ctrl.populate_dense_xi` is set, a serial EL dense pass is appended and replaces the - parallel `odet`, so `u_store` / `du_store` / `xi_s_store` come back in the axis basis that - PerturbedEquilibrium requires. Δ' is computed from the parallel BVP either way and is - bit-identical between the two. See the `populate_dense_xi` entry in the - [`ForceFreeStatesControl`](@ref) docstring for the cost trade-off. +- `odet.u_store` holds chunk-endpoint Riccati states, not dense Euler-Lagrange ξ, and + `odet.u_store_el_basis` stays `false`: this integrator never claims the EL basis, so + PerturbedEquilibrium and the HDF5 forward-integration ξ datasets require the forward path. **Bidirectional integration for large-N accuracy:** The crossing chunk (nearest to each rational surface singL[j]) is integrated *backward* @@ -1644,16 +1541,15 @@ LU solve in `apply_propagator_inverse!`. This follows the same principle as STRI (Glasser 2018 Phys. Plasmas 25, 032501). The all-forward path had ~10% energy error for the DIIID-like example (N=26, n=1); bidirectional reduces this to within 2%. """ -function parallel_eulerlagrange_integration( +function riccati_eulerlagrange_integration( ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal ) odet = _initialize_parallel_odet(ctrl, equil, ffit, intr) chunks, propagators, odet_proxies = _setup_parallel_chunks_and_proxies(odet, ctrl, intr) - bvp_threads = max(1, min(Threads.nthreads(), ctrl.parallel_threads)) - _log_parallel_start(ctrl, odet, equil, chunks, bvp_threads) + _log_parallel_start(ctrl, odet, equil, chunks) - _run_parallel_bvp_phase!(propagators, chunks, ctrl, equil, ffit, intr, odet_proxies, bvp_threads) + _run_parallel_bvp_phase!(propagators, chunks, ctrl, equil, ffit, intr, odet_proxies) # Harvest solver-step counts accumulated thread-locally in each proxy during the BVP phase. # The outer re-integration below uses riccati_integrate_chunk!, which counts via its callback. @@ -1676,11 +1572,6 @@ function parallel_eulerlagrange_integration( odet.nzero = evaluate_stability_criterion!(odet, equil.profiles) transform_u!(odet, intr) # no-op when ifix=0 (no Gaussian reduction) - # Replace BVP `odet` with a dense serial-EL pass so HDF5 `ForwardIntegration/xi_*` carries - # valid DCON ξ in axis basis for PerturbedEquilibrium. Skipped when force_termination=true. - if ctrl.populate_dense_xi && !ctrl.force_termination - odet = _populate_dense_xi_via_serial_el!(odet, ctrl, equil, ffit, intr) - end return odet, propagators, chunks, S_at_surface_left end @@ -1722,33 +1613,25 @@ end function _log_parallel_start(ctrl::ForceFreeStatesControl, odet::OdeState, equil::Equilibrium.PlasmaEquilibrium, - chunks::Vector{IntegrationChunk}, bvp_threads::Int) + chunks::Vector{IntegrationChunk}) ctrl.verbose || return @info " ψ = $((@sprintf "%.3f" odet.psifac)), q = $((@sprintf "%.3f" equil.profiles.q_spline(odet.psifac)))" - @info " Parallel FM: $(length(chunks)) chunks, $bvp_threads BVP thread$(bvp_threads == 1 ? "" : "s") (julia_nthreads=$(Threads.nthreads()), ctrl.parallel_threads=$(ctrl.parallel_threads))" + @info " Riccati FM: $(length(chunks)) chunks over $(Threads.nthreads()) thread$(Threads.nthreads() == 1 ? "" : "s")" end -# Integrate each chunk's FM propagator from identity IC. Serial when bvp_threads == 1 -# (bit-deterministic; ~20% slower than 2-thread but immune to thread- -# schedule sensitivity). Parallel uses :static scheduler so Threads.threadid() returns a -# stable index into odet_proxies. If a parallel run ever diverges on a delicate equilibrium, -# drop to parallel_threads = 1 rather than use_parallel = false — the latter is silently wrong. +# Integrate each chunk's FM propagator from identity IC across whatever threads `julia -t` +# provides. The :static scheduler makes Threads.threadid() a stable index into odet_proxies. +# Each chunk is independent (identity IC, no accumulated state), so the result does not +# depend on how chunks are distributed across threads. function _run_parallel_bvp_phase!(propagators::Vector{ChunkPropagator}, chunks::Vector{IntegrationChunk}, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::ForceFreeStatesInternal, - odet_proxies::Vector{OdeState}, bvp_threads::Int) - if bvp_threads == 1 - for i in eachindex(chunks) - integrate_propagator_chunk!(propagators[i], chunks[i], ctrl, equil, ffit, intr, - odet_proxies[1]) - end - else - Threads.@threads :static for i in eachindex(chunks) - integrate_propagator_chunk!(propagators[i], chunks[i], ctrl, equil, ffit, intr, - odet_proxies[Threads.threadid()]) - end + odet_proxies::Vector{OdeState}) + Threads.@threads :static for i in eachindex(chunks) + integrate_propagator_chunk!(propagators[i], chunks[i], ctrl, equil, ffit, intr, + odet_proxies[Threads.threadid()]) end end @@ -1789,9 +1672,7 @@ function _assemble_propagators_serially!(odet::OdeState, propagators::Vector{Chu last_crossing_step = odet.step - 1 else # Save non-crossing end-of-chunk state. These columns are FM/Riccati chunk - # endpoints, not the Euler-Lagrange state — when ctrl.populate_dense_xi=true the - # entire odet is replaced by a serial-EL pass at the end of - # parallel_eulerlagrange_integration. + # endpoints, not the Euler-Lagrange state, so the odet never claims the EL basis. odet.u_store_el_basis = false if odet.step >= size(odet.u_store, 4) resize_storage!(odet) @@ -1896,85 +1777,3 @@ function _handle_edge_dW_scan!(odet::OdeState, chunks::Vector{IntegrationChunk}, end return chunks, propagators end - -""" - _populate_dense_xi_via_serial_el!(odet, ctrl, equil, ffit, intr) -> fresh_odet - -Replace the propagator-BVP's `odet` with a fresh serial-EL `odet` that has -dense `u_store` / `du_store` populated in axis basis (the PerturbedEquilibrium -convention). The caller's `odet` is fully replaced by the fresh one because -`free_run` / `normalize_eigenfunctions!` downstream use `odet.u[:,:,1,end]` to -normalize `odet.u_store`, so both must be in the same basis. The parallel BVP -results that survive -downstream are stored in `intr` (psilim/qlim, sing[*].delta_prime, …) and in -the externally-returned `propagators` / `chunks` / `S_at_surface_left` — -none of those live on `odet`, so replacing `odet` is safe. - -The dense pass uses the **serial EL path** (`sing_der!` with standard -`integrator_callback!`, Gaussian reduction, and `transform_u!`) so that -`u_store` is in the axis basis — the only convention the PerturbedEquilibrium -/ FieldReconstruction downstream code is known to consume correctly. - -We do save and restore the `intr.psilim` / `intr.qlim` / `intr.sing[*]` fields -that the parallel BVP populated, because the dense EL pass would otherwise -overwrite them (its standard `cross_ideal_singular_surf!` runs unconditionally -and does NOT populate `delta_prime`; we keep the parallel pass's values -which `compute_delta_prime_matrix!` uses). - -Called from `parallel_eulerlagrange_integration` when -`ctrl.populate_dense_xi = true`. Approximate cost: one serial -EL integration on top of the parallel BVP phase. Required to make -`use_parallel = true` produce DCON eigenfunctions usable by the -PerturbedEquilibrium downstream pipeline. -""" -function _populate_dense_xi_via_serial_el!( - odet::OdeState, ctrl::ForceFreeStatesControl, - equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, - intr::ForceFreeStatesInternal -) - msing = intr.msing - - # Preserve parallel-BVP state on intr/odet that the serial-EL pass would otherwise - # overwrite. PE downstream (SingularCoupling.jl) is calibrated against the (S, I) - # Riccati gauge of `ca_l`/`ca_r`, so keeping the parallel-BVP values is critical. - saved = ( - psilim = intr.psilim, - qlim = intr.qlim, - ca_l = copy(odet.ca_l), - ca_r = copy(odet.ca_r), - sing_state = [( - delta_prime = copy(intr.sing[s].delta_prime), - delta_prime_col = copy(intr.sing[s].delta_prime_col), - ua_left = copy(intr.sing[s].ua_left), - psi_ua_left = intr.sing[s].psi_ua_left, - ) for s in 1:msing], - ) - - if ctrl.verbose - @info " S → ξ: serial EL dense pass for HDF5 ForceFreeStates/Solutions/ForwardIntegration/xi_*" - end - - # Run the serial branch but suppress logging - fresh_odet, _, _, _ = serial_eulerlagrange_integration(ctrl, equil, ffit, intr; verbose=false) - - # Restore BVP-result fields on `intr`. - intr.psilim = saved.psilim - intr.qlim = saved.qlim - for s in 1:msing - intr.sing[s].delta_prime = saved.sing_state[s].delta_prime - intr.sing[s].delta_prime_col = saved.sing_state[s].delta_prime_col - intr.sing[s].ua_left = saved.sing_state[s].ua_left - intr.sing[s].psi_ua_left = saved.sing_state[s].psi_ua_left - end - - # Restore the parallel BVP's Riccati-gauge `ca_l` / `ca_r` onto the - # fresh EL odet — these feed PE's `SingularCoupling.jl` which is - # written against the (S, I) Riccati convention. - fresh_odet.ca_l .= saved.ca_l - fresh_odet.ca_r .= saved.ca_r - - # Return the fresh serial-EL odet (self-consistent for ξ-function - # storage in axis basis; `ca_l`/`ca_r` carry the parallel-BVP - # Riccati-gauge values needed by PE downstream). - return fresh_odet -end diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 519dd952c..27ba49f43 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -70,7 +70,8 @@ using .ForceFreeStates: find_kinetic_singular_surfaces! using .ForceFreeStates: eulerlagrange_integration, free_run, normalize_eigenfunctions! using .ForceFreeStates: galerkin_solve, write_galerkin!, GalerkinResult, gal_matched_odestate -const _DEPRECATED_FFS_KEYS = ("mer_flag", "force_wv_symmetry", "ode_flag", "cyl_flag", "mat_flag") +const _DEPRECATED_FFS_KEYS = ("mer_flag", "force_wv_symmetry", "ode_flag", "cyl_flag", "mat_flag", + "use_riccati", "use_parallel", "parallel_threads", "populate_dense_xi") const _DEPRECATED_EQUIL_KEYS = ("power_bp", "power_b", "power_r", "power_rc") # Drop deprecated keys from a parsed gpec.toml section so legacy files keep parsing @@ -577,9 +578,9 @@ function main_from_inputs( ) pe_intr = PerturbedEquilibrium.PerturbedEquilibriumInternal(; dir_path=intr.dir_path) - # DRIVEN (RPEC): feed the coil-matched gal solution to PE instead of the shooting solution. + # DRIVEN (RPEC): feed the coil-matched gal solution to PE instead of the forward solution. # The matched OdeState is in the identity-at-edge basis; build_flux_matrix rederives the edge BC - # from u_store[:,:,1,step], so PE consumes it unchanged. The shooting odet is left untouched for + # from u_store[:,:,1,step], so PE consumes it unchanged. The forward odet is left untouched for # the Force-Free States HDF5 output. pe_odet = odet if ctrl.gal_flag && ctrl.gal_match_flag && gal_data !== nothing && gal_data.match !== nothing diff --git a/src/PerturbedEquilibrium/SingularCoupling.jl b/src/PerturbedEquilibrium/SingularCoupling.jl index 8a2c12bd5..dd24dbfb1 100644 --- a/src/PerturbedEquilibrium/SingularCoupling.jl +++ b/src/PerturbedEquilibrium/SingularCoupling.jl @@ -56,14 +56,14 @@ end Evaluate the `resnum` row of Ξ_ψ and Ξ′_ψ at `psi` from the stored ODE solution: cubic Hermite for the value, chord slope across the bracketing nodes for the derivative. -Least accurate method, kept for solution paths outside the serial EL integrator +Least accurate method, kept for solution paths outside the forward EL integrator (gal-matched, Riccati) whose stored derivatives cover only Ξ′. """ function _chord_solution_at(psi::Float64, resnum::Int, odet::OdeState, nstep::Int) isempty(odet.du_store) && error( "_chord_solution_at: no derivative store. The solution is in a basis " * - "the Euler-Lagrange kernel cannot be re-applied to (sparse parallel path); " * - "set populate_dense_xi = true for PerturbedEquilibrium runs." + "the Euler-Lagrange kernel cannot be re-applied to (sparse Riccati path); " * + "dense Ξ′ requires the Forward integrator." ) il, ir, _ = _psi_bracket(odet.psi_store, psi, nstep) psi_a, psi_b = odet.psi_store[il], odet.psi_store[ir] diff --git a/src/Tearing/Dispersion/CoupledFullMatch.jl b/src/Tearing/Dispersion/CoupledFullMatch.jl index 8bf5fb97c..9b0f136c9 100644 --- a/src/Tearing/Dispersion/CoupledFullMatch.jl +++ b/src/Tearing/Dispersion/CoupledFullMatch.jl @@ -94,7 +94,7 @@ end Construct the 4m × 4m dispersion matrix driver. `dp_raw` must be the 2m × 2m matrix in side-major ordering (the `intr.delta_prime_raw` field populated by `ForceFreeStates.compute_delta_prime_matrix!` on the -`use_parallel=true` path). `rotation[k]` is the per-surface rotation +Riccati path). `rotation[k]` is the per-surface rotation frequency; it shifts the per-surface inner Q argument by `i·ntor·rotation[k]`. Default zero rotation matches the static-equilibrium case. diff --git a/test/runtests_eulerlagrange.jl b/test/runtests_eulerlagrange.jl index 7adb75203..6547352e6 100644 --- a/test/runtests_eulerlagrange.jl +++ b/test/runtests_eulerlagrange.jl @@ -432,7 +432,7 @@ end example_dir = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example") inputs = TOML.parsefile(joinpath(example_dir, "gpec.toml")) inputs["ForceFreeStates"]["verbose"] = false - inputs["ForceFreeStates"]["use_parallel"] = false + inputs["ForceFreeStates"]["integrator"] = "forward" inputs["ForceFreeStates"]["write_outputs_to_HDF5"] = false intr = FFS.ForceFreeStatesInternal(; dir_path=example_dir) ctrl = FFS.ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...) diff --git a/test/runtests_parallel_integration.jl b/test/runtests_parallel_integration.jl index d23a00790..1557587ad 100644 --- a/test/runtests_parallel_integration.jl +++ b/test/runtests_parallel_integration.jl @@ -1,7 +1,7 @@ using LinearAlgebra using TOML -@testset "Parallel FM Integration Tests" begin +@testset "Riccati FM Integration Tests" begin @testset "ChunkPropagator identity on trivial interval" begin # Integrating over a zero-width interval should give the identity propagator. @@ -106,7 +106,7 @@ using TOML @testset "balance_integration_chunks produces target count" begin # Verify that balance_integration_chunks creates at least - # max(2*msing+3, 4*nthreads) chunks from a small set of base chunks. + # max(2*msing+3, 8*(msing+1)+msing) chunks from a small set of base chunks. ex = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example") inputs = TOML.parsefile(joinpath(ex, "gpec.toml")) inputs["ForceFreeStates"]["verbose"] = false @@ -134,9 +134,11 @@ using TOML base_chunks = GeneralizedPerturbedEquilibrium.ForceFreeStates.chunk_el_integration_bounds(odet, ctrl, intr) balanced = GeneralizedPerturbedEquilibrium.ForceFreeStates.balance_integration_chunks(base_chunks, ctrl, intr) - # Must mirror balance_integration_chunks' internal target_n formula - # (src/ForceFreeStates/EulerLagrange.jl). Keep this in sync. - target_n = max(2 * intr.msing + 3, 4 * Threads.nthreads(), 8 * (intr.msing + 1) + intr.msing) + # Must mirror balance_integration_chunks' internal target_n formula for nchunks = 0 + # (src/ForceFreeStates/EulerLagrange.jl). Keep this in sync. The formula reads only + # intr.msing — no thread count enters it, which is what makes Riccati outputs + # independent of how many threads `julia -t` provides. + target_n = max(2 * intr.msing + 3, 8 * (intr.msing + 1) + intr.msing) # After balancing, chunk count equals target_n: the while-loop adds exactly one # chunk per iteration (a bisection split) and exits when length(result) >= target_n, @@ -165,6 +167,27 @@ using TOML n_crossings_base = count(c -> c.needs_crossing, base_chunks) n_crossings_bal = count(c -> c.needs_crossing, balanced) @test n_crossings_bal == n_crossings_base + + # Chunking is a pure function of (chunks, ctrl, intr): repeated calls on the same + # inputs give bit-identical boundaries. Together with the thread-free target_n + # formula above, this is what guarantees Riccati results do not move with `julia -t`. + balanced_again = GeneralizedPerturbedEquilibrium.ForceFreeStates.balance_integration_chunks(base_chunks, ctrl, intr) + @test length(balanced_again) == length(balanced) + @test all(balanced_again[i].psi_start == balanced[i].psi_start for i in eachindex(balanced)) + @test all(balanced_again[i].psi_end == balanced[i].psi_end for i in eachindex(balanced)) + + # An explicit nchunks steers the target count directly. + ctrl_more = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; + (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])..., nchunks=target_n + 7) + balanced_more = GeneralizedPerturbedEquilibrium.ForceFreeStates.balance_integration_chunks(base_chunks, ctrl_more, intr) + @test length(balanced_more) == target_n + 7 + + # An nchunks below the singular-surface floor is clamped up, with a warning. + min_chunks = 2 * intr.msing + 3 + ctrl_few = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; + (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])..., nchunks=1) + balanced_few = @test_logs (:warn,) match_mode=:any GeneralizedPerturbedEquilibrium.ForceFreeStates.balance_integration_chunks(base_chunks, ctrl_few, intr) + @test length(balanced_few) == max(min_chunks, length(base_chunks)) end @testset "chunk_el_integration_bounds direction field — bidirectional mode" begin @@ -222,19 +245,19 @@ using TOML end end - @testset "Parallel FM integration matches standard ODE — Solovev example" begin - # Run standard and parallel FM integrations on the Solovev regression test. + @testset "Riccati FM integration matches forward ODE — Solovev example" begin + # Run forward and Riccati FM integrations on the Solovev regression test. # The energy eigenvalue et[1] should match to within 2%. # - # Bidirectional FM integration (crossing chunks integrated backward) is the - # default for use_parallel=true. It keeps FM propagators well-conditioned for - # both small-N (Solovev N=8, tested here) and large-N (DIIID N=26, tested below). + # Bidirectional FM integration (crossing chunks integrated backward) is what the + # Riccati path uses. It keeps FM propagators well-conditioned for both small-N + # (Solovev N=8, tested here) and large-N (DIIID N=26, tested below). ex = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example") - function run_solovev(use_parallel) + function run_solovev(integrator) inputs = TOML.parsefile(joinpath(ex, "gpec.toml")) inputs["ForceFreeStates"]["verbose"] = false - inputs["ForceFreeStates"]["use_parallel"] = use_parallel + inputs["ForceFreeStates"]["integrator"] = integrator intr = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesInternal(; dir_path=ex) ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...) @@ -258,11 +281,11 @@ using TOML return real(vac.et[1]), intr end - et_std, intr_std = run_solovev(false) - et_par, intr_par = run_solovev(true) + et_fwd, intr_fwd = run_solovev("forward") + et_ric, intr_ric = run_solovev("riccati") # Energy eigenvalue matches to 2% - @test isapprox(et_par, et_std; rtol=0.02) + @test isapprox(et_ric, et_fwd; rtol=0.02) # Per-surface Δ' assertions were removed: per-surface Δ' is a stub calculation # left in the code for future work but no longer reported, output, or tested. # The STRIDE BVP Δ' matrix (`SingularSurfaces/delta_prime_matrix`) is the canonical @@ -270,8 +293,8 @@ using TOML # values; Solovev is near marginal stability and BVP Δ' is pathological there. end - @testset "Parallel FM integration matches standard ODE — DIIID-like example (large N)" begin - # Run standard and parallel FM integrations on the DIIID-like example (N≈26 modes). + @testset "Riccati FM integration matches forward ODE — DIIID-like example (large N)" begin + # Run forward and Riccati FM integrations on the DIIID-like example (N≈26 modes). # Before bidirectional integration, the all-forward FM propagators were ill-conditioned # for large N, producing ~10% energy error. Bidirectional integration (backward crossing # chunks + forward intermediate chunks) restores accuracy to within 2%. @@ -279,10 +302,10 @@ using TOML # This is the key regression test for the bidirectional parallel FM fix. ex = joinpath(@__DIR__, "..", "examples", "DIIID-like_ideal_example") - function run_diiid(use_parallel) + function run_diiid(integrator) inputs = TOML.parsefile(joinpath(ex, "gpec.toml")) inputs["ForceFreeStates"]["verbose"] = false - inputs["ForceFreeStates"]["use_parallel"] = use_parallel + inputs["ForceFreeStates"]["integrator"] = integrator inputs["ForceFreeStates"]["write_outputs_to_HDF5"] = false intr = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesInternal(; dir_path=ex) ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; @@ -317,18 +340,18 @@ using TOML return real(vac.et[1]), intr end - et_par, intr_par = run_diiid(true) + et_ric, intr_ric = run_diiid("riccati") - # Parallel FM et[1] regression — pinned tightly, NOT bracketed. et[1] is grid- and + # Riccati FM et[1] regression — pinned tightly, NOT bracketed. et[1] is grid- and # equilibrium-sensitive (auto-mpsi gives a spurious value; a wrong grid/Ip shifts it), so - # a loose bracket would mask exactly that accuracy regression. The parallel-path value is + # a loose bracket would mask exactly that accuracy regression. The Riccati-path value is # deterministic and reproducible. - @test isapprox(et_par, 0.800637; rtol=2e-2) + @test isapprox(et_ric, 0.800637; rtol=2e-2) # Per-surface Δ' assertions removed (stub calculation; see Solovev testset # comment above). BVP Δ' matrix regression for DIIID-like is in the # `delta_prime_matrix — STRIDE BVP DIIID-like regression (large N)` testset. - # No explicit parallel-vs-standard cross-path check here: the two paths share the + # No explicit Riccati-vs-forward cross-path check here: the two paths share the # equilibrium grid (so a cross-path comparison is blind to grid/accuracy regressions), # and their agreement is already verified on the lighter Solovev case above. The tight # absolute pin above is the guard for grid/equilibrium regressions on this case. @@ -386,28 +409,15 @@ using TOML # physically meaningful. BVP Δ' regression is concentrated on the DIIID-like # fixture below (intrinsically stable, well-conditioned BVP Δ'). - @testset "ξ functions bit-identical between use_parallel modes (populate_dense_xi)" begin - # When `ctrl.use_parallel = true` and `ctrl.populate_dense_xi = true` - # (default), `parallel_eulerlagrange_integration` appends a serial - # Euler-Lagrange pass and returns that fresh `odet` instead of the - # propagator-BVP one. That dense pass invokes the SAME - # `eulerlagrange_integration` code path the serial `use_parallel = false` - # benchmark goes through with the SAME `(ctrl, equil, ffit, intr)` - # inputs (BVP-only state on `intr` saved/restored across the pass), so - # the resulting `psi_store` / `q_store` / `u_store` / `du_store` / - # `crit_store` arrays must be bit-identical to a standalone serial run. - # This is a strong correctness guarantee that the dense pass does NOT - # perturb the DCON eigenfunction calculation in any way — exactly what - # downstream PerturbedEquilibrium / FieldReconstruction needs. - # - # Run on both the small-N Solovev case and the large-N DIIID-like case - # to catch any (m, IC, ψ)-dependent regression. - - function run_and_capture(example_dir, use_parallel; populate_dense_xi=true) + @testset "Riccati leaves a sparse u_store in the Riccati basis" begin + # The Riccati path stores only chunk-endpoint snapshots, so u_store / psi_store have + # strictly fewer entries than the forward path's dense saved steps, and the stored + # state is not in the Euler-Lagrange axis basis. Downstream ξ consumers must therefore + # be fed by the forward integrator; this test pins that contract. + function run_and_capture(example_dir, integrator) inputs = TOML.parsefile(joinpath(example_dir, "gpec.toml")) inputs["ForceFreeStates"]["verbose"] = false - inputs["ForceFreeStates"]["use_parallel"] = use_parallel - inputs["ForceFreeStates"]["populate_dense_xi"] = populate_dense_xi + inputs["ForceFreeStates"]["integrator"] = integrator inputs["ForceFreeStates"]["write_outputs_to_HDF5"] = false intr = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesInternal(; dir_path=example_dir) ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; @@ -433,61 +443,23 @@ using TOML return odet end - # Compare the storage arrays that downstream code reads. All values - # must be EXACTLY equal (no tolerance — the dense pass calls the same - # ODE solver with the same inputs as the standalone serial path, so - # any nonzero difference indicates a real regression in the dense-pass - # machinery). - function assert_bit_identical(odet_a, odet_b) - @test odet_a.step == odet_b.step - @test odet_a.nzero == odet_b.nzero - @test length(odet_a.psi_store) == length(odet_b.psi_store) - @test length(odet_a.q_store) == length(odet_b.q_store) - @test size(odet_a.u_store) == size(odet_b.u_store) - @test size(odet_a.du_store) == size(odet_b.du_store) - @test size(odet_a.xi_s_store) == size(odet_b.xi_s_store) - @test odet_a.du_store_populated == odet_b.du_store_populated - @test maximum(abs.(odet_a.psi_store .- odet_b.psi_store)) == 0.0 - @test maximum(abs.(odet_a.q_store .- odet_b.q_store)) == 0.0 - @test maximum(abs.(odet_a.u_store .- odet_b.u_store)) == 0.0 - @test maximum(abs.(odet_a.du_store .- odet_b.du_store)) == 0.0 - @test maximum(abs.(odet_a.xi_s_store .- odet_b.xi_s_store)) == 0.0 - @test maximum(abs.(odet_a.crit_store .- odet_b.crit_store)) == 0.0 - end - - @testset "Solovev (small N)" begin - ex = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example") - odet_std = run_and_capture(ex, false) - odet_par = run_and_capture(ex, true; populate_dense_xi=true) - assert_bit_identical(odet_std, odet_par) - end - - @testset "DIIID-like (large N)" begin - ex = joinpath(@__DIR__, "..", "examples", "DIIID-like_ideal_example") - odet_std = run_and_capture(ex, false) - odet_par = run_and_capture(ex, true; populate_dense_xi=true) - assert_bit_identical(odet_std, odet_par) - end - - @testset "populate_dense_xi=false leaves sparse u_store (control)" begin - # Sanity-check the opposite mode: with populate_dense_xi=false, the - # parallel BVP path stores only chunk-endpoint Riccati snapshots, - # so u_store / du_store / psi_store have strictly fewer entries - # than the serial path. Catching this guarantees the bit-identical - # test above is meaningful — it's NOT trivially passing because - # both modes accidentally produce the same sparse data. - ex = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example") - odet_std = run_and_capture(ex, false) - odet_sparse = run_and_capture(ex, true; populate_dense_xi=false) - @test odet_sparse.step < odet_std.step - @test length(odet_sparse.psi_store) < length(odet_std.psi_store) - # The sparse solution is in the Riccati basis, so the derivative stores cannot be - # materialized from it and stay empty rather than holding unusable values. - @test !odet_sparse.u_store_el_basis - @test !odet_sparse.du_store_populated - @test isempty(odet_sparse.du_store) - @test isempty(odet_sparse.xi_s_store) - end + ex = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example") + odet_fwd = run_and_capture(ex, "forward") + odet_ric = run_and_capture(ex, "riccati") + + @test odet_ric.step < odet_fwd.step + @test length(odet_ric.psi_store) < length(odet_fwd.psi_store) + + # The forward path returns dense ξ in the axis basis, ready for PerturbedEquilibrium. + @test odet_fwd.u_store_el_basis + @test odet_fwd.du_store_populated + + # The Riccati solution is in the Riccati basis, so the derivative stores cannot be + # materialized from it and stay empty rather than holding unusable values. + @test !odet_ric.u_store_el_basis + @test !odet_ric.du_store_populated + @test isempty(odet_ric.du_store) + @test isempty(odet_ric.xi_s_store) end @testset "delta_prime_matrix — STRIDE BVP DIIID-like regression (large N)" begin @@ -498,7 +470,7 @@ using TOML ex = joinpath(@__DIR__, "..", "examples", "DIIID-like_ideal_example") inputs = TOML.parsefile(joinpath(ex, "gpec.toml")) inputs["ForceFreeStates"]["verbose"] = false - inputs["ForceFreeStates"]["use_parallel"] = true + inputs["ForceFreeStates"]["integrator"] = "riccati" inputs["ForceFreeStates"]["write_outputs_to_HDF5"] = false intr = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesInternal(; dir_path=ex) ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; diff --git a/test/runtests_rerun_from_h5.jl b/test/runtests_rerun_from_h5.jl index 2e364f2df..dc242426f 100644 --- a/test/runtests_rerun_from_h5.jl +++ b/test/runtests_rerun_from_h5.jl @@ -94,6 +94,23 @@ end end end + # A gpec.h5 written before the integrator refactor embeds retired keys such as + # `use_parallel` in its gpec_toml_raw, so a replay must warn and drop them rather than + # failing on an unknown keyword. The override flag stands in for that stale blob. + @testset "retired ForceFreeStates keys are warned and ignored on replay" begin + mktempdir() do out_dir + inputs, = GeneralizedPerturbedEquilibrium.build_inputs_from_h5([source_h5, "--output-dir", out_dir, "--override", "ForceFreeStates.use_parallel=true"]) + @test inputs["ForceFreeStates"]["use_parallel"] == true + + @test_logs (:warn,) GeneralizedPerturbedEquilibrium._drop_deprecated_keys!(inputs["ForceFreeStates"], GeneralizedPerturbedEquilibrium._DEPRECATED_FFS_KEYS, "ForceFreeStates") + @test !haskey(inputs["ForceFreeStates"], "use_parallel") + + # With the key gone the control struct builds again, on the default integrator. + ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...) + @test ctrl.integrator in ("forward", "riccati") + end + end + @testset "rerun refuses to overwrite its own source" begin # Default output is `_rerun.h5`, so clobbering only happens if the user # routes the rerun back onto the source file; exercise that guard. diff --git a/test/runtests_riccati.jl b/test/runtests_riccati.jl index 57b345702..afd084056 100644 --- a/test/runtests_riccati.jl +++ b/test/runtests_riccati.jl @@ -93,7 +93,7 @@ end # # Integration runs: # intr_ric / odet_ric — Riccati path (shared by most tests) - # intr_std / odet_std — Standard path (energy comparison only) + # intr_fwd / odet_fwd — Forward path (energy comparison only) ex = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example") inputs = TOML.parsefile(joinpath(ex, "gpec.toml")) @@ -110,9 +110,10 @@ end ffit = FFS.make_matrix(equil, intr_tmp, metric) N = intr_tmp.numpert_total - # Riccati integration + # Riccati integration. The driver returns (odet, propagators, chunks, S_at_surface_left); + # only odet is used here. intr_ric = make_solovev_intr(inputs, ctrl, equil, ex) - odet_ric = FFS.riccati_eulerlagrange_integration(ctrl, equil, ffit, intr_ric) + odet_ric, _, _, _ = FFS.riccati_eulerlagrange_integration(ctrl, equil, ffit, intr_ric) # Save inline Δ' values before any test that calls compute_delta_prime_from_ca! # (which overwrites intr_ric.sing[s].delta_prime) @@ -121,23 +122,23 @@ end vac_ric = FFS.free_run(odet_ric, ctrl, equil, ffit, intr_ric) et_ric = real(vac_ric.et[1]) - # Standard integration (needed only for energy comparison). eulerlagrange_integration - # returns (odet, propagators, chunks, S_at_surface_left); only odet is used here. - intr_std = make_solovev_intr(inputs, ctrl, equil, ex) - odet_std, _, _, _ = FFS.eulerlagrange_integration(ctrl, equil, ffit, intr_std) - vac_std = FFS.free_run(odet_std, ctrl, equil, ffit, intr_std) - et_std = real(vac_std.et[1]) + # Forward integration (needed only for energy comparison). + intr_fwd = make_solovev_intr(inputs, ctrl, equil, ex) + odet_fwd, _, _, _ = FFS.forward_eulerlagrange_integration(ctrl, equil, ffit, intr_fwd) + vac_fwd = FFS.free_run(odet_fwd, ctrl, equil, ffit, intr_fwd) + et_fwd = real(vac_fwd.et[1]) # ───────────────────────────────────────────────────────────────────────── - @testset "Riccati integration matches standard ODE — Solovev example" begin - # PR description claims Solovev energy eigenvalue error 0.006 % vs standard path. - # Tightened to rtol=1e-4 (matches the PR's headline claim within ≈2×). A regression - # of the Riccati/renormalization algorithm to ~1 % error would fail here loudly. - @test isapprox(et_ric, et_std; rtol=1e-4) + @testset "Riccati integration matches forward ODE — Solovev example" begin + # The two formalisms solve the same system, so the leading energy eigenvalue must + # agree closely; the measured Solovev disagreement is ≈5e-5. A regression of the + # Riccati/renormalization algorithm to ~1 % error would fail here loudly. + @test isapprox(et_ric, et_fwd; rtol=1e-4) - # Riccati uses no more than 2x as many steps as standard - @test odet_ric.step <= 2 * odet_std.step + # Riccati stores chunk endpoints, the forward path stores every saved ODE step, + # so the Riccati store is always the sparser of the two. + @test odet_ric.step <= odet_fwd.step end # Note: a Solovev per-surface Δ' regression testset previously lived here, @@ -148,8 +149,8 @@ end # STRIDE BVP Δ' matrix (see runtests_parallel_integration.jl). @testset "Riccati end state has U₂ ≈ I" begin - # After riccati_eulerlagrange_integration, odet.u[:,:,2] should be identity - # (canonical Riccati convention after final renorm) + # The outer-region re-integration that closes riccati_eulerlagrange_integration + # leaves odet.u[:,:,2] as the identity (canonical Riccati convention after final renorm) @test odet_ric.u[:, :, 2] ≈ I(N) rtol=1e-10 end diff --git a/test/test_data/regression_solovev_ideal_example/gpec.toml b/test/test_data/regression_solovev_ideal_example/gpec.toml index a16ab6898..6b9d537e6 100644 --- a/test/test_data/regression_solovev_ideal_example/gpec.toml +++ b/test/test_data/regression_solovev_ideal_example/gpec.toml @@ -45,10 +45,8 @@ eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced ucrit = 1e3 # Column-norm threshold that triggers solution renormalization -# Δ' BVP + parallel integration (see ForceFreeStatesControl docstring for details) -use_parallel = true # Run parallel FM-propagator BVP path (unlocks SingularSurfaces/delta_prime_matrix) -parallel_threads = 2 # BVP thread cap (1 = serial/bit-deterministic; 2 ≈ +20% speedup; ≥3 saturates) -populate_dense_xi = true # Append serial-EL pass so dense ξ is stored — REQUIRED with a [PerturbedEquilibrium] section +# Integrator selection and domain truncation (see ForceFreeStatesControl docstring for details) +integrator = "forward" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic) or "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix) set_psilim_via_dmlim = false # FALSE for limited/analytical equilibria — rationals sparse, dmlim would chop too much edge dmlim = 0.2 # Truncate integration at (last_rational_q + dmlim)/n (used when set_psilim_via_dmlim = true) diff --git a/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml b/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml index 0f0bc5c47..62e8569fb 100644 --- a/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml +++ b/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml @@ -45,10 +45,8 @@ eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced ucrit = 1e3 # Column-norm threshold that triggers solution renormalization -# Δ' BVP + parallel integration (see ForceFreeStatesControl docstring for details) -use_parallel = true # Run parallel FM-propagator BVP path (multi-n Δ' matrix has open issues — sing_lim! warns and skips — but ξ and energies are valid) -parallel_threads = 2 # BVP thread cap (1 = serial/bit-deterministic; 2 ≈ +20% speedup; ≥3 saturates) -populate_dense_xi = true # Append serial-EL pass so dense ξ is stored — REQUIRED with a [PerturbedEquilibrium] section +# Integrator selection and domain truncation (see ForceFreeStatesControl docstring for details) +integrator = "forward" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic) or "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix) set_psilim_via_dmlim = false # FALSE for multi-n — dmlim truncation is ambiguous when n varies (sing_lim! skips anyway) dmlim = 0.2 # Truncate integration at (last_rational_q + dmlim)/n (used when set_psilim_via_dmlim = true) diff --git a/test/test_data/regression_solovev_kinetic_calculated/gpec.toml b/test/test_data/regression_solovev_kinetic_calculated/gpec.toml index 5b87267a9..cca19bfff 100644 --- a/test/test_data/regression_solovev_kinetic_calculated/gpec.toml +++ b/test/test_data/regression_solovev_kinetic_calculated/gpec.toml @@ -42,6 +42,7 @@ mthvac = 64 # Number of points used in splines over poloidal a kinetic_source = "calculated" # Kinetic matrix source: "fixed" test matrices, or "calculated" from the kinetic NTV model kinetic_factor = 1.0 # Scaling of kinetic matrices (0 = ideal path; >0 enables kinetic mode; 1.0 = full strength) +integrator = "forward" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic) or "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix) eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced ucrit = 1e3 # Column-norm threshold that triggers solution renormalization diff --git a/test/test_data/regression_solovev_kinetic_example/gpec.toml b/test/test_data/regression_solovev_kinetic_example/gpec.toml index 1c020a086..d8349716b 100644 --- a/test/test_data/regression_solovev_kinetic_example/gpec.toml +++ b/test/test_data/regression_solovev_kinetic_example/gpec.toml @@ -45,10 +45,8 @@ eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced ucrit = 1e3 # Column-norm threshold that triggers solution renormalization -# Δ' BVP + parallel integration (see ForceFreeStatesControl docstring for details) -use_parallel = true # Run parallel FM-propagator BVP path (unlocks SingularSurfaces/delta_prime_matrix) -parallel_threads = 2 # BVP thread cap (1 = serial/bit-deterministic; 2 ≈ +20% speedup; ≥3 saturates) -populate_dense_xi = true # Append serial-EL pass so dense ξ is stored — REQUIRED with a [PerturbedEquilibrium] section +# Integrator selection and domain truncation (see ForceFreeStatesControl docstring for details) +integrator = "forward" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic) or "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix) set_psilim_via_dmlim = false # FALSE for limited/analytical equilibria — rationals sparse, dmlim would chop too much edge dmlim = 0.2 # Truncate integration at (last_rational_q + dmlim)/n (used when set_psilim_via_dmlim = true) diff --git a/test/test_data/regression_solovev_kinetic_multi_n/gpec.toml b/test/test_data/regression_solovev_kinetic_multi_n/gpec.toml index b11caab8b..4527c7826 100644 --- a/test/test_data/regression_solovev_kinetic_multi_n/gpec.toml +++ b/test/test_data/regression_solovev_kinetic_multi_n/gpec.toml @@ -45,10 +45,8 @@ eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced ucrit = 1e3 # Column-norm threshold that triggers solution renormalization -# Δ' BVP + parallel integration (see ForceFreeStatesControl docstring for details) -use_parallel = true # Run parallel FM-propagator BVP path (multi-n Δ' matrix has open issues — sing_lim! warns and skips — but ξ and energies are valid) -parallel_threads = 2 # BVP thread cap (1 = serial/bit-deterministic; 2 ≈ +20% speedup; ≥3 saturates) -populate_dense_xi = true # Append serial-EL pass so dense ξ is stored — REQUIRED with a [PerturbedEquilibrium] section +# Integrator selection and domain truncation (see ForceFreeStatesControl docstring for details) +integrator = "forward" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic) or "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix) set_psilim_via_dmlim = false # FALSE for multi-n — dmlim truncation is ambiguous when n varies (sing_lim! skips anyway) dmlim = 0.2 # Truncate integration at (last_rational_q + dmlim)/n (used when set_psilim_via_dmlim = true) diff --git a/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml b/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml index ca6eeb003..ff4bc38ae 100644 --- a/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml +++ b/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml @@ -42,6 +42,7 @@ mthvac = 64 # Number of points used in splines over poloidal a kinetic_source = "calculated" # Kinetic matrix source — exercises KineticForces.compute_calculated_kinetic_matrices callback with real physics kinetic_factor = 1.0 # Full-strength kinetic matrices (the "calculated" path is the real physics; no perturbation scaling) +integrator = "forward" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic) or "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix) eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced ucrit = 1e3 # Maximum fraction of solutions allowed before re-normalized From f8996d4f0d9a7261bc48a9fe2389028e703009dd Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Fri, 14 Aug 2026 23:08:10 -0400 Subject: [PATCH 2/3] FFS - REFACTOR - Extract local stability into standalone LocalStability module --- REFACTOR_PLAN.md | 6 +- docs/development/architecture.md | 14 ++-- docs/src/ballooning.md | 6 ++ docs/src/stability.md | 2 +- .../analyze_example.jl | 15 ++-- src/Analysis/ForceFreeStates.jl | 4 +- src/ForceFreeStates/ForceFreeStates.jl | 1 - src/GeneralizedPerturbedEquilibrium.jl | 10 ++- .../Ballooning.jl | 70 +++++++++---------- src/LocalStability/LocalStability.jl | 17 +++++ test/runtests.jl | 1 + test/runtests_resist_eval.jl | 3 +- 12 files changed, 92 insertions(+), 57 deletions(-) rename src/{ForceFreeStates => LocalStability}/Ballooning.jl (96%) create mode 100644 src/LocalStability/LocalStability.jl diff --git a/REFACTOR_PLAN.md b/REFACTOR_PLAN.md index 0f5bdc00b..de545c427 100644 --- a/REFACTOR_PLAN.md +++ b/REFACTOR_PLAN.md @@ -670,7 +670,11 @@ manual smoke: run the 4-line UX from the Context section in a REPL against from the §3 spec, both improvements: the new Δ′ example references the DIIID geqdsk by relative path instead of copying it, and the TOML sweep covered six regression fixtures (two more had landed on develop since the plan was written), all `forward`. -- [ ] PR 2 — `refactor/local-stability-module` +- [ ] PR 2 — `refactor/local-stability-module` — **implemented, in review.** One delta + from the §4 spec: the signature change also required updating two call-site groups the + section did not list — `examples/DIIID-like_ideal_example/analyze_example.jl` (five + ballooning entry points) and two docstring cross-references in + `src/Analysis/ForceFreeStates.jl`. - [ ] PR 3 — `refactor/forcefreestates-result` - [ ] PR 4 — `refactor/staged-main` - [ ] PR 5 — `feature/solve-api` diff --git a/docs/development/architecture.md b/docs/development/architecture.md index ef756c9c0..a6ddef8d9 100644 --- a/docs/development/architecture.md +++ b/docs/development/architecture.md @@ -12,7 +12,7 @@ This workflow is reflected in the modular structure and data flow. ## Module Structure -GPEC consists of **seven main modules** organized in `src/`: +GPEC consists of **eight main modules** organized in `src/`: ### Foundation Modules @@ -67,18 +67,23 @@ GPEC consists of **seven main modules** organized in `src/`: - `Fourfit.jl` - Fourier fitting routines - `FixedBoundaryStability.jl` - Fixed boundary analysis - `Free.jl` - Free boundary stability - - `Ballooning.jl` - Local stability scan: Mercier D_I, resistive interchange D_R, and high-n ballooning Δ' (s–α). Replaces the former standalone `Mercier.jl`. - Status: Stable, core DCON functionality implemented +6. **LocalStability** (`src/LocalStability/`) - Local high-n stability + - `Ballooning.jl` - Local stability scan: Mercier D_I, resistive interchange D_R, and high-n ballooning Δ' (s–α). Replaces the former standalone `Mercier.jl`. + - Depends only on Equilibrium (plus math libraries); carries no stability-solver state + - Main entry points: `compute_local_stability`, `ballooning_alpha_boundary` + - Status: Stable + ### Perturbed Equilibrium Modules -6. **ForcingTerms** (`src/ForcingTerms/`) - External field specification +7. **ForcingTerms** (`src/ForcingTerms/`) - External field specification - Handles external magnetic field perturbations (coils, RMP, etc.) - Supports ASCII and HDF5 forcing data formats - `ForcingMode` data structure specifies amplitude and phase for each (m,n) component - Status: Complete and functional -7. **PerturbedEquilibrium** (`src/PerturbedEquilibrium/`) - **GPEC-style plasma response** +8. **PerturbedEquilibrium** (`src/PerturbedEquilibrium/`) - **GPEC-style plasma response** - Computes plasma response to external forcing - Calculates singular coupling metrics at rational surfaces - Key files: @@ -191,6 +196,7 @@ GeneralizedPerturbedEquilibrium ├── Utilities (shared tools) │ └── FourierTransforms ├── Equilibrium (uses Splines) +├── LocalStability (uses Equilibrium) ├── Vacuum (uses Splines, Equilibrium, Utilities) ├── ForcingTerms (data I/O) ├── ForceFreeStates (uses Equilibrium, Vacuum, Splines) diff --git a/docs/src/ballooning.md b/docs/src/ballooning.md index 72697da20..8ef6cd10b 100644 --- a/docs/src/ballooning.md +++ b/docs/src/ballooning.md @@ -911,3 +911,9 @@ The local-stability output now stores ballooning ``\Delta'`` in the fourth `locstab_fs` entry. In the HDF5 output this is written as `LocalStability/ballooning_Delta_prime`, distinct from the tearing ``\Delta'`` outputs under `SingularSurfaces/` and `PerturbedEquilibrium/SingularCoupling/`. + +## API + +```@autodocs +Modules = [GeneralizedPerturbedEquilibrium.LocalStability] +``` diff --git a/docs/src/stability.md b/docs/src/stability.md index 0e8b72880..dadf13762 100644 --- a/docs/src/stability.md +++ b/docs/src/stability.md @@ -270,7 +270,7 @@ The Galerkin Δ′ solver (`src/ForceFreeStates/Galerkin/`) is documented separa ```@autodocs Modules = [GeneralizedPerturbedEquilibrium.ForceFreeStates] -Pages = ["ForceFreeStates.jl", "ForceFreeStatesStructs.jl", "Ballooning.jl", "Resist.jl", "EulerLagrange.jl", "Sing.jl", "Fourfit.jl", "Kinetic.jl", "FixedBoundaryStability.jl", "Utils.jl", "Free.jl", "Riccati.jl"] +Pages = ["ForceFreeStates.jl", "ForceFreeStatesStructs.jl", "Resist.jl", "EulerLagrange.jl", "Sing.jl", "Fourfit.jl", "Kinetic.jl", "FixedBoundaryStability.jl", "Utils.jl", "Free.jl", "Riccati.jl"] ``` ## Example usage diff --git a/examples/DIIID-like_ideal_example/analyze_example.jl b/examples/DIIID-like_ideal_example/analyze_example.jl index 319a7063d..715e3ec76 100644 --- a/examples/DIIID-like_ideal_example/analyze_example.jl +++ b/examples/DIIID-like_ideal_example/analyze_example.jl @@ -2,7 +2,7 @@ using Pkg; Pkg.activate(joinpath(@__DIR__, "../..")) using GeneralizedPerturbedEquilibrium, Plots using GeneralizedPerturbedEquilibrium: Analysis -using GeneralizedPerturbedEquilibrium: Equilibrium, ForceFreeStates +using GeneralizedPerturbedEquilibrium: Equilibrium, ForceFreeStates, LocalStability using LaTeXStrings, TOML isinteractive() ? plotlyjs() : gr() @@ -29,7 +29,7 @@ Plots.savefig(p_baloo, baloo_path) # ---------------------------------------------------------------------- # Local stability: s-alpha profiles, D_I / ballooning Δ', and 2D s-alpha # maps (ported from the former Bal_salpha_delta_di_summary notebook). -# These recompute directly from the equilibrium via the Ballooning.jl helpers +# These recompute directly from the equilibrium via the LocalStability helpers # rather than reading gpec.h5, so the s-alpha scan can perturb (p', q'). # ---------------------------------------------------------------------- @@ -66,7 +66,7 @@ for geqdsk_file in geqdsk_cases alpha_profile = fill(NaN, length(psi_norm)) for i in eachindex(psi_norm) try - ref = ForceFreeStates.salpha_reference(i, equil) + ref = LocalStability.salpha_reference(i, equil) s_profile[i] = ref.s_ref alpha_profile[i] = ref.alpha_ref catch err @@ -107,8 +107,7 @@ for geqdsk_file in geqdsk_cases # D_I (Mercier) and ballooning Δ' profiles locstab_fs = zeros(length(psi_norm), 5) - ctrl = ForceFreeStates.ForceFreeStatesControl(; verbose=false) - ForceFreeStates.compute_ballooning_stability!(ctrl, locstab_fs, equil; compute_delta_prime=true) + LocalStability.compute_ballooning_stability!(locstab_fs, equil; compute_delta_prime=true) delta_prime = Vector(locstab_fs[:, 4]) di_profile = fill(NaN, length(psi_norm)) @@ -152,7 +151,7 @@ for geqdsk_file in geqdsk_cases psi_idx_scan = argmin(abs.(psi_norm .- psi_target)) s_scales = collect(range(-5.0, 5.0; length=30)) alpha_scales = collect(range(-5.0, 5.0; length=30)) - scan = ForceFreeStates.scan_delta_prime_map(psi_idx_scan, equil; theta_k=0.0, s_scales=s_scales, alpha_scales=alpha_scales) + scan = LocalStability.scan_delta_prime_map(psi_idx_scan, equil; theta_k=0.0, s_scales=s_scales, alpha_scales=alpha_scales) scan_qprime = scan.dqdpsi_ref .* scan.s_scales scan_pprime = scan.pprime_ref .* scan.alpha_scales @@ -210,8 +209,8 @@ for geqdsk_file in geqdsk_cases Plots.savefig(p_zero, zero_path) # BALOO-style α vs ψ stability boundaries (1st and 2nd) over the Δ'(ψ, α) map - bnd = ForceFreeStates.ballooning_alpha_boundaries(ctrl, equil) - dpmap = ForceFreeStates.ballooning_delta_prime_map(ctrl, equil) + bnd = LocalStability.ballooning_alpha_boundaries(equil) + dpmap = LocalStability.ballooning_delta_prime_map(equil) p_alpha_bnd = Analysis.ForceFreeStates.plot_ballooning_alpha_boundaries(bnd, dpmap) display(p_alpha_bnd) alpha_bnd_path = joinpath(@__DIR__, "alpha_stability_boundary_$(case_label).png") diff --git a/src/Analysis/ForceFreeStates.jl b/src/Analysis/ForceFreeStates.jl index c92d6e698..e4ce95103 100644 --- a/src/Analysis/ForceFreeStates.jl +++ b/src/Analysis/ForceFreeStates.jl @@ -584,7 +584,7 @@ arrays leave natural gaps over always-stable surfaces without requiring explicit ### Arguments - `bnd`: NamedTuple with fields `psi`, `alpha`, `alpha_critical1`, `alpha_critical2` - (as returned by `ForceFreeStates.ballooning_alpha_boundaries`) + (as returned by `LocalStability.ballooning_alpha_boundaries`) ### Keyword arguments @@ -617,7 +617,7 @@ end plot_ballooning_alpha_boundaries(bnd, dpmap; save_path=nothing, psi_min=0.0) Same diagram drawn over a heatmap of the signed Δ' from -`ForceFreeStates.ballooning_delta_prime_map`: each surface's Δ'(α) is oriented by the +`LocalStability.ballooning_delta_prime_map`: each surface's Δ'(α) is oriented by the sign of its α=0 (stable) value so that positive is stable everywhere, regridded from its native physical α = α_ref*scale onto a shared uniform α axis, and shown with the Δ'=0 contour, the extracted boundaries, and the scan cap `max_alpha_scale*α_exp`. diff --git a/src/ForceFreeStates/ForceFreeStates.jl b/src/ForceFreeStates/ForceFreeStates.jl index 4e6d3c0f1..a88f9e9b8 100644 --- a/src/ForceFreeStates/ForceFreeStates.jl +++ b/src/ForceFreeStates/ForceFreeStates.jl @@ -24,7 +24,6 @@ import StaticArrays: @MMatrix # Include all necessary files include("ForceFreeStatesStructs.jl") -include("Ballooning.jl") include("Resist.jl") include("EulerLagrange.jl") include("Sing.jl") diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 27ba49f43..cd7497dc0 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -20,6 +20,11 @@ include("Equilibrium/Equilibrium.jl") import .Equilibrium as Equilibrium export Equilibrium +# Local high-n stability (Mercier, resistive interchange, ballooning Δ'); depends only on Equilibrium. +include("LocalStability/LocalStability.jl") +import .LocalStability as LocalStability +export LocalStability + include("Vacuum/Vacuum.jl") import .Vacuum as Vacuum export Vacuum @@ -64,7 +69,6 @@ include("Rerun.jl") # Import ForceFreeStates types and functions needed for main using .ForceFreeStates: ForceFreeStatesInternal, ForceFreeStatesControl, DebugSettings, FreeBoundaryResult, OdeState, FourFitVars using .ForceFreeStates: sing_lim!, sing_min!, sing_find!, resist_eval_all!, resist_geometry, ResistGeometry -using .ForceFreeStates: compute_local_stability, compute_ballooning_stability!, ballooning_alpha_boundary, ballooning_alpha_boundaries using .ForceFreeStates: make_metric, make_matrix, make_kinetic_matrix using .ForceFreeStates: find_kinetic_singular_surfaces! using .ForceFreeStates: eulerlagrange_integration, free_run, normalize_eigenfunctions! @@ -336,9 +340,9 @@ function main_from_inputs( locstab = nothing ballooning_boundary = (psi=Float64[], alpha=Float64[], alpha_critical=Float64[]) if ctrl.local_stability_flag - locstab = compute_local_stability(ctrl, equil) + locstab = LocalStability.compute_local_stability(equil; verbose=ctrl.verbose) # First ballooning stability boundary (α vs ψ_N) for BALOO-style diagnostics. - ballooning_boundary = ballooning_alpha_boundary(ctrl, equil) + ballooning_boundary = LocalStability.ballooning_alpha_boundary(equil; verbose=ctrl.verbose) end # Find all singular surfaces in the equilibrium diff --git a/src/ForceFreeStates/Ballooning.jl b/src/LocalStability/Ballooning.jl similarity index 96% rename from src/ForceFreeStates/Ballooning.jl rename to src/LocalStability/Ballooning.jl index 51da221be..95e7c0a08 100644 --- a/src/ForceFreeStates/Ballooning.jl +++ b/src/LocalStability/Ballooning.jl @@ -2,8 +2,6 @@ # Main Driver for Ballooning Stability Analysis # Computes ballooning stability criterion over all flux surfaces # ====================================================================== -using LinearAlgebra -using StaticArrays: SVector const BALLOONING_THETA_MAX_CAP = 16.5 const BALLOONING_THETA_SCALE_MULTIPLIER = 10.0 @@ -17,7 +15,7 @@ _in_ballooning_scan_window(psi::Float64, psi_edge::Float64) = BALLOONING_SCAN_PSI_WINDOW[1] <= psi <= min(BALLOONING_SCAN_PSI_WINDOW[2], psi_edge) """ - compute_ballooning_stability!(ctrl, locstab_fs, plasma_eq) + compute_ballooning_stability!(locstab_fs, plasma_eq) Main driver routine for local high-n stability analysis. Iterates over all magnetic flux surfaces, prepares ballooning coefficients, stores `det(d0bar)` @@ -26,9 +24,9 @@ to compute Delta Prime. ## Arguments - - `ctrl::ForceFreeStatesControl`: Control parameters for the analysis. - `locstab_fs::Matrix{Float64}`: Local stability matrix to store results (modified in place). - `plasma_eq::Equilibrium.PlasmaEquilibrium`: Plasma equilibrium data. + - `verbose::Bool`: Print progress messages. This function modifies `locstab_fs` in place with: @@ -38,14 +36,14 @@ This function modifies `locstab_fs` in place with: - Column 4: Delta Prime (Δ') """ function compute_ballooning_stability!( - ctrl::ForceFreeStatesControl, locstab_fs::Matrix{Float64}, plasma_eq::Equilibrium.PlasmaEquilibrium; theta_k::Float64=0.0, - compute_delta_prime::Bool=true + compute_delta_prime::Bool=true, + verbose::Bool=false ) - if ctrl.verbose + if verbose println("Evaluating local high-n ballooning stability...") end @@ -72,22 +70,22 @@ function compute_ballooning_stability!( end end - if ctrl.verbose + if verbose println("Ballooning analysis complete.") end end """ - compute_local_stability(ctrl, plasma_eq) -> CubicSeriesInterpolant + compute_local_stability(plasma_eq; verbose=false) -> CubicSeriesInterpolant Local stability profile spline over `plasma_eq.profiles.xs`, with the columns filled by [`compute_ballooning_stability!`](@ref): 1 = `D_I·ψ`, 2 = `D_R·ψ`, 4 = ballooning `Δ'`. """ -function compute_local_stability(ctrl::ForceFreeStatesControl, plasma_eq::Equilibrium.PlasmaEquilibrium) +function compute_local_stability(plasma_eq::Equilibrium.PlasmaEquilibrium; verbose::Bool=false) xs = plasma_eq.profiles.xs locstab_fs = zeros(Float64, length(xs), 5) - compute_ballooning_stability!(ctrl, locstab_fs, plasma_eq) + compute_ballooning_stability!(locstab_fs, plasma_eq; verbose=verbose) return cubic_interp(xs, Series(locstab_fs); extrap=ExtendExtrap()) end @@ -669,10 +667,10 @@ ballooning evaluator, and returns delta-prime and `det(d0bar)` `Di` maps. function scan_delta_prime_map( psi_idx::Int, plasma_eq::Equilibrium.PlasmaEquilibrium; - ctrl::ForceFreeStatesControl=ForceFreeStatesControl(; verbose=false), theta_k::Float64=0.0, s_scales::AbstractVector{<:Real}, - alpha_scales::AbstractVector{<:Real} + alpha_scales::AbstractVector{<:Real}, + verbose::Bool=false ) ref = salpha_reference(psi_idx, plasma_eq) @@ -697,7 +695,7 @@ function scan_delta_prime_map( theta_k=theta_k ) catch err - if ctrl.verbose + if verbose @warn "s-alpha scan point failed" psi_idx is ia err end nothing @@ -841,7 +839,7 @@ function critical_ballooning_alpha( end """ - ballooning_alpha_boundary(ctrl, plasma_eq; theta_k=0.0, n_scan=24) + ballooning_alpha_boundary(plasma_eq; theta_k=0.0, n_scan=24, verbose=false) Profile driver for the BALOO-style ballooning stability diagram. Loops over flux surfaces returning the experimental pressure gradient `alpha` (from @@ -861,10 +859,10 @@ far-edge surfaces dominate the scan cost. Per-surface failures, skipped surfaces surfaces with no boundary within range are returned as `NaN`. """ function ballooning_alpha_boundary( - ctrl::ForceFreeStatesControl, plasma_eq::Equilibrium.PlasmaEquilibrium; theta_k::Float64=0.0, - n_scan::Int=24 + n_scan::Int=24, + verbose::Bool=false ) xs = plasma_eq.profiles.xs npsi = length(xs) @@ -878,7 +876,7 @@ function ballooning_alpha_boundary( alpha[i] = salpha_reference(i, plasma_eq).alpha_ref alpha_critical[i] = critical_ballooning_alpha(i, plasma_eq; theta_k=theta_k, n_scan=n_scan).alpha_crit catch err - if ctrl.verbose + if verbose @warn "ballooning alpha boundary failed" psi_idx = i exception = err end end @@ -910,7 +908,7 @@ function second_critical_ballooning_alpha( end """ - ballooning_alpha_boundaries(ctrl, plasma_eq; theta_k=0.0, max_alpha_scale=8.0, n_scan=24) + ballooning_alpha_boundaries(plasma_eq; theta_k=0.0, max_alpha_scale=8.0, n_scan=24, verbose=false) Profile driver returning the experimental pressure gradient `alpha`, the first stability boundary `alpha_critical1` (lowest Δ' zero), and the second stability boundary @@ -924,11 +922,11 @@ smooth boundary curves. Surfaces outside `BALLOONING_SCAN_PSI_WINDOW` are skippe (returned as `NaN`), like per-surface failures. """ function ballooning_alpha_boundaries( - ctrl::ForceFreeStatesControl, plasma_eq::Equilibrium.PlasmaEquilibrium; theta_k::Float64=0.0, max_alpha_scale::Float64=8.0, - n_scan::Int=24 + n_scan::Int=24, + verbose::Bool=false ) xs = plasma_eq.profiles.xs npsi = length(xs) @@ -946,7 +944,7 @@ function ballooning_alpha_boundaries( alpha_critical1[i] = cr.alphas[1] length(cr.alphas) >= 2 && (alpha_critical2[i] = cr.alphas[2]) catch err - if ctrl.verbose + if verbose @warn "ballooning alpha boundaries failed" psi_idx=i exception=err end end @@ -1002,7 +1000,7 @@ function ballooning_qprime_crossings( end """ - ballooning_qprime_boundaries(ctrl, plasma_eq; theta_k=0.0, min_qprime_scale=-2.0, max_qprime_scale=4.0, n_scan=24) + ballooning_qprime_boundaries(plasma_eq; theta_k=0.0, min_qprime_scale=-2.0, max_qprime_scale=4.0, n_scan=24, verbose=false) Profile driver returning the experimental shear profile `qprime` (`dq/dpsi_norm`), the critical shear `qprime_critical1` below which the surface is ballooning unstable at the @@ -1012,12 +1010,12 @@ normalized flux `psi`, from [`ballooning_qprime_crossings`](@ref) at each surfac Arrays contain `NaN` where no boundary exists. """ function ballooning_qprime_boundaries( - ctrl::ForceFreeStatesControl, plasma_eq::Equilibrium.PlasmaEquilibrium; theta_k::Float64=0.0, min_qprime_scale::Float64=-2.0, max_qprime_scale::Float64=4.0, - n_scan::Int=24 + n_scan::Int=24, + verbose::Bool=false ) xs = plasma_eq.profiles.xs npsi = length(xs) @@ -1037,7 +1035,7 @@ function ballooning_qprime_boundaries( qprime_critical1[i] = cr.qprimes[1] length(cr.qprimes) >= 2 && (qprime_critical2[i] = cr.qprimes[2]) catch err - if ctrl.verbose + if verbose @warn "ballooning qprime boundaries failed" psi_idx=i exception=err end end @@ -1047,7 +1045,7 @@ function ballooning_qprime_boundaries( end """ - ballooning_delta_prime_map(ctrl, plasma_eq; theta_k=0.0, max_alpha_scale=8.0, n_alpha=61, max_surfaces=40) + ballooning_delta_prime_map(plasma_eq; theta_k=0.0, max_alpha_scale=8.0, n_alpha=61, max_surfaces=40, verbose=false) Map the ballooning Δ' over the (ψ_N, α) plane at fixed magnetic shear: at each flux surface, scan the pressure-gradient scaling from 0 to `max_alpha_scale` times the @@ -1058,12 +1056,12 @@ Returns a NamedTuple `(psi, alpha_scales, alpha_ref, delta_prime)` where `delta_ is `[n_surf × n_alpha]` with `NaN` on failed surfaces; physical `α = alpha_ref * scale`. """ function ballooning_delta_prime_map( - ctrl::ForceFreeStatesControl, plasma_eq::Equilibrium.PlasmaEquilibrium; theta_k::Float64=0.0, max_alpha_scale::Float64=8.0, n_alpha::Int=61, - max_surfaces::Int=40 + max_surfaces::Int=40, + verbose::Bool=false ) xs = plasma_eq.profiles.xs idx_all = [i for i in eachindex(xs) if xs[i] <= 1.0] @@ -1077,11 +1075,11 @@ function ballooning_delta_prime_map( for (k, i) in enumerate(idx) try - res = scan_delta_prime_map(i, plasma_eq; ctrl=ctrl, theta_k=theta_k, s_scales=[1.0], alpha_scales=alpha_scales) + res = scan_delta_prime_map(i, plasma_eq; theta_k=theta_k, s_scales=[1.0], alpha_scales=alpha_scales, verbose=verbose) alpha_ref[k] = res.reference.alpha_ref delta_prime[k, :] .= vec(res.delta_prime) catch err - if ctrl.verbose + if verbose @warn "ballooning delta-prime map failed" psi_idx=i exception=err end end @@ -1091,7 +1089,7 @@ function ballooning_delta_prime_map( end """ - ballooning_qprime_delta_prime_map(ctrl, plasma_eq; theta_k=0.0, min_qprime_scale=-2.0, max_qprime_scale=4.0, n_qprime=61, max_surfaces=40) + ballooning_qprime_delta_prime_map(plasma_eq; theta_k=0.0, min_qprime_scale=-2.0, max_qprime_scale=4.0, n_qprime=61, max_surfaces=40, verbose=false) Map the ballooning Δ' over the (ψ_N, q') plane at the fixed experimental pressure gradient: at each flux surface, scan the magnetic shear scaling from `min_qprime_scale` @@ -1104,13 +1102,13 @@ Returns a NamedTuple `(psi, qprime_scales, qprime_ref, delta_prime)` where `q' = qprime_ref * scale` in `dq/dpsi_norm` units. """ function ballooning_qprime_delta_prime_map( - ctrl::ForceFreeStatesControl, plasma_eq::Equilibrium.PlasmaEquilibrium; theta_k::Float64=0.0, min_qprime_scale::Float64=-2.0, max_qprime_scale::Float64=4.0, n_qprime::Int=61, - max_surfaces::Int=40 + max_surfaces::Int=40, + verbose::Bool=false ) xs = plasma_eq.profiles.xs idx_all = [i for i in eachindex(xs) if xs[i] <= 1.0] @@ -1124,11 +1122,11 @@ function ballooning_qprime_delta_prime_map( for (k, i) in enumerate(idx) try - res = scan_delta_prime_map(i, plasma_eq; ctrl=ctrl, theta_k=theta_k, s_scales=qprime_scales, alpha_scales=[1.0]) + res = scan_delta_prime_map(i, plasma_eq; theta_k=theta_k, s_scales=qprime_scales, alpha_scales=[1.0], verbose=verbose) qprime_ref[k] = res.reference.qprime_norm_ref delta_prime[k, :] .= vec(res.delta_prime) catch err - if ctrl.verbose + if verbose @warn "ballooning qprime delta-prime map failed" psi_idx=i exception=err end end diff --git a/src/LocalStability/LocalStability.jl b/src/LocalStability/LocalStability.jl new file mode 100644 index 000000000..eb96780e2 --- /dev/null +++ b/src/LocalStability/LocalStability.jl @@ -0,0 +1,17 @@ +module LocalStability + +# Local high-n stability: Mercier D_I, resistive interchange D_R, and the ballooning Δ' +# scans. Depends only on Equilibrium (plus math libraries) — no stability-solver state. +using LinearAlgebra +using FFTW +using OrdinaryDiffEq +using FastInterpolations +using StaticArrays: SVector + +import ..Equilibrium + +include("Ballooning.jl") + +export compute_local_stability, compute_ballooning_stability!, ballooning_alpha_boundary, ballooning_alpha_boundaries + +end diff --git a/test/runtests.jl b/test/runtests.jl index 76eca69d1..7e0e46d5d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,6 +3,7 @@ using Pkg using GeneralizedPerturbedEquilibrium.Vacuum using GeneralizedPerturbedEquilibrium.Equilibrium using GeneralizedPerturbedEquilibrium.ForceFreeStates +using GeneralizedPerturbedEquilibrium.LocalStability using GeneralizedPerturbedEquilibrium.ForcingTerms using GeneralizedPerturbedEquilibrium.PerturbedEquilibrium using GeneralizedPerturbedEquilibrium.Utilities diff --git a/test/runtests_resist_eval.jl b/test/runtests_resist_eval.jl index 58aee3d9b..8e334a67a 100644 --- a/test/runtests_resist_eval.jl +++ b/test/runtests_resist_eval.jl @@ -3,6 +3,7 @@ using GeneralizedPerturbedEquilibrium.Equilibrium using GeneralizedPerturbedEquilibrium.ForceFreeStates using GeneralizedPerturbedEquilibrium.ForceFreeStates: SingType, ResistGeometry + using GeneralizedPerturbedEquilibrium.LocalStability using GeneralizedPerturbedEquilibrium.Utilities using GeneralizedPerturbedEquilibrium.InnerLayer using GeneralizedPerturbedEquilibrium.Tearing: build_ggj_inputs @@ -44,7 +45,7 @@ # scan stores in locstab[:,1]. Build it on the radial grid, interpolate # to a few surface ψ values, and check against the GGJ reconstruction. xs = equil.profiles.xs - di_ref = Float64[ForceFreeStates.prepare_ballooning_coefficients(i, equil).di for i in eachindex(xs)] + di_ref = Float64[LocalStability.prepare_ballooning_coefficients(i, equil).di for i in eachindex(xs)] di_spline = cubic_interp(xs, di_ref) dq = deriv_view(equil.profiles.q_spline, 1) From 480d195486e0303be47ba18909abd07a174cc428 Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Sat, 15 Aug 2026 15:09:20 -0400 Subject: [PATCH 3/3] DOCS - FIX - Correct integrator solver claims to Vern9 and cite DCON/STRIDE papers --- docs/src/stability.md | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/docs/src/stability.md b/docs/src/stability.md index abce58c67..9d74a6460 100644 --- a/docs/src/stability.md +++ b/docs/src/stability.md @@ -52,8 +52,9 @@ driver returns dense displacement profiles, the Riccati driver returns the inter ### Forward integration -`forward_eulerlagrange_integration` is the baseline driver. It integrates the EL ODE directly -in ``(U_1, U_2)`` using Tsit5 with adaptive step control. Near each rational surface the +`forward_eulerlagrange_integration` is the baseline driver — our implementation of the +standard DCON radial integration [Glasser 2016]. It integrates the EL ODE directly +in ``(U_1, U_2)`` using the adaptive 9th-order `Vern9` solver. Near each rational surface the columns of ``U_2`` that correspond to resonant modes are zeroed via Gaussian reduction (GR), keeping the solution bounded; `transform_u!` undoes the reduction at the end, so `u_store` comes back dense in the axis (Euler-Lagrange) basis. This is the reference path for @@ -68,7 +69,9 @@ integrator = "forward" ### Riccati integration -`riccati_eulerlagrange_integration` (the default) decomposes the radial domain into +`riccati_eulerlagrange_integration` (the default) is our implementation of the STRIDE +approach [Glasser 2018b], built on the dual Riccati reformulation [Glasser 2018a]. It +decomposes the radial domain into independent chunks, integrates each chunk's fundamental-matrix (FM) propagator in parallel using `Threads.@threads`, then multiplies the propagators in order and applies each singular-surface crossing serially. It is the only driver that produces the inter-surface @@ -88,7 +91,7 @@ w = Q - \bar{K}S. than integrating the quadratic Riccati ODE directly (which blows up when ``|S|`` is large), the code integrates the linear EL system with `sing_der!` as the RHS and recovers ``S = U_1 U_2^{-1}`` via periodic renormalization — an approach that is mathematically -equivalent to O(Δψ) but uses the ODE solver's full 5th-order accuracy. Renormalization is +equivalent to O(Δψ) but uses `Vern9`'s full 9th-order accuracy. Renormalization is triggered whenever ``\max(|U_1|)`` or ``\max(|U_2|)`` exceeds the threshold `ucrit`, and is forced at the end of each chunk. At singular surface crossings, `riccati_cross_ideal_singular_surf!` applies the small-asymptotic matching directly in column