diff --git a/REFACTOR_PLAN.md b/REFACTOR_PLAN.md deleted file mode 100644 index b842ec216..000000000 --- a/REFACTOR_PLAN.md +++ /dev/null @@ -1,1117 +0,0 @@ -# GPEC Julia — ForceFreeStates modularization & `solve(eq, integrator)` API -## Complete multi-PR implementation plan - -> **NOTE FOR ALL DEVELOPERS (read this first).** -> This document is the agreed, in-progress plan for a refactor of the -> ForceFreeStates ↔ PerturbedEquilibrium interface and the top-level driver, delivered -> as THREE pull requests: #381 (integrator unification), #387 (LocalStability), and one -> combined "interface PR" whose three commits carry what were originally planned as -> PRs 3-5 (the stack was collapsed once it became clear reviews would batch at the end). It is -> committed directly to `develop` (deliberately, as documentation only — no code -> 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` 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 -> not re-litigate them in PR review unless you find a factual error. -> - This file is temporary: it gets checked off as PRs merge and is deleted once the -> sequence is implemented and vetted. - ---- - -## 1. Context - -GPEC's pipeline runs through one ~520-line monolith (`main_from_inputs`, -`src/GeneralizedPerturbedEquilibrium.jl:164`) interleaving equilibrium setup, mode -resolution, singular-surface handling, local stability, matrix assembly, three -integrator code paths, the Δ′ BVP, the Galerkin solve, HDF5 writing, and the -PE → KineticForces → SLAYER stages, communicating through ~8 loose objects. - -Target UX: - -```julia -eq = PlasmaEquilibrium("input.geqdsk"; jac_type="hamada") -ffs = solve(eq, Riccati(); nn=1, delta_mlow=8, delta_mhigh=8, vac_flag=true) -rmp = RMPField("coils.dat") -pe = perturbed_equilibrium(ffs, rmp) -# calculate_quantities(...) is OUT OF SCOPE (later deliverable) -``` - -### Decision record (settled — do NOT re-litigate) - -| # | Decision | -|---|---| -| D1 | Three integrators = three formalisms: **Forward** (serial EL; rename all misuses of "shooting"), **Riccati** (the STRIDE FM-chunk driver currently behind `use_parallel`), **Galerkin** (RDCON; becomes fully standalone). | -| D2 | Riccati uses whatever threads `julia -t` provides. Its ONLY tunable is **number of chunks** (`nchunks`). `parallel_threads` is deleted. Outputs must be independent of thread count ⇒ the auto chunk count derives from problem structure only, never `Threads.nthreads()`. | -| D3 | **No merging of two integration results.** `populate_dense_xi` + `_populate_dense_xi_via_serial_el!` + the standalone serial-Riccati path (`riccati_eulerlagrange_integration`) are deleted FIRST (PR 1). Riccati-fed PE warn-and-skips profile-based outputs PERMANENTLY (D14: riccati never produces full profiles); the separate `delta_mn` work (not in this plan) restores the resonant-coupling outputs — not the profile-based ones — from `delta_coil` + surface asymptotics. | -| D4 | Kinetic (`kinetic_factor > 0`) is Forward-only. `solve`/driver raises a clear error for Riccati+kinetic and Galerkin+kinetic. | -| D5 | One result struct **`ForceFreeStatesResult`**; optional fields are `Union{Nothing,T}`; consumers use a `require(...)` helper → `@warn` + skip. | -| D6 | Local stability (Ballooning.jl) → new top-level module **`LocalStability`**, depending only on Equilibrium (+ math deps). Only ctrl dependency is `verbose` → kwarg. | -| 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. 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). | -| D13 | Inner-layer matching runs INSIDE `solve` (a `ForceFreeStatesResult` is always a closed basis). `result.solution` holds THE solve's ξ solution product — a thin `SolutionProfiles` interchange type — whenever one exists: forward always; galerkin when matched (built directly from the match — the `gal_matched_odestate` OdeState shim is DELETED); riccati permanently `nothing` — STRIDE matching yields rational-surface data (`bpen`, `delta_mn`), never profiles (D14). Closure is explicit and universal: `result.closure ∈ (:ideal, :matched)` and `result.bpen` (msing × numpert_total; zeros under ideal closure) are always present — the landing pad for any matching implementation. No transitional arbitration API: additive gal is removed in the SAME PR that introduces the result (PR 3), so one run has at most one solution and nothing like `pe_solution` is ever needed. Matching config is integrator-agnostic: a `ResistiveMatch` object (swappable `InnerLayer` model + per-surface `eta/rho/rotation`, `gamma`, `ideal`) passed as a `match=` kwarg to `solve` (PR 5). STRIDE-side matching is a future PR; until then `match` with `Riccati()` errors "not yet implemented". The `gal_*` matching TOML keys are renamed/re-homed by that future PR, not by this stack. | -| D14 | Same physics ⇒ same field, same type, across integrators, organized by the three-class taxonomy in §9 (control surface / full profiles / rational-surface resonant data). Riccati NEVER produces full ξ/ξ′ profiles — `result.solution` is permanently `nothing` for it. The next-cycle work adds `delta_mn` to riccati AND galerkin: a bpen-like matrix encoding the jump in the pitch-resonant derivative of the solution at each rational surface, from outer-solution asymptotics (for riccati: recoverable from `delta_coil`); it yields the perturbed current and shielded resonant flux, and is what PE resonant coupling consumes from a Riccati run (class 2, not class 1). Forward `delta_mn` is NOT planned — no concrete route has been identified and there may be none. There is ONE Δ′/matching data type, unified IN THIS PR: `delta_prime` carries Δ′ matrix, raw D′, `delta_coil`, and the PEST-3 blocks, produced by riccati and galerkin alike — galerkin already computes the same physics content, today under `galerkin.*` fields and different HDF5 names; its Δ′ payload merges into `delta_prime` (fields a formalism doesn't produce stay empty/`nothing`). Control-surface energies (`wp`, `free_boundary`) target all three integrators (galerkin pending its δW implementation). SLAYER consumes the unified `delta_prime`, so riccati- and galerkin-fed SLAYER both work (this PR). SLAYER + GGJ behind one abstract inner-layer interface is a later pass. | - -### Verified code facts the workers must not re-derive - -- Dispatch today: `eulerlagrange_integration` (`src/ForceFreeStates/EulerLagrange.jl:151`): - `use_parallel` → `parallel_eulerlagrange_integration` (Riccati.jl:1647; returns - `(odet, propagators, chunks, S_at_surface_left)`), `use_riccati` → - `riccati_eulerlagrange_integration` (Riccati.jl:1312; being deleted), else - `serial_eulerlagrange_integration` (EulerLagrange.jl:172). -- The Δ′ BVP (`compute_delta_prime_matrix!`, Riccati.jl:274) is called once from - `src/GeneralizedPerturbedEquilibrium.jl:470-478`, only when propagators exist. - Active assembly = `_assemble_bvp_S_axis` (Riccati S states); `_assemble_bvp_FM_axis` - is a never-used fallback. `_solve_bvp_edge_coil` fills `intr.delta_coil` when - S-axis && `wv !== nothing`. -- Serial-Riccati `u_store` is NOT usable as ξ (renorm right-multiplications never - recorded/undone); it also leaves `u_store_el_basis == true` (foot-gun; dies with the - path). -- `populate_dense_xi` = re-run `serial_eulerlagrange_integration` and splice - (`_populate_dense_xi_via_serial_el!`, Riccati.jl:1930). The "Riccati-gauge ca needed - by SingularCoupling" comment there is STALE — **PE never reads `ca_l/ca_r`**. -- `balance_integration_chunks` (EulerLagrange.jl:79) sizes chunks with - `target_n = max(2*msing+3, 4*effective_threads, 8*(msing+1)+msing)` — the middle - term must go (D2). -- PE reads of OdeState: `u_store` (BOTH components), `du_store` (dense), `xi_s_store`, - `psi_store`, `q_store`, `step`, `du_store_populated`. NOT `ca_l/ca_r`, `crit_store`, - `edge_scan`. PE calls `materialize_derivative_stores!` itself - (`src/PerturbedEquilibrium/PerturbedEquilibrium.jl:88`) and discards the Bool. -- PE reads of `ForceFreeStatesInternal`: `nlow/nhigh/mlow/mhigh/mpert/npert/ - numpert_total`, `psilim`, `qlim`, `msing`, `sing[s].psifac/.q/.q1` only. PE - recomputes its own Green's functions via `Vacuum.compute_vacuum_response`. -- PE reads `metric.fourier_coeffs` only; `ffit.amats/bmats/cmats`, - `ffit.fmats_lower/kmats`, `ffit.kinetic_populated`, and calls - `ForceFreeStates.el_derivatives!`. -- PE sub-calc order/prereqs: `compute_plasma_response!` needs `wt0` + dense stores + - ffit A/B/C + `metric.fourier_coeffs`; `compute_singular_coupling_metrics!` needs - `wt0` + `intr.plasma_response` (from the response step) + boundary `u_store` + - Ξ/Ξ′ near surfaces (+ optional `inner_bpen`, same identity-at-edge basis). -- SLAYER (`src/Tearing/Runner/run_slayer.jl:367`) reads exactly `ffs_intr.sing` and - `ffs_intr.delta_prime_matrix` (+ `equil`, `dir_path` kwarg). -- Standalone vacuum wv: `VacuumInput(equil, ψ, mthvac, nzvac, mrange, nrange)` - (`src/Vacuum/DataTypes.jl:64`) → `compute_vacuum_response(inputs, wall).wv` — no - ODE state. `free_run` applies singfac scaling in place (`Free.jl:86-88`); - `galerkin_solve` consumes the ALREADY singfac-scaled wv and multiplies by `psio²` - (`Galerkin/GalerkinSolve.jl:124-129`). -- `EquilibriumConfig` is `@kwdef` (`src/Equilibrium/EquilibriumTypes.jl:43`) — - keyword construction works today; the Dict constructor is a filter/warn wrapper. -- Writers: FFS `write_outputs_to_HDF5` defined `GeneralizedPerturbedEquilibrium.jl:700`, - called once at `:490`. PE writer `src/PerturbedEquilibrium/Utils.jl:103`, called once - at `:618`. `write_imas` (`:1020`) reads `result.free_energies.et/.n_tor_idx` and - `result.intr.npert/.nlow`; tested in `test/runtests_imas.jl:122,152,175,183,191`. -- Rerun path: `build_inputs_from_h5` (`src/Rerun.jl:199`) returns a 7-tuple funneled - into `main_from_inputs`; it never reads the FFS flags by name (opaque dict). -- Deprecation machinery: `_drop_deprecated_keys!` + `_DEPRECATED_FFS_KEYS` / - `_DEPRECATED_EQUIL_KEYS` (`GeneralizedPerturbedEquilibrium.jl:73-86`), applied at - `:126`, `:184`, and `src/Rerun.jl:269`. Pre-commit hook `toml-no-deprecated-keys` - mirrors these lists — **update the hook regex whenever the lists change**. -- Tests: `test/runtests.jl:21-50` is a hard-coded include list (no globbing). - `use_riccati` appears in NO test and NO TOML. `riccati_eulerlagrange_integration` - called directly at `test/runtests_riccati.jl:115`. `use_parallel` toggles at - `test/runtests_eulerlagrange.jl:435`, `test/runtests_parallel_integration.jl:234-501`. - `populate_dense_xi` testsets at `test/runtests_parallel_integration.jl:389-490`. -- Docs: FFS `@autodocs` at `docs/src/stability.md:273-276` (Pages list includes - `Ballooning.jl`); `docs/src/ballooning.md` has NO autodocs block; - `docs/make.jl:46` `checkdocs=:exports`; nav at `docs/make.jl:27-45`. - `docs/development/architecture.md` module list is stale and needs updating anyway. -- Deps: CommonSolve NOT in `[deps]` (indirect in Manifest — add to `[deps]`+`[compat]`). - `using OrdinaryDiffEq` (`src/ForceFreeStates/ForceFreeStates.jl:8`) already brings - `solve` (== `CommonSolve.solve`) unqualified into FFS scope — new methods MUST be - defined via `import CommonSolve: solve` (adding methods to the same generic; the - existing unqualified `solve(prob, Vern9(); ...)` calls keep working). -- "shooting" rename scope: `EulerLagrange.jl:166` docstring; - `GeneralizedPerturbedEquilibrium.jl:580,582` comments; `Galerkin/GalerkinMatch.jl:240` - docstring; benchmarks labels (`benchmarks/compare_jbgradpsi_m2.jl`, - `scan_resistivity_m2.jl`, `scan_rotation_m2.jl`). **Do NOT rename**: the GGJ - inner-layer `:shooting` backend (`src/InnerLayer/GGJ/Shooting.jl`, `:ggj_shooting` - in `src/Tearing/Runner/Control.jl`), the ballooning-doc "shooting boundary" - (`docs/src/ballooning.md:847`), and the BVP shooting-propagator names - (`uShootR/uShootL`, `_build_S_axis_shooting_propagators`) — those are correct - shooting-method/STRIDE terminology. - ---- - -## 2. PR sequence overview - -Branch from `develop`, PR back into `develop`. **Every PR requires third-party human -review before merge — non-negotiable.** Run the regression harness once per PR and -report the table (differences are expected and get accepted knowingly; see D10). -All code must be JuliaFormatter-clean per `.JuliaFormatter.toml` before commit. - -| PR | Branch | Content | -|----|--------|---------| -| #381 | `refactor/riccati-unification` | Delete serial-Riccati + `populate_dense_xi` + `parallel_threads`; `integrator=` ctrl key; `nchunks` knob; thread-independent chunking; shooting→forward rename | -| #387 | `refactor/local-stability-module` | Extract Ballooning.jl → `LocalStability` module; drop ctrl dependency (stacked on #381) | -| interface PR | `refactor/forcefreestates-result` | ONE PR, three slice-pure commits: **(a)** §5 `ForceFreeStatesResult` + warn-and-skip consumers + standalone Galerkin; **(b)** §6 staged `main`; **(c)** §7 `solve` API (stacked on #387) | - -Commit discipline for the interface PR: commit boundaries now do the job PR boundaries -did — keep each commit slice-pure (fixes amend into the right slice before review -starts; ordinary follow-up commits after). Per-slice numerical isolation stays -verifiable via the harness with commit SHAs as refs. - ---- - -## 3. PR 1 — `refactor/riccati-unification` - -### 3.1 Control struct (`src/ForceFreeStates/ForceFreeStatesStructs.jl`) - -- DELETE fields + docstring entries: `use_riccati` (:297), `use_parallel` (:298), - `parallel_threads` (:290, docstring :258), `populate_dense_xi` (:299, docstring :259). -- ADD fields: - - `integrator::String = "riccati"` — `"forward" | "riccati" | "galerkin"` is - validated at dispatch (`"galerkin"` only becomes legal in PR 4; until then it - errors with "not yet a standalone integrator — use gal_flag"). - - `nchunks::Int = 0` — Riccati chunk-count target; `0` = auto (structure-derived). -- Validation (where `ctrl` is constructed is a splat; add checks at the top of - `eulerlagrange_integration`): error if `integrator == "riccati" && kinetic_factor > 0` - ("kinetic runs require integrator=\"forward\""); error on unknown integrator string. - -### 3.2 Integration code - -- `src/ForceFreeStates/EulerLagrange.jl`: - - `eulerlagrange_integration` dispatch: `integrator=="riccati"` → - `riccati_eulerlagrange_integration` (the renamed STRIDE driver), else forward. - - RENAME `serial_eulerlagrange_integration` → `forward_eulerlagrange_integration` - (keep the `verbose` kwarg; update the "Serial shooting branch" docstring at :166). - - `balance_integration_chunks` (:79): remove the `4 * effective_threads` term and - the `ctrl.parallel_threads` read (:90). New sizing: - `target_n = ctrl.nchunks > 0 ? max(ctrl.nchunks, 2*intr.msing + 3) : max(2*intr.msing + 3, 8*(intr.msing + 1) + intr.msing)` - — with `@warn` when an explicit `nchunks` is clamped up. NO `Threads.nthreads()` - anywhere in chunk sizing. -- `src/ForceFreeStates/Riccati.jl`: - - DELETE `riccati_eulerlagrange_integration` (:1312-1397) and - `_populate_dense_xi_via_serial_el!` (:1900-1980). - - RENAME `parallel_eulerlagrange_integration` → `riccati_eulerlagrange_integration` - (name is now free; update its docstring: "the Riccati/STRIDE integrator", drop the - populate_dense_xi paragraph and the `Enable via use_parallel` line). Remove the - `ctrl.populate_dense_xi && !ctrl.force_termination` block (:1681-1683). - - Thread pool: replace `bvp_threads = max(1, min(Threads.nthreads(), ctrl.parallel_threads))` - (:1653) with `Threads.nthreads()` used directly by `_run_parallel_bvp_phase!`; - per-thread proxies keep sizing by `Threads.maxthreadid()`. - - After the parallel path, `odet.u_store_el_basis` stays `false` (already set at - :1795) — this is now the permanent contract: Riccati's odet never claims EL basis. -- `src/PerturbedEquilibrium/SingularCoupling.jl:66-69`: update the hard `error()` - message (references `populate_dense_xi`) → "dense Ξ′ requires the Forward - integrator" (message only; the structural gate arrives in PR 3). -- `src/GeneralizedPerturbedEquilibrium.jl`: comments at :580/:582 ("shooting - solution") → "forward solution". Add the four removed keys to - `_DEPRECATED_FFS_KEYS` (:73). NOTE: with `use_parallel` warn-ignored, old TOMLs and - gpec.h5 replays (whose `gpec_toml_raw` embeds old keys) fall through to the default - `integrator="riccati"` — same physics path as before, so replays stay valid. -- `src/ForceFreeStates/Galerkin/GalerkinMatch.jl:240`: docstring "shooting - integrator's" → "forward integrator's". - -### 3.3 Pre-commit hook + TOML sweep - -- Update the `toml-no-deprecated-keys` pygrep pattern in `.pre-commit-config.yaml` - to include the four new deprecated keys. -- All 12 `examples/*/gpec.toml` + 4 `test/test_data/regression_*/gpec.toml` - (canonical annotation source = `examples/DIIID-like_ideal_example/gpec.toml` per - `docs/development/toml-conventions.md`): remove `use_parallel`, `parallel_threads`, - `populate_dense_xi` lines; add `integrator = "…"` with a convention-conform comment. - Assignment: - - `integrator = "forward"` for every deck with a `[PerturbedEquilibrium]` section or - `kinetic_factor > 0`: `DIIID-like_ideal_example`, `Solovev_ideal_example`, - `Solovev_kinetic_NTV_example`, `Solovev_kinetic_calculated_example`, - `a10_kinetic_example`, and the 4 regression fixtures. - - `integrator = "riccati"` for Δ′/stability-only decks: `Solovev_ideal_example_multi_n`, - `Solovev_ideal_example_3D`, `LAR_beta_scan`, `LAR_epsilon_scan`, - `DIIID-like_SLAYER_example` (needs `delta_prime_matrix`). - - gal decks (`DIIID-like_gal_resistive*`, `LAR_*_match_test`) keep `gal_flag=true` - and use `integrator = "riccati"` (gal stays additive until PR 4). - - NEW example `examples/DIIID-like_riccati_deltaprime_example/` (copy of - DIIID-like_ideal minus `[PerturbedEquilibrium]`/`[ForcingTerms]`, with - `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 - `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`); - `benchmarks/compare_jbgradpsi_m2.jl`, `scan_resistivity_m2.jl`, - `scan_rotation_m2.jl`: label text "shooting" → "forward". - -### 3.4 Tests - -- `test/runtests_riccati.jl`: replace the direct call at :115 with the renamed driver - (`FFS.riccati_eulerlagrange_integration(ctrl, equil, ffit, intr)` now returns the - 4-tuple — destructure) or route via `ctrl` with `integrator="riccati"`. Keep the - energy-agreement assertion vs the forward path (:127). Delete the "(S, I) identity" - check tied to the deleted serial path (:151) or re-target it to the driver's - outer-region state. -- `test/runtests_parallel_integration.jl`: `use_parallel` toggles → `integrator=` - strings; DELETE the `populate_dense_xi` testsets (:389-490); keep/extend the sparse - u_store control test as "riccati leaves sparse u_store". ADD a unit test that - `balance_integration_chunks` output is identical for any `Threads.nthreads()` - (call with same inputs; assert no thread dependence — pure function now) and that - `nchunks` steering works and clamps with a warning. -- `test/runtests_eulerlagrange.jl:435`: `use_parallel=false` → `integrator="forward"`. -- `test/runtests_rerun_from_h5.jl`: fixture decks pick up new keys automatically; the - replay of PRE-refactor h5 files exercises the deprecated-key warn path — assert the - warning fires once (cheap regression for the deprecation mechanism). - -### 3.5 Docs - -- `docs/src/stability.md`: rewrite the `use_riccati`/`use_parallel` passages - (:61, :90, :243, :310) around `integrator = "forward"|"riccati"` and `nchunks`. -- `ForceFreeStatesControl` docstring: new entries for `integrator`, `nchunks`. - -### 3.6 Verification - -1. `julia --project=. test/runtests.jl test/runtests_riccati.jl test/runtests_parallel_integration.jl test/runtests_eulerlagrange.jl test/runtests_rerun_from_h5.jl test/runtests_fullruns.jl` -2. Full suite: `julia --project=. -e 'using Pkg; Pkg.activate("."); include("test/runtests.jl")'` -3. Regression harness: `julia --project=regression-harness regression-harness/regress.jl --cases diiid_n1,solovev_n1 --refs develop,local` — expected: forward-deck quantities unchanged vs develop where the deck previously ran `use_parallel+populate_dense_xi` (ξ was already forward-produced; Δ′ dataset disappears from forward decks — flagged, accepted); riccati decks match develop's parallel path bit-for-bit. -4. Docs build: `julia --project=. build_docs_local.jl`. - ---- - -## 4. PR 2 — `refactor/local-stability-module` - -### 4.1 Module extraction - -- `git mv src/ForceFreeStates/Ballooning.jl src/LocalStability/Ballooning.jl`; create - `src/LocalStability/LocalStability.jl`: - ```julia - module LocalStability - using LinearAlgebra, FFTW, OrdinaryDiffEq, FastInterpolations - using StaticArrays: SVector - import ..Equilibrium - include("Ballooning.jl") - export compute_local_stability, compute_ballooning_stability!, - ballooning_alpha_boundary, ballooning_alpha_boundaries - end - ``` - (Exact `using` set = what Ballooning.jl actually touches; it currently free-rides on - the FFS module imports — FFTW via `FFTW.fft/ifft`, OrdinaryDiffEq via - `ODEProblem/solve/DP5/ReturnCode`, FastInterpolations via - `cubic_interp/Series/PeriodicBC/CubicFit/ExtendExtrap/integrate/cumulative_integrate`.) -- Top module (`src/GeneralizedPerturbedEquilibrium.jl`): `include` + `import .LocalStability` - + `export LocalStability` after Equilibrium, before Vacuum. Remove the ballooning - names from the FFS import line (:67). -- Signature changes (drop the `ForceFreeStatesControl` argument everywhere; it only - supplied `verbose`): - - `compute_local_stability(plasma_eq; verbose=false)` - - `compute_ballooning_stability!(locstab_fs, plasma_eq; theta_k=0.0, compute_delta_prime=true, verbose=false)` - - `ballooning_alpha_boundary(plasma_eq; theta_k=0.0, n_scan=24, verbose=false)` - - `ballooning_alpha_boundaries`, `ballooning_qprime_boundaries`, - `ballooning_delta_prime_map`, `ballooning_qprime_delta_prime_map`, - `scan_delta_prime_map` — same pattern (`ctrl::ForceFreeStatesControl=...` kwarg in - `scan_delta_prime_map` becomes `verbose::Bool=false`). -- `src/ForceFreeStates/ForceFreeStates.jl`: remove `include("Ballooning.jl")` (:27). - FFS keeps `local_stability_flag` in its control struct for now (driver reads it); - a `[LocalStability]` TOML section is future work, out of scope. -- Driver call sites (`GeneralizedPerturbedEquilibrium.jl:338,340`): - `LocalStability.compute_local_stability(equil; verbose=ctrl.verbose)` / - `LocalStability.ballooning_alpha_boundary(equil; verbose=ctrl.verbose)`. -- Cross-check test `test/runtests_resist_eval.jl:47` - (`ForceFreeStates.prepare_ballooning_coefficients`) → `LocalStability.…`. - -### 4.2 Docs - -- `docs/src/stability.md:273-276`: remove `"Ballooning.jl"` from Pages. -- `docs/src/ballooning.md`: append an `@autodocs` block - (`Modules = [GeneralizedPerturbedEquilibrium.LocalStability]`) — required because - `checkdocs=:exports` (`docs/make.jl:46`) now sees the new exports. -- `docs/development/architecture.md`: add LocalStability to the module list and the - dependency tree (the list is stale anyway; fix minimally — add LocalStability, note - it depends only on Equilibrium). - -### 4.3 Verification - -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` (`LocalStability/*` datasets must be identical; note the #363 group name already matches the new module name). - ---- - -## 5. Interface PR, commit (a) — result struct, consumers, standalone Galerkin - -### 5.1 New file `src/ForceFreeStates/Result.jl` (included from ForceFreeStates.jl) - -Reuse existing types wholesale (`SingType`, `OdeState`, `FreeBoundaryResult`, -`GalerkinResult`, `FourFitVars`, `MetricData`, `EdgeScanState`); new types are -`DeltaPrimeData`, `SolutionProfiles`, and the result itself: - -```julia -"Δ′ outputs of the Riccati STRIDE BVP (moved off ForceFreeStatesInternal at result-build time)." -struct DeltaPrimeData - matrix::Matrix{ComplexF64} # msing×msing PEST3 Δ′ (was intr.delta_prime_matrix) - raw::Matrix{ComplexF64} # 2msing×2msing side-major D′ (was intr.delta_prime_raw) - coil::Matrix{ComplexF64} # 2msing×numpert_total edge coil response (was intr.delta_coil) -end - -"The solve's ξ solution, in the exact shape PerturbedEquilibrium consumes. Field names - mirror the OdeState store subset so PE internals change minimally." -struct SolutionProfiles - basis::Symbol # :el_axis (forward) | :gal_native (matched galerkin) - step::Int # number of stored radial nodes - psi_store::Vector{Float64} - q_store::Vector{Float64} - u_store::Array{ComplexF64,4} # (N, N, 2, step) — Ξ_ψ and conjugate momentum - du_store::Array{ComplexF64,3} # (N, N, step) dΞ_ψ/dψ, ALWAYS populated - xi_s_store::Array{ComplexF64,3} # (N, N, step) Ξ_s, ALWAYS populated -end - -struct ForceFreeStatesResult - integrator::Symbol # :forward | :riccati | :galerkin - control::ForceFreeStatesControl # provenance snapshot (carries mthvac, verbose, …) - equil::Equilibrium.PlasmaEquilibrium # possibly re-formed (two-pass) - # mode space & domain (copied out of intr — plain immutable data) - mlow::Int; mhigh::Int; mpert::Int - nlow::Int; nhigh::Int; npert::Int; numpert_total::Int - psilow::Float64; psilim::Float64; qlim::Float64; q1lim::Float64 - dir_path::String - wall_settings::Vacuum.WallShapeSettings - # assembly products (always present) - metric::MetricData - ffit::FourFitVars - surfaces::Vector{SingType} # alias of intr.sing (ua/restype/α live here) - kinetic::@NamedTuple{kmsing::Int, kinsing::Vector{SingType}, scan_psi::Vector{Float64}, scan_cond::Vector{Float64}, scan_threshold::Float64} - # closure of the basis at the rationals (D13) — ALWAYS present - closure::Symbol # :ideal (jump condition imposed) | :matched (inner layer) - bpen::Matrix{ComplexF64} # (msing × numpert_total) penetrated resonant field; zeros under :ideal - # per-integrator products (presence == capability) - solution::Union{Nothing,SolutionProfiles} # THE solve's ξ solution; nothing when none exists (riccati; unmatched gal) - diagnostics::Union{Nothing,OdeState} # the integrator's raw odet (crit, edge scan, ψ trace, ca); writer-only - wp::Union{Nothing,Matrix{ComplexF64}} # fixed-boundary plasma energy W_p at psilim; present for any EL sweep even with vac_flag=false (aliases free_boundary.wp when free_run ran) - free_boundary::Union{Nothing,FreeBoundaryResult} - delta_prime::Union{Nothing,DeltaPrimeData} - galerkin::Union{Nothing,GalerkinResult} -end -``` - -Contract (D13 — final, no transitional states): -- Forward → `solution` = `SolutionProfiles(:el_axis, …)` aliasing the odet's stores (zero - copy), `diagnostics` = the same odet, `closure = :ideal`, `bpen` = zeros. -- Riccati → `solution = nothing` PERMANENTLY (chunk-endpoint states are not a ξ solution, - and no reconstruction is planned; the future STRIDE matching populates - `closure = :matched`, `bpen`, and `delta_mn` — rational-surface data, never profiles), - `diagnostics` = its odet (ψ/q/crit/edge scan/ca are valid), `closure = :ideal`. -- Galerkin, matched → `solution` = `SolutionProfiles(:gal_native, …)` built DIRECTLY from - `GalerkinResult.match`/`solution` (drop `issing` points, analytic Ξ′, `compute_node_xi_s!` - for Ξ_s — the useful guts of the deleted `gal_matched_odestate`, minus the OdeState - costume), `diagnostics = nothing`, `closure = :matched` (`:ideal` under `gal_ideal_flag`), - `bpen = galerkin.match.bpen`. -- Galerkin, unmatched → `solution = nothing` (raw homogeneous gal columns are not a driven - response basis), `closure = :ideal`. - -There is NO `pe_solution` and NO stored-basis arbitration: additive gal is removed in this -PR (§5.3), so a run has at most one solution and PE reads `result.solution` directly. - -Helpers (same file): - -```julia -"Warn-and-skip gate: true iff the optional `field` is populated." -function require(result::ForceFreeStatesResult, field::Symbol, calc::AbstractString) - getfield(result, field) === nothing || return true - @warn "Skipping $calc: `$field` was not produced by the $(result.integrator) integrator" - return false -end - -"Specialized message for the ξ-solution gate." -require_solution(result, calc) = result.solution !== nothing ? true : - (@warn "Skipping $calc: no ξ solution — dense profiles require a Forward (or matched Galerkin) run; " * - "this result came from the $(result.integrator) integrator"; false) - -"Assemble the published result once the solve is finished." -build_result(integrator, ctrl, equil, intr, metric, ffit, odet, free_energies, gal_data) -> ForceFreeStatesResult -``` - -`build_result` responsibilities (the ONLY place with assembly logic): -- `delta_prime` from the intr Δ′ fields when non-empty; `free_boundary = free_energies`; - `galerkin = gal_data`; `diagnostics = odet`. -- Forward: call `materialize_derivative_stores!(odet, …)` HERE (moving the call out of the - writer and PE — one site, always-populated `du_store`/`xi_s_store`), then wrap the stores - in `SolutionProfiles(:el_axis, …)`. -- Matched gal: build `SolutionProfiles(:gal_native, …)` from the match (see contract above). -- `closure = (gal_data !== nothing && gal_data.match !== nothing && !ctrl.gal_ideal_flag) ? :matched : :ideal`; - `bpen` = match bpen or zeros(msing, numpert_total). -`ForceFreeStatesInternal` stays as internal scratch during the solve; it no longer -crosses module boundaries after `build_result`. - -### 5.2 Consumers - -- **PE** (`src/PerturbedEquilibrium/PerturbedEquilibrium.jl`): new signature - `compute_perturbed_equilibrium(result::ForceFreeStates.ForceFreeStatesResult, ft_ctrl, ctrl, intr)` - (drop `equil/odet/wt0/mthvac/ffs_intr/metric/ffit` — all read off `result`). - Internals: - - `initialize_mode_arrays!` reads mode fields from `result`. - - PE's working solution IS `result.solution::SolutionProfiles` (never an OdeState; PE - internals re-type from `OdeState` to `SolutionProfiles` — field names match, so the - change is annotations, not logic). No materialize call in PE: `du_store`/`xi_s_store` - arrive populated. - - Response step: `require(result, :free_boundary, "plasma response") && - require_solution(result, "plasma response")` else skip. - - Coupling step: same two gates + existing internal `plasma_response` gate. - - All `ffs_intr.X` reads → `result.X`; `wt0` → `result.free_boundary.wt0`; - `mthvac` → `result.control.mthvac`. - - `pe_intr.odet_from_gal` ↔ `result.solution.basis == :gal_native`; - `pe_intr.inner_bpen = result.bpen` (driver; the gal special-case `if` is deleted). -- **FFS HDF5 writer**: re-signature to - `write_outputs_to_HDF5(result; git_version, inputs, forcing_modes, locstab, ballooning_boundary)` - — body is today's `:700-991` with `ctrl/equil/intr/odet/free_energies/ffit/gal_data` - spelled `result.*`; every group that came from an optional field gets the existing - empty-array fallback (already the pattern for FreeBoundaryStability). Δ′ datasets - read from `result.delta_prime`. Solution-adjacent datasets split by source: - `ForwardIntegration/xi_psi|u2|dxi_psi|xi_s` from `result.solution` when - `basis == :el_axis` (empty otherwise — the gal-native solution is already persisted - under the Galerkin group); `psi|q|nstep|nstep_total|crit`, `SingularSurfaces/ca_*`, - and `EdgeScan/*` from `result.diagnostics` when present (empty otherwise). Output is - byte-identical for every forward/riccati deck; the four gal decks become gal-only - files (§5.3). **Dataset names/paths unchanged** (D11). -- **SLAYER**: `Runner.run_slayer(result, control; dir_path)` — reads - `result.surfaces`, `result.delta_prime === nothing ? empty : result.delta_prime.matrix`, - `result.equil`. Keep a thin internal method for the old `(equil, sing, dpm)` shape if - convenient; update `test/runtests_slayer_runner.jl`. -- **`write_imas`** + **`main` return value**: `main`/`main_from_inputs` return - `(; ffs::ForceFreeStatesResult, pe, slayer)` (pe/slayer possibly `nothing`). - `write_imas(dd, ret)` reads `ret.ffs.free_boundary` (skip+warn if `nothing`), - `ret.ffs.npert/.nlow`. Update `test/runtests_imas.jl` call sites. -- Kinetic-forces stage keeps consuming `pe_state` + `result` fields analogously - (`set_perturbation_data!(kf_intr, pe_state, result, …)` — mode/metric reads only). - -### 5.3 Standalone Galerkin + additive-gal removal (pulled forward from PR 4) - -Additive gal is what would force a two-solutions-per-run transitional state; it dies in -this PR so the result contract above is final from day one. - -- Factor the wv computation out of `free_run` into a shared helper in - `src/ForceFreeStates/Free.jl`: - `compute_scaled_wv(ctrl, equil, intr) -> (wv, vac)` — the `VacuumInput` + - `compute_vacuum_response` + Chance singfac scaling block (no OdeState involved). - `free_run` calls it; identical numerics by construction. -- `integrator = "galerkin"` becomes legal: the driver's gal branch skips EL integration - and `free_run` entirely; runs `sing_min!` + (when `vac_flag`) `compute_scaled_wv` + - `galerkin_solve` (+ `gal_match_rpec` via the existing flags); `build_result` fills the - gal fields per the §5.1 contract. Errors if `kinetic_factor > 0`. `npert == 1` enforced - by `galerkin_solve` already. -- DELETE: `gal_matched_odestate` (GalerkinMatch.jl) and the driver's additive-gal PE - block (`pe_odet` selection). The additive path (`gal_flag=true` alongside another - integrator) is REMOVED; `gal_flag` joins `_DEPRECATED_FFS_KEYS` + the pre-commit hook. -- RETAIN `_chord_solution_at` (SingularCoupling.jl) as an uncalled helper: re-typed to - `SolutionProfiles`, hard-error branch dropped, stub-style docstring. Kept pending the - `delta_mn` resonant-coupling design (chord-slope derivatives may be useful when PE - consumes rational-surface data instead of profiles) — do NOT re-delete as dead code. -- Gal → PE this cycle: PerturbedEquilibrium's response step requires the free-boundary - δW (`wt0`), which the Galerkin formalism does not produce — so PE warn-skips entirely - on gal results (both gates: `free_boundary` missing kills response, and coupling needs - the response). The gal-native `solution` consumer path in PE therefore stays dormant - until the gal-side δW work lands (next cycle, with the STRIDE matching); the contract - and tests are already in place for it. -- Retoml the four gal decks to `integrator = "galerkin"` (drop `gal_flag`): - `DIIID-like_gal_resistive_example`, `DIIID-like_gal_resistive_pe_example`, - `LAR_ideal_match_test`, `LAR_resistive_match_test`. Their HDF5 outputs become gal-only - (FFS-side integration/energy datasets empty) — accepted per D10; gal datasets identical - because `galerkin_solve` inputs are unchanged. `gal_*` sub-knobs stay (they become - `Galerkin(...)` / `ResistiveMatch` fields in PR 5). - -### 5.4 Tests - -- New `test/runtests_result_struct.jl` (add to `test/runtests.jl` include list): - build a Solovev case; assert Forward result has `solution.basis == :el_axis`, - populated `du_store`/`xi_s_store`, `closure == :ideal`, `iszero(bpen)`, - `delta_prime === nothing`; Riccati result has `delta_prime !== nothing`, - `solution === nothing`, `diagnostics !== nothing`; `require_solution` warns exactly - once (`@test_logs (:warn,)`) and PE skips without throwing on a Riccati result with a - `[PerturbedEquilibrium]` deck; a matched gal deck (LAR_ideal_match_test-class) yields - `solution.basis == :gal_native` and `bpen == galerkin.match.bpen` (zeros under - `gal_ideal_flag`, with `closure == :ideal` there). -- Update every test that consumed `main`'s old named-tuple return - (`runtests_fullruns.jl`, `runtests_imas.jl`, `runtests_rerun_from_h5.jl`, - `runtests_parallel_integration.jl` capture helpers). - -### 5.5 Verification - -Full suite; `runtests_fullruns.jl` (forward decks produce byte-identical HDF5 vs the -stack base, riccati decks emit empty `ForwardIntegration/xi_*` + PE-skip warnings, gal -decks become gal-only files); harness vs the stack base -(`--cases diiid_n1,diiid_n1_riccati,solovev_n1 --refs refactor/local-stability-module,local` -— tracked quantities unchanged; gal-flavored cases re-baselined); docs build. - ---- - -## 6. Interface PR, commit (b) — staged `main` (staging ONLY — gal work is in commit (a)) - -### 6.1 Stage functions (all in `src/GeneralizedPerturbedEquilibrium.jl`; `main_from_inputs` becomes ~40 lines of orchestration) - -```julia -resolve_mode_space!(intr, ctrl) # today's :187-208 n-range block -load_kinetic_context(inputs, intr, ctrl, equil) # :218-236 kf_ctrl + kinetic_profiles -maybe_reform_equilibrium(equil, eq_config, additional_input, intr, ctrl, kin) # :241-268 two-pass -snapshot_forcing_modes(inputs, path, ctrl, preloaded) # :296-313 -prepare_force_free_states!(intr, ctrl, equil) # sing_lim!/sing_find!/filter (:322-360), - # sing_min! (gal), resist_eval_all!, - # m-range (:378-396), make_metric/make_matrix/ - # make_kinetic_matrix (+kinsing finder) -run_force_free_states(ctrl, equil, ffit, intr, metric) -> ForceFreeStatesResult - # integrator dispatch + free_run + - # compute_delta_prime_matrix! + galerkin - # + build_result -run_perturbed_equilibrium(result, inputs, forcing_snapshot, preloaded_coils) -> pe_state -run_kinetic_forces(inputs, result, pe_state, kf_ctrl, kinetic_profiles) -run_slayer_stage(result, inputs, pe_file) # today's closure :512-541, un-closured -``` - -Rules: rerun (`build_inputs_from_h5` → 7-tuple) and IMAS (`dd` kwarg) entry paths -funnel into the same orchestration untouched; `force_termination` early-exits preserved -(both return the new `(; ffs, pe=nothing, slayer)` shape); the two-pass equilibrium -logic stays a pre-FFS stage but is owned by the FFS-facing function -(`maybe_reform_equilibrium` calls `ForceFreeStates.rational_psi_nodes` + -`Equilibrium.refined_psi_grid`/`setup_equilibrium` exactly as today). - -### 6.2 Tests / verification - -Pure code motion: full suite unchanged; harness vs commit (a) must be identical for ALL -cases (no re-baselining in this slice); docs build. Standalone Galerkin and the -additive-gal removal live in commit (a) (§5.3). - ---- - -## 6A. Interface PR, commit (b2) — unified Δ′/matching payload (D14) - -Galerkin computes the same Δ′ physics riccati does (Δ′ matrix, raw D′, `delta_coil`, -PEST-3 blocks), today under separate `galerkin.*` fields and different HDF5 names. This -commit merges the two payloads into the ONE `delta_prime` field so consumers never care -which formalism produced it. - -- **Inventory first (mandatory)**: enumerate every Δ′-flavored field in `GalerkinResult` - and every field in `DeltaPrimeData`, and produce the exact mapping (name, shape, - normalization, sign/side conventions) BEFORE moving anything. Do not assume the two - formalisms' arrays are layout-identical — verify shapes/conventions and document any - genuine mismatch in the type's docstring rather than silently coercing. -- **Type**: extend `DeltaPrimeData` to the union of both payloads (PEST-3 blocks join it). - Fields a formalism doesn't produce stay empty/`nothing`. `build_result` fills it from - whichever formalism ran; the Δ′ payload LEAVES the `galerkin` field, which keeps only - solver internals / FEM diagnostics / RPEC match data (post-inventory list goes in the - struct docstrings). -- **HDF5**: one set of dataset paths for Δ′ outputs regardless of formalism — the - riccati/shared paths are canonical; gal's Δ′ datasets move there (clean break per - `docs/development/hdf5-conventions.md`: update writer, readers, and harness case TOMLs - together; no legacy-path shim). Coordinate with the pending #364 reconciliation so the - paths are renamed once, not twice. -- **SLAYER**: `run_slayer` routes through the unified `delta_prime` — gal-fed SLAYER now - works. Update `runtests_slayer_runner.jl` accordingly. -- **Verification**: gal Δ′ values byte-identical to the pre-unification `galerkin.*` - datasets (only paths/fields move); riccati decks byte-identical throughout; result-struct - testsets extended for the unified field on both formalisms; gal harness cases re-baseline - (h5paths updated). - -## 7. Interface PR, commit (c) — `solve` API - -### 7.1 Dependencies - -- `Project.toml`: add `CommonSolve` to `[deps]` and `[compat]` (`"0.2"`). It is - already in the Manifest transitively — no resolver churn expected. Do NOT remove - anything from Project.toml. - -### 7.2 Integrator structs (`src/ForceFreeStates/Integrators.jl`, new file) - -```julia -abstract type AbstractIntegrator end -Base.@kwdef struct Forward <: AbstractIntegrator end -Base.@kwdef struct Riccati <: AbstractIntegrator - nchunks::Int = 0 # 0 = auto (structure-derived); threads come from julia -t -end -Base.@kwdef struct Galerkin <: AbstractIntegrator - # mirror every gal_* ctrl field with identical defaults, WITHOUT the gal_ prefix: - solver::String = "LU"; nx::Int = 256; nq::Int = 6; pfac::Float64 = 0.001 - dx0::Float64 = 5e-4; dx1::Float64 = 1e-3; dx2::Float64 = 1e-3; cutoff::Int = 10 - tol::Float64 = 1e-10; gnstep::Int = 20000; dx1dx2_flag::Bool = true - sing_order::Int = 6; sing_order_ceiling::Bool = true - rpec_flag::Bool = false; edge_onesided::Bool = false -end - -# D13: inner-layer matching config, integrator-agnostic (NOT part of any integrator struct) -Base.@kwdef struct ResistiveMatch - model = InnerLayer.GGJModel(solver=:ray) # swappable inner layer; backend knobs - # (xfac/nx/nq/cutoff/kmax ← gal_inner_*) live on the model - eta::Vector{Float64} = Float64[] # per-surface, core→edge (← gal_eta) - rho::Vector{Float64} = Float64[] # (← gal_rho) - rotation::Vector{Float64} = Float64[] # Hz; γ_s = 2πi·n·f_s (← gal_rotation) - gamma::Float64 = 5 / 3 # (← gal_gamma) - ideal::Bool = false # (← gal_ideal_flag) -end -``` - -`match !== nothing` replaces `gal_match_flag`. Inside `solve`, matching dispatches per -integrator: Galerkin → `gal_match_rpec`; Riccati → errors "not yet implemented" until -the STRIDE resonant-matching PR lands (that PR also renames/deprecates the `gal_*` -matching TOML keys — until then the TOML keys map onto `ResistiveMatch` internally). - -Mapping helpers `_integrator_symbol(alg)` and `_apply_alg!(ctrl_kwargs, alg)` -translate an alg struct into the `ForceFreeStatesControl` keyword set (pure -translation — `ForceFreeStatesControl` remains the single source of truth for the -solve; the TOML `integrator=` + flat `gal_*`/`nchunks` keys keep working unchanged). - -### 7.3 `solve` + constructors - -- In `ForceFreeStates`: `import CommonSolve: solve` (coexists with the - OrdinaryDiffEq-re-exported `solve`; same generic), then - ```julia - function solve(equil::Equilibrium.PlasmaEquilibrium, alg::AbstractIntegrator; - nn::Union{Int,UnitRange{Int}}, wall::Vacuum.WallShapeSettings=Vacuum.WallShapeSettings(), - match::Union{Nothing,ResistiveMatch}=nothing, - dir_path::String=".", kwargs...) # kwargs = any ForceFreeStatesControl field - -> ForceFreeStatesResult - ``` - Body: build `ctrl` from `alg` + kwargs (`nn_low/nn_high` from `nn`), build `intr`, - then call the PR-4 stages `resolve_mode_space!` → (two-pass reform if the equilibrium - was built with `grid_type="auto"` and not yet refined — reuse - `maybe_reform_equilibrium`) → `prepare_force_free_states!` → `run_force_free_states`. - Top module: `import CommonSolve` and `export solve` (re-export the generic), plus - `export Forward, Riccati, Galerkin, ForceFreeStatesResult`. -- `Equilibrium`: outer constructor - `PlasmaEquilibrium(path::AbstractString; eq_type::String="efit", kwargs...) = - setup_equilibrium(EquilibriumConfig(; eq_type, eq_filename=abspath(path), kwargs...))` - (the `@kwdef` config makes this a 3-liner; `sol/lar/tj` analytic types keep using - `setup_equilibrium(config, analytic_config)` directly — documented, not wrapped). -- `RMPField` (in `ForcingTerms`, exported): a lazy forcing description — - ```julia - struct RMPField - ctrl::ForcingTermsControl # format/file/machine/coil_sets_raw as today - scale::Float64 # uniform multiplier applied to loaded amplitudes/currents - end - RMPField(path::AbstractString; format=_infer_format(path), scale=1.0, kwargs...) - RMPField(coil_sets::Vector{Dict{String,Any}}; scale=1.0, kwargs...) # TOML-shaped coil blocks - ``` - Constraint (verified): ForcingTerms has no n-keyed amplitude concept — amplitude is - per-conductor currents (coil format) or per-mode `ForcingMode.amplitude` (file - formats). `scale` multiplies whichever applies at materialization. A per-n amplitude - dict is deferred (needs ForcingTerms design work; note in docstring). -- `perturbed_equilibrium(ffs::ForceFreeStatesResult, rmp::RMPField; kwargs...)` - (top module): builds `PerturbedEquilibriumControl` from kwargs + - `PerturbedEquilibriumInternal(dir_path=ffs.dir_path)`, materializes forcing modes - from `rmp` against `ffs.equil` (the logic currently inside - `compute_perturbed_equilibrium`'s loading block, `PerturbedEquilibrium.jl:92-124`), - pulls `inner_bpen` from `ffs.galerkin` when `:gal_native`, and calls - `compute_perturbed_equilibrium(ffs, ft_ctrl, pe_ctrl, pe_intr)`. The TOML driver - (`run_perturbed_equilibrium`) is rewired through this same function so there is ONE - forcing-materialization path. - -### 7.4 TOML & docs & tests - -- Finalize `_DEPRECATED_FFS_KEYS` (now includes `use_riccati, use_parallel, - parallel_threads, populate_dense_xi, gal_flag`) + pre-commit hook regex. -- Docs: new "Scripting API" page (`docs/src/api.md` or extend `workflow.md`) with the - four-line UX example; `@autodocs`/`@docs` entries for `solve`, the alg structs, - `RMPField`, `perturbed_equilibrium`, `ForceFreeStatesResult` (checkdocs=:exports - will enforce); nav entry in `docs/make.jl`. -- New `test/runtests_solve_api.jl` (added to runtests.jl list): Solovev end-to-end via - the API only — `PlasmaEquilibrium(...)`; `solve(eq, Forward(); nn=1, ...)` matches a - TOML-driven `main` run on key numbers (`free_boundary.et[1]`, `nzero`); - `solve(eq, Riccati(nchunks=40); nn=1)` produces `delta_prime` matching the TOML run; - `solve(eq, Galerkin(); nn=1)` returns a gal-only result; `perturbed_equilibrium` - round-trip on the forward result; kwarg validation errors (`Riccati` + kinetic). - -### 7.5 Verification - -Full suite; harness (all cases, `--refs develop,local`, report table); docs build; -manual smoke: run the 4-line UX from the Context section in a REPL against -`examples/DIIID-like_ideal_example` inputs. - ---- - -## 7A. Interface PR, commit (d) — ξ unification + tearing surface identity (IMPLEMENTED 2026-08-15) - -Final scope (converged with the user; supersedes the earlier "minimal transpose" reading): - -- **ξ unification (the real one)**: closed axis-to-edge ξ profiles are written from - `result.solution` into the producing formalism's Solutions group with IDENTICAL names and - (mode, solution, psi) axis order: `Solutions/ForwardIntegration/*` (unchanged) and new - `Solutions/GalerkinIntegration/{psi, q, xi_psi, dxi_psidpsi, xi_s}` (the gal grid, issing - nodes dropped — the same arrays as `result.solution`, which IS `Match/xi` repacked). The - gal closure (ideal jump or inner-layer Δ) always yields these profiles; a no-closure gal - run is Δ′-only and writes none. `Match/xi`/`Match/dxidpsi` datasets are REMOVED (they were - the profiles, mislabeled as matching diagnostics); `Match/` keeps cout/cin/Delta_r/bpen/ - rpec_eig/Inner/ only. -- **Raw outer basis demoted to debug output** (user call: solver internals, like dumping an - ODE work array): the old `GalerkinIntegration/Solution/` group is now `Basis/`, written - ONLY under the new `DebugSettings.gal_basis_output` flag ([DEBUG] deck section / `debug=` - API kwarg), transposed to the shared axis order. `ForceFreeStatesResult` now carries - `debug_settings` so the writer sees the flag. verify_gal_{solution,ideal}.jl need the flag. -- **Tearing surface identity (#388 item 2)**: `SLAYERResult` gained `rational_psi`/ - `rational_q` (aligned with `params`; empty when built from bare parameters); - `run_slayer_from_inputs` takes them as kwargs; the loose `run_slayer` fills them from - `surfaces[p.ising]`; writer emits `Tearing/PerSurface/rational_psi|rational_q` when - present + annotations. Gal-fed SLAYER output now identifies its surface subset. -- Benchmarks repointed (verify_gal_match/ideal/solution, compare_gal_vs_el, - scan_{rotation,resistivity}_m2, compare_jbgradpsi_m2 — the filtered psi grid is now - first-class so several scripts simplified); annotation tables updated (axis-order - warning dropped); hdf5-conventions.md updated; result-struct testsets assert - file == result.solution + Basis gating; slayer round-trip asserts surface identity. - -NOT done (stays on #388): item 3 (PE empty-placeholder pattern — align with #368), items -4–5 (schema-owner calls), items 6–8 (comment-audit pass). Full shared-Solutions schema for -closed profiles across formalisms (one group, grid-semantics contract) is future work with -the two-stage PE. - -## 7B. Settled design (2026-08-15): source algebra, two-stage PE, deck-as-serialization - -Discussion CLOSED with the user; decisions D15/D16 below are binding. Commit (c) is -implemented but UNCOMMITTED, so its concrete `RMPField` is REPLACED in place (no shim). - -### D15 — `RMPField` is abstract, with lazy linear algebra (lands in the (c) revision) - -- `RMPField` = the user-facing ABSTRACT supertype of every forcing source. File modes, - coil set + currents, or (future, #377) fields given on ψ=1 / an arbitrary surface via - equivalent surface currents — "they are all just external fields." Constructors on the - abstract type return concrete internal subtypes (today: one leaf wrapping - `ForcingTermsControl`; a surface-field leaf arrives with #377). -- Lazy `+`, `-`, scalar `*`: return a formal linear combination WITHOUT materializing. - Valid because PE is linear in the forcing — materialization commutes with summation. - Both current leaf kinds materialize to the same normalized `Vector{ForcingMode}` basis, - so summation = match (m,n), add amplitudes. Prefer ComplexF64 scale (coil phase - rotation is physical); scale must apply to the MATERIALIZED modes, format-independent. - -### D16 — the deck is the API, serialized (one path) - -Every TOML section corresponds 1:1 to an API object/call; the keys ARE the kwargs -(the `@kwdef` splat is the mapping). Consequences, in delivery order: - -1. **#393 (this PR)**: (c) revision per D15 + commit (d). Nothing else grows scope. - ctrl→TOML serialization explicitly deferred to step 2. -2. **Next PR: "main = 20 lines" (REORDERED ahead of the PE split, user call 2026-08-15: - close FFS completely before touching PE)** — kinetic profiles become an OPTIONAL - ATTRIBUTE OF PlasmaEquilibrium (`kinetic::Union{Nothing,KineticProfiles}`, loaded - data not file path; species/factor knobs are loader kwargs; rationale: the two-pass - grid refinement needs the profiles at equilibrium FORMATION, before any solve exists; - `solve` with kinetic_factor>0 then gates on `eq.kinetic`). COORDINATE with #367 - (struct freeze) — the field addition lands after Jake's PR. SLAYER gets an API entry - point; kinetic + SLAYER get API homes; - `main()` becomes a deck INTERPRETER (parse file → same constructors and calls a - script would make); `main_from_inputs` and the stage functions dissolve. The writer - serializes the RESOLVED ctrl structs (defaults included) into every output — same - blob for TOML and API runs — so every gpec.h5 is replayable and h5→toml regeneration - is just extracting it. Scripting users get the SAME per-section loaders main uses - (e.g. `PlasmaEquilibrium("case_dir/")` reads the `[Equilibrium]` section); no second - config system, ever. Deck completeness is automatic: the deck schema IS the struct - schema, and TOML array-of-tables (`[[ForcingTerms.source]]` with per-block scale) - serializes even the source algebra. -3. **Then: two-stage PE (stacked, AFTER FFS is closed)** — `GeneralPE = - perturbed_equilibrium(ffs)` builds the source-independent response/coupling - operators; `force(GeneralPE, fields)` (or callable `GeneralPE(fields)`) materializes - sources, applies P, computes derived quantities. Pairs with the delta_mn - resonant-coupling work (same territory, same cycle). Payoff: coil scans and - optimization reuse one GeneralPE across many cheap force() calls; a TOML deck maps - onto "GeneralPE + one force()" with no deck-format change. - -Defaults contract (established, keep): both paths splat over the same `@kwdef` struct -defaults — one defaults table. API is deliberately more explicit in two spots (no -default alg; `nn` required, `nn_low/nn_high` kwargs rejected). Deprecated deck keys -warn-and-ignore; unknown API kwargs hard-error (decks are archival, scripts fail fast). - -### Reviewer constraints from Nik (Slack, 2026-08-15 — binding on the follow-on PRs) - -- **No source-type zoo.** The common currency is the control-surface spectrum per source; - keep the concrete RMPField kinds minimal. Endpoint: at most ONE more leaf kind, ever — a - spectrum-literal ("here are control-surface modes, computed elsewhere") — and the #377 - equivalent-surface-currents solve becomes a UTILITY converting fields-on-a-surface into - that spectrum, NOT a type. External couplings (thincurr/surfmn/ferritic tools) cost GPEC - zero adapters: they produce spectra, directly or via the utility. -- **`scale` is a linear-combination weight, never a physical amplitude** (amplitudes are - ambiguous for magnetic materials, coil sets with dropouts, etc.). A degraded coil set is - `nominal - failed_coil`, not `0.9 * nominal`; material fields are computed at the - operating point by the code owning their physics, weight meaningful only for small linear - excursions. Docstrings reworded accordingly (2026-08-15, in the (c) revision). -- Nik explicitly likes the multi-shift/tilt-in-one-run capability (his bookkeeping win) — - keep it central in the two-stage-PE PR spec. - -### Plasma-response methods (Fortran resp_index — binding requirement, 2026-08-15) - -Fortran GPEC computes the plasma inductance / permeability P by FIVE selectable methods -(`plas_indmats(0:4)`, `resp_index`): j=0 = ENERGY method (wt0-based when -resp_induct_flag, else eigenmode energies et) — the Fortran default and the ONLY method -ported to Julia (`compute_plasma_response!`, Response.jl); j=1..4 = SURFACE-CURRENT -methods built from the four `kapmats`/`chpmats` variants (surface current κ and scalar -potential χ per identity-at-edge drive, gpresp_eigen → gpeq_surface at psilim) — these -need only the solutions' EDGE VALUES + vacuum Green's functions, NOT δW. Under gal_flag -Fortran computes only j=1 and forces resp_index=1: gal PE worked via surface currents -from the gal eigenfunctions. - -Consequences (correcting the earlier "PE requires δW" premise): -- The Julia gal→PE skip is a PORTING GAP artifact, not physics: the one ported method is - the one method gal cannot feed. Gal's matched solution already provides the - identity-at-edge columns the surface-current methods consume. -- REQUIREMENT for the two-stage-PE PR: preserve method multiplicity — a ResponseMethod - selection (energy | surface-current variants, the resp_index analog, as a typed - argument not a magic integer), with the surface-current port unlocking gal-fed PE - independently of the gal-δW work. The gal δW work remains scheduled for free-boundary - stability of gal runs and method-0 parity. -- gal_resistive_pe harness expectations change when either route lands. - -### North-star usage sketch (user's, verbatim intent; syntax deliberately sloppy — -### requirements catalog for the two-stage-PE PR, NOT #393 scope) - -```julia -Source_A = RMPField(coil1) -Source_B = RMPField(ferritic_material_fields_at_psi1) # needs #377 -Total_fields = Source_A + Source_B # fast: just records both sources - -GeneralPE = perturbed_equilibrium(ffs_result) -SpecificPE = force(GeneralPE, Total_fields) # Biot-Savart for A, Laplace/current-potential - # solve for B, sum on the control surface, - # apply P, derived quantities per output flags - -# Error-field sensitivity workflow: per-unit sources built by coil manipulation + algebra -PF1U_nominal = RMPField(pf1u_dat, 1) # 1 A -PF1U_shifted = shift_coil(PF1U_nominal, 1e-3) - PF1U_nominal # field per mm of shift - -# Named source SETS: force() runs per key, results in per-key (xarray-like) datasets -rmp_set = ("PF1U_shift"=PF1U_shifted, "PF1U_tilt"=PF1U_tilted, - "ferritic_welds"=surfmn_fields, "REMC"=thincurr_fields) -iter_pe = force(GeneralPE, rmp_set) - -# Keyed, labeled linear algebra on operators and results ("@" = xarray-like matmul): -overlaps_per_amp_per_mm = GeneralPE.C_xe @ iter_pe.Phi_sources_root_area_normalized - -# Collapse per-unit sources to a physical case: keyed scalar sets with wildcards, -# elementwise multiply, then sum to a single total field -tilts_shifts = ("PF1U_shift"=1.1e-3, "PF1U_tilt"=0.9e-3, "ferritic_welds"=1) -currents = ("PF1U_*"=14e3,) -total = sum(tilts_shifts * currents * rmp_set) -real_pe = force(GeneralPE, total; profile_output=true) -jbgradpsi = real_pe.Jbgradpsi -``` - -Requirements this implies for the two-stage-PE PR (catalogue, to be specced there): -named source sets with per-key PE results; coil-geometry manipulation (`shift_coil`, -tilts) composing with source algebra to build per-unit error-field bases; keyed scalar -sets with wildcard matching, elementwise `*` against source sets, `sum` collapsing to -one field; labeled (xarray-style) operator/result access so couplings contract naturally -per key; a `profile_output`-style flag family for derived profile quantities. - - -## 8. Cross-cutting execution rules (for every PR) - -1. **Never merge without third-party human review. State this in every PR body.** -2. Every commit and push requires explicit per-instance maintainer approval. -3. Commit messages: `Area - TAG - message` (e.g. `ForceFreeStates - REFACTOR - Unify Riccati integrator`), - with closed Area and TAG vocabularies per `docs/development/naming.md`. -4. JuliaFormatter-clean (margin 180, kwargs `f(x; a=1)`, no trailing whitespace, LF, - single trailing newline). TOML edits follow `docs/development/toml-conventions.md` - (header block, per-line `# description` copied from the struct docstring, - descriptions identical across files, no Fortran references). -5. No PR/issue numbers in source comments; no step-numbered comments; struct fields - documented in the struct docstring. -6. Docstrings are CommonMark — no bare `[x] (y)` bracket-paren sequences. -7. Run the regression harness before requesting review; paste the report into the PR. -8. Test files are registered in `test/runtests.jl`'s hard-coded include list. -9. Keep this `REFACTOR_PLAN.md` updated (check off completed PRs); delete it in a - final cleanup commit after PR 5 is merged and the Fortran re-comparison is done. - -## 9. Sanity map: capability targets by integrator - -This is the TARGET matrix (D14): outputs representing the same physics are unified across -integrators — one field, one data type, regardless of which formalism produced it. Outputs -fall into three physics classes: - -- **Control surface**: quantities on the plasma boundary (`wp`, `free_boundary` energies). - Every integrator can supply these (gal pending its δW implementation). -- **In-plasma class 1 — full profiles**: ξ/ξ′ (or equivalent) across the volume - (`solution`), used to construct spectral, full-volume perturbed equilibria. - Forward and matched-Galerkin only; Riccati will NEVER produce these. -- **In-plasma class 2 — rational-surface resonant data**: quantities AT the rational - surfaces that quantify island-opening drive: `bpen`, and (future) `delta_mn` — the - matrix encoding the jump in the pitch-resonant derivative of the solution at each - rational surface, from outer-solution asymptotics (for Riccati: recoverable from - `delta_coil`). `delta_mn` yields the perturbed current and the shielded resonant flux, - and is what PE's resonant coupling will consume — no full profiles required. - -Legend: ✅ implemented · 🔜 target pending the named follow-on work · ❌ never · — N/A. - -| Output | Forward | Riccati | Galerkin | -|---|---|---|---| -| `wp` (control surface) | ✅ | ✅ | 🔜 gal δW work | -| `free_boundary` energies (control surface) | ✅ | ✅ | 🔜 gal δW work | -| `solution` — full ξ/ξ′ profiles (class 1) | ✅ `:el_axis` | ❌ (class 2 covers resonant coupling) | ✅ `:gal_native` | -| `closure` / `bpen` (class 2; always present, zeros under `:ideal`) | ✅ `:ideal` | ✅ `:ideal` (🔜 `:matched` with STRIDE matching) | ✅ `:ideal` or `:matched` | -| `delta_mn` (class 2; resonant-derivative jump) | ❌ not planned (no concrete route identified; may not exist) | 🔜 next-week work, from `delta_coil` | 🔜 next-week work | -| `delta_prime` — ONE unified type: Δ′ matrix, raw D′, `delta_coil`, PEST-3 blocks | — | ✅ | ✅ (PEST-3 blocks persisted; riccati recovers them via `pest3_decompose`) | -| raw integrator odet (`diagnostics`: crit, nzero, edge scan, ca) | ✅ | ✅ | — (no radial ODE sweep) | -| kinetic (`kinetic_factor>0`) | ✅ | error | error | -| SLAYER inputs (surfaces + Δ′ matrix) | surfaces only (diag fallback) | ✅ | ✅ via unified `delta_prime` | - -SLAYER is an inner-layer consumer: SLAYER + GGJ should eventually sit behind one abstract -inner-layer interface (same family as the `ResistiveMatch` models, D13). Later pass, not this one. - -## 10. Progress - -### Live status (updated 2026-08-15 — read this first when resuming) - -- **#381 and #387 MERGED into develop** (a0c270f8, 2026-08-15): riccati unification + - LocalStability module are in. Branches deleted; #393 auto-retargeted to develop and - shows MERGEABLE. -- **Interface PR = #393** (`refactor/forcefreestates-result`, worktree `../result-pr3`, - DRAFT, base = develop): - - Commit (a) = 8f8e1645, done: result struct + SolutionProfiles + closure/bpen/wp + - standalone Galerkin + additive-gal removal. Verified: 82/82 result-struct tests, - 357/357 across six files, forward byte-identity (145 datasets), gal-group equivalence - (LAR_ideal_match_test, 12+16 datasets) — all vs f8996d4f, i.e. PRE-#364 base. - - Commit (b) committed: staged-main decomposition - per §6. Verified pure motion — normalized diffs of every stage body vs its old inline - block are character-identical (only function-boundary lines differ); both - force_termination early-exits preserved; one inert reorder (local stability hoisted - ahead of sing_lim!/sing_find!; it reads only equil). Gates: 82/82 result-struct - tests; fresh byte-identity of the coarsened Solovev fixture vs the commit (a) - artifact, 143/143 compared datasets identical (145 total incl. git_version + toml - blob). Review protocol for motion commits: read resulting functions top-down + - behavioral gates, NOT the raw diff; locally use `git diff --color-moved=dimmed-zebra - --color-moved-ws=allow-indentation-change --histogram`. - - Commit (c) implemented, reviewed, and REVISED per D15 (not yet committed): solve API - per §7, then RMPField reworked in place — now an ABSTRACT type with RMPSource leaf - (ComplexF64 scale) and RMPFieldSum lazy linear combinations (+, -, scalar *; flattened - term list); sum materialization evaluates each leaf via a scratch - PerturbedEquilibriumInternal and merges amplitudes per (n,m), sorted; - compute_perturbed_equilibrium accepts Union{ForcingTermsControl,RMPField}; algebra - tests added (type-level testset + one PE call asserting 3A-A == 2A); api.md gained a - Combining-forcing-sources section. THEN materialization made PURE (user request, fewer - !-functions for multithreading): materialize_forcing_modes(ffs, forcing; dir_path, - preloaded_coil_sets, verbose) -> (modes, coil_sets), three dispatch methods, no - mutation; the preload guard + state writes live ONLY in compute_perturbed_equilibrium - (double-apply bugs structurally impossible); driver pre-materialize call deleted. - scale reworded everywhere per Nik: linear-combination WEIGHT, never physical amplitude - (dropout example: nominal - failed_coil, not 0.9*nominal). Final gates: 70/70 solve - API + 17/17 fullruns after the refactor; docs build clean. - THEN problem-type form added (user design call): EulerLagrangeProblem(equil; nn, wall, - match, dir_path, debug, ctrl kwargs) names WHAT is solved (SciML problem/alg split — - PlasmaEquilibrium hosts many future problems, so solve(eq, alg) alone was namespace- - greedy); solve(prob, alg) is canonical, solve(eq, alg; kwargs...) retained as sugar - forwarding to it; nn_low/nn_high rejection lives in the problem constructor. Name - chosen over StabilityProblem because kinetic runs make stability an imprecise label. - Deviations recorded: `solve` lives in the TOP module (prepare_force_free_states! - needs the KineticForces callback; FFS cannot import KineticForces — same CommonSolve - generic, so ForceFreeStates.solve still resolves); ResistiveMatch is a plain config - mapping 1:1 onto gal_* keys (forces gal_rpec_flag=true); solve mirrors TOML side - effects (HDF5 write, local stability); forcing materialization unified in - PerturbedEquilibrium.materialize_forcing_modes! and the TOML driver rewired through - perturbed_equilibrium (ONE forcing path). Verified: 59/59 solve-api + 114/114 - result-struct + 17/17 fullruns (agent + independent rerun), TOML byte-identity - 207/207 datasets after the rewiring, docs build exit 0. - FOUND pre-existing bug (filed as #396, cross-linked from #377): TOML file-forcing - never applies convert_forcing_normalization! (snapshot preloads raw modes; the - isempty guard skips the convert branch) — factor 16.85 on Solovev amplitude-linear - PE outputs; present since the forcing-snapshot PR; NOT fixed here (needs a design - decision re: replay double-conversion; fixing moves TOML outputs). - - Commit (d) IMPLEMENTED by the coordinator directly (not yet committed; §7A has the - full final scope): ξ unification (closed gal profiles in the shared Solutions layout - from result.solution; raw basis debug-gated as Basis/), Tearing/PerSurface - rational_psi/rational_q. Gates GREEN: 133/133 result-struct (file == result.solution, - Basis gating, Match/xi absent), 73/73 slayer (surface identity), 17/17 fullruns, - 6/6 + 14/14 h5-schema (metadata contract on all new/moved datasets), forward fixture - byte-identical 137/137 vs pre-(d) tree. - - Commit (b2) implemented and reviewed (not yet committed; §6A, D14): `DeltaPrimeData` - (now in ForceFreeStatesStructs.jl for include order) carries matrix/raw/coil + gal-only - A/B/Gamma; `galerkin_solve` returns `(GalerkinResult, DeltaPrimeData)`; canonical HDF5 - paths `SingularSurfaces/{Delta_prime_matrix,Delta_prime_raw,Delta_coil,pest3_*}` written - once from `result.delta_prime`; `GalerkinDeltaPrime/` group deleted (per-surface - identifiers moved to `GalerkinIntegration/`); gal-fed SLAYER works. Convention gate - verified (PEST-3 combinations term-identical). Found+fixed pre-existing bug: old gal - `Delta_prime_raw` dataset was (2msing+mpert)×2msing with coil rows duplicated inside. - Verified: 114/114 result-struct, 71/71 slayer (independently rerun), gal Δ′ values - byte-identical under new paths (147/147 common), forward deck untouched (138/138), - benchmarks/ readers repointed. Harness gal_resistive_diiid triage CLOSED: the "3 - changed" rows were the invoking repo's renamed case TOML reading develop's RICCATI - datasets (the additive deck writes both formalisms, and riccati's datasets sit at - exactly the new canonical names) against local's GAL datasets — cross-formalism - apples-to-oranges, not numerical movement. Fresh dual-run proved gal==gal bit-for-bit - (leading raw block isequal, pest3 diag ratio 1.0, coil isequal). Action: re-baseline - the case once; harness cross-ref comparisons spanning the rename boundary are - confounded for this case and should not be repeated. - Also per D14: riccati will NEVER produce full ξ profiles — next-cycle work is the - `delta_mn` rational-surface matrix (from `delta_coil` asymptotics) for PE resonant - coupling, not profile reconstruction. -- **#364 reconciliation DONE** (merge commit b803788e in result-pr3): develop merged - bottom-up (#381 ← develop, #387 ← #381, result-pr3 ← #381-combined). The FFS-writer - conflict resolved as our-structure + #364's literature dataset names; two scope bugs - in auto-merged #364 machinery fixed (`write_root_attrs!` and `apply_main_h5_metadata!` - referenced the deleted `intr` local); `dVdpsi_spline` kwarg threaded through - `run_kinetic_forces`; `diiid_n1_riccati.toml` h5paths renamed (10 paths); stale - `LocalStability/di|dr` docstring in Ballooning.jl fixed (stale on develop too). - Post-merge smoke: 82/82 result-struct + 66/66 slayer. - STILL OWED: fresh byte-identity + gal-equivalence re-runs vs the post-merge base, full - suite, docs build, and one harness re-baseline. -- Standing decisions in force: `_chord_solution_at` retained as uncalled helper (§5.3 — - do not re-delete); gal→PE warn-skips this cycle (no gal δW yet); matching work lands - in a new `Matching/` directory (§ follow-on); directory reorg is a separate post-#367 - post-formatter-PR pure-move PR — never folded into feature commits; comment-audit PRs - follow the #354 pattern, separate from moves. -- Process rules (unchanged): ask before EVERY commit and EVERY push; no formatter ever; - slice-pure commits; third-party human review before ANY merge — non-negotiable. - - -- [ ] 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`. -- [ ] Interface PR (`refactor/forcefreestates-result`) — three commits: (a) §5, (b) §6, (c) §7. - Commit (a) — **implemented (re-sliced §5), reviewed.** - Carries the pivot: no transitional API. `SolutionProfiles` is the one solution slot, - `closure`/`bpen` are unconditional on the result, standalone Galerkin and additive-gal - removal are pulled forward from PR 4, and `pe_solution` / `gal_matched_odestate` are - deleted rather than deferred. Deltas from the §5 spec: - 1. §5 did not say how the ForceFreeStates kernels PE calls keep working once - `ForceFreeStatesInternal` stops crossing the module boundary. Added an abstract - `ModeSpace` supertype (`ForceFreeStatesStructs.jl`) that both - `ForceFreeStatesInternal` and `ForceFreeStatesResult` subtype, and relaxed the - mode-space-only kernels to it: `el_derivatives!`, `materialize_derivative_stores!`, - `build_kinetic_metric_matrices`. - 2. `ForceFreeStatesResult` is parameterized on the equilibrium and `FourFitVars` types - (both are themselves parametric), so `result.equil` / `result.ffit` stay concretely - typed instead of becoming inference barriers on the PE hot paths. - 3. Two call sites outside `src/` consumed `main`'s old named tuple and are updated: - `benchmarks/benchmark_diiid_ideal_ntv_torque.jl` and - `examples/DIIID-like_ideal_example_IMAS/run_imas_example.jl`. - 4. Of the tests §5.4 lists for update, only `runtests_imas.jl` needed it — - `runtests_fullruns.jl`, `runtests_rerun_from_h5.jl` and - `runtests_parallel_integration.jl` never read `main`'s return value (the last drives - the low-level API directly and is unaffected). Coverage was added instead to - `runtests_slayer_runner.jl` (result-facing `run_slayer` dispatch) and - `runtests_imas.jl` (the `free_boundary === nothing` warn-and-skip). - 5. `_chord_solution_at` (PerturbedEquilibrium/SingularCoupling.jl) is deleted: with - `SolutionProfiles.du_store` populated by contract, its `!du_store_populated` branch is - unreachable. The gal-native / ideal-EL / kinetic branches are unchanged. - 6. The `integrator` TOML description changed in all 21 decks that carry the key (the - three-way value list), per the identical-descriptions rule in - `docs/development/toml-conventions.md`. - - Accepted output changes (D10), all spec'd in §5.1/§5.3/§5.4: - - Riccati decks write `ForceFreeStates/Solutions/ForwardIntegration/xi_psi` and `u2` - empty instead of the sparse chunk-endpoint snapshots (`dxi_psi`/`xi_s` were already - empty there). No harness case tracks those datasets. - - The four gal decks become gal-only files: their Galerkin datasets are unchanged, and - the FFS-side integration/energy datasets that the removed additive Riccati run used to - produce are now empty or absent. Verified dataset by dataset (§5.5 gate c). - - Observation for a later PR, not changed here: `result.bpen` has `msing` rows counted - from `intr.sing` under `:ideal` closure but from the Galerkin surface set under - `:matched`. The two can differ when `sing_min!` raises `psilow`. This reproduces the - pre-pivot behavior exactly (the driver previously assigned `gal_data.match.bpen` - directly, and `SingularCoupling` guards with `s <= size(inner_bpen, 1)`), so it is a - pre-existing row-alignment wart, not a regression. - - - [ ] Commit (b) — staged `main` (§6) - - [ ] Commit (c) — `solve` API (§7) -- [ ] Fortran re-comparison of all important quantities -- [ ] Delete this file diff --git a/benchmarks/scan_resistivity_m2.jl b/benchmarks/scan_resistivity_m2.jl deleted file mode 100644 index 00b7c0949..000000000 --- a/benchmarks/scan_resistivity_m2.jl +++ /dev/null @@ -1,57 +0,0 @@ -# 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 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] - -using HDF5, Plots, Printf, TOML - -outpng = length(ARGS) >= 1 ? ARGS[1] : joinpath(@__DIR__, "scan_resistivity_m2.png") -mtarget = length(ARGS) >= 2 ? parse(Int, ARGS[2]) : 2 - -to_c(a) = eltype(a) <: Complex ? ComplexF64.(a) : map(x -> ComplexF64(x.re, x.im), a) - -# read m=target area-normalized b^ψ on the run's PE grid -function read_m2(h5; gal::Bool) - h5open(h5) do f - pa = to_c(read(f["PerturbedEquilibrium/Response/psi_area"])) # [npsi, mpert] - col = mtarget - read(f["Info/mlow"]) + 1 - psi = gal ? read(f["ForceFreeStates/Solutions/GalerkinIntegration/psi"]) : - read(f["ForceFreeStates/Solutions/ForwardIntegration/psi"]) - (psi, pa[:, col]) - end -end - -# collect scan dirs that finished -scandirs = filter(d -> isfile(joinpath(d, "gpec.h5")), - ["/tmp/etascan_0p01", "/tmp/etascan_0p1", "/tmp/etascan_1", "/tmp/etascan_10", "/tmp/etascan_100"]) -etas = [TOML.parsefile(joinpath(d, "gpec.toml"))["ForceFreeStates"]["gal_eta"][1] for d in scandirs] -ord = sortperm(etas) -scandirs, etas = scandirs[ord], etas[ord] -eta_ref = 8e-8 -@printf("%d scan runs: η = %s (×η_ref where η_ref=%.0e)\n", length(etas), - join((@sprintf("%.0e", e) for e in etas), ", "), eta_ref) - -# rational surface for m=target -sing_psi, sing_m = h5open(joinpath(scandirs[1], "gpec.h5")) do f - (read(f["ForceFreeStates/Solutions/GalerkinIntegration/rational_psi"]), read(f["ForceFreeStates/Solutions/GalerkinIntegration/rational_m"])) -end -psi_res = mtarget in sing_m ? sing_psi[findfirst(==(mtarget), sing_m)] : NaN - -cols = cgrad(:viridis, max(length(etas), 2); categorical=true) -plt = plot(; size=(1000, 640), xlabel="ψ_N", ylabel="|b^ψ / ⟨J·|∇ψ|⟩| (area-normalized)", - title="m=$mtarget perturbed normal field — resistivity scan (resistive gal-matched PE)", - legend=:topleft, left_margin=13Plots.mm, bottom_margin=5Plots.mm, right_margin=4Plots.mm) -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 -# 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") - -savefig(plt, outpng) -println("saved: ", abspath(outpng)) diff --git a/benchmarks/scan_rotation_m2.jl b/benchmarks/scan_rotation_m2.jl deleted file mode 100644 index 2664a45c2..000000000 --- a/benchmarks/scan_rotation_m2.jl +++ /dev/null @@ -1,51 +0,0 @@ -# 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 -# 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 - -outpng = length(ARGS) >= 1 ? ARGS[1] : joinpath(@__DIR__, "scan_rotation_m2.png") -mtarget = length(ARGS) >= 2 ? parse(Int, ARGS[2]) : 2 - -to_c(a) = eltype(a) <: Complex ? ComplexF64.(a) : map(x -> ComplexF64(x.re, x.im), a) - -function read_m2(h5; gal::Bool) - h5open(h5) do f - pa = to_c(read(f["PerturbedEquilibrium/Response/psi_area"])) - col = mtarget - read(f["Info/mlow"]) + 1 - psi = gal ? read(f["ForceFreeStates/Solutions/GalerkinIntegration/psi"]) : - read(f["ForceFreeStates/Solutions/ForwardIntegration/psi"]) - (psi, pa[:, col]) - end -end - -scandirs = filter(d -> isfile(joinpath(d, "gpec.h5")), - ["/tmp/rotscan_1", "/tmp/rotscan_2", "/tmp/rotscan_4", "/tmp/rotscan_8", "/tmp/rotscan_16"]) -rots = [TOML.parsefile(joinpath(d, "gpec.toml"))["ForceFreeStates"]["gal_rotation"][1] for d in scandirs] -ord = sortperm(rots) -scandirs, rots = scandirs[ord], rots[ord] -@printf("%d scan runs: rotation f = %s Hz (η fixed = 8e-8)\n", length(rots), join((@sprintf("%g", r) for r in rots), ", ")) - -sing_psi, sing_m = h5open(joinpath(scandirs[1], "gpec.h5")) do f - (read(f["ForceFreeStates/Solutions/GalerkinIntegration/rational_psi"]), read(f["ForceFreeStates/Solutions/GalerkinIntegration/rational_m"])) -end -psi_res = mtarget in sing_m ? sing_psi[findfirst(==(mtarget), sing_m)] : NaN - -cols = cgrad(:plasma, max(length(rots), 2); categorical=true) -plt = plot(; size=(1000, 640), xlabel="ψ_N", ylabel="|b^ψ / ⟨J·|∇ψ|⟩| (area-normalized)", - title="m=$mtarget perturbed normal field — rotation scan (η=8e-8, resistive gal-matched PE)", - legend=:topleft, left_margin=13Plots.mm, bottom_margin=5Plots.mm, right_margin=4Plots.mm) -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/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") - -savefig(plt, outpng) -println("saved: ", abspath(outpng)) diff --git a/benchmarks/verify_gal_ideal.jl b/benchmarks/verify_gal_ideal.jl index b4711f3a4..4eb981a07 100644 --- a/benchmarks/verify_gal_ideal.jl +++ b/benchmarks/verify_gal_ideal.jl @@ -9,7 +9,7 @@ h5path = length(ARGS) >= 1 ? ARGS[1] : "/tmp/gal_ideal_test/gpec.h5" to_c(a) = eltype(a) <: Complex ? ComplexF64.(a) : map(x -> ComplexF64(x.re, x.im), a) cout, deltar, mxi, mdxi, sols, sols_d, iss, sing_psi = h5open(h5path) do f - (to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/cout"])), to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/Delta_r"])), + (to_c(read(f["SingularSurfaces/Match/cout"])), to_c(read(f["SingularSurfaces/Match/Delta_r"])), to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/xi_psi"])), to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/dxi_psidpsi"])), to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Basis/xi_psi"])), to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Basis/dxi_psidpsi"])), Bool.(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Basis/is_rational"])), diff --git a/benchmarks/verify_gal_match.jl b/benchmarks/verify_gal_match.jl index 51dcdd868..507f5bfaa 100644 --- a/benchmarks/verify_gal_match.jl +++ b/benchmarks/verify_gal_match.jl @@ -1,5 +1,5 @@ # Piece 2 verification: RPEC outer↔inner matched solution (closed profiles under -# GalerkinIntegration/xi_psi + matching diagnostics under GalerkinIntegration/Match/*). +# GalerkinIntegration/xi_psi + matching diagnostics under SingularSurfaces/Match/*). # 1. linear-solve residual ‖mat·cof − rmat‖/‖rmat‖ # 2. matched ξ / ξ′ finiteness # 3. edge column == identity basis: each coil drive j must give ξ_edge = e_j (the j-th harmonic), @@ -13,9 +13,9 @@ h5path = length(ARGS) >= 1 ? ARGS[1] : "examples/DIIID-like_gal_resistive_exampl xi, dxi, cout, cin, deltar, eig, resid, sing_psi = h5open(h5path) do f (read(f["ForceFreeStates/Solutions/GalerkinIntegration/xi_psi"]), read(f["ForceFreeStates/Solutions/GalerkinIntegration/dxi_psidpsi"]), - read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/cout"]), read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/cin"]), - read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/Delta_r"]), read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/rpec_eig"]), - read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/residual"]), read(f["ForceFreeStates/Solutions/GalerkinIntegration/rational_psi"])) + read(f["SingularSurfaces/Match/cout"]), read(f["SingularSurfaces/Match/cin"]), + read(f["SingularSurfaces/Match/Delta_r"]), read(f["SingularSurfaces/Match/rpec_eig"]), + read(f["SingularSurfaces/Match/residual"]), read(f["ForceFreeStates/Solutions/GalerkinIntegration/rational_psi"])) end # HDF5 stores ComplexF64 as a compound (re,im); convert if needed to_c(a) = eltype(a) <: Complex ? a : map(x -> ComplexF64(x.re, x.im), a) diff --git a/docs/development/architecture.md b/docs/development/architecture.md index 37aad79b3..059539a8b 100644 --- a/docs/development/architecture.md +++ b/docs/development/architecture.md @@ -63,7 +63,7 @@ Splines are provided by the external `FastInterpolations` package rather than by - `Surfaces/` - Singular-surface finding, Frobenius asymptotics, and GGJ coefficients - `Riccati/` - Chunked fundamental-matrix (STRIDE) driver and Δ' boundary-value problem - `Galerkin/` - RDCON outer-region singular Galerkin Δ' solver - - `Matching/` - Outer↔inner resistive matching (`DeltaPrimeData`, `resonant_match_rpec`) + - `Matching/` - Outer↔inner resistive matching (`DeltaPrimeData`, `MatchProblem`, `MatchResult`) - `Fourfit.jl` - Fourier fitting routines (`MatrixSplines`) - `FixedBoundaryStability.jl` - Fixed boundary analysis - `Free.jl` - Free boundary stability diff --git a/docs/development/hdf5-conventions.md b/docs/development/hdf5-conventions.md index c4631c984..85cc1b709 100644 --- a/docs/development/hdf5-conventions.md +++ b/docs/development/hdf5-conventions.md @@ -40,9 +40,9 @@ Top level (11 groups): | `Info/` | Run metadata: `git_version`, mode-number ranges (`mpert`, `mlow`, …, `mn_index`), `psilim`, `qlim`, `Runtimes/` (per-stage wall-clock seconds) | | `Input/` | Rerun snapshot: `gpec_toml_raw`, `RawInputs/{Equilibrium, ForcingTerms, Coils/}` | | `Equilibrium/` | Scalars (`beta_N`, `q_axis`, `q_95`, `I_p`, …) plus `Profiles/` (1-D on `psi`: 2piF, mu0p, dVdpsi, q) and `Geometry/` (2-D on `psi`×`theta`: rcoords, offset, nu, jac) | -| `ForceFreeStates/` | `Solutions/ForwardIntegration/` (u-solutions), `Solutions/GalerkinIntegration/` (closed ξ profiles in the shared layout, `Match/` diagnostics, the gal surface list, debug-gated `Basis/`), `EulerLagrangeMatrices/{Ideal,Kinetic}`, `FreeBoundaryStability/`, `EdgeScan/` | +| `ForceFreeStates/` | `Solutions/ForwardIntegration/` (u-solutions), `Solutions/GalerkinIntegration/` (closed ξ profiles in the shared layout, the gal surface list, debug-gated `Basis/`), `EulerLagrangeMatrices/{Ideal,Kinetic}`, `FreeBoundaryStability/`, `EdgeScan/` | | `LocalStability/` | Mercier `D_I`, resistive interchange `D_R`, `ballooning_Delta_prime` on `psi`; the ballooning α boundary on `ballooning_psi` | -| `SingularSurfaces/` | Per-rational-surface data: `rational_psi`/`rational_q`/`rational_m`/`rational_n`, GGJ coefficients, `Delta_prime_matrix`/`Delta_prime_raw`/`Delta_coil`/`pest3_A`/`pest3_B`/`pest3_Gamma` (Riccati or Galerkin alike), `Kinetic/` | +| `SingularSurfaces/` | Per-rational-surface data: `rational_psi`/`rational_q`/`rational_m`/`rational_n`, GGJ coefficients, `Delta_prime_matrix`/`Delta_prime_raw`/`Delta_coil`/`pest3_A`/`pest3_B`/`pest3_Gamma` (Riccati or Galerkin alike), the inner-layer `Match/` diagnostics, `Kinetic/` | | `PerturbedEquilibrium/` | `ForcingModes/`, `Response/`, `ResponseMatrices/`, `SingularCoupling/`, `Energies/`, control-surface spectra | | `KineticForces/` | `/` (torque/energy profiles, `EnergyIntegrals/`, `KineticMatrices/`); multi-ion runs add `PerSpecies///` with the same per-method layout, summing to the top-level total | | `ErrorFields/` | `CoilSensitivities/` (per-coil-set control-surface spectra and their rigid shift/tilt derivatives, `DominantMode/` full-window projection); `MonteCarlo/` (intrinsic and corrected `\|δ\|` histograms over the sampled tolerances, per batch and averaged); `Risk/` (threshold density, `P(lock\|δ)`, locking probabilities, `ToleranceScan/`); `NTV/` (correction-coil overlap and NTV torque couplings per kAt) | @@ -53,7 +53,7 @@ Reserved (documented, not yet written): `ForceFreeStates/Solutions/RiccatiIntegr ## Metadata contract (self-describing datasets) -Every dataset outside `Input/` (raw snapshot) and `GalerkinIntegration/Match/` (debug-only) must answer "what is this, in what units, plotted against what" without opening the source — enforced by `test/runtests_h5_schema.jl`: +Every dataset outside `Input/` (raw snapshot) and `SingularSurfaces/Match/` (debug-only) must answer "what is this, in what units, plotted against what" without opening the source — enforced by `test/runtests_h5_schema.jl`: - **`long_name`** — plain-text physics description. - **`units`** — SI string (`"T"`, `"Wb/rad"`, `"A"`, `"m"`, `"J"`, `"N*m"`, `"Hz"`, `"Ohm*m"`); `"1"` for dimensionless (CF convention). Normalized quantities state the normalization in `long_name` (e.g. the power-normalized stability energies are per unit ⟨|ξ|²⟩, not joules). diff --git a/docs/src/api.md b/docs/src/api.md index b2406c53d..dae66d632 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -78,18 +78,56 @@ Which products each formalism can supply differs; a result carries `nothing` in its integrator does not produce and consumers warn and skip rather than erroring. See the [Stability Analysis](stability.md) page for the result struct and its capability gates. -Inner-layer matching is requested with the integrator-agnostic `match` keyword, which closes -the basis with a resistive layer solution instead of the ideal jump condition: +## Inner-layer matching + +Inner-layer matching is its own problem, posed on a finished solve: a `MatchProblem` holds +the outer Δ′ the solve published plus the per-surface layer parameters, and the inner-layer +model passed to `solve` computes the layer response at the prescribed rotation. Solving it +returns a new result with the closure changed from `:ideal` to `:matched` and the +eigenfunctions replaced — the expensive outer solve is reused, so layer-parameter scans +cost one cheap match solve per point: + +```julia +ffs = solve(eq, Galerkin(; rpec_flag=true, cut_solution=true); nn=1) + +matched = solve(MatchProblem(ffs; eta=[1e-6, 2e-6], rho=[1e-7, 1e-7], rotation=[0.0, 0.0]), GGJModel()) +@assert matched.closure === :matched + +# A rotation scan reuses the one outer solve: +bpens = [solve(MatchProblem(ffs; eta=[1e-6, 2e-6], rho=[1e-7, 1e-7], rotation=[f, f]), GGJModel()).bpen + for f in 0.0:50.0:500.0] +``` + +The per-surface η/ρ/rotation can also be derived from the kinetic profiles attached to the +equilibrium (`layer_parameters`), with the explicit vectors as overrides. The problem needs +a Δ′ payload with coil-response columns, so it accepts Galerkin (`rpec_flag=true`) and +Riccati (`vac_flag=true`) results; a Riccati-fed match fills `bpen` and the resonant data but keeps +`solution === nothing` (no outer basis is retained). Only a closure-capable model is +accepted — `GGJModel()` today; `SLAYERModel()` is slab-only and drives the free-eigenvalue tearing +solve instead. + +## Tearing stability + +The free-eigenvalue tearing solve is the second flavor of inner-layer matching: instead of +prescribing the layer rotation, a `TearingProblem` holds the outer Δ′ fixed and root-finds +the growth rate where the inner-layer response matches it. The same model slot applies — +`SLAYERModel()` is the slab layer that exists for exactly this problem. `GGJModel()` also runs +through the scan, but GGJ growth-rate extraction is not implemented yet: its γ are placeholders. ```julia -ffs = solve(eq, Galerkin(); nn=1, - match=ResistiveMatch(; eta=[1e-6, 2e-6], rho=[1e-7, 1e-7], rotation=[0.0, 0.0])) -@assert ffs.closure === :matched +attach_kinetic_profiles!(eq, "kin.h5") # n, T, ω for the layer parameters +ffs = solve(eq, Riccati(); nn=1, vac_flag=true) # Δ′ matrix for the dispersion +tear = solve(TearingProblem(ffs; coupling_mode=:coupled), SLAYERModel()) +tear.gamma_Hz, tear.rational_q # root-found rates per surface ``` -Only the Galerkin formalism implements the match today; requesting one from `Forward` or -`Riccati` errors. Kinetic runs (`kinetic_factor > 0`) need the `[KineticForces]` profiles and -remain TOML-driven. +Keyword arguments of `TearingProblem` are the `[SLAYER]` deck section's procedure knobs +(scan mode and Q-domain, coupling mode, critical-Δ convention, extraction filters); +kinetic profiles come from `profile_file` or, when none is named, from the profiles +attached to the equilibrium. + +Kinetic runs (`kinetic_factor > 0`) need the `[KineticForces]` profiles and remain +TOML-driven. ## Entry points diff --git a/docs/src/galerkin.md b/docs/src/galerkin.md index 04c4aba22..413628f77 100644 --- a/docs/src/galerkin.md +++ b/docs/src/galerkin.md @@ -17,8 +17,7 @@ whichever formalism produced it. These are the outer-region inputs to resistive Select it with `integrator = "galerkin"` in `[ForceFreeStates]`. It replaces the radial ODE integration rather than supplementing it: the run computes its own vacuum response at the control surface (when `vac_flag`) and produces no free-boundary energies or ODE trace. -Setting `gal_match_flag` additionally matches the inner layer, giving a driven ξ solution that -`PerturbedEquilibrium` consumes in place of a forward solution. +Setting `gal_match_flag` additionally matches the inner layer, giving a driven ξ solution. The implementation lives in `src/ForceFreeStates/Galerkin/`: @@ -28,14 +27,13 @@ The implementation lives in `src/ForceFreeStates/Galerkin/`: | `GalerkinGrid.jl` | Packed grid construction and local→global DOF mapping | | `GalerkinAssembly.jl` | Element-level assembly: Hermite basis, Gauss-Lobatto stiffness, resonant and extension cells, boundary conditions | | `GalerkinSolution.jl` | Reconstruct ξ(ψ) and analytic ξ′(ψ) on the gal-native grid | -| `GalerkinMatch.jl` | DRIVEN/RPEC outer↔inner asymptotic matching, whose matched solution PerturbedEquilibrium consumes | | `GalerkinSolve.jl` | Top-level driver `galerkin_solve`, banded solve, Δ′ extraction, PEST-3 blocks, HDF5 output | ## API Reference ```@autodocs Modules = [GeneralizedPerturbedEquilibrium.ForceFreeStates] -Pages = ["Galerkin/GalerkinStructs.jl", "Galerkin/GalerkinGrid.jl", "Galerkin/GalerkinAssembly.jl", "Galerkin/GalerkinSolution.jl", "Galerkin/GalerkinMatch.jl", "Galerkin/GalerkinSolve.jl"] +Pages = ["Galerkin/GalerkinStructs.jl", "Galerkin/GalerkinGrid.jl", "Galerkin/GalerkinAssembly.jl", "Galerkin/GalerkinSolution.jl", "Galerkin/GalerkinSolve.jl"] ``` ## See also diff --git a/docs/src/stability.md b/docs/src/stability.md index c61ca716e..b987cab1a 100644 --- a/docs/src/stability.md +++ b/docs/src/stability.md @@ -144,8 +144,8 @@ zeroing vs GR), not from ODE tolerance; it is present at every thread count. the outer region is discretized on packed Hermite-cubic elements and solved as one global banded system, giving the RDCON resistive ``\Delta'`` matrix and the PEST-3 matching blocks. It computes its own vacuum response and returns no free-boundary energies, no ODE trace, and no fixed-boundary `crit` scan. With -`gal_match_flag` it also matches the inner layer, producing a driven ``\xi`` solution that -`PerturbedEquilibrium` consumes. Kinetic runs are not supported. See +`gal_match_flag` it also matches the inner layer, producing a driven ``\xi`` solution. +Kinetic runs are not supported. See `docs/src/galerkin.md` for the solver and its `gal_*` knobs. Enable with: @@ -292,7 +292,7 @@ The Galerkin Δ′ solver (`src/ForceFreeStates/Galerkin/`) is documented separa ```@autodocs Modules = [GeneralizedPerturbedEquilibrium.ForceFreeStates] -Pages = ["ForceFreeStates.jl", "CoreTypes.jl", "Surfaces/Types.jl", "Riccati/Types.jl", "Matching/DeltaPrime.jl", "Result.jl", "Surfaces/Resist.jl", "Surfaces/ResistEval.jl", "Matching/ResonantMatch.jl", "EulerLagrange.jl", "Surfaces/Finding.jl", "Surfaces/Asymptotics.jl", "Fourfit.jl", "Kinetic.jl", "FixedBoundaryStability.jl", "Utils.jl", "Free.jl", "Riccati/Propagators.jl", "Riccati/Crossings.jl", "Riccati/DeltaPrimeBVP.jl", "Riccati/Driver.jl"] +Pages = ["ForceFreeStates.jl", "CoreTypes.jl", "Surfaces/Types.jl", "Riccati/Types.jl", "Matching/DeltaPrime.jl", "Result.jl", "Surfaces/ResistEval.jl", "Matching/ResonantMatch.jl", "Matching/LayerParameters.jl", "Matching/MatchProblem.jl", "EulerLagrange.jl", "Surfaces/Finding.jl", "Surfaces/Asymptotics.jl", "Fourfit.jl", "Kinetic.jl", "FixedBoundaryStability.jl", "Utils.jl", "Free.jl", "Riccati/Propagators.jl", "Riccati/Crossings.jl", "Riccati/DeltaPrimeBVP.jl", "Riccati/Driver.jl"] ``` ## Example usage diff --git a/examples/DIIID-like_SLAYER_example/gpec.toml b/examples/DIIID-like_SLAYER_example/gpec.toml index a8a99f4c7..b28b9bc35 100644 --- a/examples/DIIID-like_SLAYER_example/gpec.toml +++ b/examples/DIIID-like_SLAYER_example/gpec.toml @@ -67,7 +67,7 @@ dmlim = 0.2 # Truncate integration at (last_rational_q + dmli # invalid, no Re/Im contour crossing polishes to a true zero, and the validity # gate correctly reports them as no_root rather than a spurious growth rate. enabled = true # Run the SLAYER tearing-mode analysis -inner_model = "slayer_fitzpatrick" # Inner-layer Δ(Q) model: "slayer_fitzpatrick", "ggj_shooting", or "ggj_galerkin" +inner_model = "slayer_fitzpatrick" # Inner-layer Δ(Q) model: "slayer_fitzpatrick", "ggj_ray", "ggj_shooting", or "ggj_galerkin" scan_mode = "amr" # Q-plane scan strategy: "amr" (adaptive refinement) or "brute_force" coupling_mode = "uncoupled" # "uncoupled" (per-surface) or "coupled" (multi-surface determinant) dc_type = "none" # Critical-Δ offset selector: "none", "lar", "rfitzp", or "toroidal" diff --git a/regression-harness/cases/gal_resistive_diiid.toml b/regression-harness/cases/gal_resistive_diiid.toml index b13388d8a..42bd21f87 100644 --- a/regression-harness/cases/gal_resistive_diiid.toml +++ b/regression-harness/cases/gal_resistive_diiid.toml @@ -77,9 +77,9 @@ order = 31 # the (grid- and ULP-fragile) near-singular points — an ill-posed benchmark. The producer is validated # separately by benchmarks/verify_gal_solution.jl (finite-difference self-consistency). -# Inner-layer matching data Δ(Q) per surface (resist_eval geometry + GGJ inner solver). (msing × 2) +# Inner-layer matching data Δ(Q) per surface (surface geometry + GGJ inner solver). (msing × 2) [quantities.gal_match_deltar_norm] -h5path = "ForceFreeStates/Solutions/GalerkinIntegration/Match/Delta_r" +h5path = "SingularSurfaces/Match/Delta_r" type = "complex_matrix" extract = "norm" label = "||gal inner-layer Δ||" @@ -88,7 +88,7 @@ order = 41 # Outer-region matched coefficients cout (2·msing × mcoil) — the 4·msing matching assembly + solve. [quantities.gal_match_cout_norm] -h5path = "ForceFreeStates/Solutions/GalerkinIntegration/Match/cout" +h5path = "SingularSurfaces/Match/cout" type = "complex_matrix" extract = "norm" label = "||gal match cout||" @@ -97,7 +97,7 @@ order = 42 # Matching linear-solve residual ‖mat·cof − rmat‖/‖rmat‖ — health check (should stay ~machine eps). [quantities.gal_match_residual] -h5path = "ForceFreeStates/Solutions/GalerkinIntegration/Match/residual" +h5path = "SingularSurfaces/Match/residual" type = "real_scalar" extract = "value" label = "gal match residual" diff --git a/src/Equilibrium/Equilibrium.jl b/src/Equilibrium/Equilibrium.jl index 9022cbd37..53d177c4c 100644 --- a/src/Equilibrium/Equilibrium.jl +++ b/src/Equilibrium/Equilibrium.jl @@ -27,7 +27,7 @@ include("KineticProfiles.jl") # --- Expose types and functions to the user --- export setup_equilibrium, EquilibriumConfig, PlasmaEquilibrium, EquilibriumParameters, ProfileSplines, GeometryProfileSplines, compute_geometry_profiles, - KineticProfileSplines, load_kinetic_profiles, shift_exb_rotation, + KineticProfileSplines, load_kinetic_profiles, shift_exb_rotation, attach_kinetic_profiles!, KineticProfileData, read_kinetic_file, write_kinetic_h5 export flux_surface_metric, flux_surface_area export wants_two_pass, refined_psi_grid, merge_mandatory_nodes, bracket_mandatory_nodes, enforce_min_spacing, implied_knot_count diff --git a/src/Equilibrium/EquilibriumTypes.jl b/src/Equilibrium/EquilibriumTypes.jl index ac5a212fc..52149c9bb 100644 --- a/src/Equilibrium/EquilibriumTypes.jl +++ b/src/Equilibrium/EquilibriumTypes.jl @@ -893,6 +893,8 @@ This object provides a complete representation of the processed plasma equilibri - `ingest::EquilibriumIngest`: raw arrays forwarded from the equilibrium input for the `gpec.h5` rerun snapshot — a [`DirectIngest`](@ref)/[`InverseIngest`](@ref) for file-based equilibria, or `nothing` for analytic ones (regenerated from their TOML section on replay) + - `kinetic::Union{Nothing,KineticProfileSplines}`: kinetic profiles attached with + [`attach_kinetic_profiles!`](@ref), or `nothing`. """ mutable struct PlasmaEquilibrium{P<:ProfileSplines,G<:GeometryProfileSplines,I2D<:FastInterpolations.CubicInterpolantND} config::EquilibriumConfig @@ -920,13 +922,15 @@ mutable struct PlasmaEquilibrium{P<:ProfileSplines,G<:GeometryProfileSplines,I2D psio::Float64 ingest::EquilibriumIngest + kinetic::Union{Nothing,KineticProfileSplines} end -# Solvers build the equilibrium before setup_equilibrium forwards eq_input.ingest, so allow -# construction without it; ingest defaults to nothing and is assigned post-construction. +# Solvers build the equilibrium before setup_equilibrium forwards eq_input.ingest and before +# any kinetic attach, so allow construction without either; both default to nothing and are +# assigned post-construction. PlasmaEquilibrium(config, params, profiles, geometry, rzphi_xs, rzphi_ys, rzphi_rsquared, rzphi_offset, rzphi_nu, rzphi_jac, eqfun_B, eqfun_metric1, eqfun_metric2, ro, zo, psio) = PlasmaEquilibrium(config, params, profiles, geometry, rzphi_xs, rzphi_ys, rzphi_rsquared, rzphi_offset, rzphi_nu, rzphi_jac, - eqfun_B, eqfun_metric1, eqfun_metric2, ro, zo, psio, nothing) + eqfun_B, eqfun_metric1, eqfun_metric2, ro, zo, psio, nothing, nothing) diff --git a/src/Equilibrium/KineticProfiles.jl b/src/Equilibrium/KineticProfiles.jl index 934256554..85da09f30 100644 --- a/src/Equilibrium/KineticProfiles.jl +++ b/src/Equilibrium/KineticProfiles.jl @@ -578,3 +578,13 @@ function shift_exb_rotation(kp::KineticProfileSplines, Δω::Real) return KineticProfileSplines(xs, sample(kp.ni_spline), sample(kp.ne_spline), sample(kp.Ti_spline), sample(kp.Te_spline), sample(kp.omegaE_spline) .+ Float64(Δω), sample(kp.loglam_spline), sample(kp.nui_spline), sample(kp.nue_spline), sample(kp.zeff_spline)) end + +""" + attach_kinetic_profiles!(equil, kinetic_file; kwargs...) -> equil + +Load `kinetic_file` into `equil.kinetic`; keywords go to [`load_kinetic_profiles`](@ref). +""" +function attach_kinetic_profiles!(equil, kinetic_file::AbstractString; kwargs...) + equil.kinetic = load_kinetic_profiles(kinetic_file; chi1=2π * equil.psio, kwargs...) + return equil +end diff --git a/src/ForceFreeStates/CoreTypes.jl b/src/ForceFreeStates/CoreTypes.jl index 5544b0e60..cf0442a60 100644 --- a/src/ForceFreeStates/CoreTypes.jl +++ b/src/ForceFreeStates/CoreTypes.jl @@ -207,19 +207,20 @@ gpec.toml. gal_sing_order_ceiling::Bool = true # auto-raise order by ceil(2·Re(α)) per surface (high Mercier index) gal_rpec_flag::Bool = false # append mpert coil-response columns to the Δ′ solve (RDCON rpec_flag): unit boundary sources whose plasma response is recorded; needed for the driven (resistive perturbed-equilibrium) Δ_gw gal_edge_onesided::Bool = false # pack the two end intervals one-sided toward their single rational end (vs the Fortran symmetric "both" pack); avoids the fine edge cell that inflates cond(A). Default false = faithful to gal.f. + gal_cut_solution::Bool = false # also keep the cut solution for a later match's composite inner profiles; implied by gal_match_flag # --- DRIVEN (RPEC) outer↔inner asymptotic matching (rmatch match_rpec port) --- gal_match_flag::Bool = false # enable the RPEC inner-layer matching: solve the coil-driven matched ξ(ψ) from the gal Δ′ + the inner-layer Δ(Q). Requires gal_rpec_flag=true. gal_ideal_flag::Bool = false # within the match, build the IDEAL solution: skip the inner-layer Δ, use bare coil columns (cout=0). Mirrors Fortran rmatch coil%ideal_flag (the EL reference). eta/rho/rotation ignored. gal_inner_solver::String = "ray" # inner-layer Δ backend for the match: "ray" (rotated-contour collocation, certified Δ at the optimal θ = arg(Q)/4; robust for |Q| ≳ 1) or "galerkin" (Hermite-cubic inps; drifts for |Q| ≳ 1) - # --- Inner-layer "galerkin" backend knobs (used only when gal_inner_solver = "galerkin") --- - gal_inner_xfac::Float64 = 10.0 # asymptotic-matching radius multiplier (inps_xfac: xmax × 10) - gal_inner_nx::Int = 1280 # inner-layer grid cells (128 · xfac in the reference) - gal_inner_nq::Int = 5 # quadrature order per cell - gal_inner_cutoff::Int = 5 # cells carrying the large solution as driving term - gal_inner_kmax::Int = 8 # large-x asymptotic series order (↔ order_pow) + # --- Inner-layer "galerkin" backend overrides (gal_inner_solver = "galerkin"; unset keeps the backend default) --- + gal_inner_xfac::Union{Nothing,Float64} = nothing # asymptotic-matching radius multiplier + gal_inner_nx::Union{Nothing,Int} = nothing # inner-layer grid cells + gal_inner_nq::Union{Nothing,Int} = nothing # quadrature order per cell + gal_inner_cutoff::Union{Nothing,Int} = nothing # cells carrying the large solution as driving term + gal_inner_kmax::Union{Nothing,Int} = nothing # large-x asymptotic series order gal_eta::Vector{Float64} = Float64[] # per-surface resistivity η (length msing, core→edge); Fortran rmatch `eta` gal_rho::Vector{Float64} = Float64[] # per-surface mass density ρ [kg/m³] (length msing, core→edge); Fortran rmatch `massden` gal_rotation::Vector{Float64} = Float64[] # per-surface rotation frequency f [Hz] (length msing, core→edge); forced eigenvalue γ_s = 2πi·n·f. Fortran rmatch `rotation` - gal_gamma::Float64 = 5 / 3 # ratio of specific heats Γ for the resistive-layer coefficients (resist_eval G term) + gal_gamma::Float64 = 5 / 3 # ratio of specific heats Γ for the resistive-layer coefficients (ggj_parameters G term) frobenius_psi_max::Float64 = 0.0 end diff --git a/src/ForceFreeStates/ForceFreeStates.jl b/src/ForceFreeStates/ForceFreeStates.jl index 1dc93e776..dbaf5c62a 100644 --- a/src/ForceFreeStates/ForceFreeStates.jl +++ b/src/ForceFreeStates/ForceFreeStates.jl @@ -13,6 +13,7 @@ using AdaptiveArrayPools using Roots using FastGaussQuadrature: gausslobatto using QuadGK: quadgk, quadgk! +import CommonSolve import ..Equilibrium import ..Utilities @@ -34,10 +35,10 @@ include("EulerLagrange.jl") # Singular-surface machinery: finding/filtering, Frobenius asymptotics, GGJ coefficients include("Surfaces/Finding.jl") include("Surfaces/Asymptotics.jl") -include("Surfaces/Resist.jl") include("Surfaces/ResistEval.jl") # Outer<->inner resistive matching +include("Matching/LayerParameters.jl") include("Matching/ResonantMatch.jl") include("FixedKineticMatrices.jl") @@ -57,7 +58,6 @@ include("Galerkin/GalerkinStructs.jl") include("Galerkin/GalerkinGrid.jl") include("Galerkin/GalerkinAssembly.jl") include("Galerkin/GalerkinSolution.jl") -include("Galerkin/GalerkinMatch.jl") include("Galerkin/GalerkinSolve.jl") # Scripting-API integrator selectors: pure configuration translated onto ForceFreeStatesControl. @@ -66,6 +66,10 @@ include("Integrators.jl") # The published solve product; last, so it can name every type the stages above define. include("Result.jl") +# Post-solve inner-layer matching: consumes and produces ForceFreeStatesResult, so it loads +# after the result machinery. +include("Matching/MatchProblem.jl") + # These are used for various small tolerances and root finders throughout ForceFreeStates global eps = 1e-10 global itmax = 50 diff --git a/src/ForceFreeStates/Galerkin/GalerkinMatch.jl b/src/ForceFreeStates/Galerkin/GalerkinMatch.jl deleted file mode 100644 index 7a6f2d978..000000000 --- a/src/ForceFreeStates/Galerkin/GalerkinMatch.jl +++ /dev/null @@ -1,234 +0,0 @@ -# GalerkinMatch.jl -# -# DRIVEN (RPEC) outer↔inner asymptotic matching. Port of rmatch `match_rpec` (match.f) and the -# outer total-solution construction of `match_output_solution` (match.f). -# -# Per rational surface: the forced eigenvalue γ_s = 2πi·n·f_s (gal_rotation is the rotation f in Hz, so γ -# is the layer Doppler frequency), the inner-layer matching data Δ(Q) via the Galerkin GGJ solver -# (InnerLayer.solve_inner), then the 4·msing matching -# system mat·[cout;cin] = rmat is assembled (rmat = −transpose of the gal Δ′ coil block) and solved for -# the per-coil outer/inner coefficients. The matched outer solution for coil drive j is -# ξ_j = Σ_{isol=1}^{2 msing} cout[isol,j]·sols[:,:,isol] + sols[:,:,2 msing+j] -# and identically for ξ′ with sols_deriv. Stacking over j gives the matched fundamental matrix in the -# identity-at-edge basis (coil column j has ξ_edge = e_j). - -""" - gal_match_rpec(ctrl, equil, intr, gal_result, dp) -> GalMatchResult - -Solve the coil-driven RPEC matching from the outer Δ′ payload `dp` and the per-surface inner-layer Δ. -Requires `gal_result.solution` (reconstructed outer ξ/ξ′) and the rpec coil block `dp.coil`. - -The resistive path uses per-surface inputs `ctrl.gal_eta`/`gal_rho`/`gal_rotation` (length `msing`) + -`ctrl.gal_gamma`. When `ctrl.gal_ideal_flag`, the inner layer is skipped and `cout=0`, so the matched -solution is the bare ideal coil column (Fortran rmatch `coil%ideal_flag`) — the DCON/EL reference. -""" -function gal_match_rpec(ctrl::ForceFreeStatesControl, equil, intr::ForceFreeStatesInternal, - gal_result::GalerkinResult, dp::DeltaPrimeData) - - msing = gal_result.msing - mpert = intr.numpert_total - nn = intr.nlow - mcoil = mpert # rpec: one coil column per poloidal harmonic - - gal_result.solution !== nothing || - error("gal_match_rpec: gal_result.solution is missing (need the gal-reconstructed ξ/ξ′)") - isempty(dp.coil) && - error("gal_match_rpec: the coil-response block is empty — gal_rpec_flag must be true for RPEC matching") - - if ctrl.gal_ideal_flag - # Ideal limit (Fortran rmatch coil%ideal_flag, match.f): skip the inner layer entirely - # and set the resistive plasma combination to zero. The gal outer solve is already fully ideal, so - # the shared construction below collapses to the bare ideal coil column sols(:,:,csol). - cout = zeros(ComplexF64, 2msing, mcoil) - cin = zeros(ComplexF64, 2msing, mcoil) - deltar = zeros(ComplexF64, msing, 2) - bpen = zeros(ComplexF64, msing, mcoil) - inner_psi = Vector{Float64}[] # no inner layer in the ideal limit - inner_xi = Matrix{ComplexF64}[] - inner_b = Matrix{ComplexF64}[] - inner_params = InnerLayer.GGJParameters[] # inner layer skipped in the ideal limit - rpec_eig = zeros(ComplexF64, msing) - residual = 0.0 - else - for (name, v) in (("gal_eta", ctrl.gal_eta), ("gal_rho", ctrl.gal_rho), ("gal_rotation", ctrl.gal_rotation)) - length(v) == msing || error("gal_match_rpec: $name has length $(length(v)), expected msing=$msing (one value per surface, core→edge)") - end - - # Re-derive the gal singular-surface set (same helper as galerkin_solve) for the SingType objects - # resist_eval needs. - sings, _, _ = gal_resonant_surfaces(intr, equil) - length(sings) == msing || - error("gal_match_rpec: re-derived $(length(sings)) surfaces, expected msing=$msing") - ctrl.gal_inner_solver in ("ray", "galerkin") || - error("gal_match_rpec: gal_inner_solver = \"$(ctrl.gal_inner_solver)\" (expected \"ray\" or \"galerkin\")") - - # --- inner-layer matching data Δ(Q) per surface (deltac_run; match.f) --- - # solve_inner_profile returns the same Δ as solve_inner plus the reconstructed inner-layer field, - # so the layer-center value (penetrated field) comes for free from the single matching solve. - deltar = zeros(ComplexF64, msing, 2) - rpec_eig = zeros(ComplexF64, msing) - # Per-surface layer-center field weights pen[i,k] = scale·Ψ_k(0), parity k=1,2 (match.f intotsol_b). - chi1 = 2π * equil.psio - pen = zeros(ComplexF64, msing, 2) - # Per-surface inner-layer ξ_ψ building blocks for the matched solution (match.f intotsol, deltac comp 2): - # the ψ grid ψ_s ± X·x0/v1 (left reversed then right) and the resc-scaled odd/even parity profiles. - inner_psi = Vector{Vector{Float64}}(undef, msing) - inner_odd = Vector{Vector{ComplexF64}}(undef, msing) # Ξ₁, antisymmetric across ψ_s - inner_even = Vector{Vector{ComplexF64}}(undef, msing) # Ξ₂, symmetric across ψ_s - inner_bodd = Vector{Vector{ComplexF64}}(undef, msing) # b^ψ₁ = scale·resc·Ψ₁, symmetric (Ψ′₁(0)=0) - inner_beven = Vector{Vector{ComplexF64}}(undef, msing) # b^ψ₂ = scale·resc·Ψ₂, antisymmetric (Ψ₂(0)=0) - inner_params = Vector{InnerLayer.GGJParameters}(undef, msing) - for i in 1:msing - params = resist_eval(sings[i], equil, intr; eta=ctrl.gal_eta[i], rho=ctrl.gal_rho[i], - gamma=ctrl.gal_gamma, ising=i) - inner_params[i] = params - γ = 2π * im * nn * ctrl.gal_rotation[i] # forced eigenvalue; gal_rotation is f [Hz], γ = 2πi·n·f - rpec_eig[i] = γ - # gal_inner_solver = "ray" (default) uses the rotated-ray backend (certified optimal-θ Δ - # + θ=0 real-axis profile re-solve with runtime certificate — see the InnerLayer - # solve_inner_profile(::GGJModel{:ray}, ...) docstring); "galerkin" uses the Hermite-FEM - # real-axis solve with the ctrl.gal_inner_* knobs (defaults match the Fortran rmatch - # deltac/inps reference). - inner = if ctrl.gal_inner_solver == "ray" - InnerLayer.solve_inner_profile(InnerLayer.GGJModel(; solver=:ray), params, γ) - else - InnerLayer.solve_inner_profile(InnerLayer.GGJModel(; solver=:galerkin), params, γ; - xfac=ctrl.gal_inner_xfac, nx=ctrl.gal_inner_nx, nq=ctrl.gal_inner_nq, cutoff=ctrl.gal_inner_cutoff, kmax=ctrl.gal_inner_kmax) - end - deltar[i, 1] = inner.Δ[1] - deltar[i, 2] = inner.Δ[2] - profΨ, profΞ, xg = inner.Ψ, inner.Ξ, inner.x - # Amplitude rescale: inner profiles are normalized in X = v1·δψ/x0 (inner_psi below); - # inner.rescale converts a big-branch amplitude to the outer δψ-normalization. - resc = inner.rescale - # b-field scaling, derived from the code's own outer convention (SingularCoupling): - # b_m = 2πi·χ₁·(m−nq)·ξ_m, m−nq = −n·q′·δψ, δψ = dψdx·X, resc·Ξ(X) ↔ ξ_m, - # and the far-field identity Ψ = X·Ξ (GWP2016 Eq. 16; Ψ is the normal-field variable, Eq. A17): - # b_m = −2πi·χ₁·n·q′·dψdx·resc·Ψ. - scale = -2π * chi1 * im * nn * sings[i].q1 * inner.dψdx - pen[i, 1] = scale * profΨ[1, 1] * resc # layer center X=0, parity 1 (Ψ(0)≠0) - pen[i, 2] = scale * profΨ[1, 2] * resc # parity 2 (Ψ(0)=0 ⇒ ~0) - xvar = xg .* inner.dψdx # inner X → ψ-distance (deltac.f:1822) - inner_psi[i] = vcat(reverse(sings[i].psifac .- xvar), sings[i].psifac .+ xvar) - inner_odd[i] = resc .* vcat(reverse(.-profΞ[:, 1]), profΞ[:, 1]) # comp 2, parity 1 (odd: −left,+right) - inner_even[i] = resc .* vcat(reverse(profΞ[:, 2]), profΞ[:, 2]) # comp 2, parity 2 (even) - # b^ψ profiles on the same two-sided grid: Ψ is the normal-field variable (GWP2016 A17), - # b_m(δψ) = scale·resc·Ψ(X) throughout the layer (→ outer frozen-in relation via Ψ = XΞ). - inner_bodd[i] = (scale * resc) .* vcat(reverse(profΨ[:, 1]), profΨ[:, 1]) # parity 1: Ψ even - inner_beven[i] = (scale * resc) .* vcat(reverse(.-profΨ[:, 2]), profΨ[:, 2]) # parity 2: Ψ odd - end - - # --- assemble the 4·msing matching system (match.f) --- - mat = zeros(ComplexF64, 4msing, 4msing) - rmat = zeros(ComplexF64, 4msing, mcoil) - @views mat[(2msing+1):4msing, 1:2msing] .= transpose(dp.raw) # Δ_out - @views rmat[(2msing+1):4msing, :] .= .-dp.coil # −Δ_coil source (already surface-side × edge mode) - for ising in 1:msing - idx1 = 2ising - 1 - idx2 = 2ising - idx3 = idx1 + 2msing - idx4 = idx2 + 2msing - delta1 = deltar[ising, 1] - delta2 = deltar[ising, 2] - mat[idx1, idx1] = 1 - mat[idx2, idx2] = 1 - mat[idx1, idx3] = -1 - mat[idx1, idx4] = 1 - mat[idx2, idx3] = -1 - mat[idx2, idx4] = -1 - # inner-layer Δ block signs per match.f match_rpec - mat[idx3, idx3] = -delta1 - mat[idx3, idx4] = delta2 - mat[idx4, idx3] = -delta1 - mat[idx4, idx4] = -delta2 - end - - # --- solve mat·cof = rmat for the outer/inner coefficients (match.f) --- - cof = mat \ rmat - residual = norm(mat * cof - rmat) / max(norm(rmat), 1e-300) - cout = cof[1:2msing, :] - cin = cof[(2msing+1):4msing, :] - - # Inner-layer penetrated (reconnected) resonant field at each rational surface, read off the GGJ - # inner solution at the layer center exactly as Fortran match_output_solution builds intotsol_b - # (match.f) — cusp-free, fit-free. bpen[i,j] = pen₁(i)·cin[2i,j] + pen₂(i)·cin[2i-1,j]. - bpen = zeros(ComplexF64, msing, mcoil) - for i in 1:msing, j in 1:mcoil - bpen[i, j] = pen[i, 1] * cin[2i, j] + pen[i, 2] * cin[2i-1, j] - end - - # Matched inner-layer ξ_ψ(ψ) per surface, per coil drive (match.f intotsol): odd parity weighted by - # cin[2i] (cofin(2·ising)), even parity by cin[2i-1] (cofin(2·ising-1)). - inner_xi = [inner_odd[i] * transpose(cin[2i, :]) .+ inner_even[i] * transpose(cin[2i-1, :]) for i in 1:msing] - # Matched inner-layer b^ψ(ψ) per surface, per coil drive. - inner_b = [inner_bodd[i] * transpose(cin[2i, :]) .+ inner_beven[i] * transpose(cin[2i-1, :]) for i in 1:msing] - end - - # --- matched outer ξ/ξ′ per coil drive (match.f); ideal: cout=0 ⇒ bare coil column --- - sols = gal_result.solution.xi # (mpert, ngrid, nsol) - sols_deriv = gal_result.solution.xi_deriv - ngrid = size(sols, 2) - xi = zeros(ComplexF64, mpert, ngrid, mcoil) - xi_deriv = zeros(ComplexF64, mpert, ngrid, mcoil) - for j in 1:mcoil - csol = 2msing + j # this coil's particular-solution column - @views xi[:, :, j] .= sols[:, :, csol] - @views xi_deriv[:, :, j] .= sols_deriv[:, :, csol] - for isol in 1:2msing - @views xi[:, :, j] .+= cout[isol, j] .* sols[:, :, isol] - @views xi_deriv[:, :, j] .+= cout[isol, j] .* sols_deriv[:, :, isol] - end - end - - # --- composite inner-region solution: cut outer background + layer (match.f intotsol/intotsol_b) --- - # The layer solution alone carries only the resonant content it resolves; the smooth background - # removed by the cut has to be added back for the inner and outer solutions to overlap in the - # matching region. Without this the inner profile does not graft onto the outer eigenfunction. - if !isempty(inner_psi) - sols_cut = gal_result.solution.xi_cut - isempty(sols_cut) && error("gal_match_rpec: solution.xi_cut is empty — the cut solution is " * - "required to build the composite inner-region solution") - cut_range = gal_result.solution.cut_range - keep = .!gal_result.solution.issing - psi_keep = gal_result.solution.psi[keep] - chi1_c = 2π * equil.psio - for i in 1:msing - m_res = round(Int, nn * sings[i].q) - ires = m_res - intr.mlow + 1 - # Clip the layer to the window where the cut is active; outside it the cut removes - # nothing and the composite is undefined (Fortran writes no points there). - lo, hi = cut_range[i, 1], cut_range[i, 2] - inside = findall(p -> lo <= p <= hi, inner_psi[i]) - if isempty(inside) - @warn "gal_match_rpec: inner layer at ψ=$(round(sings[i].psifac, digits=5)) lies outside the " * - "resonant/extension cells (ψ ∈ [$lo, $hi]); raise gal_dx1/gal_dx2 to overlap the layer" surface = i - continue - end - length(inside) < length(inner_psi[i]) && @info "gal_match_rpec: surface $i inner region clipped to the " * - "cut window ($(length(inside)) of $(length(inner_psi[i])) points)" - inner_psi[i] = inner_psi[i][inside] - inner_xi[i] = inner_xi[i][inside, :] - inner_b[i] = inner_b[i][inside, :] - - # cut outer background for this surface's resonant harmonic, per coil drive - cutmn = Matrix{ComplexF64}(undef, length(psi_keep), mcoil) - for j in 1:mcoil - @views cutmn[:, j] .= sols_cut[ires, keep, 2msing+j] - for isol in 1:2msing - @views cutmn[:, j] .+= cout[isol, j] .* sols_cut[ires, keep, isol] - end - end - itp = cubic_interp(psi_keep, Series(cutmn); extrap=ExtendExtrap()) - buf = Vector{ComplexF64}(undef, mcoil) - hint = Ref(1) - for (ip, psi_p) in enumerate(inner_psi[i]) - itp(buf, psi_p; hint=hint) - singfac = m_res - nn * equil.profiles.q_spline(psi_p) - @views inner_xi[i][ip, :] .+= buf - @views inner_b[i][ip, :] .+= (2π * im * chi1_c * singfac) .* buf - end - end - end - - return GalMatchResult(cout, cin, xi, xi_deriv, deltar, bpen, inner_psi, inner_xi, inner_b, inner_params, rpec_eig, residual) -end diff --git a/src/ForceFreeStates/Galerkin/GalerkinSolution.jl b/src/ForceFreeStates/Galerkin/GalerkinSolution.jl index 6ad3a5a16..4d16b1eca 100644 --- a/src/ForceFreeStates/Galerkin/GalerkinSolution.jl +++ b/src/ForceFreeStates/Galerkin/GalerkinSolution.jl @@ -12,7 +12,7 @@ # performed. The `cut=true` path swaps the resonant Frobenius series for its leading-order term # (`sing_get_ua_gal_cut`); the resulting cut solution supplies the regular background that the # resistive inner-layer solution is added to when forming the composite solution at a rational -# surface (see `gal_match_rpec`). +# surface (see MatchProblem). # Sampling points per cell, matching Fortran interp_np_res / interp_np (gal.f): coarse in regular cells, # dense in resonant/extension cells to resolve the near-singular asymptotic series. @@ -152,8 +152,7 @@ it for every solution column. Port of `gal_output_solution` (gal.f); the binary/ and the `b_flag` ξ→b^ψ conversion (off by default) are not ported — we keep ξ itself. When `delta` (the `(nsol, 2·msing)` small-solution coefficient matrix) is supplied, the cut solution -`xi_cut` is evaluated on the same grid; it is the background of the composite inner-region solution -built by `gal_match_rpec`. +`xi_cut` is evaluated on the same grid, as the background of MatchProblem's composite inner solution. """ function gal_output_solution(ws::GalWorkspace, asymps::Vector{GalSingAsymp}, sings::Vector{SingType}, intr::ForceFreeStatesInternal, profiles, psihigh::Float64; diff --git a/src/ForceFreeStates/Galerkin/GalerkinSolve.jl b/src/ForceFreeStates/Galerkin/GalerkinSolve.jl index a06e37e82..881361095 100644 --- a/src/ForceFreeStates/Galerkin/GalerkinSolve.jl +++ b/src/ForceFreeStates/Galerkin/GalerkinSolve.jl @@ -3,7 +3,7 @@ # Top-level driver for the RDCON outer-region singular Galerkin Δ′ solve, plus the per-cell assembly # orchestration, the banded solve, Δ′ extraction, PEST-3 blocks, and HDF5 output. # Ports gal_make_arrays (gal.f), gal_solve (gal.f), and gal_write_pest3_data -# (gal.f). The DRIVEN/RPEC inner-layer matching is wired in via gal_match_rpec (GalerkinMatch.jl). +# (gal.f). Inner-layer matching is applied afterwards by MatchProblem. """ gal_make_arrays!(ws, ctrl, equil, mats, intr, asymps, sings, nn, wv_edge) @@ -41,7 +41,7 @@ end """Empty `GalerkinResult` for a domain with no resonant surfaces.""" function empty_galerkin_result() - return GalerkinResult(0, Float64[], Float64[], Int[], Int[], Float64[], ComplexF64[], nothing, nothing) + return GalerkinResult(0, Float64[], Float64[], Int[], Int[], Float64[], ComplexF64[], nothing) end """ @@ -179,28 +179,17 @@ function galerkin_solve(ctrl::ForceFreeStatesControl, equil, mats::MatrixSplines dp_coil = ncoil > 0 ? permutedims(delta[(2*msing+1):(2*msing+ncoil), :]) : Matrix{ComplexF64}(undef, 0, 0) dp = DeltaPrimeData(Deltap, dp_raw, dp_coil, Ap, Bp, Gammap) - # Reconstruct ξ(ψ) AND analytic ξ′(ψ) on the gal-native grid (gal_output_solution). + # Reconstruct ξ(ψ) AND analytic ξ′(ψ) on the gal-native grid (gal_output_solution); the cut + # solution only when a later match needs the composite inner profiles. ctrl.verbose && @info "Reconstructing outer-region ξ and analytic ξ′ on the gal grid" solution = gal_output_solution(ws, asymps, sings, intr, equil.profiles, psihigh; - delta=(ctrl.gal_match_flag ? delta : nothing)) + delta=((ctrl.gal_cut_solution || ctrl.gal_match_flag) ? delta : nothing)) sing_psi = [s.psifac for s in sings] sing_q = [s.q for s in sings] sing_m = [s.m[1] for s in sings] sing_n = [s.n[1] for s in sings] - result = GalerkinResult(msing, sing_psi, sing_q, sing_m, sing_n, di, alpha, solution, nothing) - - # DRIVEN (RPEC) outer↔inner matching: build the coil-driven matched ξ/ξ′ (gal_match_rpec). - ctrl.gal_match_flag || return result, dp - ctrl.gal_rpec_flag || error("galerkin_solve: gal_match_flag=true requires gal_rpec_flag=true") - ctrl.verbose && @info( - ctrl.gal_ideal_flag ? - "RPEC matching: IDEAL solution (inner layer skipped, bare coil columns)" : - "RPEC matching: inner-layer Δ(Q) + outer↔inner solve for the coil-driven ξ" - ) - match = gal_match_rpec(ctrl, equil, intr, result, dp) - ctrl.gal_ideal_flag || (ctrl.verbose && @info "RPEC matching: linear-solve residual = $(match.residual)") - return GalerkinResult(msing, sing_psi, sing_q, sing_m, sing_n, di, alpha, solution, match), dp + return GalerkinResult(msing, sing_psi, sing_q, sing_m, sing_n, di, alpha, solution), dp end """ @@ -253,7 +242,7 @@ end write_galerkin!(out_h5, result::GalerkinResult; basis_output=false) Write the Galerkin solver outputs into the open HDF5 file, under -`ForceFreeStates/Solutions/GalerkinIntegration/`: the matching diagnostics and the surface list +`ForceFreeStates/Solutions/GalerkinIntegration/`: the surface list the solve ran over (a subset of `SingularSurfaces/` when the domain or the m-band excludes rationals). The Δ′/PEST-3 matrices and the closed ξ profiles are NOT written here — both go to formalism-independent homes from the driver writer (`SingularSurfaces/` off `result.delta_prime`; @@ -284,36 +273,11 @@ function write_galerkin!(out_h5, result::GalerkinResult; basis_output::Bool=fals isempty(sol.xi_cut) || (out_h5["$gal/Basis/xi_psi_cut"] = permutedims(sol.xi_cut, (1, 3, 2))) isempty(sol.cut_range) || (out_h5["$gal/Basis/cut_range"] = sol.cut_range) end - if result.match !== nothing - m = result.match - out_h5["$gal/Match/cout"] = m.cout - out_h5["$gal/Match/cin"] = m.cin - out_h5["$gal/Match/Delta_r"] = m.deltar - out_h5["$gal/Match/bpen"] = m.bpen - out_h5["$gal/Match/rpec_eig"] = m.rpec_eig - # Per-surface inner-layer ξ_ψ(ψ) (match.f intotsol); ragged grids → one dataset pair per surface. - for i in eachindex(m.inner_psi) - out_h5["$gal/Match/Inner/psi_$i"] = m.inner_psi[i] - out_h5["$gal/Match/Inner/xi_$i"] = m.inner_xi[i] - out_h5["$gal/Match/Inner/b_$i"] = m.inner_b[i] - end - out_h5["$gal/Match/residual"] = m.residual - if !isempty(m.inner_params) - for f in (:E, :F, :G, :H, :K, :M) - out_h5["$gal/Match/InnerParams/$(f)"] = [getfield(pp, f) for pp in m.inner_params] - end - # Literature names, matching the Tearing PerSurface mapping for the same fields. - out_h5["$gal/Match/InnerParams/tau_A"] = [pp.taua for pp in m.inner_params] - out_h5["$gal/Match/InnerParams/tau_R"] = [pp.taur for pp in m.inner_params] - out_h5["$gal/Match/InnerParams/dVdpsi"] = [pp.v1 for pp in m.inner_params] - end - end annotate_galerkin!(out_h5) return nothing end -# Metadata tables for the Galerkin outputs (Match/** is debug-only and exempt from the -# metadata contract; see docs/development/hdf5-conventions.md). +# Metadata tables for the Galerkin outputs (see docs/development/hdf5-conventions.md). const GALERKIN_H5_ANNOTATIONS = [ "ForceFreeStates/Solutions/GalerkinIntegration/rational_count" => (; long_name="number of rational (singular) surfaces in the Galerkin solve"), "ForceFreeStates/Solutions/GalerkinIntegration/psi" => (; long_name="normalized poloidal flux ψ_N grid of the closed Galerkin solution", scale="psi_gal"), diff --git a/src/ForceFreeStates/Galerkin/GalerkinStructs.jl b/src/ForceFreeStates/Galerkin/GalerkinStructs.jl index f4e79deb0..dbf0c8a28 100644 --- a/src/ForceFreeStates/Galerkin/GalerkinStructs.jl +++ b/src/ForceFreeStates/Galerkin/GalerkinStructs.jl @@ -145,46 +145,6 @@ struct GalerkinSolution cut_range::Matrix{Float64} end -""" - GalMatchResult - -Coil-driven RPEC matched solution from the outer↔inner asymptotic matching (Fortran rmatch -`match_rpec`). Populated by `gal_match_rpec` (GalerkinMatch.jl) when `gal_match_flag`. - - - `cout::Matrix{ComplexF64}` — `(2·msing, mcoil)` outer-region plasma-solution coefficients. - - `cin::Matrix{ComplexF64}` — `(2·msing, mcoil)` inner-region coefficients. - - `xi::Array{ComplexF64,3}`, `xi_deriv::Array{ComplexF64,3}` — `(mpert, ngrid, mcoil)` matched ξ(ψ) and - analytic ξ′(ψ) on the gal grid, one column per coil drive (identity-at-edge basis). - - `deltar::Matrix{ComplexF64}` — `(msing, 2)` inner-layer matching data `(Δ₁, Δ₂)` per surface. - - `bpen::Matrix{ComplexF64}` — `(msing, mcoil)` inner-layer penetrated (reconnected) resonant field at - each rational surface, one column per coil drive. Read off the GGJ inner solution at the layer center - (X=0) exactly as Fortran `match_output_solution` builds `intotsol_b` (match.f) — cusp-free, fit-free. - Zero for the ideal branch (`gal_ideal_flag`), where the inner layer is skipped. - - `inner_psi::Vector{Vector{Float64}}` — per surface, the inner-layer ψ grid `ψ_s ± X·x0/v1` (left wing - reversed then right wing, so ψ ascends through `ψ_s`). Empty in the ideal branch. - - `inner_xi::Vector{Matrix{ComplexF64}}` — per surface, the inner-layer displacement `ξ_ψ(ψ)` on - `inner_psi`, `(length(inner_psi[s]), mcoil)`, one column per coil drive. This is Fortran `match_solution`'s - `intotsol` (deltac component 2): `resc·(±Ξ₁·cin[2s] + Ξ₂·cin[2s-1])`, odd parity `Ξ₁` antisymmetric across - `ψ_s`, even parity `Ξ₂` symmetric, `resc=(v1/x0)^(1/2+p1)`. The raw inner-layer contribution only (no - regular outer background added), so it shares the singular asymptote with the outer near `ψ_s`. - - `rpec_eig::Vector{ComplexF64}` — forced eigenvalues `γ_s = 2πi·n·f_s` per surface. - - `residual::Float64` — relative linear-solve residual `‖mat·cof − rmat‖/‖rmat‖`. -""" -struct GalMatchResult - cout::Matrix{ComplexF64} - cin::Matrix{ComplexF64} - xi::Array{ComplexF64,3} - xi_deriv::Array{ComplexF64,3} - deltar::Matrix{ComplexF64} - bpen::Matrix{ComplexF64} - inner_psi::Vector{Vector{Float64}} - inner_xi::Vector{Matrix{ComplexF64}} - inner_b::Vector{Matrix{ComplexF64}} # per surface, matched inner-layer b^ψ(ψ) on inner_psi - inner_params::Vector{InnerLayer.GGJParameters} # per-surface layer coefficients from resist_eval - rpec_eig::Vector{ComplexF64} - residual::Float64 -end - """ GalerkinResult @@ -203,7 +163,6 @@ published as `ForceFreeStatesResult.delta_prime`. - `di::Vector{Float64}`, `alpha::Vector{ComplexF64}` — Mercier index and exponent per surface. - `solution::Union{Nothing,GalerkinSolution}` — reconstructed radial ξ(ψ) and analytic ξ′(ψ) on the gal-native grid; `nothing` if no resonant surfaces. - - `match::Union{Nothing,GalMatchResult}` — RPEC matched solution (`gal_match_flag`); `nothing` otherwise. """ struct GalerkinResult msing::Int @@ -214,5 +173,4 @@ struct GalerkinResult di::Vector{Float64} alpha::Vector{ComplexF64} solution::Union{Nothing,GalerkinSolution} - match::Union{Nothing,GalMatchResult} end diff --git a/src/ForceFreeStates/Integrators.jl b/src/ForceFreeStates/Integrators.jl index 8b123b2d0..7ef4567e1 100644 --- a/src/ForceFreeStates/Integrators.jl +++ b/src/ForceFreeStates/Integrators.jl @@ -59,8 +59,9 @@ matching `gal_*` control key without the prefix. - `dx1dx2_flag::Bool` - Enable the special dx1/dx2 treatment of resonant and extension elements. - `sing_order::Int` - Base power-series order for the singular asymptotics. - `sing_order_ceiling::Bool` - Auto-raise the order per surface for a high Mercier index. - - `rpec_flag::Bool` - Append the mpert coil-response columns to the Δ′ solve. Forced on when a [`ResistiveMatch`](@ref) is requested. + - `rpec_flag::Bool` - Append the mpert coil-response columns to the Δ′ solve. Required for a later `MatchProblem` solve on the result. - `edge_onesided::Bool` - Pack the two end intervals one-sided toward their single rational end instead of the Fortran symmetric pack. + - `cut_solution::Bool` - Also keep the cut solution, for a later `MatchProblem`'s composite inner profiles. """ @kwdef struct Galerkin <: AbstractIntegrator solver::String = "LU" @@ -78,46 +79,7 @@ matching `gal_*` control key without the prefix. sing_order_ceiling::Bool = true rpec_flag::Bool = false edge_onesided::Bool = false -end - -""" - ResistiveMatch(; eta=[], rho=[], rotation=[], gamma=5/3, ideal=false, inner_solver="ray", ...) - -Inner-layer matching configuration, passed to `solve` as `match=` and independent of the -integrator that produced the outer solution. Requesting a match closes the basis with a -resistive inner-layer solution instead of the ideal jump condition, so the result carries -`closure = :matched` and a non-zero `bpen`. - -Only [`Galerkin`](@ref) implements the match today; a `Riccati` or `Forward` solve with -`match` set errors. The per-surface vectors are ordered core to edge and must have one -entry per matched rational surface. - -## Fields - - - `eta::Vector{Float64}` - Per-surface resistivity η. - - `rho::Vector{Float64}` - Per-surface mass density ρ in kg/m³. - - `rotation::Vector{Float64}` - Per-surface rotation frequency f in Hz; the forced eigenvalue is γ_s = 2πi·n·f. - - `gamma::Float64` - Ratio of specific heats Γ in the resistive-layer coefficients. - - `ideal::Bool` - Build the ideal (perfectly shielded) matched solution: skip the inner layer and use the bare coil columns. `eta`, `rho` and `rotation` are then unread. - - `inner_solver::String` - Inner-layer Δ backend, `"ray"` (rotated-contour collocation) or `"galerkin"` (Hermite-cubic elements). - - `inner_xfac::Float64` - Asymptotic-matching radius multiplier of the `"galerkin"` backend. - - `inner_nx::Int` - Grid cells of the `"galerkin"` backend. - - `inner_nq::Int` - Quadrature order per cell of the `"galerkin"` backend. - - `inner_cutoff::Int` - Cells carrying the large solution as driving term in the `"galerkin"` backend. - - `inner_kmax::Int` - Large-x asymptotic series order of the `"galerkin"` backend. -""" -@kwdef struct ResistiveMatch - eta::Vector{Float64} = Float64[] - rho::Vector{Float64} = Float64[] - rotation::Vector{Float64} = Float64[] - gamma::Float64 = 5 / 3 - ideal::Bool = false - inner_solver::String = "ray" - inner_xfac::Float64 = 10.0 - inner_nx::Int = 1280 - inner_nq::Int = 5 - inner_cutoff::Int = 5 - inner_kmax::Int = 8 + cut_solution::Bool = false end """ @@ -165,26 +127,3 @@ function _apply_alg!(kwargs::Dict{Symbol,Any}, alg::Galerkin) return kwargs end -""" - _apply_match!(kwargs, match, alg) -> kwargs - -Translate a [`ResistiveMatch`](@ref) into the `gal_*` matching keywords, or error for an -integrator whose resonant matching is not implemented yet. `nothing` leaves `kwargs` alone, -which is the ideal-closure default. -""" -_apply_match!(kwargs::Dict{Symbol,Any}, ::Nothing, ::AbstractIntegrator) = kwargs - -function _apply_match!(kwargs::Dict{Symbol,Any}, ::ResistiveMatch, alg::AbstractIntegrator) - return error("resonant matching for this integrator is not yet implemented (requested with $(nameof(typeof(alg))))") -end - -function _apply_match!(kwargs::Dict{Symbol,Any}, match::ResistiveMatch, ::Galerkin) - kwargs[:gal_match_flag] = true - # The match consumes the coil-response columns, so it implies the rpec solve. - kwargs[:gal_rpec_flag] = true - for name in fieldnames(ResistiveMatch) - key = name === :ideal ? :gal_ideal_flag : Symbol(:gal_, name) - _set_ctrl!(kwargs, key, getfield(match, name), match) - end - return kwargs -end diff --git a/src/ForceFreeStates/Matching/LayerParameters.jl b/src/ForceFreeStates/Matching/LayerParameters.jl new file mode 100644 index 000000000..18c029a4b --- /dev/null +++ b/src/ForceFreeStates/Matching/LayerParameters.jl @@ -0,0 +1,77 @@ +# LayerParameters.jl +# +# Per-surface η, ρ and rotation for the GGJ inner layer, from kinetic profiles. + +using ..Utilities.PhysicalConstants: M_P, E_CHG +using ..Utilities: KineticProfiles +using ..Utilities.NeoclassicalResistivity: NeoResistivityModel, SpitzerModel, SauterNeoModel, + coulomb_log_e, eta_spitzer, nu_star_e, eta_neoclassical + +""" + layer_parameters(surfaces, equil; profiles=equil.kinetic, eta=nothing, rho=nothing, + rotation=nothing, mu_i=2.0, zeff=1.0, resistivity_model=SauterNeoModel(), + lnLambda_form=:nrl) -> (; eta, rho, rotation) + +Resistivity η in Ω·m, mass density ρ = μᵢ·m_p·n_e in kg/m³ and E×B rotation frequency +f = ω_E/2π in Hz at each surface, from `profiles`: the attached kinetic profiles, or a +`KineticProfiles` table. `eta`, `rho` or `rotation` vectors, one value per surface, override +the derived values. The default η is Sauter neoclassical; `resistivity_model=SpitzerModel()` +gives Spitzer. +""" +function layer_parameters( + surfaces::AbstractVector, + equil; + profiles=equil === nothing ? nothing : equil.kinetic, + eta::Union{Nothing,AbstractVector{<:Real}}=nothing, + rho::Union{Nothing,AbstractVector{<:Real}}=nothing, + rotation::Union{Nothing,AbstractVector{<:Real}}=nothing, + mu_i::Real=2.0, + zeff::Real=1.0, + resistivity_model::NeoResistivityModel=SauterNeoModel(), + lnLambda_form::Symbol=:nrl +) + msing = length(surfaces) + for (name, v) in (("eta", eta), ("rho", rho), ("rotation", rotation)) + v === nothing || length(v) == msing || + error("layer_parameters: $name has length $(length(v)), expected one value per surface (msing=$msing, core to edge)") + end + + needs_derivation = eta === nothing || rho === nothing || rotation === nothing + if needs_derivation + profiles === nothing && + error("layer_parameters: deriving η/ρ/rotation needs kinetic profiles on the equilibrium; attach them " * + "before solving (or to ffs.equil after), or pass all three of eta, rho, rotation") + + eta_out = Vector{Float64}(undef, msing) + rho_out = Vector{Float64}(undef, msing) + rot_out = Vector{Float64}(undef, msing) + for (k, sing) in enumerate(surfaces) + n_e, t_e, omega_E = _layer_profile(profiles, sing.psifac) + if eta === nothing + lnLamb = coulomb_log_e(n_e, t_e; form=lnLambda_form) + if resistivity_model isa SpitzerModel + eta_out[k] = eta_spitzer(n_e, t_e, zeff; lnLamb=lnLamb) + else + rg = sing.restype + rg === nothing && + error("layer_parameters: surface $k has restype = nothing — the neoclassical resistivity " * + "closure needs the trapped fraction from resist_eval_all!; run it first or use SpitzerModel()") + nuestar = nu_star_e(n_e, t_e, rg.R_major, rg.eps_local, sing.q, zeff; lnLamb=lnLamb) + eta_out[k] = eta_neoclassical(resistivity_model, n_e, t_e, zeff, rg.f_trap, nuestar; lnLamb=lnLamb) + end + end + rho_out[k] = mu_i * M_P * n_e + rot_out[k] = omega_E / (2π) + end + end + + return ( + eta=eta === nothing ? eta_out : collect(Float64, eta), + rho=rho === nothing ? rho_out : collect(Float64, rho), + rotation=rotation === nothing ? rot_out : collect(Float64, rotation) + ) +end + +# n_e in m⁻³, T_e in eV, and the E×B frequency at ψ from either profile container. +_layer_profile(kp::Equilibrium.KineticProfileSplines, ψ) = (kp.ne_spline(ψ), kp.Te_spline(ψ) / E_CHG, kp.omegaE_spline(ψ)) +_layer_profile(kp::KineticProfiles, ψ) = (p = kp(ψ); (p.n_e, p.T_e, p.omega)) diff --git a/src/ForceFreeStates/Matching/MatchProblem.jl b/src/ForceFreeStates/Matching/MatchProblem.jl new file mode 100644 index 000000000..8d882e6e6 --- /dev/null +++ b/src/ForceFreeStates/Matching/MatchProblem.jl @@ -0,0 +1,262 @@ +# MatchProblem.jl +# +# Driven inner-layer matching on a finished solve (rmatch match_rpec + match_output_solution, match.f). + +""" + MatchProblem(ffs; eta=nothing, rho=nothing, rotation=nothing, gamma=5/3, ideal=false, + mu_i=2.0, zeff=1.0, resistivity_model=SauterNeoModel(), lnLambda_form=:nrl) + +Driven inner-layer matching on a finished force-free-states solve `ffs`: each rational surface +is matched at the forced eigenvalue γ = 2πi·n·f, with f the plasma rotation frequency there. +`solve(prob, GGJModel())` returns a copy of `ffs` with the matched solution and the penetrated +field `bpen`. `ffs` needs the Δ′ coil block, from `Galerkin(; rpec_flag=true)` or `Riccati()` +with `vac_flag=true`. + +# Keywords +- `eta`, `rho`, `rotation`: per-surface resistivity in Ω·m, mass density in kg/m³ and rotation + frequency in Hz, core to edge. Any left `nothing` is derived from `ffs.equil.kinetic` by + [`layer_parameters`](@ref) with `mu_i`, `zeff`, `resistivity_model` and `lnLambda_form`. +- `gamma`: ratio of specific heats. +- `ideal`: skip the inner layer and return the perfectly shielding solution. + +Fields: `ffs`, the matched `surfaces`, and the resolved `eta`, `rho`, `rotation`, `gamma`, `ideal`. + +Ref: Fortran rmatch `match_rpec` (match.f). +""" +struct MatchProblem{R<:ForceFreeStatesResult} + ffs::R + surfaces::Vector{SingType} + eta::Vector{Float64} + rho::Vector{Float64} + rotation::Vector{Float64} + gamma::Float64 + ideal::Bool +end + +function MatchProblem( + ffs::ForceFreeStatesResult; + eta::Union{Nothing,AbstractVector{<:Real}}=nothing, + rho::Union{Nothing,AbstractVector{<:Real}}=nothing, + rotation::Union{Nothing,AbstractVector{<:Real}}=nothing, + gamma::Real=5 / 3, + ideal::Bool=false, + mu_i::Real=2.0, + zeff::Real=1.0, + resistivity_model::NeoResistivityModel=SauterNeoModel(), + lnLambda_form::Symbol=:nrl +) + ffs.npert == 1 || error("MatchProblem: matching is single-n; the result spans n = $(ffs.nlow):$(ffs.nhigh)") + ffs.closure === :matched && error("MatchProblem: the result is already matched; pose the match on the outer solve") + dp = ffs.delta_prime + dp === nothing && + error("MatchProblem: the $(ffs.integrator) result has no Δ′; solve with Galerkin(; rpec_flag=true), or Riccati() with vac_flag=true") + isempty(dp.coil) && + error("MatchProblem: the Δ′ coil block is empty; solve with Galerkin(; rpec_flag=true), or Riccati() with vac_flag=true") + + sings = _matched_surfaces(ffs) + size(dp.raw, 1) == 2 * length(sings) || + error("MatchProblem: Δ′ raw block is $(size(dp.raw, 1))×$(size(dp.raw, 2)) but the result carries $(length(sings)) matchable surfaces") + + if ideal + msing = length(sings) + return MatchProblem(ffs, sings, zeros(msing), zeros(msing), zeros(msing), Float64(gamma), true) + end + params = layer_parameters(sings, ffs.equil; eta=eta, rho=rho, rotation=rotation, + mu_i=mu_i, zeff=zeff, resistivity_model=resistivity_model, lnLambda_form=lnLambda_form) + return MatchProblem(ffs, sings, params.eta, params.rho, params.rotation, Float64(gamma), false) +end + +# The surfaces the Δ′ rows are indexed by (see DeltaPrimeData): all of them for Riccati, the +# Galerkin solve's own subset otherwise. +_matched_surfaces(ffs::ForceFreeStatesResult) = + ffs.galerkin === nothing ? copy(ffs.surfaces) : [s for s in ffs.surfaces if s.psifac in ffs.galerkin.sing_psi] + +""" + closure_capable(model) -> Bool + +Whether `model` can solve a [`MatchProblem`](@ref): `GGJModel` with the `:ray` or `:galerkin` +backend can (they reconstruct the layer profiles); `SLAYERModel` (slab, tearing parity only) +cannot. +""" +closure_capable(::InnerLayer.GGJModel{S}) where {S} = S in (:ray, :galerkin) +closure_capable(::InnerLayer.SLAYERModel) = false + +""" + solve(prob::MatchProblem, model) -> ForceFreeStatesResult + +Match `prob.ffs` to the inner layer of `model`. The result has `closure = :matched` (`:ideal` +with `ideal=true`), `bpen`, the [`MatchResult`](@ref) in `result.match` and, for a Galerkin +solve, the matched ξ. A `:matched` result drops the ideal δW (`wp`, `free_boundary`), so +PerturbedEquilibrium skips its response rather than mixing ideal and matched physics. +""" +function CommonSolve.solve(prob::MatchProblem, model::InnerLayer.InnerLayerModel) + closure_capable(model) || + error("$(nameof(typeof(model))) cannot solve a MatchProblem (use GGJModel(); SLAYERModel is for TearingProblem)") + match = _compute_match(prob, model) + return _matched_result(prob.ffs, match, prob.ideal) +end + +# Inner-layer Δ per surface, the matching system, then the matched outer and inner profiles. +function _compute_match(prob::MatchProblem, model::InnerLayer.GGJModel) + ffs = prob.ffs + dp = ffs.delta_prime + sings = prob.surfaces + msing = length(sings) + mcoil = size(dp.coil, 2) + nn = ffs.nlow + equil = ffs.equil + chi1 = 2π * equil.psio + + gal_sol = ffs.galerkin === nothing ? nothing : ffs.galerkin.solution + + if prob.ideal + # Ideal limit (rmatch coil%ideal_flag): no inner layer, cout = 0. + cout = zeros(ComplexF64, 2msing, mcoil) + cin = zeros(ComplexF64, 2msing, mcoil) + deltar = zeros(ComplexF64, msing, 2) + bpen = zeros(ComplexF64, msing, mcoil) + inner_psi = Vector{Float64}[] + inner_xi = Matrix{ComplexF64}[] + inner_b = Matrix{ComplexF64}[] + inner_params = InnerLayer.GGJParameters[] + rpec_eig = zeros(ComplexF64, msing) + residual = 0.0 + else + # Inner-layer Δ and profiles per surface (match.f deltac_run, intotsol). + deltar = zeros(ComplexF64, msing, 2) + rpec_eig = zeros(ComplexF64, msing) + # Layer-center field weight per parity (match.f intotsol_b). + pen = zeros(ComplexF64, msing, 2) + inner_psi = Vector{Vector{Float64}}(undef, msing) # ψ_s ± X·dψdx, ascending + inner_odd = Vector{Vector{ComplexF64}}(undef, msing) # Ξ₁, odd about ψ_s + inner_even = Vector{Vector{ComplexF64}}(undef, msing) # Ξ₂, even + inner_bodd = Vector{Vector{ComplexF64}}(undef, msing) # b^ψ from Ψ₁, even + inner_beven = Vector{Vector{ComplexF64}}(undef, msing) # b^ψ from Ψ₂, odd + inner_params = Vector{InnerLayer.GGJParameters}(undef, msing) + for i in 1:msing + params = ggj_parameters(sings[i], equil; eta=prob.eta[i], rho=prob.rho[i], + gamma=prob.gamma, ising=i) + inner_params[i] = params + γ = 2π * im * nn * prob.rotation[i] # forced eigenvalue, f in Hz + rpec_eig[i] = γ + inner = InnerLayer.solve_inner_profile(model, params, γ) + deltar[i, 1] = inner.Δ[1] + deltar[i, 2] = inner.Δ[2] + profΨ, profΞ, xg = inner.Ψ, inner.Ξ, inner.x + resc = inner.rescale # inner big-branch amplitude → outer δψ normalization + # b_m = 2πi·χ₁·(m−nq)·ξ_m with m−nq = −n·q′·dψdx·X and Ψ = XΞ (GWP2016 Eq. 16) ⇒ b_m = scale·resc·Ψ. + scale = -2π * chi1 * im * nn * sings[i].q1 * inner.dψdx + pen[i, 1] = scale * profΨ[1, 1] * resc # X = 0, parity 1 + pen[i, 2] = scale * profΨ[1, 2] * resc # parity 2 (Ψ(0) = 0) + xvar = xg .* inner.dψdx # X → ψ distance + inner_psi[i] = vcat(reverse(sings[i].psifac .- xvar), sings[i].psifac .+ xvar) + inner_odd[i] = resc .* vcat(reverse(.-profΞ[:, 1]), profΞ[:, 1]) + inner_even[i] = resc .* vcat(reverse(profΞ[:, 2]), profΞ[:, 2]) + inner_bodd[i] = (scale * resc) .* vcat(reverse(profΨ[:, 1]), profΨ[:, 1]) + inner_beven[i] = (scale * resc) .* vcat(reverse(.-profΨ[:, 2]), profΨ[:, 2]) + end + + cout, cin, residual = _match_system(dp.raw, dp.coil, deltar) + + # Penetrated field at each layer center (match.f intotsol_b). + bpen = zeros(ComplexF64, msing, mcoil) + for i in 1:msing, j in 1:mcoil + bpen[i, j] = pen[i, 1] * cin[2i, j] + pen[i, 2] * cin[2i-1, j] + end + + # Matched inner ξ_ψ and b^ψ per coil: odd parity × cin[2i], even × cin[2i-1] (match.f intotsol). + inner_xi = [inner_odd[i] * transpose(cin[2i, :]) .+ inner_even[i] * transpose(cin[2i-1, :]) for i in 1:msing] + inner_b = [inner_bodd[i] * transpose(cin[2i, :]) .+ inner_beven[i] * transpose(cin[2i-1, :]) for i in 1:msing] + end + + reconnected_flux = Matrix{ComplexF64}(dp.coil .+ transpose(dp.raw) * cout) + + # Matched outer ξ per coil: coil column + Σ cout·plasma columns (needs the Galerkin basis). + if gal_sol === nothing + xi = Array{ComplexF64,3}(undef, 0, 0, 0) + xi_deriv = Array{ComplexF64,3}(undef, 0, 0, 0) + empty!(inner_psi); empty!(inner_xi); empty!(inner_b) + else + sols = gal_sol.xi # (mpert, ngrid, nsol) + sols_deriv = gal_sol.xi_deriv + mpert = size(sols, 1) + ngrid = size(sols, 2) + xi = zeros(ComplexF64, mpert, ngrid, mcoil) + xi_deriv = zeros(ComplexF64, mpert, ngrid, mcoil) + for j in 1:mcoil + csol = 2msing + j + @views xi[:, :, j] .= sols[:, :, csol] + @views xi_deriv[:, :, j] .= sols_deriv[:, :, csol] + for isol in 1:2msing + @views xi[:, :, j] .+= cout[isol, j] .* sols[:, :, isol] + @views xi_deriv[:, :, j] .+= cout[isol, j] .* sols_deriv[:, :, isol] + end + end + + # Composite inner solution: add back the outer background the cut removed (match.f intotsol). + if !isempty(inner_psi) + sols_cut = gal_sol.xi_cut + if isempty(sols_cut) + @warn "MatchProblem: the solve retained no cut solution (set cut_solution=true on Galerkin), " * + "so the composite inner-region profiles are skipped; bpen and the matched outer ξ are unaffected" + empty!(inner_psi); empty!(inner_xi); empty!(inner_b) + else + cut_range = gal_sol.cut_range + keep = .!gal_sol.issing + psi_keep = gal_sol.psi[keep] + for i in 1:msing + m_res = round(Int, nn * sings[i].q) + ires = m_res - ffs.mlow + 1 + # Keep only points inside the cut window. + lo, hi = cut_range[i, 1], cut_range[i, 2] + inside = findall(p -> lo <= p <= hi, inner_psi[i]) + if isempty(inside) + @warn "MatchProblem: inner layer at ψ=$(round(sings[i].psifac, digits=5)) lies outside the " * + "resonant/extension cells (ψ ∈ [$lo, $hi]); raise gal_dx1/gal_dx2 to overlap the layer" surface = i + continue + end + length(inside) < length(inner_psi[i]) && @info "MatchProblem: surface $i inner region clipped to the " * + "cut window ($(length(inside)) of $(length(inner_psi[i])) points)" + inner_psi[i] = inner_psi[i][inside] + inner_xi[i] = inner_xi[i][inside, :] + inner_b[i] = inner_b[i][inside, :] + + # Cut outer background of the resonant harmonic, per coil. + cutmn = Matrix{ComplexF64}(undef, length(psi_keep), mcoil) + for j in 1:mcoil + @views cutmn[:, j] .= sols_cut[ires, keep, 2msing+j] + for isol in 1:2msing + @views cutmn[:, j] .+= cout[isol, j] .* sols_cut[ires, keep, isol] + end + end + itp = cubic_interp(psi_keep, Series(cutmn); extrap=ExtendExtrap()) + buf = Vector{ComplexF64}(undef, mcoil) + hint = Ref(1) + for (ip, psi_p) in enumerate(inner_psi[i]) + itp(buf, psi_p; hint=hint) + singfac = m_res - nn * equil.profiles.q_spline(psi_p) + @views inner_xi[i][ip, :] .+= buf + @views inner_b[i][ip, :] .+= (2π * im * chi1 * singfac) .* buf + end + end + end + end + end + + return MatchResult(cout, cin, deltar, rpec_eig, residual, bpen, reconnected_flux, + xi, xi_deriv, inner_psi, inner_xi, inner_b, inner_params) +end + +# Copy of ffs with the match applied; the ideal reference keeps closure :ideal, its zero bpen and δW. +function _matched_result(ffs::ForceFreeStatesResult, match::MatchResult, ideal::Bool) + gal = ffs.galerkin + solution = (gal !== nothing && gal.solution !== nothing && !isempty(match.xi)) ? + _matched_gal_profiles(gal, match, ffs.mats, ffs) : ffs.solution + ideal && return _with(ffs; match=match, solution=solution) + # The ideal δW does not describe the matched plasma; PE gates on its absence. + return _with(ffs; closure=:matched, bpen=match.bpen, match=match, solution=solution, + wp=nothing, free_boundary=nothing) +end + +# A copy of the immutable `x` with the named fields replaced. +_with(x; fields...) = typeof(x).name.wrapper((get(fields, f, getfield(x, f)) for f in fieldnames(typeof(x)))...) diff --git a/src/ForceFreeStates/Matching/ResonantMatch.jl b/src/ForceFreeStates/Matching/ResonantMatch.jl index 4f3fd3056..354fc4b15 100644 --- a/src/ForceFreeStates/Matching/ResonantMatch.jl +++ b/src/ForceFreeStates/Matching/ResonantMatch.jl @@ -1,80 +1,96 @@ -# Outer<->inner resistive match, Wang et al. 2020 (PoP 27, 122509) Eq. 11: -# C = -(Δ_out - Δ_in(i2πf))^{-1} Δ_coil -# Raw STRIDE outer Δ' + raw coil drive matched to the GGJ inner layer (resist_eval -> solve_inner). -struct ResonantMatchResult - cout::Matrix{ComplexF64} # outer coeffs (2msing × ncoil) - cin::Matrix{ComplexF64} # inner coeffs (2msing × ncoil) - deltar::Matrix{ComplexF64} # inner-layer Δ per surface (msing × 2) - rpec_eig::Vector{ComplexF64} # forced eigenvalue γ_s = 2πi·n·f - reconnected_flux::Matrix{ComplexF64} # reconnected resonant flux (2msing × ncoil) - bpen::Matrix{ComplexF64} # area-weighted penetrated field (msing × ncoil); empty until Stage 2 - residual::Float64 -end - -function resonant_match_rpec(delta_out_raw::AbstractMatrix, delta_coil_raw::AbstractMatrix, - sings::Vector{SingType}, equil::Equilibrium.PlasmaEquilibrium, - intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl) +# ResonantMatch.jl +# +# Matching result and the 4·msing outer↔inner system (rmatch match_rpec; Wang et al. 2020 PoP 27 122509, Eq. 11). - msing = size(delta_out_raw, 1) ÷ 2 - ncoil = size(delta_coil_raw, 2) - nn = intr.nlow - empty_bpen = Matrix{ComplexF64}(undef, 0, 0) +""" + MatchResult - size(delta_out_raw) == (2msing, 2msing) || error("delta_out_raw $(size(delta_out_raw)) != (2msing,2msing)") - length(sings) == msing || error("sings $(length(sings)) != msing $msing") - size(delta_coil_raw, 1) == 2msing || error("delta_coil_raw rows $(size(delta_coil_raw,1)) != 2msing") - - if ctrl.gal_ideal_flag # ideal limit: no inner layer, no reconnection - return ResonantMatchResult(zeros(ComplexF64,2msing,ncoil), zeros(ComplexF64,2msing,ncoil), - zeros(ComplexF64,msing,2), zeros(ComplexF64,msing), Matrix{ComplexF64}(delta_coil_raw), empty_bpen, 0.0) - end - for (nm,v) in (("gal_eta",ctrl.gal_eta),("gal_rho",ctrl.gal_rho),("gal_rotation",ctrl.gal_rotation)) - length(v) == msing || error("$nm length $(length(v)) != msing $msing") - end +Matching coefficients and matched profiles, found in `result.match`. The profile +fields are empty when the solve kept no Galerkin basis; the inner fields are empty with `ideal`. - chi1 = 2π * equil.psio - deltar = zeros(ComplexF64, msing, 2) - rpec_eig = zeros(ComplexF64, msing) - # Layer-center (X=0) penetrated-field weights pen[i,k] = scale·Ψ_k(0)·rescale (match.f intotsol_b); - # solve_inner_profile returns the same Δ as solve_inner plus the inner-layer field needed for pen. - pen = zeros(ComplexF64, msing, 2) - for i in 1:msing - params = resist_eval(sings[i], equil, intr; eta=ctrl.gal_eta[i], rho=ctrl.gal_rho[i], gamma=ctrl.gal_gamma, ising=i) - γ = 2π*im*nn*ctrl.gal_rotation[i] - rpec_eig[i] = γ - inner = InnerLayer.solve_inner_profile(InnerLayer.GGJModel(; solver=:galerkin), params, γ; - xfac=ctrl.gal_inner_xfac, nx=ctrl.gal_inner_nx, nq=ctrl.gal_inner_nq, cutoff=ctrl.gal_inner_cutoff, kmax=ctrl.gal_inner_kmax) - deltar[i,1] = inner.Δ[1]; deltar[i,2] = inner.Δ[2] - scale = -2π * chi1 * im * nn * sings[i].q1 * inner.dψdx # b_m = −2πi·χ₁·n·q′·dψdx·rescale·Ψ (GalerkinMatch.jl) - pen[i,1] = scale * inner.Ψ[1,1] * inner.rescale # parity 1 (Ψ(0)≠0) - pen[i,2] = scale * inner.Ψ[1,2] * inner.rescale # parity 2 (Ψ(0)=0 ⇒ ~0) - end +## Fields + - `cout`, `cin`: `(2msing, ncoil)` outer and inner coefficients. + - `deltar`: `(msing, 2)` inner-layer Δ per surface and parity. + - `rpec_eig`: forced eigenvalue γ = 2πi·n·f per surface. + - `residual`: relative residual of the linear solve. + - `bpen`: `(msing, ncoil)` penetrated resonant field at each layer center. + - `reconnected_flux`: `(2msing, ncoil)` `Δ_coil + Δ_outᵀ·cout`. + - `xi`, `xi_deriv`: `(mpert, ngrid, ncoil)` matched outer ξ and ξ′, one column per coil. + - `inner_psi`, `inner_xi`, `inner_b`: per surface, the inner ψ grid and the composite ξ_ψ and b^ψ on it. + - `inner_params`: per-surface `GGJParameters` used. +""" +struct MatchResult + cout::Matrix{ComplexF64} + cin::Matrix{ComplexF64} + deltar::Matrix{ComplexF64} + rpec_eig::Vector{ComplexF64} + residual::Float64 + bpen::Matrix{ComplexF64} + reconnected_flux::Matrix{ComplexF64} + xi::Array{ComplexF64,3} + xi_deriv::Array{ComplexF64,3} + inner_psi::Vector{Vector{Float64}} + inner_xi::Vector{Matrix{ComplexF64}} + inner_b::Vector{Matrix{ComplexF64}} + inner_params::Vector{InnerLayer.GGJParameters} +end - mat = zeros(ComplexF64, 4msing, 4msing) +# Solve the 4·msing system mat·[cout; cin] = rmat (rmatch match_rpec). +function _match_system(dp_raw::AbstractMatrix, dp_coil::AbstractMatrix, deltar::AbstractMatrix) + msing = size(dp_raw, 1) ÷ 2 + ncoil = size(dp_coil, 2) + mat = zeros(ComplexF64, 4msing, 4msing) rmat = zeros(ComplexF64, 4msing, ncoil) - @views mat[2msing+1:4msing, 1:2msing] .= transpose(delta_out_raw) - @views rmat[2msing+1:4msing, :] .= .-delta_coil_raw - for i in 1:msing - a=2i-1; b=2i; c=a+2msing; d=b+2msing - d1=deltar[i,1]; d2=deltar[i,2] - mat[a,a]=1; mat[b,b]=1 - mat[a,c]=-1; mat[a,d]=1 - mat[b,c]=-1; mat[b,d]=-1 - mat[c,c]=-d1; mat[c,d]=d2 - mat[d,c]=-d1; mat[d,d]=-d2 + @views mat[(2msing+1):4msing, 1:2msing] .= transpose(dp_raw) # Δ_out + @views rmat[(2msing+1):4msing, :] .= .-dp_coil # −Δ_coil + for ising in 1:msing + idx1 = 2ising - 1 + idx2 = 2ising + idx3 = idx1 + 2msing + idx4 = idx2 + 2msing + delta1 = deltar[ising, 1] + delta2 = deltar[ising, 2] + mat[idx1, idx1] = 1 + mat[idx2, idx2] = 1 + mat[idx1, idx3] = -1 + mat[idx1, idx4] = 1 + mat[idx2, idx3] = -1 + mat[idx2, idx4] = -1 + # inner-layer Δ block (match.f signs) + mat[idx3, idx3] = -delta1 + mat[idx3, idx4] = delta2 + mat[idx4, idx3] = -delta1 + mat[idx4, idx4] = -delta2 end cof = mat \ rmat - residual = norm(mat*cof - rmat) / max(norm(rmat), 1e-300) - cout = cof[1:2msing, :] - cin = cof[2msing+1:4msing, :] - - reconnected_flux = delta_coil_raw .+ transpose(delta_out_raw)*cout - # Inner-layer penetrated (reconnected) resonant field per surface — ONE quantity per surface, read off - # the inner solution at the layer center (match.f intotsol_b; GalerkinMatch.jl): bpen[i,j] = pen₁(i)·cin[2i,j] + pen₂(i)·cin[2i-1,j]. - bpen = zeros(ComplexF64, msing, ncoil) - for i in 1:msing, j in 1:ncoil - bpen[i,j] = pen[i,1]*cin[2i,j] + pen[i,2]*cin[2i-1,j] + residual = norm(mat * cof - rmat) / max(norm(rmat), 1e-300) + return cof[1:2msing, :], cof[(2msing+1):4msing, :], residual +end + +# Write `m` under SingularSurfaces/Match/ (debug-only: exempt from the metadata contract). +function write_match!(out_h5, m::MatchResult) + g = "SingularSurfaces/Match" + out_h5["$g/cout"] = m.cout + out_h5["$g/cin"] = m.cin + out_h5["$g/Delta_r"] = m.deltar + out_h5["$g/bpen"] = m.bpen + out_h5["$g/rpec_eig"] = m.rpec_eig + # Per-surface inner-layer ξ_ψ(ψ) (match.f intotsol); ragged grids → one dataset pair per surface. + for i in eachindex(m.inner_psi) + out_h5["$g/Inner/psi_$i"] = m.inner_psi[i] + out_h5["$g/Inner/xi_$i"] = m.inner_xi[i] + out_h5["$g/Inner/b_$i"] = m.inner_b[i] + end + out_h5["$g/residual"] = m.residual + if !isempty(m.inner_params) + for f in (:E, :F, :G, :H, :K, :M) + out_h5["$g/InnerParams/$(f)"] = [getfield(pp, f) for pp in m.inner_params] + end + # Literature names, matching the Tearing PerSurface mapping for the same fields. + out_h5["$g/InnerParams/tau_A"] = [pp.taua for pp in m.inner_params] + out_h5["$g/InnerParams/tau_R"] = [pp.taur for pp in m.inner_params] + out_h5["$g/InnerParams/dVdpsi"] = [pp.v1 for pp in m.inner_params] end - return ResonantMatchResult(cout, cin, deltar, rpec_eig, reconnected_flux, bpen, residual) + return nothing end diff --git a/src/ForceFreeStates/Result.jl b/src/ForceFreeStates/Result.jl index da8fe3648..721ec5615 100644 --- a/src/ForceFreeStates/Result.jl +++ b/src/ForceFreeStates/Result.jl @@ -58,9 +58,11 @@ so bpen and closure are always present. - `surfaces::Vector{SingType}` - Ideal singular surfaces in the integration domain, with asymptotic bases and GGJ coefficients. - `kinetic::NamedTuple` - Kinetic singular-surface scan (`kmsing`, `kinsing`, `scan_psi`, `scan_cond`, `scan_threshold`); empty unless the finder ran. - `closure::Symbol` - How the basis is closed at the rationals: `:ideal` (the ideal jump - condition was imposed) or `:matched` (an inner-layer solution was matched in). + condition was imposed) or `:matched` (an inner-layer solution was matched in). A + `:matched` result carries no ideal δW: `wp` and `free_boundary` are `nothing`. - `bpen::Matrix{ComplexF64}` - `(msing × numpert_total)` penetrated resonant field from the inner layer, per surface and driving mode; all zeros under `:ideal` closure. + - `match::Union{Nothing,MatchResult}` - The [`MatchResult`](@ref) of a `MatchProblem` solve; `nothing` otherwise. - `solution::Union{Nothing,SolutionProfiles}` - THE solve's ξ solution; `nothing` when the formalism produced none (Riccati, and Galerkin without a match). - `diagnostics::Union{Nothing,OdeState}` - The integrator's raw ODE state (ψ trace, `crit`, @@ -73,8 +75,7 @@ so bpen and closure are always present. - `delta_prime::Union{Nothing,DeltaPrimeData}` - Δ′/outer-region matching payload from whichever formalism ran (Riccati BVP or Galerkin); `nothing` when neither produced one. - `galerkin::Union{Nothing,GalerkinResult}` - RDCON Galerkin solver internals and FEM - diagnostics, including the RPEC inner-layer match when requested. Its Δ′ payload lives - in `delta_prime`, not here. + diagnostics. Its Δ′ payload lives in `delta_prime`, not here. """ struct ForceFreeStatesResult{E<:Equilibrium.PlasmaEquilibrium,F<:MatrixSplines} <: ModeSpace integrator::Symbol @@ -106,6 +107,7 @@ struct ForceFreeStatesResult{E<:Equilibrium.PlasmaEquilibrium,F<:MatrixSplines} # Closure of the basis at the rationals, always present. closure::Symbol bpen::Matrix{ComplexF64} + match::Union{Nothing,MatchResult} # Per-formalism products; presence is the capability signal. solution::Union{Nothing,SolutionProfiles} @@ -124,7 +126,10 @@ that produced `result`. Warns naming the calculation being skipped otherwise. """ function require(result::ForceFreeStatesResult, field::Symbol, calc::AbstractString) getfield(result, field) === nothing || return true - @warn "Skipping $calc: `$field` was not produced by the $(result.integrator) integrator" + why = result.closure === :matched && field in (:wp, :free_boundary) ? + "a :matched result carries no δW (the resistive δW is not implemented)" : + "`$field` was not produced by the $(result.integrator) integrator" + @warn "Skipping $calc: $why" return false end @@ -147,9 +152,8 @@ end # analytic Galerkin derivative rather than a differenced value spline, and # Ξ_s = −A⁻¹(B·Ξ′ + C·Ξ) is the same outer ideal-MHD relation `sing_der!` uses. The grid runs # inner→edge, so the last node is the control surface and carries the edge boundary condition. -function _matched_gal_profiles(gal_result::GalerkinResult, mats::MatrixSplines, intr::ModeSpace) +function _matched_gal_profiles(gal_result::GalerkinResult, m::MatchResult, mats::MatrixSplines, intr::ModeSpace) sol = gal_result.solution - m = gal_result.match npert = intr.numpert_total keep = .!sol.issing @@ -183,9 +187,9 @@ end Assemble the published result once the solve is finished — the one place that decides what a formalism's raw output means downstream. Beyond materializing the forward path's derivative -stores, packing the matched Galerkin solution, and forming the fixed-boundary `W_p` when the -free-boundary stage did not run, every field is copied or aliased from what the stages -already produced. +stores and forming the fixed-boundary `W_p` when the free-boundary stage did not run, every +field is copied or aliased from what the stages already produced. The integrator always +publishes the ideal closure; a `MatchProblem` solve transforms the result afterwards. `odet` is the integrator's ODE state (`nothing` for Galerkin); `gal_data`/`gal_dp` the Galerkin solver internals and its Δ′ payload (`nothing` otherwise). The two formalisms are never both @@ -203,16 +207,12 @@ function build_result( gal_data::Union{Nothing,GalerkinResult}, gal_dp::Union{Nothing,DeltaPrimeData} ) - matched = gal_data !== nothing && gal_data.match !== nothing - # The forward sweep is the only formalism whose stores need materializing; doing it here # keeps `SolutionProfiles.du_store`/`xi_s_store` populated by construction. solution = if integrator === :forward && odet !== nothing materialize_derivative_stores!(odet, equil, mats, intr) SolutionProfiles(:el_axis, odet.step, odet.psi_store, odet.q_store, odet.u_store, odet.du_store, odet.xi_s_store) - elseif matched - _matched_gal_profiles(gal_data, mats, intr) else nothing end @@ -229,11 +229,8 @@ function build_result( kinetic = (kmsing=intr.kmsing, kinsing=intr.kinsing, scan_psi=intr.kinsing_scan_psi, scan_cond=intr.kinsing_scan_cond, scan_threshold=intr.kinsing_scan_threshold) - # The ideal-flag match deliberately skips the inner-layer Δ, so its basis is ideal-closed - # and carries no penetrated field. - closure = (matched && !ctrl.gal_ideal_flag) ? :matched : :ideal - bpen = closure === :matched ? gal_data.match.bpen : - zeros(ComplexF64, intr.msing, intr.numpert_total) + closure = :ideal + bpen = zeros(ComplexF64, intr.msing, intr.numpert_total) # Fixed-boundary plasma energy matrix at the edge; free of any vacuum dependence, so a # vac_flag=false run still publishes its energy product. Aliases free_run's when it ran. @@ -250,7 +247,7 @@ function build_result( intr.mlow, intr.mhigh, intr.mpert, intr.nlow, intr.nhigh, intr.npert, intr.numpert_total, intr.psilow, intr.psilim, intr.qlim, intr.q1lim, intr.dir_path, intr.wall_settings, intr.debug_settings, metric, mats, intr.sing, kinetic, - closure, bpen, + closure, bpen, nothing, solution, odet, wp, free_energies, delta_prime, gal_data ) end diff --git a/src/ForceFreeStates/Surfaces/Resist.jl b/src/ForceFreeStates/Surfaces/Resist.jl deleted file mode 100644 index fd4a4b46a..000000000 --- a/src/ForceFreeStates/Surfaces/Resist.jl +++ /dev/null @@ -1,126 +0,0 @@ -# Resist.jl -# -# Per-rational-surface resistive inner-layer parameters (Glasser–Greene–Johnson). -# Julia port of the RDCON `resist_eval` (resist.f): evaluates the flux-surface-averaged -# curvature/mass coefficients E, F, H, M, G, K and the local Alfvén / resistive time scales -# τ_A, τ_R at a singular surface. -# -# Returns an `InnerLayer.GGJParameters` directly (the inner-layer matching input type), so a -# resistive ForceFreeStates run can hand its per-surface coefficients straight to the GGJ inner -# solver. InnerLayer is included before ForceFreeStates precisely so this type is available here. -# -# Resistivity η and mass density ρ are REQUIRED inputs (per surface) — there is no built-in -# resistivity/density model. We deliberately do NOT port the Fortran `resist.f` Spitzer block -# (hardcoded ne/Te → η, ρ); the caller supplies η and ρ from whatever transport/profile model is -# appropriate. τ_A and τ_R are then formed from physically separable pieces: a purely geometric -# factor times the explicit η and ρ (no RMATCH-style divide-out / re-multiply). -# -# Differences from the Fortran reference (deliberate): -# - Rational surfaces lie off the equilibrium grid (q = m/n between nodes), so the geometry -# is evaluated with the callable bicubic interpolants at the exact (ψ_s, θ) — using -# DerivOp(0,1) for θ-derivatives — rather than the nodal-grid access mercier_scan! uses. -# - η and ρ are caller-supplied (no Spitzer model, no hardcoded ne/Te/mi defaults). -# -# Reference: Glasser, Greene & Johnson, Phys. Fluids 18, 875 (1975); Glasser PoP 23, 072505 (2016). - -""" - resist_eval(sing, equil, intr; eta, rho, gamma, ising=0) -> InnerLayer.GGJParameters - -Compute the GGJ resistive inner-layer parameters at the rational surface `sing`. -Port of RDCON `resist_eval` (resist.f), but with η and ρ supplied by the caller (no -built-in Spitzer/density model). Returns an `InnerLayer.GGJParameters` carrying the curvature/mass -coefficients E, F, G, H, K, M, the local time scales τ_A (`taua`), τ_R (`taur`), `v1` (V′ normalized -by total volume), and the surface index `ising`. - -# Required keyword arguments - - `eta::Real` — plasma resistivity η at this surface [Ω·m] - - `rho::Real` — mass density ρ at this surface [kg/m³] - - `gamma::Real` — ratio of specific heats (enters G; physical thermodynamic value, e.g. 5/3) - -# Physics -The six flux-surface averages (Fortran `avg(1:6)`) are -`⟨B²/|∇ψ|²⟩, ⟨1/|∇ψ|²⟩, ⟨1/B²⟩, ⟨1/(B²|∇ψ|²)⟩, ⟨B²⟩, ⟨|∇ψ|²/B²⟩`, each weighted by `J/V′` -and integrated over θ ∈ [0,1) — exactly the resist.f integrands. The timescales are -`τ_A = √(ρ·M·μ₀) / |2π n q₁ χ₁ / V′|` and `τ_R = (⟨B²/|∇ψ|²⟩/⟨B²⟩)·μ₀/η`, with η and ρ entering -once, explicitly. -""" -function resist_eval(sing::SingType, equil::Equilibrium.PlasmaEquilibrium, - intr::ForceFreeStatesInternal; eta::Real, rho::Real, gamma::Real, ising::Int=0) - - profiles = equil.profiles - psifac = sing.psifac - nn = sing.n[1] - - # --- 1D surface quantities at the (off-grid) rational ψ --- - hint = Ref(1) - twopif = profiles.F_spline(psifac; hint=hint) - p = profiles.P_spline(psifac; hint=hint) - p1 = profiles.P_deriv(psifac; hint=hint) - v1 = profiles.dVdpsi_spline(psifac; hint=hint) - v2 = profiles.dVdpsi_deriv(psifac; hint=hint) - q = sing.q - q1 = sing.q1 - chi1 = 2π * equil.psio - - # --- θ-integrate the six geometric averages (resist.f) --- - nth = length(equil.rzphi_ys) - ff = zeros(nth, 6) - hint2d = (Ref(1), Ref(1)) - for itheta in 1:nth - theta = equil.rzphi_ys[itheta] - f1 = equil.rzphi_rsquared((psifac, theta); hint=hint2d) - f2 = equil.rzphi_offset((psifac, theta); hint=hint2d) - jac = equil.rzphi_jac((psifac, theta); hint=hint2d) - fy1 = equil.rzphi_rsquared((psifac, theta); deriv=DerivOp(0, 1), hint=hint2d) - fy2 = equil.rzphi_offset((psifac, theta); deriv=DerivOp(0, 1), hint=hint2d) - fy3 = equil.rzphi_nu((psifac, theta); deriv=DerivOp(0, 1), hint=hint2d) - - rfac = sqrt(f1) - eta_geo = 2π * (theta + f2) # geometric poloidal angle (local; not the resistivity η) - r = equil.ro + rfac * cos(eta_geo) - - v21 = fy1 / (2.0 * rfac * jac) - v22 = (1.0 + fy2) * 2π * rfac / jac - v23 = fy3 * r / jac - v33 = 2π * r / jac - bsq = chi1^2 * (v21^2 + v22^2 + (v23 + q * v33)^2) - dpsisq = (2π * r)^2 * (v21^2 + v22^2) - - ff[itheta, 1] = bsq / dpsisq - ff[itheta, 2] = 1.0 / dpsisq - ff[itheta, 3] = 1.0 / bsq - ff[itheta, 4] = 1.0 / (bsq * dpsisq) - ff[itheta, 5] = bsq - ff[itheta, 6] = dpsisq / bsq - @views ff[itheta, :] .*= jac / v1 - end - # Snap the repeated endpoint exactly equal to the start - @views ff[end, :] .= ff[1, :] - itp = cubic_interp(equil.rzphi_ys, Series(ff); bc=PeriodicBC()) - avg = FastInterpolations.integrate(itp) - - # --- curvature terms E, F, H (resist.f) --- - E = p1 * v1 / (q1 * chi1^2)^2 * avg[1] * (twopif * q1 * chi1 / avg[5] - v2) - F = (p1 * v1 / (q1 * chi1^2))^2 * (avg[1] * avg[3] + (twopif / chi1)^2 * (avg[1] * avg[4] - avg[2]^2)) - H = twopif * p1 * v1 / (q1 * chi1^3) * (avg[2] - avg[1] / avg[5]) - - # --- mass factor M and derived G, K (resist.f) --- - M = avg[1] * (avg[6] + (twopif / chi1)^2 * (avg[3] - 1 / avg[5])) - G = avg[5] / (M * gamma * p) - K = (q1 * chi1^2 / (p1 * v1))^2 * avg[5] / (M * avg[1]) - - # --- Alfvén and resistive times from caller-supplied η, ρ (resist.f) --- - # τ_A = √(ρ·M·μ₀) / |2π n q₁ χ₁ / V′| — geometric factor × √ρ. - # τ_R = (⟨B²/|∇ψ|²⟩ / ⟨B²⟩) · μ₀ / η — geometric factor × 1/η. - # η, ρ enter once, explicitly (no RMATCH divide-out / re-multiply). - mu0 = 4π * 1e-7 - geom_taua = sqrt(M * mu0) / abs(2π * nn * q1 * chi1 / v1) - taua = sqrt(rho) * geom_taua - geom_taur = avg[1] / avg[5] * mu0 - taur = geom_taur / eta - - v1norm = equil.params.volume === nothing ? v1 : v1 / equil.params.volume - - return InnerLayer.GGJParameters(; E=E, F=F, G=G, H=H, K=K, M=M, - taua=taua, taur=taur, v1=v1norm, ising=ising) -end diff --git a/src/ForceFreeStates/Surfaces/ResistEval.jl b/src/ForceFreeStates/Surfaces/ResistEval.jl index 09bd6dd80..fb51feb40 100644 --- a/src/ForceFreeStates/Surfaces/ResistEval.jl +++ b/src/ForceFreeStates/Surfaces/ResistEval.jl @@ -7,10 +7,8 @@ # # Port of Fortran RDCON `resist_eval` (geometric part only). # Unlike the Fortran, this routine produces *only* the pure-equilibrium -# quantities; kinetic timescales (τ_A, τ_R) are built on top in the -# downstream `build_ggj_inputs` helper using the same KineticProfiles that -# feed SLAYER, rather than Fortran's hardcoded `ne=1e14, te=3e3` -# parameter defaults. +# quantities; `ggj_parameters` forms τ_A / τ_R on top from caller-supplied +# η and ρ, rather than Fortran's hardcoded `ne=1e14, te=3e3` defaults. # # The 6 theta-integrands match the Fortran layout: # 1: B² / |∇ψ|² @@ -198,6 +196,25 @@ function resist_geometry(equil::Equilibrium.PlasmaEquilibrium, ) end +""" + ggj_parameters(sing, equil; eta, rho, gamma=5/3, ising=0) -> InnerLayer.GGJParameters + +GGJ layer parameters at `sing` from its `restype` geometry and the local η (Ω·m) and ρ (kg/m³). +""" +function ggj_parameters(sing::SingType, equil::Equilibrium.PlasmaEquilibrium; + eta::Real, rho::Real, gamma::Real=5 / 3, ising::Int=0) + rg = sing.restype + rg === nothing && throw(ArgumentError("ggj_parameters: the surface at ψ=$(sing.psifac) has restype = nothing; run resist_eval_all! first")) + equil.params.volume === nothing && throw(ArgumentError("ggj_parameters: equil.params.volume is nothing")) + MU_0 = Utilities.PhysicalConstants.MU_0 + chi1 = 2π * equil.psio + taua = sqrt(rho * rg.M * MU_0) / abs(2π * Int(sing.n[1]) * sing.q1 * chi1 / rg.v1_local) + taur = (rg.avg_bsq_over_dpsisq / rg.avg_bsq) * MU_0 / eta + G = rg.avg_bsq / (rg.M * gamma * rg.p_local) + return InnerLayer.GGJParameters(; E=rg.E, F=rg.F, G=G, H=rg.H, K=rg.K, M=rg.M, + taua=taua, taur=taur, v1=rg.v1_local / equil.params.volume, ising=ising) +end + """ resist_eval_all!(intr::ForceFreeStatesInternal, equil; gamma=5/3) diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 23de76dad..dd655ff3b 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -81,12 +81,15 @@ using .ForceFreeStates: sing_lim!, sing_min!, sing_find!, remove_singular_surfs! using .ForceFreeStates: make_metric, build_matrix_splines, build_kinetic_matrix_splines using .ForceFreeStates: find_kinetic_singular_surfaces! using .ForceFreeStates: eulerlagrange_integration, free_run, normalize_eigenfunctions! -using .ForceFreeStates: galerkin_solve, write_galerkin! +using .ForceFreeStates: galerkin_solve, write_galerkin!, write_match! # Scripting-API surface: the integrator selectors, the published result, the equilibrium # constructor and the forcing description, re-exported so a user needs one `using`. -using .ForceFreeStates: AbstractIntegrator, Forward, Riccati, Galerkin, ResistiveMatch -using .Equilibrium: PlasmaEquilibrium +using .ForceFreeStates: AbstractIntegrator, Forward, Riccati, Galerkin +using .ForceFreeStates: MatchProblem, layer_parameters +using .InnerLayer: GGJModel, SLAYERModel +using .Tearing.Runner: TearingProblem +using .Equilibrium: PlasmaEquilibrium, attach_kinetic_profiles! using .ForcingTerms: RMPField const _DEPRECATED_FFS_KEYS = ("mer_flag", "force_wv_symmetry", "ode_flag", "cyl_flag", "mat_flag", "reform_eq_with_psilim", @@ -219,8 +222,8 @@ function main_from_inputs( equil = Equilibrium.setup_equilibrium(eq_config, additional_input) - kf_ctrl, kinetic_profiles, kf_species = load_kinetic_context(inputs, intr, ctrl, equil) - equil = maybe_reform_equilibrium(equil, eq_config, additional_input, intr, ctrl, kinetic_profiles) + kf_ctrl, kf_species = load_kinetic_context(inputs, intr, ctrl, equil) + equil = maybe_reform_equilibrium(equil, eq_config, additional_input, intr, ctrl) equil_dt = time() - equil_start push!(runtimes, "equilibrium" => equil_dt) @@ -253,7 +256,7 @@ function main_from_inputs( ffs_start = time() locstab, ballooning_boundary = run_local_stability(ctrl, equil) - metric, mats = prepare_force_free_states!(intr, ctrl, equil, kf_ctrl, kinetic_profiles; species=kf_species) + metric, mats = prepare_force_free_states!(intr, ctrl, equil, kf_ctrl; species=kf_species) ffs_result = run_force_free_states(ctrl, equil, mats, intr, metric; runtimes=runtimes) if ctrl.write_outputs_to_HDF5 @@ -309,7 +312,7 @@ function main_from_inputs( end kf_start = time() - run_kinetic_forces(inputs, ffs_result, pe_state, kf_ctrl, kinetic_profiles, kf_species) + run_kinetic_forces(inputs, ffs_result, pe_state, kf_ctrl, kf_species) if "KineticForces" in keys(inputs) push!(runtimes, "kinetic_forces" => time() - kf_start) # Mirrors the write gate inside `run_kinetic_forces` (needs a PE state to contract against). @@ -319,7 +322,7 @@ function main_from_inputs( end ef_start = time() - error_fields = run_error_fields(inputs, ffs_result, pe_state, preloaded_coil_sets, kf_ctrl, kinetic_profiles) + error_fields = run_error_fields(inputs, ffs_result, pe_state, preloaded_coil_sets, kf_ctrl) coil_sensitivities, monte_carlo, locking_risk, efc_couplings = error_fields === nothing ? (nothing, nothing, nothing, nothing) : error_fields if "ErrorFields" in keys(inputs) push!(runtimes, "error_fields" => time() - ef_start) @@ -391,10 +394,10 @@ function resolve_mode_space!(intr::ForceFreeStatesInternal, ctrl::ForceFreeState end """ - load_kinetic_context(inputs, intr, ctrl, equil) -> (kf_ctrl, kinetic_profiles, species) + load_kinetic_context(inputs, intr, ctrl, equil) -> (kf_ctrl, species) -Build the KineticForces control and load the kinetic profiles once for the whole run. -`kinetic_profiles` is `nothing` when no stage asks for them. +Build the KineticForces control, attach the kinetic profiles to `equil` when a stage needs +them, and resolve the NTV species (`nothing` unless the deck requests them). """ function load_kinetic_context( inputs::Dict{String,Any}, @@ -404,8 +407,8 @@ function load_kinetic_context( ) # The profiles are reused by the grid refinement, the stability kinetic callback (via # `calculated_cb`), and the post-PE torque diagnostics block. The `"fixed"` kinetic source - # path in stability does not need kinetic_profiles, but the post-PE block always does, so we - # load whenever a [KineticForces] section is present or the stability path requests the + # path in stability does not need them, but the post-PE block always does, so we attach + # whenever a [KineticForces] section is present or the stability path requests the # calculated source. psio is invariant across grid re-formation. kf_ctrl = haskey(inputs, "KineticForces") ? @@ -413,19 +416,18 @@ function load_kinetic_context( (Symbol(k) => v for (k, v) in inputs["KineticForces"])...) : KineticForces.KineticForcesControl() - kinetic_profiles = nothing species = nothing needs_kinetic_profiles = haskey(inputs, "KineticForces") || (ctrl.kinetic_factor > 0 && ctrl.kinetic_source == "calculated") if needs_kinetic_profiles kinetic_file = joinpath(intr.dir_path, kf_ctrl.kinetic_file) - kinetic_profiles = Equilibrium.load_kinetic_profiles( - kinetic_file; - zi=kf_ctrl.zi, zimp=kf_ctrl.zimp, - mi=kf_ctrl.mi, mimp=kf_ctrl.mimp, - density_factor=kf_ctrl.density_factor, temperature_factor=kf_ctrl.temperature_factor, - ExB_rotation_factor=kf_ctrl.ExB_rotation_factor, toroidal_rotation_factor=kf_ctrl.toroidal_rotation_factor, - chi1=2π * equil.psio) + if equil.kinetic === nothing + Equilibrium.attach_kinetic_profiles!(equil, kinetic_file; + zi=kf_ctrl.zi, zimp=kf_ctrl.zimp, + mi=kf_ctrl.mi, mimp=kf_ctrl.mimp, + density_factor=kf_ctrl.density_factor, temperature_factor=kf_ctrl.temperature_factor, + ExB_rotation_factor=kf_ctrl.ExB_rotation_factor, toroidal_rotation_factor=kf_ctrl.toroidal_rotation_factor) + end if !isempty(kf_ctrl.ion_species) || kf_ctrl.electron speclist = isempty(kf_ctrl.ion_species) ? [KineticForces.IonSpecies(; z=kf_ctrl.zi, m=kf_ctrl.mi, fraction=1.0)] : kf_ctrl.ion_species @@ -437,25 +439,26 @@ function load_kinetic_context( end end - return kf_ctrl, kinetic_profiles, species + return kf_ctrl, species end """ - maybe_reform_equilibrium(equil, eq_config, additional_input, intr, ctrl, kinetic_profiles) -> equil + maybe_reform_equilibrium(equil, eq_config, additional_input, intr, ctrl) -> equil -Two-pass auto grid: measure the pass-1 equilibrium's curvature (profiles, geometry, kinetic -profiles), pin knots on rational surfaces, and re-form on the refined grid from the in-memory -input — no file re-read. Returns `equil` untouched when the configuration wants a single pass. +Two-pass auto grid: measure the pass-1 equilibrium's curvature (profiles, geometry, the +attached `equil.kinetic` profiles), pin knots on rational surfaces, and re-form on the refined +grid from the in-memory input — no file re-read. The kinetic attachment is carried over to the +re-formed equilibrium. Returns `equil` untouched when the configuration wants a single pass. """ function maybe_reform_equilibrium( equil::Equilibrium.PlasmaEquilibrium, eq_config::Equilibrium.EquilibriumConfig, additional_input, intr::ForceFreeStatesInternal, - ctrl::ForceFreeStatesControl, - kinetic_profiles + ctrl::ForceFreeStatesControl ) Equilibrium.wants_two_pass(eq_config) || return equil + kinetic_profiles = equil.kinetic mandatory = ForceFreeStates.rational_psi_nodes(equil; nlow=intr.nlow, nhigh=intr.nhigh) # Smallest |n| in the run sets the widest matching half-stencil dpsi = singfac_min/(n_min·|q′|), @@ -477,6 +480,8 @@ function maybe_reform_equilibrium( nothing # fall back to re-reading the input file end equil = Equilibrium.setup_equilibrium(eq_config, rerun_input; override_psi_nodes=psi_nodes) + # Carry the attachment across re-formation; psio (hence chi1) is invariant, so no re-load. + equil.kinetic = kinetic_profiles implied = Equilibrium.implied_knot_count(equil; tau=eq_config.psi_accuracy, kin=kinetic_profiles) if implied > 1.5 * (length(psi_nodes) - 1) @warn "Two-pass psi grid: refined equilibrium implies $implied knots vs $(length(psi_nodes) - 1) used — " * @@ -542,18 +547,18 @@ function run_local_stability(ctrl::ForceFreeStatesControl, equil::Equilibrium.Pl end """ - prepare_force_free_states!(intr, ctrl, equil, kf_ctrl, kinetic_profiles) -> (metric, mats) + prepare_force_free_states!(intr, ctrl, equil, kf_ctrl; species=nothing) -> (metric, mats) Set up the force-free-states solve on `intr`: integration limits, the surviving singular surfaces and their GGJ coefficients, the poloidal mode range, and the metric plus -Euler-Lagrange (and, when requested, kinetic) matrices. +Euler-Lagrange (and, when requested, kinetic) matrices. The `"calculated"` kinetic source +reads its profiles off `equil.kinetic`. """ function prepare_force_free_states!( intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, - kf_ctrl::KineticForces.KineticForcesControl, - kinetic_profiles; + kf_ctrl::KineticForces.KineticForcesControl; species=nothing ) # Determine psilim and qlim (where we will integrate to) @@ -573,10 +578,8 @@ function prepare_force_free_states!( sing_min!(intr, ctrl, equil) end - # Populate Glasser-Greene-Johnson geometric coefficients (E, F, G, H, - # K, M) for each surviving singular surface. Needed by the Julia GGJ - # inner-layer analysis; kinetic timescales (τ_A, τ_R) are layered on - # top by `build_ggj_inputs` using the same kinetic profiles as SLAYER. + # Populate Glasser-Greene-Johnson geometric coefficients (E, F, G, H, K, M) for each + # surviving singular surface; `ggj_parameters` layers τ_A / τ_R on top. if intr.msing > 0 ForceFreeStates.resist_eval_all!(intr, equil) end @@ -631,7 +634,7 @@ function prepare_force_free_states!( calculated_cb = (c, e, i, m, f) -> KineticForces.compute_calculated_kinetic_matrices( c, e, i, m, f; - kf_ctrl=kf_ctrl, kinetic_profiles=kinetic_profiles, species=species) + kf_ctrl=kf_ctrl, kinetic_profiles=equil.kinetic, species=species) mats = build_kinetic_matrix_splines(ctrl, equil, mats, intr, metric; calculated_source=calculated_cb) @@ -723,31 +726,58 @@ function run_force_free_states( end # Publish the solve: from here on the downstream stages read the result, never `intr`. - return build_result(Symbol(ctrl.integrator), ctrl, equil, intr, metric, mats, odet, free_energies, gal_data, gal_dp) + result = build_result(Symbol(ctrl.integrator), ctrl, equil, intr, metric, mats, odet, free_energies, gal_data, gal_dp) + + # Deck-driven matching, through the same MatchProblem as the API. + if ctrl.integrator == "galerkin" && ctrl.gal_match_flag + ctrl.gal_rpec_flag || error("gal_match_flag=true requires gal_rpec_flag=true") + ctrl.gal_ideal_flag || !(isempty(ctrl.gal_eta) || isempty(ctrl.gal_rho) || isempty(ctrl.gal_rotation)) || + error("gal_match_flag needs gal_eta, gal_rho and gal_rotation (one value per matched surface, core to edge)") + ctrl.gal_inner_solver in ("ray", "galerkin") || + error("gal_inner_solver = \"$(ctrl.gal_inner_solver)\" (expected \"ray\" or \"galerkin\")") + ctrl.verbose && @info( + ctrl.gal_ideal_flag ? + "RPEC matching: IDEAL solution (inner layer skipped, bare coil columns)" : + "RPEC matching: inner-layer Δ(Q) + outer↔inner solve for the coil-driven ξ" + ) + # The gal_inner_* knobs override the "galerkin" backend's defaults; unset keys keep them. + model = if ctrl.gal_inner_solver == "galerkin" + knobs = (xfac=ctrl.gal_inner_xfac, nx=ctrl.gal_inner_nx, nq=ctrl.gal_inner_nq, + cutoff=ctrl.gal_inner_cutoff, kmax=ctrl.gal_inner_kmax) + InnerLayer.GGJModel(; solver=:galerkin, (k => v for (k, v) in pairs(knobs) if v !== nothing)...) + else + InnerLayer.GGJModel(; solver=:ray) + end + prob = ForceFreeStates.MatchProblem(result; eta=ctrl.gal_eta, rho=ctrl.gal_rho, rotation=ctrl.gal_rotation, + gamma=ctrl.gal_gamma, ideal=ctrl.gal_ideal_flag) + result = solve(prob, model) + ctrl.gal_ideal_flag || (ctrl.verbose && @info "RPEC matching: linear-solve residual = $(result.match.residual)") + end + + return result end """ - EulerLagrangeProblem(equil; nn, wall=Vacuum.WallShapeSettings(), match=nothing, + EulerLagrangeProblem(equil; nn, wall=Vacuum.WallShapeSettings(), dir_path=".", debug=DebugSettings(), kwargs...) The perturbed-plasma Euler-Lagrange problem posed on an equilibrium: the extremization of the perturbed potential energy whose solutions are the force-free (and, via the TOML path, kinetic) perturbed states. This is the WHAT of a stability solve; the integrator passed to [`solve`](@ref) is the HOW. A `PlasmaEquilibrium` hosts many possible problems — this type -names this one, so `solve` stays unambiguous as other problem classes appear. +names this one, so `solve` stays unambiguous as other problem classes appear. Inner-layer +matching is a separate problem posed on the finished solve: see `MatchProblem`. -`nn` is the toroidal mode number or range. `wall` is the vacuum wall shape, `match` an -optional [`ResistiveMatch`](@ref) closing the basis with an inner-layer solution instead of -the ideal jump, `dir_path` the working directory outputs are written to, and `debug` the -diagnostic dump settings of the DEBUG deck section. Any remaining keyword is a -`ForceFreeStatesControl` field, so the TOML keys and the problem keywords are the same -knobs. `nn_low`/`nn_high` are rejected — they come from `nn`. +`nn` is the toroidal mode number or range. `wall` is the vacuum wall shape, `dir_path` the +working directory outputs are written to, and `debug` the diagnostic dump settings of the +DEBUG deck section. Any remaining keyword is a `ForceFreeStatesControl` field, so the TOML +keys and the problem keywords are the same knobs. `nn_low`/`nn_high` are rejected — they +come from `nn`. ## Fields - `equil::Equilibrium.PlasmaEquilibrium` - The equilibrium the problem is posed on. - `wall::Vacuum.WallShapeSettings` - Vacuum wall shape for the free-boundary energies. - - `match::Union{Nothing,ForceFreeStates.ResistiveMatch}` - Optional inner-layer closure. - `dir_path::String` - Working directory for outputs. - `debug::DebugSettings` - Diagnostic dump settings. - `ctrl_kwargs::Dict{Symbol,Any}` - `ForceFreeStatesControl` keywords, `nn` already folded in. @@ -755,7 +785,6 @@ knobs. `nn_low`/`nn_high` are rejected — they come from `nn`. struct EulerLagrangeProblem equil::Equilibrium.PlasmaEquilibrium wall::Vacuum.WallShapeSettings - match::Union{Nothing,ForceFreeStates.ResistiveMatch} dir_path::String debug::DebugSettings ctrl_kwargs::Dict{Symbol,Any} @@ -765,7 +794,6 @@ function EulerLagrangeProblem( equil::Equilibrium.PlasmaEquilibrium; nn::Union{Int,AbstractUnitRange{Int}}, wall::Vacuum.WallShapeSettings=Vacuum.WallShapeSettings(), - match::Union{Nothing,ForceFreeStates.ResistiveMatch}=nothing, dir_path::AbstractString=".", debug::DebugSettings=DebugSettings(), kwargs... @@ -775,7 +803,7 @@ function EulerLagrangeProblem( error("the toroidal mode range comes from the `nn` keyword; drop nn_low/nn_high") ctrl_kwargs[:nn_low] = first(nn) ctrl_kwargs[:nn_high] = last(nn) - return EulerLagrangeProblem(equil, wall, match, String(dir_path), debug, ctrl_kwargs) + return EulerLagrangeProblem(equil, wall, String(dir_path), debug, ctrl_kwargs) end """ @@ -788,7 +816,7 @@ Solve the perturbed-plasma [`EulerLagrangeProblem`](@ref) with the formalism `al a `gpec.toml` run of `main` does and produces the same result object. The second form is sugar building the problem from an equilibrium and the problem keywords in one call. -Knobs owned by `alg` or `match` are rejected as `ForceFreeStatesControl` keywords. Kinetic +Knobs owned by `alg` are rejected as `ForceFreeStatesControl` keywords. Kinetic runs are TOML-driven this cycle: `kinetic_factor > 0` needs the `[KineticForces]` profiles and errors here. @@ -805,7 +833,6 @@ function solve(prob::EulerLagrangeProblem, alg::ForceFreeStates.AbstractIntegrat equil = prob.equil ctrl_kwargs = copy(prob.ctrl_kwargs) ForceFreeStates._apply_alg!(ctrl_kwargs, alg) - ForceFreeStates._apply_match!(ctrl_kwargs, prob.match, alg) ctrl = ForceFreeStatesControl(; ctrl_kwargs...) ctrl.kinetic_factor > 0 && @@ -825,11 +852,11 @@ function solve(prob::EulerLagrangeProblem, alg::ForceFreeStates.AbstractIntegrat @warn "Two-pass auto grid needs the equilibrium's raw ingest, which analytic and IMAS equilibria do not carry; " * "solving on the single-pass grid. Set mpsi explicitly to choose the grid." else - equil = maybe_reform_equilibrium(equil, equil.config, nothing, intr, ctrl, nothing) + equil = maybe_reform_equilibrium(equil, equil.config, nothing, intr, ctrl) end locstab, ballooning_boundary = run_local_stability(ctrl, equil) - metric, mats = prepare_force_free_states!(intr, ctrl, equil, kf_ctrl, nothing) + metric, mats = prepare_force_free_states!(intr, ctrl, equil, kf_ctrl) result = run_force_free_states(ctrl, equil, mats, intr, metric) if ctrl.write_outputs_to_HDF5 @@ -954,18 +981,18 @@ function perturbed_equilibrium( end """ - run_kinetic_forces(inputs, result, pe_state, kf_ctrl, kinetic_profiles) + run_kinetic_forces(inputs, result, pe_state, kf_ctrl, species=nothing) Compute and write the neoclassical toroidal viscosity torque diagnostics when the deck carries a -`[KineticForces]` section. No-op when the perturbed-equilibrium state the operators contract -against is missing. +`[KineticForces]` section, reading the kinetic profiles off `result.equil.kinetic` (per-species +profiles ride the `species` resolution instead). No-op when the perturbed-equilibrium state the +operators contract against is missing. """ function run_kinetic_forces( inputs::Dict{String,Any}, result::ForceFreeStatesResult, pe_state, kf_ctrl::KineticForces.KineticForcesControl, - kinetic_profiles, species=nothing ) # ---------------------------------------------------------------- @@ -981,13 +1008,12 @@ function run_kinetic_forces( if pe_state === nothing @info "Skipping NTV torque diagnostics: no perturbed-equilibrium data (e.g. kinetic_source=\"calculated\")." else - # kf_ctrl and kinetic_profiles were loaded once before the equilibrium was re-formed. kf_intr = KineticForces.KineticForcesInternal(result.equil; verbose=kf_ctrl.verbose) KineticForces.set_perturbation_data!(kf_intr, pe_state, result, result.equil, result.metric) if species === nothing kf_state = KineticForces.KineticForcesState() - KineticForces.compute_torque_all_methods!(kf_state, kf_intr, kf_ctrl, result.equil, kinetic_profiles) + KineticForces.compute_torque_all_methods!(kf_state, kf_intr, kf_ctrl, result.equil, result.equil.kinetic) if kf_ctrl.write_outputs_to_HDF5 h5open(joinpath(result.dir_path, kf_ctrl.HDF5_filename), "cw") do h5file KineticForces.write_to_hdf5!(h5file, kf_state; @@ -1043,7 +1069,7 @@ function forcing_terms_control(inputs::Dict{String,Any}) end """ - run_error_fields(inputs, result, pe_state, preloaded_coil_sets, kf_ctrl, kinetic_profiles) -> (sensitivities, monte_carlo, risk, couplings) or nothing + run_error_fields(inputs, result, pe_state, preloaded_coil_sets, kf_ctrl) -> (sensitivities, monte_carlo, risk, couplings) or nothing Linearize every coil set's resonant drive with respect to its rigid shifts and tilts and write `ErrorFields/CoilSensitivities/` when the deck carries an `[ErrorFields]` section. When the @@ -1052,7 +1078,7 @@ Carlo on the full-window dominant mode with the `[ErrorFields.MonteCarlo]` setti `ErrorFields/MonteCarlo/`; with an `[ErrorFields.scenario]` table as well, evaluate the locking risk (and the tolerance scan when `[ErrorFields.Risk]` names `scan_scales`) and write `ErrorFields/Risk/`; with `[ErrorFields.NTV]` naming correction arrays, evaluate their overlap and NTV -torque couplings per kilo-ampere-turn (needs the kinetic context) and write `ErrorFields/NTV/`. +torque couplings per kilo-ampere-turn (needs `result.equil.kinetic`) and write `ErrorFields/NTV/`. Stages not requested return `nothing`. Needs the perturbed-equilibrium state's singular-coupling matrix and coil-format forcing, and errors otherwise: a deck asking for error-field sensitivities without them is a misconfiguration, not a @@ -1063,8 +1089,7 @@ function run_error_fields( result::ForceFreeStatesResult, pe_state, preloaded_coil_sets::Union{Nothing,Vector{ForcingTerms.CoilSet}}, - kf_ctrl::KineticForces.KineticForcesControl=KineticForces.KineticForcesControl(), - kinetic_profiles=nothing + kf_ctrl::KineticForces.KineticForcesControl=KineticForces.KineticForcesControl() ) ("ErrorFields" in keys(inputs)) || return nothing @@ -1144,7 +1169,7 @@ function run_error_fields( missing = setdiff(ntv_ctrl.efc_coils, [cs.name for cs in efc_sets]) isempty(missing) || error("[ErrorFields.NTV] efc_coils not among the run's coil sets: $(join(missing, ", "))") ntv_start = time() - couplings = efc_couplings(result, efc_sets, rc, dom, cfg, kf_ctrl, kinetic_profiles; method=ntv_ctrl.method, verbose=ntv_ctrl.verbose, + couplings = efc_couplings(result, efc_sets, rc, dom, cfg, kf_ctrl, result.equil.kinetic; method=ntv_ctrl.method, verbose=ntv_ctrl.verbose, rotation_scan=ntv_ctrl.rotation_scan, scan_points=ntv_ctrl.rotation_scan_points, scan_max_points=ntv_ctrl.rotation_scan_max_points, scan_tolerance=ntv_ctrl.rotation_scan_tolerance, span_factor=ntv_ctrl.rotation_span_factor, offset_factor=ntv_ctrl.rotation_offset_factor) @info "NTV couplings of $(length(couplings)) correction arrays in $(@sprintf("%.1f", time() - ntv_start)) s" @@ -1336,8 +1361,8 @@ function run_slayer_stage(result::ForceFreeStatesResult, inputs::Dict{String,Any slayer_ctrl.enabled || return nothing @info "\n SLAYER\n$_SECTION" slayer_start = time() - slayer_result = Runner.run_slayer(result, slayer_ctrl; - dir_path=result.dir_path) + model = Runner._build_inner_model(slayer_ctrl.inner_model) + slayer_result = solve(Runner.TearingProblem(result, slayer_ctrl), model) slayer_dt = time() - slayer_start runtimes === nothing || push!(runtimes, "tearing" => slayer_dt) @info "SLAYER completed in $(@sprintf("%.3f", slayer_dt)) s" @@ -1565,7 +1590,7 @@ function write_outputs_to_HDF5( # (populated by ForceFreeStates.resist_eval_all! after sing_find!). # Both kinetic-free (E, F, G, H, K, M) and geometry-only # (avg_bsq_over_dpsisq, avg_bsq) quantities are written so - # downstream consumers (Tearing.InnerLayer.GGJ.build_ggj_inputs) + # downstream consumers (ForceFreeStates.ggj_parameters) # can reconstruct τ_A / τ_R from any kinetic-profile source. if all(s -> s.restype !== nothing, result.surfaces) out_h5["SingularSurfaces/E"] = [s.restype.E for s in result.surfaces] @@ -1611,6 +1636,7 @@ function write_outputs_to_HDF5( dp.B === nothing || (out_h5["SingularSurfaces/pest3_B"] = dp.B) dp.Gamma === nothing || (out_h5["SingularSurfaces/pest3_Gamma"] = dp.Gamma) end + result.match === nothing || write_match!(out_h5, result.match) # Write kinetic singular surface data (det(F̄) near-zeros) and the cond(F̄) scan # used to find them. Populated only when kinetic crossings were searched for. @@ -1776,7 +1802,8 @@ end export main, write_imas export solve, perturbed_equilibrium -export PlasmaEquilibrium, EulerLagrangeProblem, Forward, Riccati, Galerkin, ResistiveMatch, ForceFreeStatesResult, RMPField +export PlasmaEquilibrium, attach_kinetic_profiles!, EulerLagrangeProblem, Forward, Riccati, Galerkin, ForceFreeStatesResult, RMPField +export MatchProblem, TearingProblem, GGJModel, SLAYERModel, layer_parameters include("Precompile.jl") diff --git a/src/InnerLayer/GGJ/GGJ.jl b/src/InnerLayer/GGJ/GGJ.jl index 80a4bd617..cdd6e37d6 100644 --- a/src/InnerLayer/GGJ/GGJ.jl +++ b/src/InnerLayer/GGJ/GGJ.jl @@ -43,16 +43,36 @@ import ..InnerLayerModel, ..InnerLayerResponse, ..InnerLayerParameters import ..solve_inner, ..solve_inner_profile """ - GGJModel{S} <: InnerLayerModel + GGJModel{S,O} <: InnerLayerModel + GGJModel(; solver=:ray, options...) Glasser–Greene–Johnson resistive inner-layer model. `S` selects the solver backend: `:ray` (default; robust at large |Q| on/near the imaginary axis), `:galerkin` (real-axis Hermite FEM; degrades for |Q| ≳ 1), or `:shooting` -(|Q| ≪ 1 only). The backends take different numerical-knob keywords. +(|Q| ≪ 1 only). Backend keywords given at construction, e.g. +`GGJModel(; solver=:galerkin, nx=1280)`, apply to every solve. + +## Fields + + - `options::NamedTuple` - Backend keywords forwarded to every solve; call-site keywords win. """ -struct GGJModel{S} <: InnerLayerModel end +struct GGJModel{S,O<:NamedTuple} <: InnerLayerModel + options::O +end + +GGJModel(; solver::Symbol=:ray, options...) = GGJModel{solver,typeof(values(options))}(values(options)) +GGJModel{S}() where {S} = GGJModel{S,@NamedTuple{}}((;)) -GGJModel(; solver::Symbol=:ray) = GGJModel{solver}() +# The backends implement the option-free model; a model carrying options forwards them. +const _BareGGJ{S} = GGJModel{S,@NamedTuple{}} +function solve_inner(m::GGJModel{S}, params, γ::Number; kwargs...) where {S} + m isa _BareGGJ && throw(MethodError(solve_inner, (m, params, γ))) + return solve_inner(GGJModel{S}(), params, γ; m.options..., kwargs...) +end +function solve_inner_profile(m::GGJModel{S}, params, γ::Number; kwargs...) where {S} + m isa _BareGGJ && throw(MethodError(solve_inner_profile, (m, params, γ))) + return solve_inner_profile(GGJModel{S}(), params, γ; m.options..., kwargs...) +end include("GGJParameters.jl") include("InnerAsymptotics.jl") diff --git a/src/InnerLayer/GGJ/Galerkin.jl b/src/InnerLayer/GGJ/Galerkin.jl index 0db8811d0..a6c287438 100644 --- a/src/InnerLayer/GGJ/Galerkin.jl +++ b/src/InnerLayer/GGJ/Galerkin.jl @@ -873,7 +873,7 @@ Hermite-FEM implementation of the [`solve_inner_profile`](@ref) interface: real-axis solve, so `Δ` and the profiles come from the same solution. Same numerics/kwargs as `solve_inner(GGJModel(; solver=:galerkin), ...)`. """ -function solve_inner_profile(::GGJModel{:galerkin}, params::GGJParameters, γ::Number; kwargs...) +function solve_inner_profile(::_BareGGJ{:galerkin}, params::GGJParameters, γ::Number; kwargs...) Δ, _, prof, _ = solve_inner_profile(params, γ; kwargs...) return (; Δ=Δ, x=prof.x, Ψ=prof.Ψ, Ξ=prof.Ξ, _profile_conversions(params)...) end @@ -894,7 +894,7 @@ Returns the parity-projected matching data (GWP2016 Eqs. 34–35) with the respectively — no parity swap (see the boundary-condition block above for the parity derivation). """ -function solve_inner(::GGJModel{:galerkin}, params::GGJParameters, γ::Number; +function solve_inner(::_BareGGJ{:galerkin}, params::GGJParameters, γ::Number; kmax::Int=8, nx::Int=512, nq::Int=4, pfac::Float64=1.0, cutoff::Int=5, xfac::Float64=1.0, tol_res::Float64=1e-5) Q = inner_Q(params, γ) diff --git a/src/InnerLayer/GGJ/Ray.jl b/src/InnerLayer/GGJ/Ray.jl index cb188f4aa..8440432f6 100644 --- a/src/InnerLayer/GGJ/Ray.jl +++ b/src/InnerLayer/GGJ/Ray.jl @@ -933,7 +933,7 @@ channels are swapped into the named fields of [`InnerLayerResponse`](@ref). Preferred for |Q| ≳ 1 and near the imaginary axis; use `solve_ray` directly when the full [`RaySolveResult`](@ref) is wanted. """ -function solve_inner(::GGJModel{:ray}, params::GGJParameters, γ::Number; kwargs...) +function solve_inner(::_BareGGJ{:ray}, params::GGJParameters, γ::Number; kwargs...) res = solve_ray(params, inner_Q(params, γ); kwargs...) return InnerLayerResponse(res.Δ[2], res.Δ[1]) end @@ -954,7 +954,7 @@ is returned as `certΔ` and warns above `certify_rtol`. `npc` sets the profile points per mesh cell; extra keywords forward to both [`solve_ray`](@ref) calls (θ is fixed by the method — do not pass it). """ -function solve_inner_profile(::GGJModel{:ray}, params::GGJParameters, γ::Number; +function solve_inner_profile(::_BareGGJ{:ray}, params::GGJParameters, γ::Number; npc::Int=8, certify_rtol::Float64=1e-3, kwargs...) Q = inner_Q(params, γ) rr = solve_ray(params, Q; kwargs...) # certified Δ (optimal θ) diff --git a/src/InnerLayer/GGJ/Shooting.jl b/src/InnerLayer/GGJ/Shooting.jl index 50e7c403a..b83b2c369 100644 --- a/src/InnerLayer/GGJ/Shooting.jl +++ b/src/InnerLayer/GGJ/Shooting.jl @@ -357,7 +357,7 @@ Tolerances `reltol`/`abstol` are the integrator tolerances; `rtol_origin` controls the truncation error of the origin Frobenius series and the choice of `tmin`. """ -function solve_inner(::GGJModel{:shooting}, params::GGJParameters, γ::Number; +function solve_inner(::_BareGGJ{:shooting}, params::GGJParameters, γ::Number; reltol::Float64=1e-6, abstol::Float64=1e-6, rtol_origin::Float64=1e-6, nps::Int=8, fmax::Float64=1.0, solver=Tsit5()) @@ -381,5 +381,5 @@ function solve_inner(::GGJModel{:shooting}, params::GGJParameters, γ::Number; return InnerLayerResponse(Δ_rescaled[2], Δ_rescaled[1]) end -solve_inner(::GGJModel{:shooting}, params::GGJParameters, γ::Real; kwargs...) = +solve_inner(::_BareGGJ{:shooting}, params::GGJParameters, γ::Real; kwargs...) = solve_inner(GGJModel{:shooting}(), params, ComplexF64(γ); kwargs...) diff --git a/src/InnerLayer/SLAYER/LayerInputs.jl b/src/InnerLayer/SLAYER/LayerInputs.jl index 252e81ae9..b2fb418c0 100644 --- a/src/InnerLayer/SLAYER/LayerInputs.jl +++ b/src/InnerLayer/SLAYER/LayerInputs.jl @@ -292,7 +292,7 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; prof = profiles(psi) # Take ω_*e, ω_*i from the spline derivatives, or from `profiles` when the caller - # supplies them directly. `run_slayer` supplies zeros, so the latter is a library path. + # supplies them directly. The tearing solve supplies zeros, so the latter is a library path. ω_e_use, ω_i_use = compute_omega_star ? _omega_star_at(psi, n_res) : (prof.omega_e, prof.omega_i) # Pull geometric trapped-fraction inputs from ResistGeometry when diff --git a/src/Tearing/LayerInputs.jl b/src/Tearing/LayerInputs.jl deleted file mode 100644 index 3c1e2cc35..000000000 --- a/src/Tearing/LayerInputs.jl +++ /dev/null @@ -1,135 +0,0 @@ -# LayerInputs.jl (GGJ) -# -# Build per-surface `GGJParameters` from a solved `PlasmaEquilibrium`, the -# `SingType` rational-surface list (each carrying a populated -# `restype::ResistGeometry` from `ForceFreeStates.resist_eval_all!`), and a -# `KineticProfiles` object — the same three ingredients `build_slayer_inputs` -# consumes. Produces the (E, F, G, H, K, τ_A, τ_R) tuple that GGJ's -# `solve_inner` needs, with τ_A / τ_R built from kinetic profiles using the -# same Spitzer resistivity and mass-density formulas SLAYER uses. -# -# Deliberately does *not* use any hardcoded `ne`/`te` defaults. The kinetic -# content enters through `profiles` alone; this keeps GGJ and SLAYER using -# bit-identical plasma inputs when both are driven by the same -# `KineticProfiles`. - -using ..Utilities: KineticProfiles -using ..Utilities.PhysicalConstants: MU_0, M_E, M_P, E_CHG, EPS_0 -using ..Utilities.NeoclassicalResistivity -using ..Utilities.NeoclassicalResistivity: NeoResistivityModel, SpitzerModel, - SauterNeoModel, RedlNeoModel, - coulomb_log_e, eta_spitzer, nu_star_e, eta_neoclassical -using ..ForceFreeStates: ResistGeometry -using ..InnerLayer.GGJ: GGJParameters - -""" - build_ggj_inputs(equil, sings, profiles; mu_i=2.0, zeff=1.0, - v1_scale=1.0, - resistivity_model::NeoResistivityModel=SpitzerModel(), - lnLambda_form::Symbol=:nrl) -> Vector{GGJParameters} - -Construct a `GGJParameters` for each rational surface in `sings`. Each -surface's geometric coefficients (E, F, G, H, K, M) come from the -`sing.restype::ResistGeometry` populated by `resist_eval_all!`. Kinetic -timescales are derived from the `KineticProfiles` at `sing.psifac`: - -``` -ρ(ψ) = μ_i · m_p · n_e(ψ) -η(ψ) = eta_neoclassical(model, n_e, T_e, Z_eff, f_t, ν*_e) [Ω·m] -τ_A = √(ρ · M · μ_0) / |2π · n · q' · χ₁ / V'| [Alfvén time] -τ_R = (⟨B²/|∇ψ|²⟩ / ⟨B²⟩) · μ_0 / η [resistive diffusion] -``` - -The mode number `n` is taken from `sings[k].n[1]` (first resonant mode at -the surface). `χ₁ = 2π · psio`. The `v1_scale` kwarg is an optional -multiplicative factor on `V'` in the τ_A denominator (a `v1 / volume` -normalization option); default `1.0` means use the raw `V'`. - -# Resistivity model - -`resistivity_model` selects the η closure: - - - `SpitzerModel()` (default) — Sauter 1999 Eq. 18a (Zeff-aware Spitzer), - with the NRL Coulomb log. - - `SauterNeoModel()` — multiplies by Sauter 1999 F_33 using f_t and ν*_e - from the surface's `ResistGeometry`. Produces the physically-correct - trapped-particle-corrected η for H-mode tearing stability. - - `RedlNeoModel()` — Redl 2021 F_33 (improved high-ν* fit). - -`lnLambda_form` selects `:nrl` (default), `:sauter`, or `:wesson`. - -Throws if any surface's `restype` is still `nothing` — call -`ForceFreeStates.resist_eval_all!(intr, equil)` first. -""" -function build_ggj_inputs(equil, sings, profiles::KineticProfiles; - mu_i::Real=2.0, zeff::Real=1.0, - v1_scale::Real=1.0, - resistivity_model::NeoResistivityModel=SpitzerModel(), - lnLambda_form::Symbol=:nrl) - psio = equil.psio - chi1 = 2π * psio - - out = Vector{GGJParameters}(undef, length(sings)) - for (k, sing) in enumerate(sings) - rg = sing.restype - rg === nothing && - throw( - ArgumentError( - "build_ggj_inputs: surface $k has " * - "restype = nothing. Call " * - "ForceFreeStates.resist_eval_all!(intr, equil) " * - "after sing_find! to populate it." - ) - ) - rg isa ResistGeometry || - throw(ArgumentError("build_ggj_inputs: surface $k has " * - "restype of unexpected type $(typeof(rg)).")) - - # Kinetic profiles at this surface - prof = profiles(sing.psifac) - n_e = prof.n_e # [m⁻³] - t_e = prof.T_e # [eV] - - # Shared Coulomb log and resistivity closure (identical to SLAYER - # when the same resistivity_model is selected). - lnLamb = coulomb_log_e(n_e, t_e; form=lnLambda_form) - if resistivity_model isa SpitzerModel - eta_use = eta_spitzer(n_e, t_e, zeff; lnLamb=lnLamb) - else - nuestar = nu_star_e(n_e, t_e, rg.R_major, rg.eps_local, - sing.q, zeff; lnLamb=lnLamb) - eta_use = eta_neoclassical(resistivity_model, n_e, t_e, zeff, - rg.f_trap, nuestar; lnLamb=lnLamb) - end - rho = mu_i * M_P * n_e - - # Alfvén time at the rational surface - n_tor = Int(sing.n[1]) - v1 = rg.v1_local * v1_scale - taua = sqrt(rho * rg.M * MU_0) / - abs(2π * n_tor * sing.q1 * chi1 / v1) - - # Resistive diffusion time - taur = (rg.avg_bsq_over_dpsisq / rg.avg_bsq) * MU_0 / eta_use - - # dV/dψ normalized by total plasma volume (`v1 = v1/volume`). This is - # the `v1` consumed by `rescale_delta` as v1^(2p1); NOT the raw V' - # used in τ_A above. - equil.params.volume === nothing && - throw( - ArgumentError( - "build_ggj_inputs: equil.params.volume " * - "is nothing. Ensure the equilibrium " * - "solver populated the total plasma " * - "volume before building GGJ inputs." - ) - ) - v1_norm = rg.v1_local / equil.params.volume - - out[k] = GGJParameters(; - E=rg.E, F=rg.F, G=rg.G, H=rg.H, K=rg.K, M=rg.M, - taua=taua, taur=taur, v1=v1_norm, ising=k - ) - end - return out -end diff --git a/src/Tearing/Runner/Control.jl b/src/Tearing/Runner/Control.jl index 1a12ff7e6..5fa570ad7 100644 --- a/src/Tearing/Runner/Control.jl +++ b/src/Tearing/Runner/Control.jl @@ -16,8 +16,8 @@ constructor. # Core toggles - `enabled` -- run the analysis at all (default `false`) - - `inner_model` -- `:slayer_fitzpatrick` (default), `:ggj_shooting`, or - `:ggj_galerkin` + - `inner_model` -- `:slayer_fitzpatrick` (default), `:ggj_ray`, + `:ggj_shooting`, or `:ggj_galerkin` - `scan_mode` -- `:amr` (default) or `:brute_force` - `coupling_mode` -- `:uncoupled` (default, per-surface) or `:coupled` (multi-surface determinant) @@ -157,7 +157,7 @@ there is one consistent interface for resistive and kinetic profiles. store_scan::Bool = false end -const _VALID_INNER_MODELS = (:slayer_fitzpatrick, :ggj_shooting, :ggj_galerkin) +const _VALID_INNER_MODELS = (:slayer_fitzpatrick, :ggj_ray, :ggj_shooting, :ggj_galerkin) const _VALID_SCAN_MODES = (:amr, :brute_force) const _VALID_COUPLING_MODES = (:uncoupled, :coupled) const _VALID_DC_TYPES = (:none, :lar, :rfitzp, :toroidal) diff --git a/src/Tearing/Runner/Result.jl b/src/Tearing/Runner/Result.jl index 60a728869..6cd862815 100644 --- a/src/Tearing/Runner/Result.jl +++ b/src/Tearing/Runner/Result.jl @@ -7,7 +7,7 @@ """ SLAYERResult -Output of `run_slayer`. Carries both summary eigenvalues (ω_Hz, γ_Hz) and +Output of the tearing solve. Carries both summary eigenvalues (ω_Hz, γ_Hz) and full diagnostic detail (valid roots, poles, filtered roots, contours) for downstream inspection and HDF5 output. diff --git a/src/Tearing/Runner/Runner.jl b/src/Tearing/Runner/Runner.jl index 919088065..dc172869b 100644 --- a/src/Tearing/Runner/Runner.jl +++ b/src/Tearing/Runner/Runner.jl @@ -27,16 +27,17 @@ using LinearAlgebra using Statistics: mean, median using HDF5 +import CommonSolve using FastInterpolations: cubic_interp using ..Utilities using ..Utilities: KineticProfiles using ...Equilibrium: read_kinetic_file, KineticProfileData +using ...ForceFreeStates: ForceFreeStatesResult, layer_parameters, ggj_parameters using ..InnerLayer using ..InnerLayer: InnerLayerParameters, InnerLayerResponse, solve_inner, SLAYERModel, SLAYERParameters, build_slayer_inputs, GGJModel, GGJParameters, LayerWidths, slayer_layer_thickness -import ..build_ggj_inputs # defined at the Tearing level (needs ForceFreeStates) using ..Dispersion using ..Dispersion: SurfaceCoupling, surface_coupling, MultiSurfaceCoupling, multi_surface_coupling, @@ -48,11 +49,12 @@ using ..Dispersion: SurfaceCoupling, surface_coupling, include("Control.jl") include("Result.jl") include("run_slayer.jl") +include("TearingProblem.jl") include("HDF5Output.jl") export SLAYERControl, slayer_control_from_toml, validate export SLAYERResult, empty_slayer_result -export run_slayer, run_slayer_from_inputs, ggj_inner_deltas +export TearingProblem, run_slayer_from_inputs, ggj_inner_deltas export write_slayer_hdf5! end # module Runner diff --git a/src/Tearing/Runner/TearingProblem.jl b/src/Tearing/Runner/TearingProblem.jl new file mode 100644 index 000000000..4723defcb --- /dev/null +++ b/src/Tearing/Runner/TearingProblem.jl @@ -0,0 +1,136 @@ +# TearingProblem.jl +# +# Tearing growth rates from a finished solve: root-find Δ_inner(Q) = Δ′_outer per surface. + +""" + TearingProblem(ffs; kwargs...) + +Tearing-mode growth rates from a finished force-free-states solve: hold the outer Δ′ fixed +and find the growth rate at which the inner layer matches it. Solve with +`solve(prob, SLAYERModel())` or `solve(prob, GGJModel())`. Keywords are +[`SLAYERControl`](@ref) fields, the `[SLAYER]` deck knobs. Kinetic profiles come from +`profile_file` if set, else from `ffs.equil.kinetic`. Without a full Δ′ matrix the per-surface +Δ′ stubs are used, with a warning. + +Fields: `ffs`, `control::SLAYERControl`. +""" +# `ffs` is duck-typed (equil, surfaces, delta_prime, dir_path) so tests can use stand-ins. +struct TearingProblem{R} + ffs::R + control::SLAYERControl +end + +function TearingProblem(ffs::ForceFreeStatesResult; kwargs...) + haskey(kwargs, :inner_model) && error("TearingProblem: the inner-layer model is the solve argument; drop inner_model") + return TearingProblem(ffs, SLAYERControl(; enabled=true, kwargs...)) +end + +# Attached profiles as a KineticProfiles table (eV; ω* recomputed downstream; no χ, so the scalar fallbacks apply). +function _profiles_from_equilibrium(kp) + xs = kp.xs + E_CHG = Utilities.PhysicalConstants.E_CHG + return (profiles=KineticProfiles(; + psi=xs, + n_e=[kp.ne_spline(ψ) for ψ in xs], + T_e=[kp.Te_spline(ψ) / E_CHG for ψ in xs], + T_i=[kp.Ti_spline(ψ) / E_CHG for ψ in xs], + omega=[kp.omegaE_spline(ψ) for ψ in xs], + omega_e=zeros(length(xs)), + omega_i=zeros(length(xs))), + chi_perp=nothing, chi_tor=nothing) +end + +function _tearing_params(::GGJModel, equil, surfaces, loaded, control) + lp = layer_parameters(surfaces, equil; profiles=loaded.profiles, + mu_i=control.mu_i, + zeff=control.zeff, + resistivity_model=_build_resistivity_model(control.resistivity_model), + lnLambda_form=control.lnLambda_form) + return [ggj_parameters(s, equil; eta=lp.eta[k], rho=lp.rho[k], ising=k) for (k, s) in enumerate(surfaces)] +end + +function _tearing_params(::SLAYERModel, equil, surfaces, loaded, control) + # `equil.config.b0exp` is a NORMALIZATION (commonly exactly 1.0), not the toroidal + # field, so substituting it here silently ran the layer physics at B_T = 1 T. Pass the + # control value through instead: `nothing` makes build_slayer_inputs compute the + # physical B_T = F(psi)/(2*pi*R_0) per surface from the equilibrium's F-spline, which is + # what its docstring already prescribes. + bt = control.bt + # χ⊥/χ_φ from the kinetic file when present, else the scalar fallbacks. + chi_perp = loaded.chi_perp === nothing ? control.chi_perp : loaded.chi_perp + chi_tor = loaded.chi_tor === nothing ? control.chi_tor : loaded.chi_tor + (loaded.chi_perp === nothing || loaded.chi_tor === nothing) && @warn( + "SLAYER: no usable chi_e/chi_phi profile(s) (dataset absent or all-zero); " * + "using the scalar chi_perp/chi_tor fallback for the missing one(s).") + return build_slayer_inputs(equil, surfaces, loaded.profiles; + bt=bt, + mu_i=control.mu_i, + zeff=control.zeff, + chi_perp=chi_perp, + chi_tor=chi_tor, + dr_val=control.dr_val, + dgeo_val=control.dgeo_val, + dc_type=control.dc_type, + theta=control.theta_sample, + resistivity_model=_build_resistivity_model(control.resistivity_model), + lnLambda_form=control.lnLambda_form) +end + +""" + solve(prob::TearingProblem, model) -> SLAYERResult + +Tearing growth rates per surface with the inner layer of `model`. +""" +function CommonSolve.solve(prob::TearingProblem, model::InnerLayer.InnerLayerModel) + control = prob.control + ffs = prob.ffs + equil = ffs.equil + surfaces = ffs.surfaces + + validate(control) + control.enabled || return empty_slayer_result(control) + isempty(surfaces) && return empty_slayer_result(control) + + loaded = if !isempty(control.profile_file) + _load_profiles(control, ffs.dir_path) + elseif equil.kinetic !== nothing + _profiles_from_equilibrium(equil.kinetic) + else + error("TearingProblem: no kinetic profiles — set profile_file, or attach profiles " * + "to the equilibrium with attach_kinetic_profiles!(equil, file)") + end + + params = _tearing_params(model, equil, surfaces, loaded, control) + + # Δ' matrix: prefer the full inter-surface matrix; fall back to a + # diagonal built from each SingType's scalar delta_prime. + dpm = ffs.delta_prime === nothing ? Matrix{ComplexF64}(undef, 0, 0) : ffs.delta_prime.matrix + dp = if !isempty(dpm) && size(dpm) == (length(params), length(params)) + Matrix{ComplexF64}(dpm) + else + # The full Δ' matrix is unavailable (e.g. the parallel-FM stage that + # populates it was not run). The scalar-diagonal fallback uses + # `sing.delta_prime`, which is a coarse per-surface stub; surfaces + # with no entry default to Δ'=0, giving γ computed from zero drive. + n_missing = count(s -> isempty(s.delta_prime), surfaces) + @warn( + "SLAYER: delta_prime_matrix is empty or wrong-sized " * + "($(size(dpm)) vs " * + "($(length(params)),$(length(params)))); falling back to the " * + "diagonal `sing.delta_prime` stub. Growth rates use a coarse " * + "per-surface Δ' and may be unreliable" * + (n_missing > 0 ? "; $n_missing surface(s) have NO Δ' entry and " * + "default to Δ'=0 (zero tearing drive)." : ".") + ) + M = zeros(ComplexF64, length(params), length(params)) + for (k, s) in enumerate(surfaces) + M[k, k] = isempty(s.delta_prime) ? 0.0 + 0im : s.delta_prime[1] + end + M + end + + rational_psi = Float64[surfaces[p.ising].psifac for p in params] + rational_q = Float64[surfaces[p.ising].q for p in params] + return run_slayer_from_inputs(model, params, dp, control; + rational_psi=rational_psi, rational_q=rational_q) +end diff --git a/src/Tearing/Runner/run_slayer.jl b/src/Tearing/Runner/run_slayer.jl index 2459087cf..a4750fa6a 100644 --- a/src/Tearing/Runner/run_slayer.jl +++ b/src/Tearing/Runner/run_slayer.jl @@ -1,17 +1,7 @@ -# Runner.jl +# run_slayer.jl # -# Top-level orchestration for the SLAYER tearing-mode analysis. Given a -# `ForceFreeStatesResult` (which supplies the equilibrium, the rational-surface -# list and the outer-region Δ' matrix) + a -# populated `SLAYERControl`, `run_slayer` loads kinetic profiles, builds -# per-surface SLAYER parameters, runs the requested scan mode, extracts -# growth rates by contour intersection, and returns a `SLAYERResult`. -# -# A secondary entry point `run_slayer_from_inputs` takes pre-built -# per-surface parameters + a Δ' matrix and bypasses the -# equilibrium-driven `build_slayer_inputs` step. This is what the test -# suite drives; it keeps the end-to-end code covered without requiring a -# full equilibrium solve in every test. +# Tearing scan core: profile loading, the slab Δ′ reference conversion, and +# `run_slayer_from_inputs` (scan + growth-rate extraction). `TearingProblem` builds its inputs. # --------------------------------------------------------------------- # Profile loading @@ -62,6 +52,8 @@ end function _build_inner_model(name::Symbol) if name === :slayer_fitzpatrick return SLAYERModel(; variant=:fitzpatrick) + elseif name === :ggj_ray + return GGJModel(; solver=:ray) elseif name === :ggj_shooting return GGJModel(; solver=:shooting) elseif name === :ggj_galerkin @@ -197,17 +189,16 @@ end # Core analysis entry point that takes pre-built parameters. # --------------------------------------------------------------------- """ - run_slayer_from_inputs(params::Vector{SLAYERParameters}, - dp_matrix::AbstractMatrix, - control::SLAYERControl) -> SLAYERResult - -Run the SLAYER tearing analysis given pre-built per-surface -`SLAYERParameters` and the outer-region Δ' matrix. Bypasses the -equilibrium-driven `build_slayer_inputs` step — use this when the -parameters are already known (e.g. in unit tests or when rebuilding -from cached HDF5 output). + run_slayer_from_inputs(model::InnerLayerModel, params, dp_matrix, control) -> SLAYERResult + +Run the tearing analysis scan core given the inner-layer dispatch `model`, the pre-built +per-surface parameters, and the outer-region Δ' matrix. Bypasses the orchestration of the +`TearingProblem` solve — use this when the parameters are already known (e.g. in unit +tests or when rebuilding from cached HDF5 output). The model is a typed object end to +end; nothing below this point reads `control.inner_model`. """ -function run_slayer_from_inputs(params::AbstractVector{<:InnerLayerParameters}, +function run_slayer_from_inputs(model::InnerLayer.InnerLayerModel, + params::AbstractVector{<:InnerLayerParameters}, dp_matrix::AbstractMatrix, control::SLAYERControl; rational_psi::Vector{Float64}=Float64[], @@ -222,18 +213,16 @@ function run_slayer_from_inputs(params::AbstractVector{<:InnerLayerParameters}, "≠ ($n, $n)")) dp = Matrix{ComplexF64}(dp_matrix) - model = _build_inner_model(control.inner_model) - # Guard: the inner-layer model and the parameter eltype must match, or # `_build_surface_coupling` throws an opaque MethodError downstream. expected_P = _is_ggj(model) ? GGJParameters : SLAYERParameters all(p -> p isa expected_P, params) || throw( ArgumentError( - "run_slayer: inner_model=$(control.inner_model) requires " * + "run_slayer: a $(typeof(model)) inner model requires " * "$(expected_P) per-surface parameters, but got eltype " * "$(eltype(params)). Build inputs with the matching builder " * - "(build_slayer_inputs for SLAYER, build_ggj_inputs for GGJ).") + "(build_slayer_inputs for SLAYER, ggj_parameters for GGJ).") ) # Slab-layer path: convert Δ' from its ψ_N reference length to the @@ -396,112 +385,3 @@ function ggj_inner_deltas(params::AbstractVector{GGJParameters}, Q::Number; end return out end - -# --------------------------------------------------------------------- -# Full pipeline: equilibrium + ForceFreeStates → parameters → analysis -# --------------------------------------------------------------------- -""" - run_slayer(result, control; dir_path="./") -> SLAYERResult - -Orchestrate the full SLAYER analysis against a `ForceFreeStates.ForceFreeStatesResult`, -reading its equilibrium, singular surfaces and Δ' matrix. Kinetic profiles are -read from `control.profile_file` (relative to `dir_path`) through the shared -`Equilibrium.read_kinetic_file` reader; when the file carries `chi_e`/`chi_phi` -profiles they set χ⊥(ψ)/χ_φ(ψ), otherwise the scalar `control.chi_perp`/ -`chi_tor` fallbacks are used. - -The toroidal field comes from `control.bt`; leaving it unset (the default) makes -`build_slayer_inputs` evaluate the physical `B_T = F(ψ)/(2π·R₀)` per surface. - -Returns an `enabled=false` `SLAYERResult` when `control.enabled` is -false. -""" -function run_slayer(result, control::SLAYERControl; dir_path::AbstractString="./") - dpm = result.delta_prime === nothing ? Matrix{ComplexF64}(undef, 0, 0) : result.delta_prime.matrix - return run_slayer(result.equil, result.surfaces, dpm, control; dir_path=dir_path) -end - -""" - run_slayer(equil, surfaces, delta_prime_matrix, control; dir_path="./") -> SLAYERResult - -Loose-argument form of [`run_slayer`](@ref), taking the equilibrium, the singular-surface -vector and the outer-region Δ' matrix directly. Per-surface parameters are built via -`build_slayer_inputs`; an empty or wrong-sized `delta_prime_matrix` falls back to a diagonal -built from the `sing.delta_prime` stubs. -""" -function run_slayer(equil, surfaces::AbstractVector, delta_prime_matrix::AbstractMatrix, - control::SLAYERControl; dir_path::AbstractString="./") - validate(control) - control.enabled || return empty_slayer_result(control) - isempty(surfaces) && return empty_slayer_result(control) - - loaded = _load_profiles(control, dir_path) - profiles = loaded.profiles - - if control.inner_model in (:ggj_shooting, :ggj_galerkin) - # GGJ γ-extraction is future work; `run_slayer_from_inputs` emits the - # warning once the model is built (so direct callers see it too). - params = build_ggj_inputs(equil, surfaces, profiles; - mu_i=control.mu_i, - zeff=control.zeff, - resistivity_model=_build_resistivity_model(control.resistivity_model), - lnLambda_form=control.lnLambda_form) - else - # `equil.config.b0exp` is a NORMALIZATION (commonly exactly 1.0), not the toroidal - # field, so substituting it here silently ran the layer physics at B_T = 1 T. Pass the - # control value through instead: `nothing` makes build_slayer_inputs compute the - # physical B_T = F(psi)/(2*pi*R_0) per surface from the equilibrium's F-spline, which is - # what its docstring already prescribes. - bt = control.bt - # χ⊥/χ_φ from the kinetic file when present, else the scalar fallbacks. - chi_perp = loaded.chi_perp === nothing ? control.chi_perp : loaded.chi_perp - chi_tor = loaded.chi_tor === nothing ? control.chi_tor : loaded.chi_tor - (loaded.chi_perp === nothing || loaded.chi_tor === nothing) && @warn( - "SLAYER: kinetic file has no usable chi_e/chi_phi profile(s) " * - "(dataset absent or all-zero); using the scalar " * - "control.chi_perp/chi_tor fallback for the missing one(s).") - params = build_slayer_inputs(equil, surfaces, profiles; - bt=bt, - mu_i=control.mu_i, - zeff=control.zeff, - chi_perp=chi_perp, - chi_tor=chi_tor, - dr_val=control.dr_val, - dgeo_val=control.dgeo_val, - dc_type=control.dc_type, - theta=control.theta_sample, - resistivity_model=_build_resistivity_model(control.resistivity_model), - lnLambda_form=control.lnLambda_form) - end - - # Δ' matrix: prefer the full parallel-FM matrix; fall back to a - # diagonal built from each SingType's scalar delta_prime. - dp = if !isempty(delta_prime_matrix) && - size(delta_prime_matrix) == (length(params), length(params)) - Matrix{ComplexF64}(delta_prime_matrix) - else - # The full Δ' matrix is unavailable (e.g. the parallel-FM stage that - # populates it was not run). The scalar-diagonal fallback uses - # `sing.delta_prime`, which is a coarse per-surface stub; surfaces - # with no entry default to Δ'=0, giving γ computed from zero drive. - n_missing = count(s -> isempty(s.delta_prime), surfaces) - @warn( - "SLAYER: delta_prime_matrix is empty or wrong-sized " * - "($(size(delta_prime_matrix)) vs " * - "($(length(params)),$(length(params)))); falling back to the " * - "diagonal `sing.delta_prime` stub. Growth rates use a coarse " * - "per-surface Δ' and may be unreliable" * - (n_missing > 0 ? "; $n_missing surface(s) have NO Δ' entry and " * - "default to Δ'=0 (zero tearing drive)." : ".") - ) - M = zeros(ComplexF64, length(params), length(params)) - for (k, s) in enumerate(surfaces) - M[k, k] = isempty(s.delta_prime) ? 0.0 + 0im : s.delta_prime[1] - end - M - end - - rational_psi = Float64[surfaces[p.ising].psifac for p in params] - rational_q = Float64[surfaces[p.ising].q for p in params] - return run_slayer_from_inputs(params, dp, control; rational_psi=rational_psi, rational_q=rational_q) -end diff --git a/src/Tearing/Tearing.jl b/src/Tearing/Tearing.jl index 745a30857..27db0129c 100644 --- a/src/Tearing/Tearing.jl +++ b/src/Tearing/Tearing.jl @@ -12,8 +12,7 @@ # `InnerLayer` itself lives at the top level (`src/InnerLayer/`) and is loaded # before `ForceFreeStates`, which depends on it for the matched-Δ′ Galerkin # solve. Tearing re-binds it here so `Dispersion` and `Runner` reach it via -# `..InnerLayer`, and owns `build_ggj_inputs`, the equilibrium/ForceFreeStates -# glue that cannot live inside `InnerLayer` without creating a dependency cycle. +# `..InnerLayer`. module Tearing @@ -21,7 +20,6 @@ using ..Utilities import ..InnerLayer as InnerLayer -include("LayerInputs.jl") include("Dispersion/Dispersion.jl") include("Runner/Runner.jl") @@ -29,6 +27,5 @@ import .Dispersion as Dispersion import .Runner as Runner export InnerLayer, Dispersion, Runner -export build_ggj_inputs end # module Tearing diff --git a/test/runtests.jl b/test/runtests.jl index e1256e56d..747db5344 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -33,6 +33,7 @@ else include("./runtests_parallel_integration.jl") include("./runtests_result_struct.jl") include("./runtests_solve_api.jl") + include("./runtests_matching_models.jl") include("./runtests_dominant_coupling.jl") include("./runtests_error_fields.jl") include("./runtests_tolerance_toml.jl") diff --git a/test/runtests_dispersion_residual.jl b/test/runtests_dispersion_residual.jl index e8af06da4..30ebf24ff 100644 --- a/test/runtests_dispersion_residual.jl +++ b/test/runtests_dispersion_residual.jl @@ -97,7 +97,7 @@ p_ggj = glasser_wang_2020_eq55() sc_ggj = surface_coupling(GGJModel(solver=:shooting), p_ggj, -1.0 + 0.0im) - @test sc_ggj isa SurfaceCoupling{GGJModel{:shooting},GGJParameters} + @test sc_ggj isa SurfaceCoupling{<:GGJModel{:shooting},GGJParameters} @test sc_ggj(1e-3 + 0.0im) isa ComplexF64 end diff --git a/test/runtests_matching_models.jl b/test/runtests_matching_models.jl new file mode 100644 index 000000000..251d4e427 --- /dev/null +++ b/test/runtests_matching_models.jl @@ -0,0 +1,91 @@ +using TOML + +# Inner-layer model capability gates and option forwarding, and the shared layer-parameter +# builder, which must reproduce the resistivity/density closures from `equil.kinetic`. +@testset "Matching models and layer parameters" begin + GPEC = GeneralizedPerturbedEquilibrium + FFS = GPEC.ForceFreeStates + NCR = GPEC.Utilities.NeoclassicalResistivity + E_CHG = GPEC.Utilities.PhysicalConstants.E_CHG + M_P = GPEC.Utilities.PhysicalConstants.M_P + + @testset "model capability gates and option forwarding" begin + IL = GPEC.InnerLayer + @test IL.GGJModel() isa IL.GGJModel{:ray} + @test FFS.closure_capable(IL.GGJModel()) + @test !FFS.closure_capable(IL.SLAYERModel()) + @test !FFS.closure_capable(IL.GGJModel(; solver=:shooting)) # no layer profiles + @test IL.GGJModel(; solver=:galerkin, nx=640).options == (nx=640,) + # Carried options reach the backend exactly as call-site keywords do. + p = IL.glasser_wang_2020_eq55() + γ = 1e-3im + @test IL.solve_inner(IL.GGJModel(; solver=:galerkin, nx=256), p, γ) == IL.solve_inner(IL.GGJModel(; solver=:galerkin), p, γ; nx=256) + @test IL.solve_inner(IL.GGJModel(; solver=:galerkin, nx=256), p, γ; nx=512) == IL.solve_inner(IL.GGJModel(; solver=:galerkin), p, γ) + end + + @testset "override precedence needs no kinetic data" begin + fake_surfaces = [(psifac=0.3,), (psifac=0.6,)] + out = FFS.layer_parameters(fake_surfaces, nothing; + eta=[1e-7, 2e-7], rho=[1e-7, 1e-7], rotation=[0.0, 100.0]) + @test out.eta == [1e-7, 2e-7] + @test out.rho == [1e-7, 1e-7] + @test out.rotation == [0.0, 100.0] + @test_throws ErrorException FFS.layer_parameters(fake_surfaces, nothing; + eta=[1e-7], rho=[1e-7, 1e-7], rotation=[0.0, 0.0]) + # Overriding η alone skips the neoclassical closure, so no surface geometry is needed. + n = 5 + table = GPEC.Utilities.KineticProfiles(; psi=collect(range(0.0, 1.0; length=n)), n_e=fill(5e19, n), + T_e=fill(1e3, n), T_i=fill(1e3, n), omega=fill(100.0, n), omega_e=zeros(n), omega_i=zeros(n)) + partial = FFS.layer_parameters(fake_surfaces, nothing; profiles=table, eta=[1e-7, 2e-7]) + @test partial.eta == [1e-7, 2e-7] && all(>(0), partial.rho) + end + + @testset "derivation from equil.kinetic matches the shared closures" begin + template = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example") + deck = TOML.parsefile(joinpath(template, "gpec.toml")) + equil = GPEC.Equilibrium.setup_equilibrium( + GPEC.Equilibrium.EquilibriumConfig(deck["Equilibrium"], template), + GPEC.Equilibrium.SolovevConfig(deck["SOL_INPUT"]) + ) + ffs_kwargs = Dict(Symbol(k) => v for (k, v) in deck["ForceFreeStates"] + if !(k in ("integrator", "nn_low", "nn_high"))) + kin_file = joinpath(@__DIR__, "..", "examples", "Solovev_kinetic_NTV_example", "kinetic.dat") + attach_kinetic_profiles!(equil, kin_file; zi=1) + + mktempdir() do dir + ffs = solve(equil, Forward(); nn=1, dir_path=dir, ffs_kwargs...) + surfaces = ffs.surfaces + @test !isempty(surfaces) + + # A derivation without kinetic data must fail loudly. + bare = GPEC.Equilibrium.setup_equilibrium( + GPEC.Equilibrium.EquilibriumConfig(deck["Equilibrium"], template), + GPEC.Equilibrium.SolovevConfig(deck["SOL_INPUT"]) + ) + @test_throws ErrorException FFS.layer_parameters(surfaces, bare) + + kp = ffs.equil.kinetic + out = FFS.layer_parameters(surfaces, ffs.equil; resistivity_model=NCR.SpitzerModel()) + for (k, sing) in enumerate(surfaces) + ψ = sing.psifac + n_e = kp.ne_spline(ψ) + t_e = kp.Te_spline(ψ) / E_CHG + lnLamb = NCR.coulomb_log_e(n_e, t_e; form=:nrl) + @test out.eta[k] ≈ NCR.eta_spitzer(n_e, t_e, 1.0; lnLamb=lnLamb) rtol = 1e-12 + @test out.rho[k] ≈ 2.0 * M_P * n_e rtol = 1e-12 + @test out.rotation[k] ≈ kp.omegaE_spline(ψ) / (2π) rtol = 1e-12 + end + + # Partial override: eta passed through verbatim, the rest still derived. + mixed = FFS.layer_parameters(surfaces, ffs.equil; eta=fill(3e-8, length(surfaces)), resistivity_model=NCR.SpitzerModel()) + @test all(mixed.eta .== 3e-8) + @test mixed.rho == out.rho + + # The default closure is Sauter neoclassical: it reads the surface's trapped fraction and stays physical. + neo = FFS.layer_parameters(surfaces, ffs.equil) + @test neo.eta == FFS.layer_parameters(surfaces, ffs.equil; resistivity_model=NCR.SauterNeoModel()).eta + @test all(isfinite, neo.eta) && all(>(0), neo.eta) + @test neo.eta != out.eta + end + end +end diff --git a/test/runtests_resist_eval.jl b/test/runtests_resist_eval.jl index 8e334a67a..2bbe431af 100644 --- a/test/runtests_resist_eval.jl +++ b/test/runtests_resist_eval.jl @@ -6,7 +6,7 @@ using GeneralizedPerturbedEquilibrium.LocalStability using GeneralizedPerturbedEquilibrium.Utilities using GeneralizedPerturbedEquilibrium.InnerLayer - using GeneralizedPerturbedEquilibrium.Tearing: build_ggj_inputs + using GeneralizedPerturbedEquilibrium.ForceFreeStates: layer_parameters, ggj_parameters using FastInterpolations using TOML @@ -17,6 +17,12 @@ sol_cfg = Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) equil = Equilibrium.setup_equilibrium(eq_cfg, sol_cfg) + # Per-surface GGJ parameters from a kinetic-profile table, as the GGJ tearing solve builds them. + function ggj_from(equil, sings, profiles; kwargs...) + lp = layer_parameters(sings, equil; profiles=profiles, kwargs...) + return [ggj_parameters(s, equil; eta=lp.eta[k], rho=lp.rho[k], ising=k) for (k, s) in enumerate(sings)] + end + @testset "resist_geometry: returns finite values with expected signs" begin # Pick a few interior surfaces; compute q1 from the equilibrium dq = deriv_view(equil.profiles.q_spline, 1) @@ -98,7 +104,7 @@ @test intr.sing[1].restype === rg_first end - @testset "build_ggj_inputs: builds GGJParameters from sings + profiles" begin + @testset "ggj_parameters: builds GGJParameters from sings + profiles" begin # Synthetic profiles psi_pts = collect(0.0:0.1:1.0) profiles = KineticProfiles(; psi=psi_pts, @@ -120,7 +126,7 @@ intr = ForceFreeStates.ForceFreeStatesInternal(; sing=[s1], msing=1) ForceFreeStates.resist_eval_all!(intr, equil) - gs = build_ggj_inputs(equil, intr.sing, profiles; mu_i=2.0, zeff=1.0) + gs = ggj_from(equil, intr.sing, profiles; mu_i=2.0, zeff=1.0) @test length(gs) == 1 @test gs[1] isa GGJParameters @@ -143,7 +149,7 @@ @test gs[1].ising == 1 end - @testset "build_ggj_inputs: errors when restype not populated" begin + @testset "ggj_parameters: errors when restype not populated" begin # Need ≥4 points for the cubic spline psi_pts = collect(0.0:0.25:1.0) n = length(psi_pts) @@ -159,7 +165,9 @@ ua_right=zeros(ComplexF64, 0, 0, 0), psi_ua_left=0.0, psi_ua_right=0.0) @test s_unpop.restype === nothing - @test_throws ArgumentError build_ggj_inputs(equil, [s_unpop], profiles) + @test_throws ArgumentError ggj_parameters(s_unpop, equil; eta=1e-7, rho=1e-7) + # The default Sauter closure needs the same geometry, so the derivation refuses too. + @test_throws ErrorException ggj_from(equil, [s_unpop], profiles) end @testset "GGJ solve_inner runs on built parameters" begin @@ -182,7 +190,7 @@ psi_ua_left=0.0, psi_ua_right=0.0) intr = ForceFreeStates.ForceFreeStatesInternal(; sing=[s1], msing=1) ForceFreeStates.resist_eval_all!(intr, equil) - gs = build_ggj_inputs(equil, intr.sing, profiles; mu_i=2.0) + gs = ggj_from(equil, intr.sing, profiles; mu_i=2.0) # Verify D_I < 0 so the GGJ shooting solver doesn't bail @test mercier_di(gs[1]) < 0 diff --git a/test/runtests_result_struct.jl b/test/runtests_result_struct.jl index f02a7ee73..e63936a73 100644 --- a/test/runtests_result_struct.jl +++ b/test/runtests_result_struct.jl @@ -144,7 +144,7 @@ using HDF5 @test ffs.integrator === :galerkin @test ffs.galerkin !== nothing - @test ffs.galerkin.match !== nothing + @test ffs.match !== nothing sol = ffs.solution @test sol isa FFS.SolutionProfiles @@ -184,8 +184,10 @@ using HDF5 @test iszero(ffs.bpen) else @test ffs.closure === :matched - @test ffs.bpen == ffs.galerkin.match.bpen + @test ffs.bpen == ffs.match.bpen @test !iszero(ffs.bpen) + # The ideal δW does not describe the matched plasma. + @test ffs.wp === nothing && ffs.free_boundary === nothing end # The closed profiles land in the shared Solutions layout: same names and @@ -196,10 +198,10 @@ using HDF5 @test read(f["$gal/xi_psi"]) == sol.u_store[:, :, 1, :] @test read(f["$gal/dxi_psidpsi"]) == sol.du_store @test read(f["$gal/xi_s"]) == sol.xi_s_store - # Matching diagnostics keep their group; the profile datasets left it, - # and the raw outer basis is debug-gated off by default. - @test haskey(f, "$gal/Match/cout") - @test !haskey(f, "$gal/Match/xi") + # Matching diagnostics sit with the other per-surface results, and the + # raw outer basis is debug-gated off by default. + @test haskey(f, "SingularSurfaces/Match/cout") + @test !haskey(f, "$gal/Match") @test !haskey(f, "$gal/Basis") @test !haskey(f, "$gal/Solution") end diff --git a/test/runtests_slayer_runner.jl b/test/runtests_slayer_runner.jl index ee114e37f..e33f6ee49 100644 --- a/test/runtests_slayer_runner.jl +++ b/test/runtests_slayer_runner.jl @@ -1,6 +1,7 @@ -@testset "Runner: Control + run_slayer + HDF5 output" begin +@testset "Runner: Control + TearingProblem + HDF5 output" begin using GeneralizedPerturbedEquilibrium using GeneralizedPerturbedEquilibrium.InnerLayer + using GeneralizedPerturbedEquilibrium: GGJModel, SLAYERModel, TearingProblem, solve using GeneralizedPerturbedEquilibrium.Dispersion using GeneralizedPerturbedEquilibrium.Runner using HDF5 @@ -94,23 +95,31 @@ @test_throws ArgumentError slayer_control_from_toml(bad) end - @testset "run_slayer: result-facing form forwards surfaces and Δ'" begin + @testset "TearingProblem solve forwards surfaces and Δ'" begin # A result with no singular surfaces short-circuits before any equilibrium access, # so a stand-in result is enough to pin the forwarding of surfaces / delta_prime. c = SLAYERControl(; enabled=true, profile_file="unused.h5") no_surfaces = (equil=nothing, surfaces=GeneralizedPerturbedEquilibrium.ForceFreeStates.SingType[], delta_prime=nothing) - r = run_slayer(no_surfaces, c) + r = solve(TearingProblem(no_surfaces, c), SLAYERModel()) @test isempty(r.params) # A disabled control never looks at the result at all. - r_off = run_slayer(no_surfaces, SLAYERControl(; enabled=false, profile_file="unused.h5")) + r_off = solve(TearingProblem(no_surfaces, SLAYERControl(; enabled=false, profile_file="unused.h5")), SLAYERModel()) @test r_off.enabled == false end + @testset "the default :ray backend drives the tearing solve" begin + @test SLAYERControl(; inner_model=:ggj_ray).inner_model === :ggj_ray + no_surfaces = (equil=nothing, surfaces=GeneralizedPerturbedEquilibrium.ForceFreeStates.SingType[], + delta_prime=nothing, dir_path=".") + r = solve(TearingProblem(no_surfaces, SLAYERControl(; enabled=true, profile_file="unused.h5")), GGJModel()) + @test isempty(r.params) + end + # Δ′ is unified across formalisms, so a Galerkin run feeds SLAYER exactly as a Riccati one # does: `result.delta_prime.matrix` is populated and already sized to the surface list, which - # is the predicate `run_slayer` uses to accept it over the per-surface diagonal stub. + # is the predicate the tearing solve uses to accept it over the per-surface diagonal stub. @testset "Galerkin-fed SLAYER: gal Δ' drives the coupled solve" begin mktempdir() do dir deck = joinpath(@__DIR__, "..", "examples", "LAR_ideal_match_test") @@ -130,7 +139,7 @@ _mk_params(; rs=0.6, lu=2.0e7, tauk=1.2e-4, m=3, ising=2)] c = SLAYERControl(; enabled=true, coupling_mode=:coupled, scan_mode=:brute_force, Q_re_range=(-1.0, 1.0), Q_im_range=(-0.5, 0.8), nre=20, nim=20, pole_threshold=1e5) - r = run_slayer_from_inputs(params, dpm, c) + r = run_slayer_from_inputs(SLAYERModel(; variant=:fitzpatrick), params, dpm, c) @test r.enabled @test r.coupled_extraction isa GrowthRateResult end @@ -140,7 +149,7 @@ c = SLAYERControl(; enabled=false) params = [_mk_params()] dp = ComplexF64[0.0 + 0im;;] # 1×1 matrix - r = run_slayer_from_inputs(params, dp, c) + r = run_slayer_from_inputs(SLAYERModel(; variant=:fitzpatrick), params, dp, c) @test r.enabled == false @test isempty(r.Q_root) @test isempty(r.params) @@ -150,7 +159,7 @@ c = SLAYERControl(; enabled=true) params = [_mk_params()] bad_dp = ComplexF64[0.0 0.0; 0.0 0.0] - @test_throws ArgumentError run_slayer_from_inputs(params, bad_dp, c) + @test_throws ArgumentError run_slayer_from_inputs(SLAYERModel(; variant=:fitzpatrick), params, bad_dp, c) end @testset "delta_prime_to_rs_reference: ψ_N → r_s conversion" begin @@ -217,7 +226,7 @@ Q_im_range=(-0.5, 0.8), nre=80, nim=80, pole_threshold=1e5) # tuned for lu^(1/3) scale - r = run_slayer_from_inputs(params, dp, c) + r = run_slayer_from_inputs(SLAYERModel(; variant=:fitzpatrick), params, dp, c) @test r.enabled @test length(r.Q_root) == 1 # single coupled eigenvalue @test abs(r.Q_root[1] - Q_target) < 2e-2 # grid-resolution limited @@ -245,7 +254,7 @@ nre=40, nim=40, pole_threshold=1e5, store_scan=true) - r = run_slayer_from_inputs(params, dp, c; + r = run_slayer_from_inputs(SLAYERModel(; variant=:fitzpatrick), params, dp, c; rational_psi=[0.45, 0.72], rational_q=[2.0, 3.0]) mktemp() do path, io diff --git a/test/runtests_solve_api.jl b/test/runtests_solve_api.jl index d04eb9065..a5b9c7077 100644 --- a/test/runtests_solve_api.jl +++ b/test/runtests_solve_api.jl @@ -183,24 +183,111 @@ using TOML @test kwargs[:gal_nx] == 64 @test kwargs[:gal_rpec_flag] - # A match implies the coil-response columns and fills the inner-layer knobs. - FFS._apply_match!(kwargs, ResistiveMatch(; eta=[1e-6], inner_solver="ray"), Galerkin()) - @test kwargs[:gal_match_flag] - @test kwargs[:gal_rpec_flag] - @test kwargs[:gal_eta] == [1e-6] - @test kwargs[:gal_inner_solver] == "ray" - @test !kwargs[:gal_ideal_flag] - # Every key the objects own is a `ForceFreeStatesControl` field. @test all(in(fieldnames(FFS.ForceFreeStatesControl)), keys(kwargs)) end @testset "rejected keyword combinations" begin - @test_throws ErrorException solve(equil, Forward(); nn=1, dir_path=".", ffs_kwargs..., kinetic_factor=0.5) - @test_throws ErrorException solve(equil, Riccati(); nn=1, dir_path=".", ffs_kwargs..., match=ResistiveMatch()) - @test_throws ErrorException solve(equil, Forward(); nn=1, dir_path=".", ffs_kwargs..., match=ResistiveMatch()) + # A calculated-source kinetic solve gates on profiles attached to the equilibrium. + @test_throws ErrorException solve(equil, Forward(); nn=1, dir_path=".", ffs_kwargs..., + kinetic_factor=0.5, kinetic_source="calculated") @test_throws ErrorException solve(equil, Forward(); nn=1, dir_path=".", ffs_kwargs..., integrator="riccati") @test_throws ErrorException solve(equil, Riccati(); nn=1, dir_path=".", ffs_kwargs..., nchunks=8) @test_throws ErrorException solve(equil, Forward(); nn=1, dir_path=".", ffs_kwargs..., nn_low=2) end + + # Path of the first field where `a` and `b` differ (recursing into structs, arrays and + # dicts), or `nothing` when they are equal everywhere. Underscore fields are private + # lazily-built caches (e.g. a spline's transpose) and are skipped. + function first_diff(a, b, path="") + typeof(a) === typeof(b) || return "$path (type)" + a isa Union{Number,AbstractString,Symbol,Nothing,Function,Type} && return isequal(a, b) ? nothing : path + a isa AbstractArray{<:Number} && return isequal(a, b) ? nothing : path + if a isa AbstractArray + size(a) == size(b) || return "$path (size)" + for i in eachindex(a) + d = first_diff(a[i], b[i], "$path[$i]") + d === nothing || return d + end + return nothing + end + if a isa AbstractDict + keys(a) == keys(b) || return "$path (keys)" + for k in keys(a) + d = first_diff(a[k], b[k], "$path[$k]") + d === nothing || return d + end + return nothing + end + for f in fieldnames(typeof(a)) + startswith(string(f), "_") && continue + isdefined(a, f) == isdefined(b, f) || return "$path.$f (definedness)" + isdefined(a, f) || continue + d = first_diff(getfield(a, f), getfield(b, f), "$path.$f") + d === nothing || return d + end + return nothing + end + + @testset "solves leave their inputs untouched" begin + snap = deepcopy(equil) + mktempdir() do dir + gal = solve(equil, Galerkin(; nx=32, rpec_flag=true, cut_solution=true); nn=1, dir_path=dir, ffs_kwargs...) + @test first_diff(equil, snap) === nothing + # Matching reuses the outer solve: repeated match solves must not alter it. + gal_snap = deepcopy(gal) + n = length(MatchProblem(gal; ideal=true).surfaces) + layer = (eta=fill(1e-7, n), rho=fill(1e-7, n)) + slow = solve(MatchProblem(gal; layer..., rotation=fill(1.0, n)), GGJModel()) + fast = solve(MatchProblem(gal; layer..., rotation=fill(10.0, n)), GGJModel()) + @test slow.bpen != fast.bpen + @test first_diff(gal, gal_snap) === nothing + end + end + + @testset "a resistive match keeps its coefficients and drops the ideal δW" begin + mktempdir() do dir + ric = solve(equil, Riccati(); nn=1, dir_path=dir, ffs_kwargs...) + @test ric.free_boundary !== nothing && !isempty(ric.delta_prime.coil) + n = length(MatchProblem(ric; ideal=true).surfaces) + matched = solve(MatchProblem(ric; eta=fill(1e-7, n), rho=fill(1e-7, n), rotation=fill(1.0, n)), GGJModel()) + @test matched.closure === :matched + @test matched.match !== nothing && matched.bpen == matched.match.bpen + @test matched.wp === nothing && matched.free_boundary === nothing + @test matched.solution === nothing # no Galerkin basis to build a matched ξ from + # A matched result is a new solution, not an outer solve to match again. + @test_throws ErrorException MatchProblem(matched; ideal=true) + end + end + + @testset "MatchProblem gates on its inputs" begin + mktempdir() do dir + # A Forward result carries no Δ′ payload, so the problem is unconstructible. + fwd = solve(equil, Forward(); nn=1, dir_path=dir, ffs_kwargs...) + @test_throws ErrorException MatchProblem(fwd; ideal=true) + # A slab model can never close a matched solution. + gal = solve(equil, Galerkin(; nx=32, rpec_flag=true); nn=1, dir_path=dir, ffs_kwargs...) + prob = MatchProblem(gal; ideal=true) + @test_throws ErrorException solve(prob, SLAYERModel()) + # The tearing model is the solve argument, never a TearingProblem keyword. + @test_throws ErrorException TearingProblem(gal; inner_model=:ggj_ray) + # The ideal reference match keeps the ideal closure and replaces the solution with + # the bare coil columns in the identity-at-edge basis. + matched = solve(prob, GGJModel()) + @test matched.closure === :ideal + @test matched.match !== nothing + # The ideal reference keeps the ideal δW (only a :matched closure drops it). + @test matched.wp === gal.wp && matched.free_boundary === gal.free_boundary + @test matched.solution !== nothing && matched.solution.basis === :gal_native + end + end + + @testset "kinetic profiles live on the equilibrium" begin + kin_file = joinpath(@__DIR__, "..", "examples", "Solovev_kinetic_NTV_example", "kinetic.dat") + eq = attach_kinetic_profiles!(deepcopy(equil), kin_file; zi=1) + @test eq.kinetic isa GPEC.Equilibrium.KineticProfileSplines + @test eq.kinetic.ni_spline(0.5) > 0 + # Kinetic solves stay TOML-driven even with profiles attached. + @test_throws ErrorException solve(eq, Forward(); nn=1, dir_path=".", ffs_kwargs..., kinetic_factor=0.5) + end end