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..de545c427 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,8 +666,15 @@ manual smoke: run the 4-line UX from the Context section in a REPL against ## 10. Progress -- [ ] PR 1 — `refactor/riccati-unification` -- [ ] PR 2 — `refactor/local-stability-module` +- [ ] 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` — **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/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 16428c468..4c4f4f9dd 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 c31452c2d..913b48354 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/rational_psi"]), read(f["SingularSurfaces/GalerkinDeltaPrime/rational_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 bd3875c23..54e8ce105 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 dcb5204c9..1854fce57 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 b8b539a18..66cca17a5 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/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/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 ae96f9aed..9fb042e7d 100644 --- a/docs/src/stability.md +++ b/docs/src/stability.md @@ -45,66 +45,73 @@ 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 — 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. 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) 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 +``\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 `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 +`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 +133,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 +231,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 +241,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 @@ -272,7 +273,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 @@ -307,8 +308,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 +327,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/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/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..18df054b9 --- /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 a0e0a8396..feebc0c7f 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..f41d9ed0a --- /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/rational_count" +type = "int_scalar" +extract = "value" +label = "# singular surfaces" +noise_threshold = 0 +order = 30 + +[quantities.sing_psi] +h5path = "SingularSurfaces/rational_psi" +type = "real_vector" +extract = "all_real" +label = "singular psi locations" +noise_threshold = 1e-8 +order = 31 + +[quantities.sing_q] +h5path = "SingularSurfaces/rational_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/q_axis" +type = "real_scalar" +extract = "value" +label = "q0" +noise_threshold = 1e-12 +order = 60 + +[quantities.q95] +h5path = "Equilibrium/q_95" +type = "real_scalar" +extract = "value" +label = "q95" +noise_threshold = 1e-12 +order = 61 + +[quantities.betan] +h5path = "Equilibrium/beta_N" +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/Analysis/ForceFreeStates.jl b/src/Analysis/ForceFreeStates.jl index f276c058a..99dd60e2b 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/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/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/ForceFreeStates/ForceFreeStatesStructs.jl b/src/ForceFreeStates/ForceFreeStatesStructs.jl index 5d7cbb920..50eed7da0 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 95f64277d..0de2b07a5 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 @@ -65,13 +70,13 @@ 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! 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 @@ -336,9 +341,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 @@ -578,9 +583,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/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/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.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_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 34ccc7256..7e7737e06 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_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) 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 f2fc6d1f0..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 cef6ad546..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