diff --git a/.gitignore b/.gitignore index 8ecc70e4..d7121fb4 100644 --- a/.gitignore +++ b/.gitignore @@ -17,6 +17,7 @@ !docs/resources/**/*.pdf *.png !docs/src/assets/**/*.png +!docs/src/islands/figures/*.png *.jld2 .gitattributes Manifest.toml diff --git a/Project.toml b/Project.toml index 31cd4deb..3d538155 100644 --- a/Project.toml +++ b/Project.toml @@ -18,6 +18,7 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" IMASdd = "c5a45a97-b3f9-491c-b9a7-aa88c3bc0067" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" +Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7" LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" @@ -48,6 +49,7 @@ ForwardDiff = "1.4.1" HDF5 = "0.17.2" IMASdd = "8" JLD2 = "0.6.3" +Krylov = "0.10.8" LaTeXStrings = "1.4.0" LinearAlgebra = "1" OrdinaryDiffEq = "6.102.0" diff --git a/benchmarks/islands/README.md b/benchmarks/islands/README.md new file mode 100644 index 00000000..fe273953 --- /dev/null +++ b/benchmarks/islands/README.md @@ -0,0 +1,23 @@ +# Islands benchmark ladder (docs/05) + +Benchmark scripts for the Islands verification ladder +(`docs/src/islands/design/05-verification.md`). Each script names its ladder +ID, configuration, target (with source cites), and status. + +**Status policy (the `[VERIFY]` discipline, `src/Islands/CLAUDE.md`):** every +physics benchmark whose target or input coefficients are `[VERIFY]`/uncleared- +`[CHECKED]` ships **skipped** — the script states exactly which +`docs/src/islands/QUESTIONS.md` entries gate it and exits without running. +Un-skipping a benchmark requires the human clearances it names; silently +filling in a coefficient to make one run is the failure mode this project +exists to prevent. The structural A-ladder (A1–A8) runs in CI via +`test/runtests_islands_*.jl`, not here. + +Figure scripts (docs/07 pipeline) live in `figures/` and read archived +benchmark data only. + +| Script | Ladder ID | Status | +|---|---|---| +| `benchmark_B2_large_w_limits.jl` | B2 | SKIPPED — gated on Q2/Q3/Q4 | +| `benchmark_B4_polarization_omegaE.jl` | B4 | SKIPPED — gated on Q2/Q3 | +| `benchmark_B5_york_thresholds.jl` | B5a/B5b/B5c | SKIPPED — gated on Q2/Q3/Q4 | diff --git a/benchmarks/islands/benchmark_B2_large_w_limits.jl b/benchmarks/islands/benchmark_B2_large_w_limits.jl new file mode 100644 index 00000000..f9123e84 --- /dev/null +++ b/benchmarks/islands/benchmark_B2_large_w_limits.jl @@ -0,0 +1,14 @@ +# benchmark_B2_large_w_limits.jl — ladder B2 (docs/05 §B) +# +# Target: large-w analytic limits — Δ_bs + Δ_cur ∝ 1/w with coefficient against +# WCHH96 Eq. (85) mapped to the island frame (Diss19 p. 86 frame caveat), and the +# Δ_pol ∝ 1/w³ tail [CHECKED: Diss19 pp. 84–86; D21 Fig. 8-class curves]. +# +# STATUS: SKIPPED — gated per docs/src/islands/QUESTIONS.md: +# Q3 the Δ moment prefactors and electron-closure constants are uncleared +# Q4 WCHH96 is not yet in the docs/08 reference library (cited only via +# transcriptions) and the ψ̃ amplitude carries an open [VERIFY] +# Q2 D7 (re-derived equation set) unratified. + +println("SKIPPED: B2 large-w limits — gated on QUESTIONS Q2, Q3, Q4") +println(" (WCHH96 not in the reference library; ψ̃ [VERIFY] open; closure constants uncleared).") diff --git a/benchmarks/islands/benchmark_B4_polarization_omegaE.jl b/benchmarks/islands/benchmark_B4_polarization_omegaE.jl new file mode 100644 index 00000000..d16f9b08 --- /dev/null +++ b/benchmarks/islands/benchmark_B4_polarization_omegaE.jl @@ -0,0 +1,16 @@ +# benchmark_B4_polarization_omegaE.jl — ladder B4 (docs/05 §B) +# +# Target: polarization structure — (i) Wilson–Connor collisionless and Smolyakov +# collisional scalings; (ii) Δ_pol ∝ ω_E² away from zero with the sign reversal +# at ω_E ≈ −0.89 ω_dia,e, reversal point insensitive to w/ρ_θi; (iii) +# torque-balance roots at discrete ω̂_E (Diss19 benchmark root ω₀ = −0.93 +# ω_dia,e) [CHECKED: D23b Fig. 8; Diss19 Fig. 4.18]. +# +# STATUS: SKIPPED — gated per docs/src/islands/QUESTIONS.md: +# Q3 the frame-convention signs (src/Islands/frames/) are NaN-gated until +# cleared — the ω_E-dependence of Δ_pol is exactly the physics the frames +# module must own before any sign-bearing benchmark runs +# Q2 D7 unratified (L0 equation set). + +println("SKIPPED: B4 polarization ω_E structure — gated on QUESTIONS Q2, Q3") +println(" (frame-convention signs NaN-gated; L0 equation set unratified).") diff --git a/benchmarks/islands/benchmark_B5_york_thresholds.jl b/benchmarks/islands/benchmark_B5_york_thresholds.jl new file mode 100644 index 00000000..1766b183 --- /dev/null +++ b/benchmarks/islands/benchmark_B5_york_thresholds.jl @@ -0,0 +1,22 @@ +# benchmark_B5_york_thresholds.jl — ladder B5a/B5b/B5c (docs/05 §B) +# +# Targets (transcribed, awaiting human sign-off — see docs/src/islands/design/05): +# B5a DK-NTM threshold w_c ≃ 2.76 ρ_θi (half-width) ≡ 8.73 ρ_bi at ε = 0.1, +# :original drift model [CHECKED: I19 Fig. 9; PRL p. 4] +# + open [VERIFY]: run collisionality (I19 §4.2 ν_★ = 0.01 vs L23 p. 82 ν_★ = 1e-3). +# B5b RDK-NTM improved-model threshold w_c ≈ 0.45 ρ_θi ≡ 1.46 ρ_bi half-width +# [CHECKED: D21 abstract + Fig. 8; D23a abstract] +# B5c kokuchou finite-ν_★ surface w_c ≈ 0.440 ρ̂_θi + 0.0178 ν_★ − 7.54e-5 +# [CHECKED: L23 Eqs. 6.3.1–6.3.2], L23-amended equation set. +# +# STATUS: SKIPPED — running this requires assigning Level-0 physics +# coefficients that are at most [CHECKED] and not human-cleared: +# Q2 ratify D7 (re-derived L0 equation set) and D8 (the B5a/b/c triangle) +# Q3 clear the L0 coefficient set (ω̂_D/L̂_B toggle, collision kernel, +# electron closure, quasineutrality closure) +# Q4 resolve the ψ̃ amplitude [VERIFY] and the B5a collisionality [VERIFY] +# per docs/src/islands/QUESTIONS.md. Do not fill in a number to make this run. + +println("SKIPPED: B5a/B5b/B5c York thresholds — gated on QUESTIONS Q2, Q3, Q4") +println(" (uncleared [CHECKED] L0 coefficients; D7/D8 unratified).") +println(" See docs/src/islands/QUESTIONS.md and docs/src/islands/design/05-verification.md.") diff --git a/benchmarks/islands/figures/README.md b/benchmarks/islands/figures/README.md new file mode 100644 index 00000000..f0b2ad83 --- /dev/null +++ b/benchmarks/islands/figures/README.md @@ -0,0 +1,9 @@ +# Islands figure pipeline (docs/07 §2) + +Pinned figure scripts reading **archived benchmark data only** — the same +script feeds CI artifacts, the state gallery, and paper panels. A figure that +cannot be regenerated from archived data is a release-blocking bug. + +Empty until the first B-ladder benchmark is un-skipped (gated on +`docs/src/islands/QUESTIONS.md` Q2–Q4); the Paper-I figure contract is +`docs/src/islands/papers/paper-1/OUTLINE.md`. diff --git a/benchmarks/islands/figures/make_structural_figures.jl b/benchmarks/islands/figures/make_structural_figures.jl new file mode 100644 index 00000000..448dfbbe --- /dev/null +++ b/benchmarks/islands/figures/make_structural_figures.jl @@ -0,0 +1,151 @@ +# make_structural_figures.jl — pinned figure script (design docs/07 §2) +# +# Regenerates the structural (A-ladder) figures embedded in +# docs/src/islands/numerics.md from the verification machinery itself — no +# archived data needed yet because everything here is deterministic M1/M2 +# structure (manufactured coefficients, no physics). Physics (B-ladder) figures +# will read archived benchmark data once QUESTIONS Q2–Q4 clear. +# +# Usage (repo root): +# env -u LD_LIBRARY_PATH julia --project=. benchmarks/islands/figures/make_structural_figures.jl +# +# Writes PNGs into docs/src/islands/figures/ (committed as docs assets so the +# Documenter CI build does not need to run this script). + +ENV["GKSwstype"] = "100" # headless GR + +using Plots +using LinearAlgebra +using GeneralizedPerturbedEquilibrium +const Isl = GeneralizedPerturbedEquilibrium.Islands +const PS = Isl.PhaseSpace +const Op = Isl.Operators +const So = Isl.Solvers +const V = Isl.Verify +const Fi = Isl.Fields + +outdir = normpath(joinpath(@__DIR__, "..", "..", "..", "docs", "src", "islands", "figures")) +mkpath(outdir) +saved = String[] +function save!(p, name) + path = joinpath(outdir, name) + savefig(p, path) + push!(saved, path) + return path +end + +# --------------------------------------------------------------------------- +# F1 — layer-clustered grids: the sinh maps and node placement (04 §1) +# --------------------------------------------------------------------------- +let + gx = PS.MappedFDGrid(33; halfwidth=6.0, clustering=2.0, order=4) + gy = PS.MappedFDGrid(25; halfwidth=4.0, clustering=2.0, center=1.0, domain=:half, order=4) + s = range(-1, 1; length=33) + p1 = plot(s, gx.nodes; lw=2, xlabel="computational coordinate s", ylabel="x", + label="x(s), β=2", legend=:topleft, title="radial map (packs x = 0)") + scatter!(p1, s, gx.nodes; ms=3, label="nodes") + p2 = plot(; xlabel="node index", ylabel="y", title="pitch grid (packs y_c = 1)", legend=:topleft) + scatter!(p2, 1:gy.n, gy.nodes; ms=4, label="y nodes, β=2") + hline!(p2, [1.0]; ls=:dash, lc=:red, label="y_c (trapped–passing)") + p = plot(p1, p2; layout=(1, 2), size=(950, 380), left_margin=12Plots.mm, bottom_margin=6Plots.mm) + save!(p, "grids_clustering.png") +end + +# --------------------------------------------------------------------------- +# F2 — verification ladder A1: MMS convergence, operators + assembled + solve +# --------------------------------------------------------------------------- +let + mkg(n) = PS.IslandGrid(; nx=n, nxi=8, ny=n, nE=3, halfwidth_x=6.0, clustering_x=1.0, + y_max=4.0, y_c=1.0, clustering_y=0.8, order=4) + ns = [9, 17, 33] + p = plot(; xscale=:log10, yscale=:log10, xlabel="radial/pitch resolution n", + ylabel="max error", legend=:bottomleft, title="MMS convergence (ladder A1)", + size=(700, 480), left_margin=12Plots.mm, bottom_margin=6Plots.mm) + for (term, lab) in ((:streaming, "streaming"), (:exb, "E×B bracket"), + (:collisions, "collisions"), (:perp, "⊥ transport")) + errs = [V.mms_operator_error(mkg(n), term) for n in ns] + plot!(p, ns, errs; marker=:circle, lw=2, label=lab) + end + errs = [V.mms_assembled_error(mkg(n)) for n in ns] + plot!(p, ns, errs; marker=:square, lw=3, lc=:black, label="assembled residual") + solve_errs = [V.solve_mms(n).err for n in (9, 17, 33)] + plot!(p, [9, 17, 33], solve_errs; marker=:diamond, lw=3, ls=:dash, label="converged solve (A1-solve)") + guide = errs[1] .* (ns[1] ./ ns) .^ 4 + plot!(p, ns, guide; ls=:dot, lc=:gray, label="4th order") + save!(p, "mms_convergence.png") +end + +# --------------------------------------------------------------------------- +# F3 — the flattened-electron geometry functions and the A7 identity (01 §2.4) +# --------------------------------------------------------------------------- +let + Ωout = 1.02:0.02:5.0 + Ωin = -0.98:0.02:0.98 + Q = [Fi.Q_omega(Ω) for Ω in Ωout] + Qin = [Fi.Q_omega(Ω) for Ω in Ωin] + h = [Fi.h_profile(Ω; prefactor=1.0) for Ω in Ωout] + p1 = plot(Ωin, Qin; lw=2, label="Q(Ω), inside", xlabel="Ω", ylabel="Q", + title="Q(Ω) = (1/2π)∮√(Ω+cos ξ) dξ", legend=:topleft) + plot!(p1, Ωout, Q; lw=2, label="Q(Ω), outside") + vline!(p1, [1.0]; ls=:dash, lc=:red, label="separatrix") + a7 = [abs(Fi.flat_average_d2h_dx2(Ω, 1.0)) for Ω in (1.2, 2.0, 3.0, 5.0)] + p2 = plot(Ωout, h; lw=2, label="h(Ω) (unit prefactor)", xlabel="Ω", ylabel="h", + title="electron profile h(Ω)", legend=:topleft) + vline!(p2, [1.0]; ls=:dash, lc=:red, label="separatrix (h ≡ 0 inside)") + annotate!(p2, 1.6, 0.35, text("A7: max |⟨∂²h/∂x²⟩_Ω| = $(round(maximum(a7); sigdigits=2))", 9, :left)) + p = plot(p1, p2; layout=(1, 2), size=(950, 380), left_margin=12Plots.mm, bottom_margin=6Plots.mm, right_margin=6Plots.mm) + save!(p, "hQ_profiles.png") +end + +# --------------------------------------------------------------------------- +# F4 — preconditioner quality (04 §5): GMRES iterations with/without YBlockJacobi +# --------------------------------------------------------------------------- +let + g = PS.IslandGrid(; nx=9, nxi=8, ny=9, nE=2, halfwidth_x=6.0, clustering_x=1.0, + y_max=4.0, y_c=1.0, clustering_y=0.8, order=4) + nx, nξ, ny, nE, nσ = PS.nnodes(g) + P = @. g.y.nodes * (4.0 - g.y.nodes) + K, = Op.conservative_pitch_operator(g.y, P, ones(ny)) + cstiff = fill(30.0, nx, nξ, nE, nσ) + shift = fill(-1.0, nx, nξ, ny, nE, nσ) + stack = Op.IslandStack((Op.PitchAngleDiffusion(K, cstiff), Op.RadiationSink(shift)), + Op.Quasineutrality(1.3)) + f0! = So.flat_residual(stack, g) + N = Op.statelength(g) + b = sin.((1:N) ./ 7) + f!(out, u) = (f0!(out, u); out .-= b; out) + pc = So.YBlockJacobi(g, (ix, iξ, iE, iσ) -> I(ny) + cstiff[ix, iξ, iE, iσ] .* K; phi_scale=-1.3) + s0 = So.newton_krylov(f!, zeros(N); rtol=1e-10, memory=300) + s1 = So.newton_krylov(f!, zeros(N); rtol=1e-10, memory=300, precond=pc) + p = bar(["unpreconditioned", "y-block Jacobi (TSVD)"], [s0.gmres_iters, s1.gmres_iters]; + ylabel="total GMRES iterations", legend=false, + title="preconditioner: stiff collisional solve", + ylims=(0, 1.25 * s0.gmres_iters), + size=(600, 430), left_margin=12Plots.mm, bottom_margin=6Plots.mm, top_margin=4Plots.mm) + annotate!(p, [(1, s0.gmres_iters + 0.06 * s0.gmres_iters, text("$(s0.gmres_iters)", 10)), + (2, s1.gmres_iters + 0.06 * s0.gmres_iters, text("$(s1.gmres_iters)", 10))]) + save!(p, "preconditioner_gmres.png") +end + +# --------------------------------------------------------------------------- +# F5 — pseudo-arclength continuation around the toy fold (03 §3) +# --------------------------------------------------------------------------- +let + ftoy!(out, u, p) = (out[1] = u[1]^2 + p; out) + pa = So.pseudo_arclength(ftoy!, [1.0], -1.0; ds=0.15, nsteps=30, rtol=1e-12, atol=1e-12) + us = [z[1] for z in pa.us] + p = plot(pa.ps, us; marker=:circle, lw=2, xlabel="parameter p", ylabel="u", + label="continuation path", legend=:topleft, + title="pseudo-arclength steps around the fold (u² + p = 0)", + size=(650, 430), left_margin=12Plots.mm, bottom_margin=6Plots.mm) + plot!(p, -1.2:0.01:0.0, sqrt.(-(-1.2:0.01:0.0)); ls=:dash, lc=:gray, label="analytic ±√(−p)") + plot!(p, -1.2:0.01:0.0, -sqrt.(-(-1.2:0.01:0.0)); ls=:dash, lc=:gray, label="") + if !isempty(pa.folds) + k = pa.folds[1] + 1 + scatter!(p, [pa.ps[k]], [us[k]]; ms=8, mc=:red, label="fold detected") + end + save!(p, "continuation_fold.png") +end + +println("Saved figures:") +foreach(f -> println(" ", f), saved) diff --git a/docs/make.jl b/docs/make.jl index 8f701374..5fddf45d 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -41,6 +41,22 @@ makedocs(; "Analysis" => "analysis.md", "Utilities" => "utilities.md" ], + "Islands" => [ + "Overview" => "islands/index.md", + "Numerics (as implemented)" => "islands/numerics.md", + "Paper I — figure contract" => "islands/papers/paper-1/OUTLINE.md", + "Design documents" => [ + "00 — Roadmap" => "islands/design/00-roadmap.md", + "01 — Level-0 physics" => "islands/design/01-physics-level0.md", + "02 — Species and EPs" => "islands/design/02-species-and-eps.md", + "03 — Architecture" => "islands/design/03-architecture.md", + "04 — Numerics" => "islands/design/04-numerics.md", + "05 — Verification ladder" => "islands/design/05-verification.md", + "06 — Autonomy and tooling" => "islands/design/06-autonomy-and-tooling.md", + "07 — Documentation and papers" => "islands/design/07-documentation-and-papers.md", + "08 — Reference library" => "islands/design/08-reference-library.md" + ] + ], "Citations" => "citations.md", "Developer Notes" => "developer_notes.md", ], diff --git a/docs/src/islands.md b/docs/src/islands.md index 3317325a..0a616786 100644 --- a/docs/src/islands.md +++ b/docs/src/islands.md @@ -3,29 +3,55 @@ The Islands module is a steady-state, multi-species drift-kinetic solver for the resonant magnetic island/layer region in tokamaks — the nonlinear analog of SLAYER, generalizing the Modified Rutherford Equation. It is under active -development. The design documents live under `islands/design/` in the -documentation source tree, and the module conventions in `src/Islands/CLAUDE.md`. +development; the module conventions live in `src/Islands/CLAUDE.md`. -## Status — milestone M1 (skeleton) +This page is the **API reference**. The narrative documentation lives in the +**Islands** section of this site: the [project overview](islands/index.md), +the equations-and-figures chapter of what is +[implemented and verified so far](islands/numerics.md), the +[Paper I figure contract](islands/papers/paper-1/OUTLINE.md), and the full +[design document set](islands/design/00-roadmap.md). -M1 lands the numerical skeleton, not the physics numbers (module `CLAUDE.md`, the -`[VERIFY]` policy): +## Status — M1 skeleton + M2 Level-0 solve machinery (structure, gated physics) + +M1 landed the numerical skeleton and M2 the Level-0 solve machinery — structure +only, no physics numbers (module `CLAUDE.md`, the `[VERIFY]` policy): - `Islands.PhaseSpace` — the `(x, ξ, λ→y, E, σ)` phase-space grids with layer-clustered mappings (design `04 §1`): Fourier spectral `∂ξ`, high-order finite-difference `∂x`/`∂y` on stretched grids, and Gauss energy quadrature. Pure numerics; no physics coefficients. - `Islands.Operators` — the `AbstractTerm` operator stack and residual assembly - (design `03 §2`) as allocation-free, AD-compatible *structural stubs*. Every - physics coefficient is a supplied data field, never a literal — literature - numbers stay `[VERIFY]`/`[CHECKED]`-gated until human-cleared. - - `Islands.Verify` — the manufactured-solution (MMS) and AD-vs-finite-difference - JVP harness backing verification ladder **A1/A2** (design `05 §A`), exercised - by `test/runtests_islands_{grids,operators}.jl`. + (design `03 §2`) as allocation-free, AD-compatible structure, including the + mimetic (exactly conservative) pitch-angle diffusion and the neoclassical- + matching far-field boundary conditions (`01 §3` — never bare Neumann). Every + physics coefficient is a supplied data field, never a literal. + - `Islands.SpeciesLists` — first-class species arrays (`02 §1`, Decision D3): + backgrounds, `Bulk`/`Trace` roles, trace-criteria checks that warn and never + silently degrade. + - `Islands.Frames` — THE frequency/frame conversion module (`01 §5`): the + conversion *forms* with every sign/normalization a NaN-gated + `FrameConvention` field until human-cleared (QUESTIONS Q3). + - `Islands.Fields` — the Level-0 quasineutrality closure structure: `Q(Ω)`, + `h(Ω)` (supplied prefactor), the coefficient-free identity + `⟨∂²h/∂x²⟩_Ω = 0` (ladder A7), and the NaN-gated `ElectronClosure` constants. + - `Islands.Moments` — `J̄_∥` assembly and the `Δ_cos`/`Δ_sin` Ampère + projections (`01 §4`) with **required, gated** prefactors; `Ω`-average and + channel-split diagnostics. + - `Islands.Solvers` — matrix-free Newton–Krylov (ForwardDiff JVP + GMRES, + Eisenstat–Walker forcing), the TSVD-regularized physics-block preconditioner + skeleton (`04 §3, §5`), tiny-grid dense debug Jacobian, and pseudo-arclength + continuation with fold detection. + - `Islands.Verify` — MMS/JVP harness (ladder A1/A2), solve-level MMS and + zero-drive configurations (A5), and the `y_c`-block conditioning monitor + (A8), exercised by `test/runtests_islands_{grids,operators,solve}.jl`. -The physics operators (drift-frequency coefficients, collision kernels, the -Δ moments and York thresholds) land in later milestones and remain gated until -their `[VERIFY]` tags are cleared. +Structural gates **A1–A5, A7 (coefficient-free part), A8** run in CI. The +physics numbers (drift-frequency coefficients, collision kernels, closure +constants, the Δ prefactors, York thresholds) remain `[VERIFY]`-gated: the +B-ladder benchmarks in `benchmarks/islands/` ship skipped, each naming the +`QUESTIONS.md` entries (Q2–Q4) whose human clearance un-gates it. The Paper-I +figure contract is `docs/src/islands/papers/paper-1/OUTLINE.md`. ## API Reference @@ -39,12 +65,42 @@ Modules = [GeneralizedPerturbedEquilibrium.Islands] Modules = [GeneralizedPerturbedEquilibrium.Islands.PhaseSpace] ``` +### Species (`Islands.SpeciesLists`) + +```@autodocs +Modules = [GeneralizedPerturbedEquilibrium.Islands.SpeciesLists] +``` + +### Frames and parameters (`Islands.Frames`) + +```@autodocs +Modules = [GeneralizedPerturbedEquilibrium.Islands.Frames] +``` + ### Operator stack (`Islands.Operators`) ```@autodocs Modules = [GeneralizedPerturbedEquilibrium.Islands.Operators] ``` +### Field-equation closure structure (`Islands.Fields`) + +```@autodocs +Modules = [GeneralizedPerturbedEquilibrium.Islands.Fields] +``` + +### Output moments (`Islands.Moments`) + +```@autodocs +Modules = [GeneralizedPerturbedEquilibrium.Islands.Moments] +``` + +### Newton–Krylov solve (`Islands.Solvers`) + +```@autodocs +Modules = [GeneralizedPerturbedEquilibrium.Islands.Solvers] +``` + ### Verification harness (`Islands.Verify`) ```@autodocs diff --git a/docs/src/islands/LOG.md b/docs/src/islands/LOG.md index 3c7dfba9..2119e8e9 100644 --- a/docs/src/islands/LOG.md +++ b/docs/src/islands/LOG.md @@ -8,6 +8,57 @@ relevant. --- +## 2026-07-08 — M2 L0 solve machinery: Newton–Krylov + moments + species/frames/fields (structure, gated physics) + +- **Contract**: `docs/src/islands/design/M2-launch-prompt.md` (interactive /goal + run). Branch `feature/islands-m2`, **PR #324** (stacked on + `feature/islands-m1`/PR #320; retargets to `feature/islands` when #320 merges). + Full suite green locally. +- **Moved**: the full L0 solve *structure*, every physics coefficient a supplied + `[VERIFY]`-gated parameter (physics-verifier: **PASS**): + - `solvers/` — matrix-free Newton–Krylov (Krylov.jl GMRES on a preallocated + ForwardDiff JVP; Eisenstat–Walker; line search; convergence on norm AND + max-norm per `04 §5`), `YBlockJacobi` physics-block preconditioner with + TSVD-regularized pencil solves (the `04 §3` y_c treatment), dense tiny-grid + debug Jacobian, pseudo-arclength continuation with fold detection (toy fold + found at step 6 of the test problem). + - `species/` (D3 plumbing), `frames/` (conversion forms, NaN-gated + `FrameConvention`), `fields/` (Q(Ω)/h(Ω) structure + NaN-gated + `ElectronClosure`), `moments/` (J̄_∥, Δ projections with required gated + prefactors, ⟨·⟩_Ω diagnostics), operators additions (mimetic + `PitchAngleDiffusion`, `FarFieldConditions` — never bare Neumann, + `weighted_moment!`). + - **Structural gates green** (67 new tests in `runtests_islands_solve.jl`): + A5 (residual exactly 0 at g≡0), solve-MMS at design order (3.98 observed, + nx 17→33), A4 (conservation ≲1e-11, entropy sign exact), A3 parity, A7 + ⟨∂²h/∂x²⟩_Ω ≈ 1e-16, A8 σ_min monitor + singular detection. Preconditioner + cuts a stiff collisional solve from 79.5 s/28 Newton to 0.6 s/7 (GMRES + 1700→200-class); all new kernels pass `--check-bounds=yes`. + - `benchmarks/islands/` created: B2/B4/B5 scripts **skipped**, each naming its + gating QUESTIONS IDs; `regression-harness` case `islands_l0_structural` + (solve-MMS err 5.254e-2, 6 Newton/1210 GMRES, A7 8.0e-17, σ_min 0.1139). + - Paper-I figure contract: `docs/src/islands/papers/paper-1/OUTLINE.md` + (claims C1–C3 green as CI artifacts; C4–C8 gated on Q2–Q4). + - **Rendered docs story** (user-flagged gap vs docs/07's M0–M1 intent): new + `docs/src/islands/numerics.md` — the equations + figures of everything as + implemented — plus the pinned figure script + (`benchmarks/islands/figures/make_structural_figures.jl`, five structural + figures committed as docs assets) and a full "Islands" site-nav section + (overview, numerics chapter, Paper-I contract, design docs 00–08). + Remaining docs/07 infra for later milestones: anchor-sync CI check, + STATE.md dashboard. +- **Physics debugging note**: the first solve-MMS attempt failed to converge — + the generic `Collisions` (a_y ∂²y) term has no y-BCs, so its BVP + discretization is unstable under refinement; the *mimetic* divergence form + (degenerate P → 0 endpoints, zero-flux built in) is the correct structure and + the far-field x-BCs are what make the advective solve well-posed. Exactly the + design's point (`01 §3`, `04 §1`). +- **Blocked**: the York gates (B5a/b/c, B2, B4) and Paper-I claims C4–C8 — all + on the human clearance queue **Q2/Q3/Q4** (unchanged). +- **Next**: human clears Q2–Q4 → thin run fills the L0 coefficients from the + D7 re-derivation and un-skips the B-ladder; independent M2+ work: kinetic- + electron toggle (E4), io/ TOML section, trace-species linear pass. + ## 2026-07-08 — M1 skeleton: phase-space grids + operator stack + MMS/AD harness - **PR**: #320 (`feature/islands-m1` → `feature/islands`); full suite green. diff --git a/docs/src/islands/figures/continuation_fold.png b/docs/src/islands/figures/continuation_fold.png new file mode 100644 index 00000000..6431920b Binary files /dev/null and b/docs/src/islands/figures/continuation_fold.png differ diff --git a/docs/src/islands/figures/grids_clustering.png b/docs/src/islands/figures/grids_clustering.png new file mode 100644 index 00000000..06c1890f Binary files /dev/null and b/docs/src/islands/figures/grids_clustering.png differ diff --git a/docs/src/islands/figures/hQ_profiles.png b/docs/src/islands/figures/hQ_profiles.png new file mode 100644 index 00000000..1c13c178 Binary files /dev/null and b/docs/src/islands/figures/hQ_profiles.png differ diff --git a/docs/src/islands/figures/mms_convergence.png b/docs/src/islands/figures/mms_convergence.png new file mode 100644 index 00000000..5b2a4dd6 Binary files /dev/null and b/docs/src/islands/figures/mms_convergence.png differ diff --git a/docs/src/islands/figures/preconditioner_gmres.png b/docs/src/islands/figures/preconditioner_gmres.png new file mode 100644 index 00000000..d885a9b6 Binary files /dev/null and b/docs/src/islands/figures/preconditioner_gmres.png differ diff --git a/docs/src/islands/numerics.md b/docs/src/islands/numerics.md new file mode 100644 index 00000000..32f9fb74 --- /dev/null +++ b/docs/src/islands/numerics.md @@ -0,0 +1,302 @@ +# Islands — numerics as implemented (M1–M2) + +This chapter documents the Islands machinery **as it exists in the code today** +— the discretization, operator stack, solver, and output-moment assembly landed +by milestones M1 and M2, with the verification evidence behind each piece. It +is the "Physics Book" companion to the aspirational [design documents](design/00-roadmap.md): +the design docs say what Islands *will* compute; this page says what is +*implemented and verified now*, equation by equation. + +!!! warning "Where the physics is (and isn't)" + Everything on this page is **structure and numerics**. Under the module's + `[VERIFY]` policy (`src/Islands/CLAUDE.md`), no physics coefficient, sign, + or normalization from the drift-kinetic literature has been assigned a + value anywhere in `src/`: every such quantity enters as a *supplied, + gated parameter* (many deliberately default to `NaN` so an un-cleared + convention poisons results rather than guessing). The human clearance + queue that un-gates the physics is `QUESTIONS.md` entries **Q2–Q4**; until + then the physics benchmarks (`benchmarks/islands/`) ship skipped by design. + +## 1. Phase space and discretization + +The Level-0 solve lives on the orbit-averaged phase space +``(x, \xi;\, y, E, \sigma)`` — radial distance from the rational surface, +helical angle, pitch ``y = \lambda B_{\max}``, energy, and the parallel-velocity +sign ``\sigma = \pm 1`` (design `03 §1`). Implemented in `Islands.PhaseSpace`: + +**Helical angle ``\xi``** — Fourier pseudo-spectral on the periodic domain. The +dense spectral derivative on an even number ``n`` of uniform nodes is + +```math +(D_1)_{jk} \;=\; \frac{(-1)^{j-k}}{2}\,\cot\!\Big(\frac{(j-k)\,h}{2}\Big), +\qquad j \neq k,\quad h = \tfrac{2\pi}{n}, +``` + +exact for bandlimited data (verified to ``6\times10^{-15}`` in the tests). + +**Radial ``x`` and pitch ``y``** — high-order finite differences on +layer-clustered grids. A uniform computational coordinate ``s \in [-1, 1]`` maps +to the physical coordinate through a monotone ``\sinh`` stretching that packs +nodes at the internal layers the drift-kinetic problem develops (`04 §1–2`): + +```math +x(s) \;=\; x_c + L\,\frac{\sinh(\beta s)}{\sinh(\beta)} +\qquad\Longrightarrow\qquad +\frac{d}{dx} = \frac{1}{x'(s)}\frac{d}{ds},\quad +\frac{d^2}{dx^2} = \frac{1}{x'(s)^2}\frac{d^2}{ds^2} - \frac{x''(s)}{x'(s)^3}\frac{d}{ds}. +``` + +The radial grid packs toward the rational surface ``x = 0``; the pitch grid +toward the trapped–passing boundary ``y_c`` — the two layers whose widths scale +as ``\nu^{1/2}`` and set the prior art's operating floor (`04 §2`). +Derivative matrices use Fornberg weights on windows of ``\mathrm{order}+d`` +points for the ``d``-th derivative, so ``D_1`` **and** ``D_2`` hold the design +order uniformly, including at boundary rows. Composite-Simpson weights on the +same nodes (pushed through the map Jacobian) give quadrature at matching order. + +**Energy ``E``** — Gauss–Laguerre nodes and weights: the Level-0 Maxwellian +weight ``\int_0^\infty f(E)\, e^{-E}\, dE = \sum_i w_i f(E_i)`` (a slowing-down +background at Level 2 changes the map, not the machinery). + +![layer-clustered grids](figures/grids_clustering.png) + +*Implementing symbols:* `PhaseSpace.FourierGrid`, `PhaseSpace.MappedFDGrid`, +`PhaseSpace.GaussGrid`, `PhaseSpace.IslandGrid`, `PhaseSpace.fd_weights`. + +## 2. The operator stack + +The unknowns are ``U = (g,\, \tilde\Phi)`` — the orbit-averaged distribution +``g(x, \xi, y, E, \sigma)`` per species and the electrostatic potential +``\tilde\Phi(x, \xi)`` — and the steady-state residual is assembled as a sum of +independent operator applications (design `03 §2`; no term inspects which +others are active, no regime branches anywhere): + +```math +R_g(U) \;=\; \sum_{\text{terms } T} T[U], +\qquad +R_\Phi(U) \;=\; M[g] - \alpha\,\tilde\Phi , +``` + +with the Level-0 term structures (each coefficient below is **supplied data**, +its physics value gated): + +| Term | Structure | Gated coefficient | +|---|---|---| +| `ParallelStreaming` | ``a_\xi\, \partial_\xi g + a_x\, \partial_x g`` | island-induced streaming frequencies | +| `MagneticDrift` | ``c_D\, \partial_\xi g`` (with the `:original`/`:improved` ``\hat L_B^{-1}`` toggle) | the precession frequency ``\hat\omega_D(y, E; \sigma)`` | +| `ExBDrift` | ``c_E \left( \partial_\xi\tilde\Phi\, \partial_x g - \partial_x\tilde\Phi\, \partial_\xi g \right)`` — the ``(x,\xi)`` Poisson bracket, the one state-nonlinear Level-0 term | the ``E\times B`` coupling | +| `PitchAngleDiffusion` | ``c\,(K g)`` along ``y`` (mimetic form, §3) | ``\hat\nu(E)`` and the pitch diffusivity profile | +| `GradientDrive` | additive source | the ``(\mathbf v_E + \mathbf v_D + \mathbf v_{\tilde\psi})\cdot\nabla F_0`` drive | +| `Quasineutrality` | ``M[g] - \alpha\tilde\Phi`` | the closure coefficient ``\alpha`` | + +Every `apply!` kernel is allocation-free (a CI regression test holds this at +**0 bytes**) and generic over the element type, so ForwardDiff dual numbers +flow through the entire stack — that is what makes the solver's Jacobian exact +(§5). + +*Implementing symbols:* `Operators.AbstractTerm`, `Operators.apply!`, +`Operators.residual!`, `Operators.IslandStack`. + +## 3. The conservative collision structure + +Bootstrap physics is unforgiving about non-conservative collision operators +(design `00`, risk register), so the pitch-angle operator is implemented in +**mimetic divergence form**. With gradient matrix ``G`` (the ``y``-grid ``D_1``), +quadrature weights ``w_q``, a supplied non-negative diffusivity profile ``P(y)`` +and measure ``w(y)``: + +```math +K \;=\; -\,W_q^{-1}\, G^{\mathsf T}\, \mathrm{diag}(P \circ w_q)\, G, +\qquad W_q = \mathrm{diag}(w \circ w_q), +``` + +which enjoys the two Level-0 conservation properties **exactly in floating +point**, not just asymptotically: + +```math +\mathbf 1^{\mathsf T} W_q K g + = -\,(G\mathbf 1)^{\mathsf T}\,\mathrm{diag}(P \circ w_q)\,(G g) = 0 +\quad\text{(particles; } G\mathbf 1 = 0\text{)}, +``` +```math +g^{\mathsf T} W_q K g + = -\,(G g)^{\mathsf T}\,\mathrm{diag}(P \circ w_q)\,(G g) \;\le\; 0 +\quad\text{(entropy sign)}. +``` + +Verified to ``10^{-14}`` (ladder **A4**). A physically-profiled ``P`` vanishes +at the pitch-domain endpoints, so zero-flux boundary behavior is built into the +operator — no artificial ``y`` boundary conditions. (This mattered in practice: +the generic ``a_y \partial_y^2`` form *without* boundary conditions is an +unstable BVP discretization under refinement; the mimetic degenerate form is +the correct structure.) + +*Implementing symbols:* `Operators.conservative_pitch_operator`, +`Operators.PitchAngleDiffusion`. + +## 4. Boundary conditions + +Far-field rows at ``|x| = L_x`` are replaced by matching conditions +``g - g_\infty`` and ``\tilde\Phi - \tilde\Phi_\infty`` (Dirichlet-type against +supplied far-field states). **Never bare Neumann** ``\partial_x g = 0``: the +prior art traced a spurious "winged" solution branch directly to Neumann +non-uniqueness (`01 §3`). The physical far field — the no-island neoclassical +solution — is gated physics, so the code takes it as supplied data; the tests +use manufactured far fields. These conditions are also what make the +first-order-in-``x`` advective solve well-posed. + +*Implementing symbols:* `Operators.FarFieldConditions`, `Operators.apply_farfield!`. + +## 5. The Newton–Krylov solve + +Decision D2: steady-state Newton–Krylov, never time-stepping and never the +sources' nested Picard loops (which, per the prior-art forensics, *never met* +their own convergence criterion in production). The pieces (design `04 §5`): + +**Exact matrix-free Jacobian.** The directional derivative comes from one dual- +number sweep of the residual — no finite-difference Jacobian, no global sparse +matrix: + +```math +J(u)\,v \;=\; \left.\frac{d}{d\varepsilon}\right|_{\varepsilon=0} F(u + \varepsilon v) +\quad\text{via ForwardDiff duals through the stack.} +``` + +**Inexact Newton with Eisenstat–Walker forcing.** Each step solves +``J\,\delta u = -F`` by GMRES only to the tolerance the outer iteration needs, + +```math +\eta_k = \gamma \left( \frac{\lVert F_k \rVert}{\lVert F_{k-1} \rVert} \right)^{2}, +``` + +with a backtracking line search on ``\lVert F \rVert``. Convergence is declared +on the norm **and** the pointwise maximum of the residual — the array-averaged +residual famously hid locally divergent regions in the prior art. + +**Physics-block preconditioning with explicit regularization.** The stiff +pitch-direction blocks are factored per pencil by SVD and truncated below +``\epsilon\,\sigma_{\max}`` — the deliberate treatment of the intrinsically +near-singular trapped–passing matching block (`04 §3`; the prior art measured +``\mathrm{rcond} \sim 10^{-16}`` there and got machine-dependent *noise, not +crashes*, under plain LU). On a collision-dominated test solve the block-Jacobi +preconditioner cuts the work by an order of magnitude: + +![preconditioner comparison](figures/preconditioner_gmres.png) + +A companion diagnostic tracks the smallest singular value of the ``y_c``-block +of the (tiny-grid, debug) dense Jacobian — ladder **A8** — so a silent +conditioning regression is *tested for*, not observed. + +*Implementing symbols:* `Solvers.newton_krylov`, `Solvers.JVPOperator`, +`Solvers.YBlockJacobi`, `Solvers.dense_jacobian`, `Verify.yc_block_sigma_min`. + +## 6. Continuation and fold detection + +Δ-surface generation, Newton globalization, and (at Level 3) the penetration +bifurcation all ride on pseudo-arclength continuation, so fold handling is in +from day one. The corrector solves the extended system + +```math +G(z) = \begin{pmatrix} F(u, p) \\ t \cdot (z - z_{\text{pred}}) \end{pmatrix} = 0, +\qquad z = (u, p), +``` + +with a secant tangent ``t`` and fold detection via sign reversal of the +tangent's parameter component ``t_p``: + +![continuation fold](figures/continuation_fold.png) + +*Implementing symbol:* `Solvers.pseudo_arclength`. + +## 7. Island geometry, the electron-closure functions, and the Δ moments + +The island flux-surface label is the pinned module convention (half-width +``w``): + +```math +\Omega(x, \xi) = \frac{2x^2}{w^2} - \cos\xi, +\qquad \Omega = -1 \text{ at the O-point},\quad \Omega = +1 \text{ at the separatrix}, +``` + +with the flux-surface average +``\langle f \rangle_\Omega = \oint f\,(\Omega + \cos\xi)^{-1/2} d\xi \,/\, +\oint (\Omega + \cos\xi)^{-1/2} d\xi`` — a *diagnostic* only, never a solve +coordinate (Decision D1). The flattened-electron closure geometry is +implemented as structure with a supplied amplitude: + +```math +Q(\Omega) = \frac{1}{2\pi} \oint \sqrt{\Omega + \cos\xi}\; d\xi, +\qquad +h(\Omega) = \Theta(\Omega - 1)\; C \int_1^{\Omega} \frac{d\Omega'}{Q(\Omega')}, +``` + +so ``h`` is exactly flat inside the separatrix. Because ``h'(\Omega) = C/Q``, +the chain rule gives the **coefficient-free consistency identity** (ladder +**A7** — the unit target that historically caught inherited bugs in this +lineage): + +```math +\Big\langle \frac{\partial^2 h}{\partial x^2} \Big\rangle_{\!\Omega} + = \frac{4}{w^2}\left[ h''(\Omega)\, \frac{Q}{Q'} \; +\; h'(\Omega) \right] + \;\equiv\; 0 , +``` + +verified to ``10^{-16}`` for arbitrary amplitude ``C``: + +![Q and h profiles](figures/hQ_profiles.png) + +The output moments are the two Ampère projections of the species-summed +parallel current ``\bar J_\parallel = \sum_j Z_j \int W_j\, g_j`` (`01 §4`): + +```math +\Delta_{\cos} = C_{\cos} \int dx \oint d\xi\; \bar J_\parallel \cos\xi, +\qquad +\Delta_{\sin} = C_{\sin} \int dx \oint d\xi\; \bar J_\parallel \sin\xi, +``` + +where the ``\xi``-projection is spectrally exact on the periodic grid and the +prefactors ``C_{\cos}, C_{\sin}`` (physically ``\mp\mu_0 R / 2\tilde\psi``) are +**required, gated arguments** — ``\tilde\psi`` carries an open `[VERIFY]` and +the sin normalization is `[DERIVED]`-unpinned (QUESTIONS Q4). The parity +structure ``\Delta_{\cos}`` even / ``\Delta_{\sin}`` odd under ``\xi``-reflection +is verified exactly (ladder **A3**). + +*Implementing symbols:* `Moments.parallel_current!`, `Moments.delta_moments`, +`Moments.omega_average`, `Fields.Q_omega`, `Fields.h_profile`, +`Fields.flat_average_d2h_dx2`. + +## 8. Verification evidence (the A-ladder, all green in CI) + +The manufactured-solution ladder verifies discretization order and the AD +plumbing simultaneously — per operator, for the assembled residual, and through +a full converged Newton solve forced by the analytic source: + +![MMS convergence](figures/mms_convergence.png) + +| Gate | Statement | Result | +|---|---|---| +| A1 | per-operator + assembled MMS at design order | 4th order (``\xi`` spectral to ``10^{-15}``) | +| A1-solve | converged Newton solve recovers the manufactured state | observed order **3.98** | +| A2 | AD Jacobian–vector product vs. central finite differences | agree to ``\sim 10^{-9}`` | +| A3 | ``\Delta_{\cos}`` even / ``\Delta_{\sin}`` odd parity | exact | +| A4 | particle conservation + entropy sign of the collision operator | exact (``10^{-14}`` / definite) | +| A5 | zero-drive null: ``g \equiv 0 \Rightarrow R = 0`` | **exactly** machine zero | +| A7 | ``\langle \partial^2 h / \partial x^2 \rangle_\Omega = 0`` | ``10^{-16}`` | +| A8 | ``y_c``-block ``\sigma_{\min}`` monitor + singular detection | active | +| — | allocation regression on every hot kernel | 0 bytes | + +Everything above regenerates from one pinned script +(`benchmarks/islands/figures/make_structural_figures.jl`) and runs in the test +suite (`test/runtests_islands_{grids,operators,solve}.jl`); the +`islands_l0_structural` regression case tracks the headline numbers across +commits. + +## 9. What comes next + +The physics story — the drift-kinetic coefficients, the York threshold triangle +(DK-NTM ``8.73\,\rho_{bi}`` → RDK-NTM ``1.46\,\rho_{bi}``), the +``\Delta_{\text{pol}}(\omega_E)`` sign-reversal curve — is written as the +[Paper I figure contract](papers/paper-1/OUTLINE.md) and un-gates claim by +claim as the human clearance queue (`QUESTIONS.md` Q2–Q4) is worked through. +The [design documents](design/00-roadmap.md) hold the full eight-milestone +program. diff --git a/docs/src/islands/papers/paper-1/OUTLINE.md b/docs/src/islands/papers/paper-1/OUTLINE.md new file mode 100644 index 00000000..bfc3704d --- /dev/null +++ b/docs/src/islands/papers/paper-1/OUTLINE.md @@ -0,0 +1,44 @@ +# Paper I — OUTLINE (the Level-0 figure contract) + +> Created at level *start* per design doc `07 §3`: claims → figures → ladder +> IDs. This outline is the figure contract — agents implementing benchmarks +> know which figures are paper figures from day one. Every claim must be backed +> by a ladder ID; claims lacking one are flagged. Status reflects the +> `[VERIFY]` gating of `docs/src/islands/QUESTIONS.md` (Q2–Q4) — **no +> submission until every tag in the paper's equation set is cleared**. + +**Working title:** Formulation and verification of a generalized drift-kinetic +solver for magnetic-island stability (Islands, Level 0). + +**Indicative venue:** Phys. Plasmas (methods paper). **Gate:** Level 0 +(design `00`), i.e. ladder A + B1/B2/B4/B5a–c green with convergence artifacts. + +## Claims and figures + +| # | Claim | Figure(s) | Ladder ID(s) | Status | +|---|---|---|---|---| +| C1 | The `(x, ξ; y, E, σ)` discretization converges at design order: spectral in `ξ`, 4th-order FD on layer-clustered `x`/`y` grids, per operator and for the assembled solve | F1: MMS convergence panels (per-operator + assembled residual + assembled solve error vs. resolution) | A1, A2 | **green** (M1/M2; CI artifacts) | +| C2 | The steady-state Newton–Krylov solve is exact-Jacobian (AD), globally convergent from zero states, and its conditioning is monitored at the trapped–passing boundary (no silent-noise regime, contra L23 §4.2) | F2: Newton/GMRES convergence histories with and without the physics-block preconditioner; `σ_min(y_c)` track | A5, A8 (+ solver gates) | **green** (M2) | +| C3 | The discretized collision operator conserves particles exactly and has definite entropy sign — bootstrap-relevant structure holds discretely, not just asymptotically | F3: conservation/entropy residuals vs. resolution and profile | A4 | **green** (L0 parts, M2) | +| C4 | In the no-island limit the solver reproduces standard local neoclassics (bootstrap current vs. Sauter/NEO) — the strongest global check of the velocity-space discretization | F4: `J_bs` vs. `ν_★` against NEO/Sauter | B1 | **gated** — needs cleared L0 coefficient set (Q2, Q3) + external NEO runs | +| C5 | At large `w` the solver recovers the analytic MRE limits: `Δ_bs + Δ_cur ∝ 1/w` (WCHH96 Eq. 85, frame-mapped) and `Δ_pol ∝ 1/w³` | F5: `Δ` channels vs. `w` with analytic asymptotes | B2 | **gated** — Q2, Q3, Q4 (WCHH96 acquisition; `ψ̃` [VERIFY]) | +| C6 | The three-code threshold triangle is reproduced in its exact configurations: DK-NTM `w_c ≃ 2.76 ρ_θi ≡ 8.73 ρ_bi` (`:original` drift model), RDK-NTM `w_c ≈ 0.45 ρ_θi ≡ 1.46 ρ_bi` (`:improved`), kokuchou's finite-`ν_★` surface `w_c(ρ̂_θi, ν_★)` — and the 8.73 → 1.46 shift is a single drift-model toggle (E1) | F6: threshold `w_c` vs. configuration (the triangle); F7: the E1 toggle scan | B5a, B5b, B5c (E1) | **gated** — Q2 (D7/D8), Q3, Q4 (B5a collisionality) | +| C7 | The polarization contribution is frame-pinned: `Δ_pol ∝ ω_E²` away from zero with sign reversal at `ω_E ≈ −0.89 ω_dia,e`, insensitive to `w/ρ_θi`; single-`ω_E` `Δ` values are misleading (surfaces over `(w, ω_E)` are the deliverable) | F8: `Δ_pol(ω_E)` reversal curve vs. D23b Fig. 8; F9: `Δ(w, ω_E)` surface | B4 | **gated** — Q2, Q3 (frame-convention signs) | +| C8 (candidate headline) | Resolution of L23's open question: the stabilizing *electron* `Δ_pol` at `ω_E = 0` — reproduced or refuted by the `ω_E` scan with kinetic vs. flattened electrons (E4) | F10: electron-channel `Δ_pol(ω_E)` decomposition | B4 + E4 | **gated** — Q2, Q3; kinetic-electron toggle is M2+/E4 work | + +## Verification-artifact rules (docs/05 reporting) + +- Every figure names its configuration (docs/03 §2) and git SHA; no benchmark + "passes" on a single grid — convergence + tolerance archived with the result. +- Threshold numbers are reported as **half-widths** with both `ρ_θi` and + `ρ_bi = ε^{1/2} ρ_θi` stated at the run's `ε` (docs/05 rule 5). +- Disagreements with published targets are triaged per the standing rule + (docs/05: our bug / their approximation / their published-equation error / + transcription error) with `[VERIFY]` resolution logged first. + +## Dependencies for un-gating (the human clearance queue) + +`QUESTIONS.md` **Q2** (ratify D7/D8), **Q3** (clear the L0 `[CHECKED]` +coefficient set), **Q4** (resolve the `ψ̃` and B5a-collisionality `[VERIFY]`s; +acquire WCHH96 + Park 2022). C1–C3 are already green as CI artifacts; C4–C8 +un-gate in that order of effort once the queue clears. diff --git a/regression-harness/cases/islands_l0_structural.toml b/regression-harness/cases/islands_l0_structural.toml new file mode 100644 index 00000000..091b4949 --- /dev/null +++ b/regression-harness/cases/islands_l0_structural.toml @@ -0,0 +1,60 @@ +# Regression case: Islands L0 solve machinery, structural (pre-physics) quantities. +# A "computed" case (kind = "computed", no example_dir) tracking the M2 structural +# gates: the solve-level MMS error and solver iteration counts (Newton-Krylov on the +# manufactured advective stack at nx = 17), the coefficient-free flattened-electron +# closure identity |_Omega| at Omega = 2, and the y_c-block smallest singular +# value (the L23 SS4.2 silent-noise tripwire). The physics B-ladder stays skipped until +# [VERIFY] clearance (docs/src/islands/QUESTIONS.md Q2-Q4); these are the numbers that +# silently move if the discretization, solver, or quadratures drift. +[case] +name = "islands_l0_structural" +description = "Islands L0 structural: solve-MMS error, solver iterations, A7 identity, y_c sigma_min" +kind = "computed" + +[quantities.solve_mms_err] +h5path = "islands/solve_mms_err" +type = "real_scalar" +extract = "value" +label = "solve-MMS max error (nx=17)" +noise_threshold = 1e-8 +order = 10 + +[quantities.newton_iters] +h5path = "islands/newton_iters" +type = "real_scalar" +extract = "value" +label = "Newton iterations" +noise_threshold = 0.0 +order = 11 + +[quantities.gmres_iters] +h5path = "islands/gmres_iters" +type = "real_scalar" +extract = "value" +label = "GMRES iterations (total)" +noise_threshold = 0.0 +order = 12 + +[quantities.a7_identity] +h5path = "islands/a7_identity" +type = "real_scalar" +extract = "value" +label = "|_Omega| (A7)" +noise_threshold = 1e-12 +order = 13 + +[quantities.yc_sigma_min] +h5path = "islands/yc_sigma_min" +type = "real_scalar" +extract = "value" +label = "y_c block sigma_min (A8)" +noise_threshold = 1e-10 +order = 14 + +[quantities.runtime] +h5path = "" +type = "runtime" +extract = "value" +label = "Runtime (s)" +noise_threshold = 0.0 +order = 90 diff --git a/regression-harness/src/runner.jl b/regression-harness/src/runner.jl index 8705cb4d..d6aa9b21 100644 --- a/regression-harness/src/runner.jl +++ b/regression-harness/src/runner.jl @@ -2,7 +2,9 @@ Runner: orchestrates checking out commits, running GPEC, and extracting results. """ -"""Read the GPEC-only runtime from the timing file written by the subprocess.""" +""" +Read the GPEC-only runtime from the timing file written by the subprocess. +""" function _read_timing_file(path::String)::Float64 if isfile(path) return parse(Float64, strip(read(path, String))) @@ -29,7 +31,7 @@ function _materialize_rundir(example_path::String, overrides::Dict{String,Any}) for (dotted, val) in overrides ks = split(dotted, ".") d = cfg - for k in ks[1:(end - 1)] + for k in ks[1:(end-1)] d = get!(d, k, Dict{String,Any}()) end d[ks[end]] = val @@ -109,28 +111,60 @@ open(ARGS[2], "w") do f end """ +# Islands L0 structural regression: the solve-level MMS (Newton-Krylov recovers +# the manufactured state), the coefficient-free A7 closure identity, and the A8 +# y_c-block conditioning monitor. Tracks the structural (pre-physics) numbers +# that would silently move if the discretization, solver, or quadratures drift; +# the physics B-ladder stays skipped until [VERIFY] clearance (QUESTIONS Q2-Q4). +const COMPUTED_ISLANDS_SCRIPT_TEMPLATE = """ +using Pkg +%INSTANTIATE% +using GeneralizedPerturbedEquilibrium +using HDF5 +const Isl = GeneralizedPerturbedEquilibrium.Islands +t_start = time() +r = Isl.Verify.solve_mms(17) +a7 = Isl.Fields.flat_average_d2h_dx2(2.0, 1.0) +grid = Isl.PhaseSpace.IslandGrid(nx=7, nxi=8, ny=7, nE=2, halfwidth_x=6.0, clustering_x=1.0, + y_max=4.0, y_c=1.0, clustering_y=0.8, order=4) +setup = Isl.Verify.zero_drive_setup(grid) +J = Isl.Solvers.dense_jacobian(setup.f, zeros(setup.N)) +mon = Isl.Verify.yc_block_sigma_min(J, grid) +elapsed = time() - t_start +h5open(ARGS[1], "w") do fid + fid["islands/solve_mms_err"] = r.err + fid["islands/newton_iters"] = Float64(r.iterations) + fid["islands/gmres_iters"] = Float64(r.gmres_iters) + fid["islands/a7_identity"] = abs(a7) + fid["islands/yc_sigma_min"] = mon.sigma_min +end +open(ARGS[2], "w") do f + println(f, elapsed) +end +""" + """ Run GPEC for a single commit/ref and case. Dispatches to run_local for the working tree or run_at_commit for a git ref. """ function run_commit(db::SQLite.DB, commit_hash::String, ref_name::String, - case_spec::CaseSpec, repo_root::String; - force::Bool=false, verbose::Bool=false, - no_instantiate::Bool=false) + case_spec::CaseSpec, repo_root::String; + force::Bool=false, verbose::Bool=false, + no_instantiate::Bool=false) if case_spec.kind == "computed" if commit_hash == LOCAL_REF return run_computed_local(db, case_spec, repo_root; - verbose=verbose, no_instantiate=no_instantiate) + verbose=verbose, no_instantiate=no_instantiate) end return run_computed_at_commit(db, commit_hash, ref_name, case_spec, repo_root; - force=force, verbose=verbose, no_instantiate=no_instantiate) + force=force, verbose=verbose, no_instantiate=no_instantiate) end if commit_hash == LOCAL_REF return run_local(db, case_spec, repo_root; - force=force, verbose=verbose, no_instantiate=no_instantiate) + force=force, verbose=verbose, no_instantiate=no_instantiate) end return run_at_commit(db, commit_hash, ref_name, case_spec, repo_root; - force=force, verbose=verbose, no_instantiate=no_instantiate) + force=force, verbose=verbose, no_instantiate=no_instantiate) end """ @@ -141,6 +175,8 @@ function _computed_script_template(case_spec::CaseSpec) return COMPUTED_GGJ_SCRIPT_TEMPLATE elseif case_spec.name == "efit_fixedbdy_separatrix" return COMPUTED_SEPARATRIX_SCRIPT_TEMPLATE + elseif case_spec.name == "islands_l0_structural" + return COMPUTED_ISLANDS_SCRIPT_TEMPLATE end error("No computed-script template registered for case '$(case_spec.name)'") end @@ -153,11 +189,11 @@ import GeneralizedPerturbedEquilibrium), reads the resulting tempfile h5 with so callers can handle store_failed_run uniformly. """ function _execute_computed(case_spec::CaseSpec, project_root::String; - verbose::Bool, no_instantiate::Bool, - stderr_buf::IO) + verbose::Bool, no_instantiate::Bool, + stderr_buf::IO) instantiate_line = no_instantiate ? "" : "Pkg.instantiate()" script_content = replace(_computed_script_template(case_spec), - "%INSTANTIATE%" => instantiate_line) + "%INSTANTIATE%" => instantiate_line) tmpscript = tempname() * ".jl" h5path = tempname() * ".h5" timingfile = tempname() * ".timing" @@ -166,8 +202,8 @@ function _execute_computed(case_spec::CaseSpec, project_root::String; if verbose run(pipeline(`julia --project=$project_root $tmpscript $h5path $timingfile`)) else - run(pipeline(`julia --project=$project_root $tmpscript $h5path $timingfile`, - stdout=devnull, stderr=stderr_buf)) + run(pipeline(`julia --project=$project_root $tmpscript $h5path $timingfile`; + stdout=devnull, stderr=stderr_buf)) end runtime_s = _read_timing_file(timingfile) if !isfile(h5path) @@ -186,18 +222,18 @@ end Run a kind="computed" case against the working tree. """ function run_computed_local(db::SQLite.DB, case_spec::CaseSpec, repo_root::String; - verbose::Bool=false, no_instantiate::Bool=false) + verbose::Bool=false, no_instantiate::Bool=false) delete_cached(db, LOCAL_REF, case_spec.name) date = Dates.format(Dates.now(), "yyyy-mm-ddTHH:MM:SS") @info "Running: $(case_spec.name) @ local (working tree, computed)" stderr_buf = IOBuffer() try extracted, runtime_s = _execute_computed(case_spec, repo_root; - verbose=verbose, - no_instantiate=no_instantiate, - stderr_buf=stderr_buf) + verbose=verbose, + no_instantiate=no_instantiate, + stderr_buf=stderr_buf) store_run(db, LOCAL_REF, "local", date, "working tree", case_spec.name, - runtime_s, extracted) + runtime_s, extracted) @info " Completed in $(round(runtime_s, digits=3))s — $(length(extracted)) quantities extracted" catch e err_msg = if e isa ProcessFailedException @@ -212,7 +248,7 @@ function run_computed_local(db::SQLite.DB, case_spec::CaseSpec, repo_root::Strin err_msg_short = length(err_msg) > 2000 ? "..." * last(err_msg, 2000) : err_msg @warn "Run failed (local computed): $(first(err_msg_short, 200))" store_failed_run(db, LOCAL_REF, "local", date, "working tree", case_spec.name, - err_msg_short) + err_msg_short) end end @@ -220,9 +256,9 @@ end Run a kind="computed" case at a specific git commit via worktree. """ function run_computed_at_commit(db::SQLite.DB, commit_hash::String, ref_name::String, - case_spec::CaseSpec, repo_root::String; - force::Bool=false, verbose::Bool=false, - no_instantiate::Bool=false) + case_spec::CaseSpec, repo_root::String; + force::Bool=false, verbose::Bool=false, + no_instantiate::Bool=false) if !force && is_cached(db, commit_hash, case_spec.name) info = get_run_info(db, commit_hash, case_spec.name) if info !== nothing @@ -243,11 +279,11 @@ function run_computed_at_commit(db::SQLite.DB, commit_hash::String, ref_name::St try worktree_path = create_worktree(commit_hash, repo_root) extracted, runtime_s = _execute_computed(case_spec, worktree_path; - verbose=verbose, - no_instantiate=no_instantiate, - stderr_buf=stderr_buf) + verbose=verbose, + no_instantiate=no_instantiate, + stderr_buf=stderr_buf) store_run(db, commit_hash, commit_info.short, commit_info.date, - commit_info.msg, case_spec.name, runtime_s, extracted) + commit_info.msg, case_spec.name, runtime_s, extracted) @info " Completed in $(round(runtime_s, digits=3))s — $(length(extracted)) quantities extracted" catch e err_msg = if e isa ProcessFailedException @@ -262,7 +298,7 @@ function run_computed_at_commit(db::SQLite.DB, commit_hash::String, ref_name::St err_msg_short = length(err_msg) > 2000 ? "..." * last(err_msg, 2000) : err_msg @warn "Run failed (computed) for $(commit_info.short): $(first(err_msg_short, 200))" store_failed_run(db, commit_hash, commit_info.short, commit_info.date, - commit_info.msg, case_spec.name, err_msg_short) + commit_info.msg, case_spec.name, err_msg_short) finally if worktree_path !== nothing remove_worktree(worktree_path, repo_root) @@ -275,8 +311,8 @@ Run GPEC in the current working tree (uncommitted changes included). Always re-runs (local results are never cached since the working tree is mutable). """ function run_local(db::SQLite.DB, case_spec::CaseSpec, repo_root::String; - force::Bool=false, verbose::Bool=false, - no_instantiate::Bool=false) + force::Bool=false, verbose::Bool=false, + no_instantiate::Bool=false) # Always delete previous local results and re-run delete_cached(db, LOCAL_REF, case_spec.name) @@ -287,7 +323,7 @@ function run_local(db::SQLite.DB, case_spec::CaseSpec, repo_root::String; if !isdir(example_path) @warn "Example directory not found: $(case_spec.example_dir)" store_failed_run(db, LOCAL_REF, "local", date, "working tree", case_spec.name, - "Example directory not found: $(case_spec.example_dir)") + "Example directory not found: $(case_spec.example_dir)") return end @@ -309,8 +345,8 @@ function run_local(db::SQLite.DB, case_spec::CaseSpec, repo_root::String; if verbose run(pipeline(`julia --project=$repo_root $tmpscript $rundir $timingfile`)) else - run(pipeline(`julia --project=$repo_root $tmpscript $rundir $timingfile`, - stdout=devnull, stderr=stderr_buf)) + run(pipeline(`julia --project=$repo_root $tmpscript $rundir $timingfile`; + stdout=devnull, stderr=stderr_buf)) end runtime_s = _read_timing_file(timingfile) @@ -318,13 +354,13 @@ function run_local(db::SQLite.DB, case_spec::CaseSpec, repo_root::String; if !isfile(h5path) @warn "gpec.h5 not produced" store_failed_run(db, LOCAL_REF, "local", date, "working tree", case_spec.name, - "gpec.h5 not produced after successful run") + "gpec.h5 not produced after successful run") return end extracted = extract_quantities(h5path, case_spec.quantities, runtime_s) store_run(db, LOCAL_REF, "local", date, "working tree", case_spec.name, - runtime_s, extracted) + runtime_s, extracted) @info " Completed in $(round(runtime_s, digits=1))s — $(length(extracted)) quantities extracted" @@ -341,7 +377,7 @@ function run_local(db::SQLite.DB, case_spec::CaseSpec, repo_root::String; err_msg_short = length(err_msg) > 2000 ? "..." * last(err_msg, 2000) : err_msg @warn "Run failed (local): $(first(err_msg_short, 200))" store_failed_run(db, LOCAL_REF, "local", date, "working tree", case_spec.name, - err_msg_short) + err_msg_short) finally if tmpscript !== nothing rm(tmpscript; force=true) @@ -360,9 +396,9 @@ Run GPEC for a specific git commit via worktree. Stores results in the database. Skips if already cached (unless force=true). """ function run_at_commit(db::SQLite.DB, commit_hash::String, ref_name::String, - case_spec::CaseSpec, repo_root::String; - force::Bool=false, verbose::Bool=false, - no_instantiate::Bool=false) + case_spec::CaseSpec, repo_root::String; + force::Bool=false, verbose::Bool=false, + no_instantiate::Bool=false) # Check cache if !force && is_cached(db, commit_hash, case_spec.name) info = get_run_info(db, commit_hash, case_spec.name) @@ -398,8 +434,8 @@ function run_at_commit(db::SQLite.DB, commit_hash::String, ref_name::String, if !isdir(example_path) @warn "Example directory not found at commit $(commit_info.short): $(case_spec.example_dir)" store_failed_run(db, commit_hash, commit_info.short, commit_info.date, - commit_info.msg, case_spec.name, - "Example directory not found: $(case_spec.example_dir)") + commit_info.msg, case_spec.name, + "Example directory not found: $(case_spec.example_dir)") return end @@ -418,8 +454,8 @@ function run_at_commit(db::SQLite.DB, commit_hash::String, ref_name::String, if verbose run(pipeline(`julia --project=$project_root $tmpscript $rundir $timingfile`)) else - run(pipeline(`julia --project=$project_root $tmpscript $rundir $timingfile`, - stdout=devnull, stderr=stderr_buf)) + run(pipeline(`julia --project=$project_root $tmpscript $rundir $timingfile`; + stdout=devnull, stderr=stderr_buf)) end runtime_s = _read_timing_file(timingfile) @@ -428,8 +464,8 @@ function run_at_commit(db::SQLite.DB, commit_hash::String, ref_name::String, if !isfile(h5path) @warn "gpec.h5 not produced at $(commit_info.short)" store_failed_run(db, commit_hash, commit_info.short, commit_info.date, - commit_info.msg, case_spec.name, - "gpec.h5 not produced after successful run") + commit_info.msg, case_spec.name, + "gpec.h5 not produced after successful run") return end @@ -438,7 +474,7 @@ function run_at_commit(db::SQLite.DB, commit_hash::String, ref_name::String, # Store in database store_run(db, commit_hash, commit_info.short, commit_info.date, - commit_info.msg, case_spec.name, runtime_s, extracted) + commit_info.msg, case_spec.name, runtime_s, extracted) @info " Completed in $(round(runtime_s, digits=1))s — $(length(extracted)) quantities extracted" @@ -456,7 +492,7 @@ function run_at_commit(db::SQLite.DB, commit_hash::String, ref_name::String, err_msg_short = length(err_msg) > 2000 ? "..." * last(err_msg, 2000) : err_msg @warn "Run failed for $(commit_info.short): $(first(err_msg_short, 200))" store_failed_run(db, commit_hash, commit_info.short, commit_info.date, - commit_info.msg, case_spec.name, err_msg_short) + commit_info.msg, case_spec.name, err_msg_short) finally # Clean up if tmpscript !== nothing diff --git a/src/Islands/Islands.jl b/src/Islands/Islands.jl index f6693c97..9fec2f9e 100644 --- a/src/Islands/Islands.jl +++ b/src/Islands/Islands.jl @@ -19,20 +19,26 @@ milestone M1 proceeds. module Islands # Submodule layout follows docs/src/islands/design/03-architecture.md §1. -# M1 lands the discretization + operator-stack skeleton + verification harness; -# the remaining submodules (species, frames, fields, moments, solvers) land as -# later milestones proceed. +# M1 landed the discretization + operator stack + MMS/AD harness; M2 lands the +# L0 solve machinery (species, frames, solve, fields, moments) as gated +# structure. Remaining design dirs (geometry/, closures/, io/) arrive with the +# milestones that need them (L2/L4/M2-io). include("phasespace/PhaseSpace.jl") # grids (x, ξ, λ→y, E, σ), layer-clustered maps +include("species/Species.jl") # Species, backgrounds, roles (D3) +include("frames/Frames.jl") # THE frequency/frame conversion module (gated forms) include("operators/Operators.jl") # the AbstractTerm stack + residual assembly -include("verify/Verify.jl") # MMS + AD-vs-FD JVP harness (ladder A1, A2) -# include("species/Species.jl") # Species, backgrounds, roles (M2+) -# include("frames/Frames.jl") # THE frequency/frame conversion module (M2+) -# include("fields/Fields.jl") # Φ̃ quasineutrality (A_∥ Ampère at L3) (M2+/L3) -# include("moments/Moments.jl") # Δ_cos, Δ_sin, profiles, channel decomps (M2) -# include("solvers/Solvers.jl") # Newton–Krylov, continuation, trace pass (M2) +include("fields/Fields.jl") # quasineutrality closure structure, h(Ω)/Q(Ω) +include("moments/Moments.jl") # J̄_∥, Δ_cos/Δ_sin projections, ⟨·⟩_Ω diagnostics +include("solvers/Solvers.jl") # Newton–Krylov, preconditioner, continuation +include("verify/Verify.jl") # MMS + AD-vs-FD JVP harness, y_c monitor import .PhaseSpace +import .SpeciesLists +import .Frames import .Operators +import .Fields +import .Moments +import .Solvers import .Verify end # module Islands diff --git a/src/Islands/fields/Fields.jl b/src/Islands/fields/Fields.jl new file mode 100644 index 00000000..a050d988 --- /dev/null +++ b/src/Islands/fields/Fields.jl @@ -0,0 +1,143 @@ +""" + Fields + +The Level-0 field-equation layer (design `01 §3`, `03 §1`): quasineutrality is +the only field equation at Level 0 (Ampère arrives at Level 3). The residual +term itself lives in `Operators.Quasineutrality`; this module provides the +closure *structure* around it: + + - the island-geometry functions `Q(Ω)`, `h(Ω)` of the flattened-electron + (WCHH96-class) closure — implemented as **structure with a supplied + prefactor** (the physics prefactor `w_ψ/2√2` is `[CHECKED]`-uncleared, + QUESTIONS Q3); + - the coefficient-free consistency identity `⟨∂²h/∂x²⟩_Ω = 0` (ladder A7 — + kokuchou's unit set, which caught inherited DK-NTM bugs); + - [`ElectronClosure`](@ref) — the gated constant set of the closure, NaN- + poisoned until human-cleared. +""" +module Fields + +import QuadGK + +export Q_omega, dQ_domega, h_profile, dh_domega, d2h_domega2 +export flat_average_d2h_dx2, ElectronClosure, is_cleared + +""" + Q_omega(Ω; rtol=1e-10) + +`Q(Ω) = (1/2π) ∮ √(Ω + cos ξ) dξ` over the region where the radicand is +positive (`01 §2.4` structure). For `Ω > 1` the integral covers the full +period; for `|Ω| < 1` it runs between the turning points `cos ξ_b = −Ω`. +""" +function Q_omega(Ω; rtol::Real=1e-10) + if Ω > 1 + val, _ = QuadGK.quadgk(ξ -> sqrt(Ω + cos(ξ)), 0.0, 2π; rtol=rtol) + return val / (2π) + elseif -1 < Ω <= 1 + ξb = acos(-Ω) + val, _ = QuadGK.quadgk(ξ -> sqrt(max(Ω + cos(ξ), 0.0)), -ξb, ξb; rtol=rtol) + return val / (2π) + else + throw(DomainError(Ω, "Q_omega needs Ω > −1")) + end +end + +""" + dQ_domega(Ω; rtol=1e-10) + +`dQ/dΩ = (1/4π) ∮ (Ω + cos ξ)^{−1/2} dξ` (differentiating `Q_omega` under the +integral; the inside-separatrix endpoint singularity is integrable). +""" +function dQ_domega(Ω; rtol::Real=1e-10) + if Ω > 1 + val, _ = QuadGK.quadgk(ξ -> 1.0 / sqrt(Ω + cos(ξ)), 0.0, 2π; rtol=rtol) + return val / (4π) + elseif -1 < Ω < 1 + ξb = acos(-Ω) + val, _ = QuadGK.quadgk(ξ -> 1.0 / sqrt(Ω + cos(ξ)), -ξb, ξb; rtol=rtol) + return val / (4π) + else + throw(DomainError(Ω, "dQ_domega needs Ω ∈ (−1, 1) ∪ (1, ∞)")) + end +end + +""" + h_profile(Ω; prefactor, rtol=1e-10) + +The flattened-electron profile function *structure* `h(Ω) = Θ(Ω − 1) · prefactor · ∫₁^Ω dΩ′/Q(Ω′)` (`01 §2.4`): exactly flat inside the separatrix, +`→ x` far outside. The physics prefactor (`w_ψ/2√2` in the sources) is +`[CHECKED]`-uncleared (QUESTIONS Q3) and therefore **supplied**; the A7 +identity below is prefactor-independent. +""" +function h_profile(Ω; prefactor, rtol::Real=1e-10) + Ω <= 1 && return zero(float(Ω)) + val, _ = QuadGK.quadgk(Ωp -> 1.0 / Q_omega(Ωp; rtol=rtol), 1.0, Ω; rtol=rtol) + return prefactor * val +end + +""" + dh_domega(Ω; prefactor, rtol=1e-10) + +`h′(Ω) = prefactor / Q(Ω)` for `Ω > 1`, zero inside. +""" +dh_domega(Ω; prefactor, rtol::Real=1e-10) = Ω > 1 ? prefactor / Q_omega(Ω; rtol=rtol) : zero(float(Ω)) + +""" + d2h_domega2(Ω; prefactor, rtol=1e-10) + +`h″(Ω) = −prefactor · Q′(Ω)/Q(Ω)²` for `Ω > 1`, zero inside. +""" +function d2h_domega2(Ω; prefactor, rtol::Real=1e-10) + Ω <= 1 && return zero(float(Ω)) + Q = Q_omega(Ω; rtol=rtol) + return -prefactor * dQ_domega(Ω; rtol=rtol) / Q^2 +end + +""" + flat_average_d2h_dx2(Ω, w; prefactor=1.0, rtol=1e-10) + +The ladder-A7 identity integrand: `⟨∂²h/∂x²⟩_Ω` assembled from the chain rule +`∂²h/∂x² = h″(Ω)(∂Ω/∂x)² + h′(Ω) ∂²Ω/∂x²` with `∂Ω/∂x = 4x/w²` on the surface +(`(∂Ω/∂x)² = 8(Ω + cos ξ)/w²`), averaged with the `(Ω + cos ξ)^{−1/2}` weight. +Analytically **exactly zero** for any prefactor — a coefficient-free +consistency check of the `Q`, `h` quadratures and the `⟨·⟩_Ω` machinery +(L23 Eq. 4.1.1-class unit target). Valid for `Ω > 1`. +""" +function flat_average_d2h_dx2(Ω, w; prefactor::Real=1.0, rtol::Real=1e-10) + Ω > 1 || throw(DomainError(Ω, "the A7 identity applies outside the separatrix (Ω > 1)")) + hp = dh_domega(Ω; prefactor=prefactor, rtol=rtol) + hpp = d2h_domega2(Ω; prefactor=prefactor, rtol=rtol) + num, _ = QuadGK.quadgk(ξ -> (hpp * 8 * (Ω + cos(ξ)) / w^2 + hp * 4 / w^2) / sqrt(Ω + cos(ξ)), 0.0, 2π; rtol=rtol) + den, _ = QuadGK.quadgk(ξ -> 1.0 / sqrt(Ω + cos(ξ)), 0.0, 2π; rtol=rtol) + return num / den +end + +""" + ElectronClosure(; k_HS=NaN, f_p=NaN, h_prefactor=NaN, C_phi=NaN) + +The gated constant set of the Level-0 flattened-electron closure (`01 §2.4`, +`§3`), all `[CHECKED]`-uncleared (QUESTIONS Q3) and **NaN-defaulted** so an +un-cleared closure poisons results instead of committing to a guess: + +## Fields + + - `k_HS` — the Hirshman–Sigmar flow coefficient (sources: `≃ −1.173`). + - `f_p` — the passing fraction form (sources: `≃ 1 − 1.46√ε`). + - `h_prefactor` — the `h(Ω)` amplitude (sources: `w_ψ/2√2`). + - `C_phi` — the quasineutrality closure coefficient (sources: `1/2L̂_{n0}`). +""" +Base.@kwdef struct ElectronClosure + k_HS::Float64 = NaN + f_p::Float64 = NaN + h_prefactor::Float64 = NaN + C_phi::Float64 = NaN +end + +""" + is_cleared(ec::ElectronClosure) + +`true` only when every gated closure constant has been assigned (no NaN). +""" +is_cleared(ec::ElectronClosure) = !(isnan(ec.k_HS) || isnan(ec.f_p) || isnan(ec.h_prefactor) || isnan(ec.C_phi)) + +end # module Fields diff --git a/src/Islands/frames/Frames.jl b/src/Islands/frames/Frames.jl new file mode 100644 index 00000000..079167d9 --- /dev/null +++ b/src/Islands/frames/Frames.jl @@ -0,0 +1,119 @@ +""" + Frames + +THE frequency/frame conversion module (design `01 §5`, module CLAUDE.md): no +other part of Islands may contain an ω sign convention. The polarization-current +sign disputes in the literature are largely frame disputes; this module owns the +conversions so they cannot be reproduced inconsistently elsewhere. + +**Milestone-M2 status: forms only, signs gated.** The frame identities of +`01 §5` are `[CHECKED: Diss19 pp. 46–48]` but not human-cleared (QUESTIONS Q3), +so every sign/normalization here is a `FrameConvention` field with a **NaN +default**: using an un-cleared convention poisons every downstream number +instead of silently committing to a guess. Only the mechanical, sign-free +bookkeeping (`frame_shift`, round-trips, parameter validation) is active. +""" +module Frames + +export Level0Parameters, FrameConvention +export frame_shift, omega_dia_form, effective_dlnn_form, is_cleared + +""" + Level0Parameters(; w_hat, omega_E_hat, epsilon, inv_Lq_hat, q_s, tau, nu_star) + +The Level-0 input parameter vector `p` (`01 §5`). Species-resolved quantities +(gradients, `η_j`, backgrounds, roles) live on the species list; this struct +carries the scalars. + +## Fields + + - `w_hat` — island **half**-width `w/ρ_θi` (half-width convention pinned + in the module CLAUDE.md; thresholds are always reported as half-widths). + - `omega_E_hat` — `ω_E/ω_dia,e` (`≡ −ω₀/ω_dia,e`; a scanned input from day one, + ordering O4). + - `epsilon` — inverse aspect ratio `r_s/R₀`. + - `inv_Lq_hat` — `L̂_q⁻¹ = (ψ_s/q) dq/dψ` at `r_s`. + - `q_s` — safety factor at the rational surface. + - `tau` — `T_e/T_i`. + - `nu_star` — per-species banana collisionality `ν_★j`, keyed by species name. +""" +struct Level0Parameters + w_hat::Float64 + omega_E_hat::Float64 + epsilon::Float64 + inv_Lq_hat::Float64 + q_s::Float64 + tau::Float64 + nu_star::Dict{Symbol,Float64} + function Level0Parameters(w_hat, omega_E_hat, epsilon, inv_Lq_hat, q_s, tau, nu_star) + w_hat > 0 || throw(ArgumentError("w_hat must be positive (half-width)")) + epsilon > 0 || throw(ArgumentError("epsilon must be positive")) + q_s > 0 || throw(ArgumentError("q_s must be positive")) + tau > 0 || throw(ArgumentError("tau must be positive")) + all(v -> v >= 0, values(nu_star)) || throw(ArgumentError("nu_star entries must be nonnegative")) + return new(w_hat, omega_E_hat, epsilon, inv_Lq_hat, q_s, tau, Dict{Symbol,Float64}(nu_star)) + end +end + +function Level0Parameters(; w_hat, omega_E_hat, epsilon, inv_Lq_hat, q_s, tau=1.0, nu_star=Dict{Symbol,Float64}()) + return Level0Parameters(w_hat, omega_E_hat, epsilon, inv_Lq_hat, q_s, tau, nu_star) +end + +""" + FrameConvention(; C_dia=NaN, sign_omega0=NaN, C_gradient_shift=NaN) + +The gated sign/normalization set of the frame identities (`01 §5`, +`[CHECKED: Diss19 pp. 46–48]`, awaiting human clearance — QUESTIONS **Q3**). +All fields default to **NaN** so an un-cleared convention poisons results +rather than silently committing to a guess; `is_cleared` tests for this. + +## Fields + + - `C_dia` — prefactor (incl. sign) of the electron diamagnetic + frequency form `ω_dia,e = C_dia · m T̂_e L̂_n⁻¹ / q_s`. + - `sign_omega0` — the island-propagation relation `ω₀ = sign_omega0 · ω_E`. + - `C_gradient_shift` — prefactor of the frame shift of the effective density + gradient, `L_n⁻¹ = L_{n0}⁻¹ (1 + C_gradient_shift · Z_j ω_E/ω_dia,e)`. +""" +Base.@kwdef struct FrameConvention + C_dia::Float64 = NaN + sign_omega0::Float64 = NaN + C_gradient_shift::Float64 = NaN +end + +""" + is_cleared(conv::FrameConvention) + +`true` only when every gated field has been assigned (no NaN). Downstream code +must check this before producing physics numbers. +""" +is_cleared(conv::FrameConvention) = !(isnan(conv.C_dia) || isnan(conv.sign_omega0) || isnan(conv.C_gradient_shift)) + +""" + frame_shift(omega, omega_E) + +The frame-invariant combination `ω − ω_E` (`01 §5`): mechanical bookkeeping, +no sign convention (the invariance is what *defines* the shift). +""" +frame_shift(omega, omega_E) = omega - omega_E + +""" + omega_dia_form(m_mode, T_hat, inv_Ln_hat, q_s, conv) + +The electron diamagnetic frequency *form* `C_dia · m T̂ L̂_n⁻¹ / q_s` +(structure of `01 §5`; value gated on `conv.C_dia`, QUESTIONS Q3). Returns NaN +until the convention is cleared. +""" +omega_dia_form(m_mode, T_hat, inv_Ln_hat, q_s, conv::FrameConvention) = conv.C_dia * m_mode * T_hat * inv_Ln_hat / q_s + +""" + effective_dlnn_form(inv_Ln0_hat, Z, omega_E_hat, conv) + +The frame shift of the effective density gradient *form* +`L̂_n⁻¹ = L̂_{n0}⁻¹ (1 + C_gradient_shift · Z ω̂_E)` (structure of `01 §5`; value +gated on `conv.C_gradient_shift`, QUESTIONS Q3). `omega_E_hat` is already +normalized to `ω_dia,e`. Returns NaN until the convention is cleared. +""" +effective_dlnn_form(inv_Ln0_hat, Z, omega_E_hat, conv::FrameConvention) = inv_Ln0_hat * (1 + conv.C_gradient_shift * Z * omega_E_hat) + +end # module Frames diff --git a/src/Islands/moments/Moments.jl b/src/Islands/moments/Moments.jl new file mode 100644 index 00000000..a9c4d71d --- /dev/null +++ b/src/Islands/moments/Moments.jl @@ -0,0 +1,133 @@ +""" + Moments + +Output-moment assembly (design `01 §4`, `03 §1`): the parallel current +`J̄_∥(x, ξ)` from species-summed velocity moments, its `cos ξ`/`sin ξ` Ampère +projections `Δ_cos`/`Δ_sin`, and the island flux-surface-average diagnostics +(`Ω` label, `⟨·⟩_Ω`, bootstrap/polarization channel split). + +**Gating:** the projection and quadrature machinery here is pure numerics. The +physics enters through (i) the per-species velocity-space weights `W_j` (the +`v̂_∥`-structure — `[VERIFY]`-gated, QUESTIONS Q3), and (ii) the `Δ` moment +prefactors (`±μ₀R/2ψ̃` with `ψ̃` under an open `[VERIFY]`, and the sin-moment +normalization still `[DERIVED]`-unpinned — QUESTIONS Q4). Both are **required +caller-supplied arguments** with no defaults; nothing here assigns them. + +The island label convention is the module-CLAUDE.md pin: `Ω = 2x²/w² − cos ξ`, +O-point `Ω = −1`, separatrix `Ω = +1`, `w` = **half**-width. +""" +module Moments + +using LinearAlgebra +import QuadGK +import ..PhaseSpace: IslandGrid, nnodes +import ..Operators: weighted_moment! +import ..SpeciesLists: Species + +export parallel_current!, delta_moments, omega_label, omega_average, channel_split + +""" + parallel_current!(Jpar, gs, species, weights, grid) + +Assemble `J̄_∥(x, ξ) = Σ_j Z_j ∫ W_j g_j` into `Jpar[ix, iξ]` (`01 §4`): one +`weighted_moment!` per species, charge-scaled and accumulated. `gs`, `species` +and `weights` are aligned vectors; each `W_j` is the supplied `(ny, nE, nσ)` +parallel-flow velocity weight (gated physics, QUESTIONS Q3). +""" +function parallel_current!(Jpar, gs::AbstractVector, species::AbstractVector{<:Species}, weights::AbstractVector, grid::IslandGrid) + length(gs) == length(species) == length(weights) || + throw(ArgumentError("gs, species, weights must be aligned (got $(length(gs)), $(length(species)), $(length(weights)))")) + fill!(Jpar, zero(eltype(Jpar))) + for (g, sp, W) in zip(gs, species, weights) + weighted_moment!(Jpar, g, W, grid; scale=sp.Z, accumulate=true) + end + return Jpar +end + +""" + delta_moments(Jpar, grid; prefactor_cos, prefactor_sin) + +The two Ampère projections of `J̄_∥` through the island (`01 §4`): + + Δ_cos = prefactor_cos · ∫dx ∮dξ J̄_∥ cos ξ + Δ_sin = prefactor_sin · ∫dx ∮dξ J̄_∥ sin ξ + +`ξ`-integration is the uniform-grid trapezoid (spectrally exact on the periodic +grid); `x`-integration uses the grid's Simpson weights. The prefactors +(`∓μ₀R/2ψ̃`, `01 §4`) are **gated** — `ψ̃` carries an open `[VERIFY]` and the +sin normalization is `[DERIVED]`-unpinned (QUESTIONS Q4) — so both are required +arguments. Returns `(Δcos, Δsin)`. +""" +function delta_moments(Jpar, grid::IslandGrid; prefactor_cos, prefactor_sin) + nx, nξ = size(Jpar) + (nx, nξ) == (grid.x.n, grid.ξ.n) || throw(ArgumentError("Jpar must be (nx, nξ) = $((grid.x.n, grid.ξ.n))")) + wξ = grid.ξ.L / grid.ξ.n + pc = zero(eltype(Jpar)) + ps = zero(eltype(Jpar)) + @inbounds for iξ in 1:nξ + c = cos(grid.ξ.nodes[iξ]) + s = sin(grid.ξ.nodes[iξ]) + for ix in 1:nx + w = grid.x.wq[ix] * wξ * Jpar[ix, iξ] + pc += w * c + ps += w * s + end + end + return (Δcos=prefactor_cos * pc, Δsin=prefactor_sin * ps) +end + +""" + omega_label(x, ξ, w) + +The island flux-surface label `Ω = 2x²/w² − cos ξ` (pinned convention, module +CLAUDE.md; O-point `Ω = −1`, separatrix `Ω = +1`, `w` = half-width). A +diagnostic label only — never a solve coordinate (Decision D1). +""" +omega_label(x, ξ, w) = 2 * x^2 / w^2 - cos(ξ) + +# x on the Ω surface at helical angle ξ (branch = ±1 selects the x-sign). +_x_on_surface(Ω, ξ, w, branch) = branch * (w / sqrt(2)) * sqrt(max(Ω + cos(ξ), 0.0)) + +""" + omega_average(f, Ω, w; branch=+1, rtol=1e-10) + +Island flux-surface average `⟨f⟩_Ω = ∮ f (Ω + cos ξ)^{−1/2} dξ / ∮ (Ω + cos ξ)^{−1/2} dξ` +(`01 §4` decomposition weight) of a callable `f(x, ξ)`: + + - `Ω > 1` (outside the separatrix): the `±x` branches are distinct surfaces; + `branch` selects one, integrating over the full `ξ ∈ [0, 2π)`. + - `−1 < Ω < 1` (inside): the surface closes through both branches between the + turning points `cos ξ_b = −Ω`; the endpoint weight singularity is integrable + and handled by the adaptive quadrature. +""" +function omega_average(f, Ω, w; branch::Int=1, rtol::Real=1e-10) + if Ω > 1 + num, _ = QuadGK.quadgk(ξ -> f(_x_on_surface(Ω, ξ, w, branch), ξ) / sqrt(Ω + cos(ξ)), 0.0, 2π; rtol=rtol) + den, _ = QuadGK.quadgk(ξ -> 1.0 / sqrt(Ω + cos(ξ)), 0.0, 2π; rtol=rtol) + return num / den + elseif -1 < Ω < 1 + ξb = acos(-Ω) + integrand(ξ) = (f(_x_on_surface(Ω, ξ, w, +1), ξ) + f(_x_on_surface(Ω, ξ, w, -1), ξ)) / sqrt(Ω + cos(ξ)) + num, _ = QuadGK.quadgk(integrand, -ξb, ξb; rtol=rtol) + den, _ = QuadGK.quadgk(ξ -> 2.0 / sqrt(Ω + cos(ξ)), -ξb, ξb; rtol=rtol) + return num / den + else + throw(DomainError(Ω, "omega_average needs Ω ∈ (−1, 1) ∪ (1, ∞) (O-point/separatrix excluded)")) + end +end + +""" + channel_split(Jfun, Ω, w; branch=+1) + +The `01 §4` bootstrap/polarization bookkeeping at one `Ω` surface: returns +`(bs, pol)` where `bs = ⟨J⟩_Ω` (the flux-surface-constant "bootstrap+curvature" +part) and `pol(x, ξ) = Jfun(x, ξ) − bs` (the piece that `Ω`-averages to zero). +Diagnostic only — the solve never separates channels; the split is approximate +bookkeeping (L23 Eq. 2.5.3 caveat). +""" +function channel_split(Jfun, Ω, w; branch::Int=1) + bs = omega_average(Jfun, Ω, w; branch=branch) + return (bs=bs, pol=(x, ξ) -> Jfun(x, ξ) - bs) +end + +end # module Moments diff --git a/src/Islands/operators/Operators.jl b/src/Islands/operators/Operators.jl index 230dad73..3a7afba5 100644 --- a/src/Islands/operators/Operators.jl +++ b/src/Islands/operators/Operators.jl @@ -36,7 +36,10 @@ import ..PhaseSpace: IslandGrid, nnodes export IslandState, IslandCache, IslandStack, AbstractTerm export ParallelStreaming, MagneticDrift, ExBDrift, Collisions, GradientDrive, PerpTransport, RadiationSink, Quasineutrality -export apply!, residual!, velocity_moment!, statelength, flatten!, unflatten! +export PitchAngleDiffusion, conservative_pitch_operator +export FarFieldConditions, apply_farfield! +export apply!, residual!, velocity_moment!, weighted_moment!, statelength, + flatten!, unflatten!, g_flat_index, Φ_flat_index # --------------------------------------------------------------------------- # State, cache @@ -161,6 +164,53 @@ struct Collisions{A} <: AbstractTerm end Collisions(a_y, b_y; model::Symbol=:pitch_angle) = Collisions(a_y, b_y, model) +""" + PitchAngleDiffusion(K, c) + +Pitch-angle collision operator in **discretely conservative (mimetic) form** +(`01 §2.3` structure; ladder A4): adds `c ⋅ (K g)` along `y`, where `K` is the +divergence-form matrix built by [`conservative_pitch_operator`](@ref) and `c` is +a supplied coefficient array over `(x, ξ, E, σ)` — **`y`-independent by +construction**, so the exact discrete particle conservation and entropy-sign +properties of `K` are preserved (the physical `ν̂(v̂)` energy dependence lives in +`c`'s `E`-dependence and is `[VERIFY]`-gated). The momentum-restoring +field-particle piece of the Level-0 operator is gated physics (QUESTIONS Q3) +and is not part of this structure. +""" +struct PitchAngleDiffusion{M<:AbstractMatrix,A} <: AbstractTerm + K::M + c::A +end + +""" + conservative_pitch_operator(ygrid, P, wmeas) + +Build the mimetic divergence-form pitch operator on `ygrid` +(`C[g] = (1/w) d/dy(P w dg/dy)` discretized as `K = −Wq⁻¹ Gᵀ diag(P .* wq) G` +with `G = ygrid.D1` and `Wq = wmeas .* ygrid.wq`), returning `(K, Wq)`. + +Because `G` differentiates constants exactly and the quadrature weights are +positive, `K` satisfies **exactly** (to machine precision, ladder A4): + + - particle conservation: `Wqᵀ (K g) = 0` for any `g` (zero-flux built in); + - entropy sign: `gᵀ diag(Wq) K g = −(Gg)ᵀ diag(P .* wq)(Gg) ≤ 0` for `P ≥ 0`. + +`P` (the pitch-space diffusivity profile, physically the `λ√(1−λB)`-structure) +and `wmeas` (the velocity-space measure) are **supplied** profiles — their +Level-0 physics forms are `[VERIFY]`-gated (QUESTIONS Q3); any positive test +profiles exercise the conservation structure. +""" +function conservative_pitch_operator(ygrid, P::AbstractVector, wmeas::AbstractVector) + length(P) == ygrid.n || throw(ArgumentError("P must have length $(ygrid.n)")) + length(wmeas) == ygrid.n || throw(ArgumentError("wmeas must have length $(ygrid.n)")) + all(>=(0), P) || throw(ArgumentError("diffusivity profile P must be nonnegative")) + all(>(0), wmeas) || throw(ArgumentError("measure weights wmeas must be positive")) + G = ygrid.D1 + Wq = wmeas .* ygrid.wq + K = -Diagonal(1.0 ./ Wq) * (G' * Diagonal(P .* ygrid.wq) * G) + return Matrix(K), Wq +end + """ GradientDrive(drive) @@ -321,6 +371,24 @@ function apply!(R::IslandState, t::Collisions, U::IslandState, grid::IslandGrid, return R end +function apply!(R::IslandState, t::PitchAngleDiffusion, U::IslandState, grid::IslandGrid, ::IslandCache) + K = t.K + c = t.c + g = U.g + nx, nξ, ny, nE, nσ = size(g) + @inbounds for iσ in 1:nσ, iE in 1:nE, iξ in 1:nξ, ix in 1:nx + cv = c[ix, iξ, iE, iσ] + for a in 1:ny + acc = zero(eltype(R.g)) + for b in 1:ny + acc += K[a, b] * g[ix, iξ, b, iE, iσ] + end + R.g[ix, iξ, a, iE, iσ] += cv * acc + end + end + return R +end + function apply!(R::IslandState, t::GradientDrive, U::IslandState, ::IslandGrid, ::IslandCache) @inbounds @. R.g += t.drive return R @@ -388,6 +456,78 @@ function velocity_moment!(M, g, grid::IslandGrid; accumulate::Bool=false) return M end +""" + weighted_moment!(M, g, W, grid; scale=1, accumulate=false) + +Accumulate the weighted velocity moment `scale · ∫dy ∫dE Σ_σ W(y, E, σ) g` into +`M[ix, iξ]` (`03 §2`, moments). `W` is a supplied `(ny, nE, nσ)` velocity-space +weight — e.g. the `v̂_∥`-structure of the parallel-flow moment, whose Level-0 +physics form is `[VERIFY]`-gated (QUESTIONS Q3). `scale` carries a per-species +factor (e.g. the charge `Z_j`). `velocity_moment!` is the `W ≡ 1` special case. +""" +function weighted_moment!(M, g, W, grid::IslandGrid; scale=1, accumulate::Bool=false) + accumulate || fill!(M, zero(eltype(M))) + wy = grid.y.wq + wE = grid.E.weights + nx, nξ, ny, nE, nσ = size(g) + size(W) == (ny, nE, nσ) || throw(ArgumentError("weight W must have size (ny, nE, nσ) = $((ny, nE, nσ))")) + @inbounds for iσ in 1:nσ, iE in 1:nE, iy in 1:ny + w = scale * wy[iy] * wE[iE] * W[iy, iE, iσ] + for iξ in 1:nξ, ix in 1:nx + M[ix, iξ] += w * g[ix, iξ, iy, iE, iσ] + end + end + return M +end + +# --------------------------------------------------------------------------- +# Far-field boundary conditions (01 §3, 04 §1): match g and Φ to supplied +# far-field states at |x| = L_x. NEVER bare Neumann — L23 traced its spurious +# "winged" solution branch to Neumann non-uniqueness; the physical far field is +# the neoclassical (no-island) solution, which is [VERIFY]-gated physics and +# therefore *supplied* here, not computed. +# --------------------------------------------------------------------------- +""" + FarFieldConditions(g_left, g_right, Φ_left, Φ_right) + +Dirichlet-type far-field matching data at the two radial boundaries (`01 §3`): +`g → g_far` and `Φ̂ → Φ̂_far` as `|x| → L_x`. The residual rows at `ix = 1, nx` +are replaced by `(U − far)` so the Newton solve pins the far field. The +neoclassical no-island `g_far` is gated physics (QUESTIONS Q3) — callers supply +it (tests use manufactured far fields). + +## Fields + + - `g_left`, `g_right` — `(nξ, ny, nE, nσ)` far-field distributions at `ix = 1`, `ix = nx`. + - `Φ_left`, `Φ_right` — length-`nξ` far-field potentials at `ix = 1`, `ix = nx`. +""" +struct FarFieldConditions{T,A4<:AbstractArray{T,4},V<:AbstractVector{T}} + g_left::A4 + g_right::A4 + Φ_left::V + Φ_right::V +end + +""" + apply_farfield!(R, U, bc, grid) + +Replace the residual rows at the radial boundaries with the far-field matching +conditions `U − far` (`01 §3`). Called after the operator stack in +[`residual!`](@ref); allocation-free and AD-transparent. +""" +function apply_farfield!(R::IslandState, U::IslandState, bc::FarFieldConditions, grid::IslandGrid) + nx, nξ, ny, nE, nσ = size(U.g) + @inbounds for iσ in 1:nσ, iE in 1:nE, iy in 1:ny, iξ in 1:nξ + R.g[1, iξ, iy, iE, iσ] = U.g[1, iξ, iy, iE, iσ] - bc.g_left[iξ, iy, iE, iσ] + R.g[nx, iξ, iy, iE, iσ] = U.g[nx, iξ, iy, iE, iσ] - bc.g_right[iξ, iy, iE, iσ] + end + @inbounds for iξ in 1:nξ + R.Φ[1, iξ] = U.Φ[1, iξ] - bc.Φ_left[iξ] + R.Φ[nx, iξ] = U.Φ[nx, iξ] - bc.Φ_right[iξ] + end + return R +end + # --------------------------------------------------------------------------- # Stack + residual assembly # --------------------------------------------------------------------------- @@ -421,6 +561,19 @@ function residual!(R::IslandState, U::IslandState, stack::IslandStack, grid::Isl return R end +""" + residual!(R, U, stack, grid, cache, bc) + +Residual assembly with far-field boundary conditions: assemble the stack, then +replace the radial-boundary rows with the `FarFieldConditions` matching +residual (`01 §3`). +""" +function residual!(R::IslandState, U::IslandState, stack::IslandStack, grid::IslandGrid, cache::IslandCache, bc::FarFieldConditions) + residual!(R, U, stack, grid, cache) + apply_farfield!(R, U, bc, grid) + return R +end + # recursive tuple walk keeps the loop type-stable (no dynamic dispatch, no alloc) @inline _apply_kinetic!(R, ::Tuple{}, U, grid, cache) = R @inline function _apply_kinetic!(R, terms::Tuple, U, grid, cache) @@ -465,4 +618,26 @@ function unflatten!(U::IslandState, v) return U end +""" + g_flat_index(grid, ix, iξ, iy, iE, iσ) + +Flat-state-vector index of `g[ix, iξ, iy, iE, iσ]` under the `flatten!` layout +(`g` block column-major, then `Φ`). Used by the `y_c`-block conditioning +monitor (ladder A8) to address pitch-pencil sub-blocks of the dense Jacobian. +""" +function g_flat_index(grid::IslandGrid, ix::Int, iξ::Int, iy::Int, iE::Int, iσ::Int) + nx, nξ, ny, nE, = nnodes(grid) + return ix + nx * ((iξ - 1) + nξ * ((iy - 1) + ny * ((iE - 1) + nE * (iσ - 1)))) +end + +""" + Φ_flat_index(grid, ix, iξ) + +Flat-state-vector index of `Φ[ix, iξ]` under the `flatten!` layout. +""" +function Φ_flat_index(grid::IslandGrid, ix::Int, iξ::Int) + nx, nξ, ny, nE, nσ = nnodes(grid) + return nx * nξ * ny * nE * nσ + ix + nx * (iξ - 1) +end + end # module Operators diff --git a/src/Islands/phasespace/PhaseSpace.jl b/src/Islands/phasespace/PhaseSpace.jl index 9a9c241e..025132fa 100644 --- a/src/Islands/phasespace/PhaseSpace.jl +++ b/src/Islands/phasespace/PhaseSpace.jl @@ -302,11 +302,13 @@ internal layers (`x = 0`, `y = y_c`) per `04 §1`. ## Fields - - `x` — radial `MappedFDGrid` (clustered at `x = 0`). - - `ξ` — helical `FourierGrid`. - - `y` — pitch `MappedFDGrid` on `[0, y_max]` (clustered at `y_c`). - - `E` — energy `GaussGrid`. - - `σ` — `[+1.0, -1.0]`. + - `x` — radial `MappedFDGrid` (clustered at `x = 0`). + - `ξ` — helical `FourierGrid`. + - `y` — pitch `MappedFDGrid` on `[0, y_max]` (clustered at `y_c`). + - `E` — energy `GaussGrid`. + - `σ` — `[+1.0, -1.0]`. + - `y_c` — trapped–passing boundary location in `y` (the layer the pitch grid + packs toward; consumed by the `y_c`-block conditioning monitor, ladder A8). """ struct IslandGrid x::MappedFDGrid @@ -314,6 +316,7 @@ struct IslandGrid y::MappedFDGrid E::GaussGrid σ::Vector{Float64} + y_c::Float64 end function IslandGrid(; nx::Int, nxi::Int, ny::Int, nE::Int, @@ -326,7 +329,7 @@ function IslandGrid(; nx::Int, nxi::Int, ny::Int, nE::Int, y = MappedFDGrid(ny; halfwidth=y_max, clustering=clustering_y, center=y_c, domain=:half, order=order) E = GaussGrid(nE; kind=energy_kind) - return IslandGrid(x, ξ, y, E, [1.0, -1.0]) + return IslandGrid(x, ξ, y, E, [1.0, -1.0], Float64(y_c)) end """ diff --git a/src/Islands/solvers/Solvers.jl b/src/Islands/solvers/Solvers.jl new file mode 100644 index 00000000..2b4799cf --- /dev/null +++ b/src/Islands/solvers/Solvers.jl @@ -0,0 +1,331 @@ +""" + Solvers + +The Islands steady-state solve (design `03 §3`, `04 §5`, Decision D2): +matrix-free **Newton–Krylov** — never initial-value time-stepping, never nested +Picard loops (kokuchou's Picard criterion was *never met* in production; the +sources' Φ-outer/ū_∥i-inner iteration is the explicit anti-pattern). + +Pieces: + + - [`flat_residual`](@ref) — wrap an operator stack (+ optional far-field BCs) + as a flat-vector residual `f!(out, u)`, generic over `eltype` with + preallocated real and dual workspaces. + - [`JVPOperator`](@ref) — the exact Jacobian–vector product via ForwardDiff + duals through the stack (`04 §5`); no global sparse Jacobian is ever formed + except in [`dense_jacobian`](@ref) tiny-grid debug mode. + - [`YBlockJacobi`](@ref) — the physics-block preconditioner skeleton: per-pitch- + pencil blocks with **TSVD-regularized** solves, the explicit `y_c`-matching + treatment of `04 §3` (GMRES must never resolve the near-singular directions + itself). + - [`newton_krylov`](@ref) — inexact Newton with Eisenstat–Walker forcing and + a backtracking line search. + - [`pseudo_arclength`](@ref) — the continuation scaffold with fold detection + from day one (`03 §3`: the Level-3 penetration bifurcation is a fold; the + solver must step around folds, not fail at them). + +Pure numerics — no physics coefficients (the stack it solves carries the gated +parameters). +""" +module Solvers + +using LinearAlgebra +import Krylov +import ForwardDiff +import ..PhaseSpace: IslandGrid, nnodes +import ..Operators: IslandState, IslandCache, IslandStack, FarFieldConditions, + residual!, flatten!, unflatten!, statelength + +export flat_residual, JVPOperator, dense_jacobian +export YBlockJacobi, newton_krylov, pseudo_arclength + +# Tag for the solver's ForwardDiff duals (one directional derivative). +struct JVPTag end +const JVPDual = ForwardDiff.Dual{JVPTag,Float64,1} + +# --------------------------------------------------------------------------- +# Flat residual wrapper +# --------------------------------------------------------------------------- +""" + flat_residual(stack, grid; bc=nothing) + +Wrap `residual!` for `stack` on `grid` (with optional `FarFieldConditions`) as +a flat-vector function `f!(out, u)` usable by [`newton_krylov`](@ref) and +[`JVPOperator`](@ref). Real (`Float64`) and dual workspaces are preallocated so +repeated evaluations allocate nothing. +""" +function flat_residual(stack::IslandStack, grid::IslandGrid; bc::Union{Nothing,FarFieldConditions}=nothing) + U = IslandState(grid) + R = IslandState(grid) + cache = IslandCache(grid) + Ud = IslandState{JVPDual}(grid) + Rd = IslandState{JVPDual}(grid) + cached = IslandCache{JVPDual}(grid) + function f!(out, u) + if eltype(u) === Float64 + unflatten!(U, u) + bc === nothing ? residual!(R, U, stack, grid, cache) : residual!(R, U, stack, grid, cache, bc) + flatten!(out, R) + else + unflatten!(Ud, u) + bc === nothing ? residual!(Rd, Ud, stack, grid, cached) : residual!(Rd, Ud, stack, grid, cached, bc) + flatten!(out, Rd) + end + return out + end + return f! +end + +# --------------------------------------------------------------------------- +# Matrix-free Jacobian-vector product +# --------------------------------------------------------------------------- +""" + JVPOperator(f!, u) + +Matrix-free linear operator for the Jacobian of `f!` at the point `u` (which is +**aliased**, so updating `u` in place moves the linearization point): `mul!(y, J, v)` computes the exact directional derivative via one dual-mode evaluation +of `f!` (`04 §5`). Satisfies the `mul!`/`size`/`eltype` interface Krylov.jl +expects. +""" +struct JVPOperator{F} + f!::F + u::Vector{Float64} + udual::Vector{JVPDual} + rdual::Vector{JVPDual} +end + +function JVPOperator(f!, u::Vector{Float64}) + n = length(u) + return JVPOperator(f!, u, Vector{JVPDual}(undef, n), Vector{JVPDual}(undef, n)) +end + +Base.size(J::JVPOperator) = (length(J.u), length(J.u)) +Base.size(J::JVPOperator, d::Int) = d <= 2 ? length(J.u) : 1 +Base.eltype(::JVPOperator) = Float64 + +function LinearAlgebra.mul!(y::AbstractVector, J::JVPOperator, v::AbstractVector) + @inbounds for i in eachindex(J.u) + J.udual[i] = JVPDual(J.u[i], ForwardDiff.Partials((v[i],))) + end + J.f!(J.rdual, J.udual) + @inbounds for i in eachindex(y) + y[i] = ForwardDiff.partials(J.rdual[i])[1] + end + return y +end + +""" + dense_jacobian(f!, u) + +Dense Jacobian of `f!` at `u` by column-wise dual sweeps — **tiny-grid debug +mode only** (`04 §5`): used for eigenvalue/conditioning diagnostics (the +`y_c`-block singular-value monitor, ladder A8), never in the production solve. +""" +function dense_jacobian(f!, u::Vector{Float64}) + n = length(u) + J = zeros(n, n) + ud = Vector{JVPDual}(undef, n) + rd = Vector{JVPDual}(undef, n) + for j in 1:n + @inbounds for i in 1:n + ud[i] = JVPDual(u[i], ForwardDiff.Partials((i == j ? 1.0 : 0.0,))) + end + f!(rd, ud) + @inbounds for i in 1:n + J[i, j] = ForwardDiff.partials(rd[i])[1] + end + end + return J +end + +# --------------------------------------------------------------------------- +# Physics-block preconditioner skeleton (04 §3, §5) +# --------------------------------------------------------------------------- +""" + YBlockJacobi(grid, block; svd_cutoff=1e-10, phi_scale=1.0) + +Block-Jacobi preconditioner over pitch pencils — the skeleton of the `04 §5` +physics-block preconditioner. `block(ix, iξ, iE, iσ)` returns the `ny × ny` +stiff-direction block for that pencil (typically the identity-plus-collisions +operator); each block is factored by SVD and **truncated below +`svd_cutoff · σ_max`** — the explicit, regularized `y_c`-matching treatment of +`04 §3` (kokuchou's rcond ≈ 1e-16 block gave machine-dependent noise under +plain LU; TSVD is the documented fix). `phi_scale` is the diagonal applied to +the `Φ` rows. Apply via `ldiv!` (pass as `N=...` with `ldiv=true` to GMRES). +""" +struct YBlockJacobi + grid::IslandGrid + pinvs::Array{Matrix{Float64},4} + phi_scale::Float64 +end + +function YBlockJacobi(grid::IslandGrid, block; svd_cutoff::Real=1e-10, phi_scale::Real=1.0) + nx, nξ, ny, nE, nσ = nnodes(grid) + pinvs = Array{Matrix{Float64},4}(undef, nx, nξ, nE, nσ) + for iσ in 1:nσ, iE in 1:nE, iξ in 1:nξ, ix in 1:nx + B = Matrix{Float64}(block(ix, iξ, iE, iσ)) + size(B) == (ny, ny) || throw(ArgumentError("block must be ny×ny = ($ny, $ny)")) + F = svd(B) + smax = F.S[1] + sinv = [s > svd_cutoff * smax ? 1.0 / s : 0.0 for s in F.S] + pinvs[ix, iξ, iE, iσ] = F.V * Diagonal(sinv) * F.U' + end + return YBlockJacobi(grid, pinvs, Float64(phi_scale)) +end + +function LinearAlgebra.ldiv!(y::AbstractVector, P::YBlockJacobi, x::AbstractVector) + nx, nξ, ny, nE, nσ = nnodes(P.grid) + ng = nx * nξ * ny * nE * nσ + xg = reshape(view(x, 1:ng), nx, nξ, ny, nE, nσ) + yg = reshape(view(y, 1:ng), nx, nξ, ny, nE, nσ) + @inbounds for iσ in 1:nσ, iE in 1:nE, iξ in 1:nξ, ix in 1:nx + Pi = P.pinvs[ix, iξ, iE, iσ] + for a in 1:ny + acc = 0.0 + for b in 1:ny + acc += Pi[a, b] * xg[ix, iξ, b, iE, iσ] + end + yg[ix, iξ, a, iE, iσ] = acc + end + end + @inbounds for i in (ng+1):length(x) + y[i] = x[i] / P.phi_scale + end + return y +end + +# --------------------------------------------------------------------------- +# Inexact Newton-Krylov +# --------------------------------------------------------------------------- +""" + newton_krylov(f!, u0; rtol=1e-8, atol=1e-10, max_iter=25, precond=nothing, + eta0=0.5, eta_max=0.9, gamma=0.9, ls_max=10, memory=30, + verbose=false) + +Matrix-free inexact Newton–Krylov solve of `f!(out, u) = 0` (design `03 §3`, +`04 §5`, Decision D2): + + - GMRES (Krylov.jl) on the ForwardDiff [`JVPOperator`](@ref), with the + **Eisenstat–Walker** (choice-2) forcing term `η_k = γ(‖F_k‖/‖F_{k−1}‖)²` so + early Newton steps do not over-solve the linear system; + - a backtracking (Armijo-on-`‖F‖`) line search; + - optional right preconditioner `precond` applied via `ldiv!` (e.g. + [`YBlockJacobi`](@ref)). + +Convergence is declared on **both** the global residual norm and its max-norm +(`04 §5`: the array-averaged residual hid locally divergent regions in the +prior art). Returns a NamedTuple `(u, converged, iterations, residual_norms, residual_max, gmres_iters)` — the total Krylov iteration count is a tracked +preconditioner-quality metric (`04 §5`), not folklore. +""" +function newton_krylov(f!, u0::AbstractVector{Float64}; + rtol::Real=1e-8, atol::Real=1e-10, max_iter::Int=25, precond=nothing, + eta0::Real=0.5, eta_max::Real=0.9, gamma::Real=0.9, ls_max::Int=10, + memory::Int=30, verbose::Bool=false) + u = collect(u0) + n = length(u) + F = similar(u) + f!(F, u) + nF = norm(F) + nF0 = max(nF, eps()) + tol = max(atol, rtol * nF0) + hist = Float64[nF] + J = JVPOperator(f!, u) # aliases u: updates move the linearization point + Ftrial = similar(u) + utrial = similar(u) + eta = float(eta0) + nF_prev = nF + k = 0 + gmres_total = 0 + while nF > tol && k < max_iter + k += 1 + k > 1 && (eta = clamp(gamma * (nF / nF_prev)^2, 1e-4, eta_max)) + rhs = -F + du, stats = + precond === nothing ? + Krylov.gmres(J, rhs; rtol=eta, atol=0.1 * atol, memory=memory, itmax=10n) : + Krylov.gmres(J, rhs; rtol=eta, atol=0.1 * atol, memory=memory, itmax=10n, N=precond, ldiv=true) + gmres_total += stats.niter + verbose && @info "newton iter $k: ‖F‖=$nF, η=$eta, gmres iters=$(stats.niter)" + # backtracking line search on ‖F‖ + λ = 1.0 + nFtrial = Inf + for _ in 1:ls_max + @. utrial = u + λ * du + f!(Ftrial, utrial) + nFtrial = norm(Ftrial) + nFtrial <= (1 - 1e-4 * λ) * nF && break + λ /= 2 + end + nF_prev = nF + if nFtrial < nF + u .= utrial + F .= Ftrial + nF = nFtrial + else + verbose && @warn "newton line search stagnated at iter $k" + break + end + push!(hist, nF) + end + return (u=u, converged=(nF <= tol), iterations=k, residual_norms=hist, residual_max=maximum(abs, F), gmres_iters=gmres_total) +end + +# --------------------------------------------------------------------------- +# Pseudo-arclength continuation scaffold (03 §3) +# --------------------------------------------------------------------------- +""" + pseudo_arclength(f!, u0, p0; ds=0.1, nsteps=10, newton_kwargs...) + +Pseudo-arclength continuation scaffold for a parameterized residual +`f!(out, u, p)` (`03 §3`): predictor along the (secant) tangent, corrector = +[`newton_krylov`](@ref) on the extended system `[F(u, p); t·(z − z_pred)]`, and +**fold detection from day one** via sign changes of the tangent's `p`-component +(the Level-3 penetration bifurcation is a fold in this curve). Returns +`(us, ps, folds)` where `folds` lists the step indices at which `dp/ds` +reversed. +""" +function pseudo_arclength(f!, u0::AbstractVector{Float64}, p0::Real; + ds::Real=0.1, nsteps::Int=10, newton_kwargs...) + n = length(u0) + # converge onto the curve at fixed p0 + sol = newton_krylov((out, u) -> f!(out, u, p0), collect(u0); newton_kwargs...) + sol.converged || error("pseudo_arclength: initial solve at p0 did not converge") + us = [copy(sol.u)] + ps = [float(p0)] + folds = Int[] + # initial tangent from du/dp: J du = -dF/dp (finite difference in p) + Fp = zeros(n) + Fm = zeros(n) + hp = max(1e-7, 1e-7 * abs(p0)) + f!(Fp, sol.u, p0 + hp) + f!(Fm, sol.u, p0 - hp) + dFdp = (Fp .- Fm) ./ (2hp) + Jop = JVPOperator((out, u) -> f!(out, u, p0), copy(sol.u)) + dudp, _ = Krylov.gmres(Jop, -dFdp; rtol=1e-10) + t = vcat(dudp, 1.0) + t ./= norm(t) + tp_prev = t[end] + for step in 1:nsteps + zpred = vcat(us[end], ps[end]) .+ ds .* t + tloc = copy(t) # freeze tangent for this corrector + function ext!(out, z) + f!(view(out, 1:n), view(z, 1:n), z[n+1]) + out[n+1] = LinearAlgebra.dot(tloc, z) - LinearAlgebra.dot(tloc, zpred) + return out + end + sol = newton_krylov(ext!, zpred; newton_kwargs...) + sol.converged || break + push!(us, sol.u[1:n]) + push!(ps, sol.u[n+1]) + # secant tangent for the next step; fold when dp/ds changes sign + t = vcat(us[end] .- us[end-1], ps[end] - ps[end-1]) + t ./= norm(t) + if t[end] * tp_prev < 0 + push!(folds, step) + end + tp_prev = t[end] + end + return (us=us, ps=ps, folds=folds) +end + +end # module Solvers diff --git a/src/Islands/species/Species.jl b/src/Islands/species/Species.jl new file mode 100644 index 00000000..f1948150 --- /dev/null +++ b/src/Islands/species/Species.jl @@ -0,0 +1,170 @@ +""" + SpeciesLists + +Species as a first-class dimension of the Islands solve (design `02 §1`, +Decision D3): every kinetic object is indexed by species, and the species list +is an array from day one even though Level-0 *physics* is single-bulk-ion. + +Pure data structures and bookkeeping — no physics coefficients live here. +Quantities follow the Level-0 normalization (`01 §5`): densities to `n_e`, +temperatures to `T_i`, masses to the reference bulk ion, gradients as inverse +scale lengths at `r_s`. +""" +module SpeciesLists + +export AbstractBackground, Maxwellian, SlowingDown +export SpeciesRole, Bulk, Trace +export Species, is_bulk, is_trace, bulk_species, trace_species +export validate_species, check_trace_criteria + +# --------------------------------------------------------------------------- +# Backgrounds +# --------------------------------------------------------------------------- +""" + AbstractBackground + +Supertype of species background distributions (`02 §1.1`). `Maxwellian` is the +Level-0 background; `SlowingDown` enters at Level 2 (energetic particles) and +changes the energy-grid weight map, not the machinery. +""" +abstract type AbstractBackground end + +""" + Maxwellian(; n, T, dlnn_dr, dlnT_dr) + +Maxwellian background at the rational surface `r_s` (`02 §1.1`), carrying the +strictly-local constant gradients of ordering O2. + +## Fields + + - `n` — density normalized to `n_e`. + - `T` — temperature normalized to `T_i`. + - `dlnn_dr` — inverse density scale length `L̂_n⁻¹` at `r_s`. + - `dlnT_dr` — inverse temperature scale length `L̂_T⁻¹` at `r_s`. +""" +Base.@kwdef struct Maxwellian <: AbstractBackground + n::Float64 + T::Float64 + dlnn_dr::Float64 = 0.0 + dlnT_dr::Float64 = 0.0 +end + +""" + SlowingDown(; S0, v_birth, v_crit, dlnS_dr) + +Slowing-down background (isotropic at Level-2 entry, `02 §3`). Present now so +the background abstraction is exercised from day one; the Level-2 physics that +consumes it is out of Level-0 scope. + +## Fields + + - `S0` — source rate normalization. + - `v_birth` — birth speed (normalized to the reference thermal speed). + - `v_crit` — critical speed (from `T_e` and composition). + - `dlnS_dr` — inverse source scale length at `r_s`. +""" +Base.@kwdef struct SlowingDown <: AbstractBackground + S0::Float64 + v_birth::Float64 + v_crit::Float64 + dlnS_dr::Float64 = 0.0 +end + +# --------------------------------------------------------------------------- +# Roles and the species struct +# --------------------------------------------------------------------------- +""" + SpeciesRole + +`Bulk` = full nonlinear participant (its `g` enters the Newton state vector). +`Trace` = linear post-processing pass with frozen bulk fields (`02 §1.2`) — +one linear solve, an additive contribution to `Δ_cos`/`Δ_sin`. +""" +@enum SpeciesRole Bulk Trace + +""" + Species(; name, Z, m, background, role, collisional_coupling=false) + +One species of the Islands solve (`02 §1.1`). + +## Fields + + - `name` — identifying symbol (unique within a list). + - `Z` — charge number (`Float64` to allow mean-charge-state impurities). + - `m` — mass ratio to the reference bulk ion. + - `background` — an `AbstractBackground` (Maxwellian at Level 0). + - `role` — `Bulk` or `Trace` (`02 §1.2`). + - `collisional_coupling` — whether the bulk feels this species' field-particle + back-reaction (`02 §1.3`; L1+ physics, plumbing present from day one). +""" +Base.@kwdef struct Species{B<:AbstractBackground} + name::Symbol + Z::Float64 + m::Float64 + background::B + role::SpeciesRole = Bulk + collisional_coupling::Bool = false +end + +is_bulk(sp::Species) = sp.role == Bulk +is_trace(sp::Species) = sp.role == Trace + +""" + bulk_species(list) / trace_species(list) + +Filter a species list by role. +""" +bulk_species(list::AbstractVector{<:Species}) = filter(is_bulk, list) +trace_species(list::AbstractVector{<:Species}) = filter(is_trace, list) + +# --------------------------------------------------------------------------- +# Validation and the trace-criteria check +# --------------------------------------------------------------------------- +""" + validate_species(list) + +Structural validation of a species list: unique names, positive `Z`-magnitude, +`m`, density and temperature, and at least one `Bulk` species. Throws +`ArgumentError` on violation; returns the list for chaining. +""" +function validate_species(list::AbstractVector{<:Species}) + isempty(list) && throw(ArgumentError("species list is empty")) + names = [sp.name for sp in list] + length(unique(names)) == length(names) || throw(ArgumentError("duplicate species names: $names")) + any(is_bulk, list) || throw(ArgumentError("species list needs at least one Bulk species")) + for sp in list + sp.m > 0 || throw(ArgumentError("species $(sp.name): mass ratio must be positive")) + abs(sp.Z) > 0 || throw(ArgumentError("species $(sp.name): charge must be nonzero")) + if sp.background isa Maxwellian + sp.background.n > 0 || throw(ArgumentError("species $(sp.name): density must be positive")) + sp.background.T > 0 || throw(ArgumentError("species $(sp.name): temperature must be positive")) + end + end + return list +end + +""" + check_trace_criteria(list; threshold=0.05) + +Check every `Trace` species against both trace criteria (`02 §1.2`): charge +trace (`n_j Z_j ≪ n_e`) and current trace (`n_j ≪ n_e`), with densities already +normalized to `n_e`. A violating species gets an `@warn` recommending `Bulk` +promotion — **warn, never silently degrade** (`02 §1.2` promotion rule). +`threshold` is a diagnostic heuristic knob for "≪", not a physics coefficient. +Returns the vector of violating species names (empty when clean). +""" +function check_trace_criteria(list::AbstractVector{<:Species}; threshold::Real=0.05) + violators = Symbol[] + for sp in trace_species(list) + sp.background isa Maxwellian || continue + n = sp.background.n + charge_frac = n * abs(sp.Z) + if charge_frac > threshold || n > threshold + push!(violators, sp.name) + @warn "species $(sp.name) violates trace criteria (n Z = $(charge_frac), n = $(n) vs threshold $(threshold)); run it as Bulk" maxlog = 1 + end + end + return violators +end + +end # module SpeciesLists diff --git a/src/Islands/verify/Verify.jl b/src/Islands/verify/Verify.jl index 4c585b04..b41f3054 100644 --- a/src/Islands/verify/Verify.jl +++ b/src/Islands/verify/Verify.jl @@ -27,15 +27,17 @@ using LinearAlgebra import ForwardDiff import ..PhaseSpace: IslandGrid, MappedFDGrid, GaussGrid, nnodes import ..Operators: IslandState, IslandCache, IslandStack, residual!, velocity_moment!, - apply!, statelength, flatten!, unflatten!, + apply!, statelength, flatten!, unflatten!, g_flat_index, ParallelStreaming, MagneticDrift, ExBDrift, Collisions, GradientDrive, - PerpTransport, Quasineutrality + PerpTransport, RadiationSink, Quasineutrality, FarFieldConditions +import ..Solvers: flat_residual, newton_krylov export manufactured_state, test_coefficients, build_stack, mms_operator_error, mms_assembled_error, estimate_order, jvp_fd_maxerror, moment_selfconvergence, - term_allocations, residual_allocations + term_allocations, residual_allocations, + yc_block_sigma_min, solve_mms, zero_drive_setup # --------------------------------------------------------------------------- # Manufactured solution: smooth, separable, ξ-periodic and bandlimited. @@ -347,4 +349,112 @@ function residual_allocations(grid::IslandGrid; coeffs=test_coefficients(grid)) return @allocated residual!(R, U, stack, grid, cache) end +# --------------------------------------------------------------------------- +# Solve-level MMS configurations (ladder A1-solve, A5) — the manufactured +# configurations live here per 04 §4 so tests and benchmark scripts share them. +# --------------------------------------------------------------------------- +""" + zero_drive_setup(grid) + +The ladder-A5 configuration: the manufactured advective stack with **all drives +off** (no `GradientDrive`, no source), so `g ≡ 0, Φ ≡ 0` is the exact solution +and the residual there is machine zero (`01 §6`). Returns `(f, N)` — the flat +residual function and state length — for the null test and the trivial-Newton +check. The `RadiationSink(−1)` entry is a **unit relaxation shift** (a +manufactured test coefficient, not physics) making the zero state locally unique. +""" +function zero_drive_setup(grid::IslandGrid) + c = test_coefficients(grid) + shift = fill(-1.0, size(c.drive)) + kin = (ParallelStreaming(c.a_xi, c.a_x), MagneticDrift(c.c_D; variant=:original), + ExBDrift(c.c_E), Collisions(c.a_y, c.b_y; model=:pitch_angle), RadiationSink(shift)) + stack = IslandStack(kin, Quasineutrality(c.α)) + return (f=flat_residual(stack, grid), N=statelength(grid)) +end + +""" + solve_mms(nx; nxi=8, ny=9, nE=2, rtol=1e-10, memory=300) + +The assembled **solve-level** MMS (the ladder-A1 extension to a full converged +Newton–Krylov solve): the advective manufactured stack (streaming + drift + +E×B + unit relaxation shift + quasineutrality) with far-field matching BCs +taken from the manufactured state, forced by the *analytic* continuous source +so the discrete solution differs from `g*` by the discretization error. The +first-order-in-`x` stack needs only the `x` far-field conditions (the pitch +operator's degenerate zero-flux structure is exercised separately, ladder A4). +Returns `(err, converged, iterations, gmres_iters)` where `err` is the max-norm +solution error against the manufactured state — refining `nx` must show the +design order. +""" +function solve_mms(nx::Int; nxi::Int=8, ny::Int=9, nE::Int=2, rtol::Real=1e-10, memory::Int=300) + grid = IslandGrid(; nx=nx, nxi=nxi, ny=ny, nE=nE, halfwidth_x=6.0, clustering_x=1.0, + y_max=4.0, y_c=1.0, clustering_y=0.8, order=4) + c = test_coefficients(grid) + shift = fill(-1.0, size(c.drive)) + kin = (ParallelStreaming(c.a_xi, c.a_x), MagneticDrift(c.c_D; variant=:original), + ExBDrift(c.c_E), RadiationSink(shift)) + stack = IslandStack(kin, Quasineutrality(c.α)) + Ustar, deriv = manufactured_state(grid) + bc = FarFieldConditions(Ustar.g[1, :, :, :, :], Ustar.g[nx, :, :, :, :], Ustar.Φ[1, :], Ustar.Φ[nx, :]) + f0! = flat_residual(stack, grid; bc=bc) + N = statelength(grid) + # continuous source: the analytic values of every active term at (g*, Φ*) + Sg = c.a_xi .* deriv.dgdξ .+ c.a_x .* deriv.dgdx .+ c.c_D .* deriv.dgdξ .+ Ustar.g + dPx = reshape(deriv.dΦdx, nx, nxi, 1, 1, 1) + dPξ = reshape(deriv.dΦdξ, nx, nxi, 1, 1, 1) + Sg .+= c.c_E .* (dPξ .* deriv.dgdx .- dPx .* deriv.dgdξ) + Sg[1, :, :, :, :] .= 0.0 + Sg[nx, :, :, :, :] .= 0.0 # BC rows carry no source + MΦ = zeros(nx, nxi) + velocity_moment!(MΦ, Ustar.g, grid) + SΦ = MΦ .- c.α .* Ustar.Φ # field rows: discretely consistent source + SΦ[1, :] .= 0.0 + SΦ[nx, :] .= 0.0 + S = zeros(N) + flatten!(S, IslandState(Sg, SΦ)) + f!(out, u) = (f0!(out, u); out .-= S; out) + sol = newton_krylov(f!, zeros(N); rtol=rtol, atol=1e-13, max_iter=30, memory=memory) + ustar = zeros(N) + flatten!(ustar, Ustar) + return (err=maximum(abs, sol.u .- ustar), converged=sol.converged, + iterations=sol.iterations, gmres_iters=sol.gmres_iters) +end + +# --------------------------------------------------------------------------- +# y_c matching-block conditioning monitor (ladder A8, 04 §3). +# --------------------------------------------------------------------------- +""" + yc_block_sigma_min(J, grid; window=3) + +Ladder-A8 conditioning monitor: the smallest singular value of the pitch-pencil +sub-block of the (tiny-grid debug) Jacobian `J` nearest the trapped–passing +boundary `y_c`, minimized over all `(x, ξ, E, σ)` pencils. The prior art's +`y_c` matching system was intrinsically near-singular (rcond ≈ 1e-16, L23 §4.2) +and produced machine-dependent **noise, not crashes** — so this must be +*tested for*, not observed. Returns `(sigma_min, pencil)` where `pencil` is the +`(ix, iξ, iE, iσ)` of the minimizing block. +""" +function yc_block_sigma_min(J::AbstractMatrix, grid::IslandGrid; window::Int=3) + nx, nξ, ny, nE, nσ = nnodes(grid) + window <= ny || throw(ArgumentError("window ($window) exceeds ny ($ny)")) + # contiguous y-index window centered on the node nearest y_c + ic = argmin(abs.(grid.y.nodes .- grid.y_c)) + lo = clamp(ic - window ÷ 2, 1, ny - window + 1) + iys = lo:(lo+window-1) + σmin = Inf + pencil = (0, 0, 0, 0) + idx = Vector{Int}(undef, window) + for iσ in 1:nσ, iE in 1:nE, iξ in 1:nξ, ix in 1:nx + for (k, iy) in enumerate(iys) + idx[k] = g_flat_index(grid, ix, iξ, iy, iE, iσ) + end + s = svdvals(J[idx, idx])[end] + if s < σmin + σmin = s + pencil = (ix, iξ, iE, iσ) + end + end + return (sigma_min=σmin, pencil=pencil) +end + end # module Verify diff --git a/test/runtests.jl b/test/runtests.jl index 7679f942..98a7bfb5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -38,4 +38,5 @@ else include("./runtests_rerun_from_h5.jl") include("./runtests_islands_grids.jl") include("./runtests_islands_operators.jl") + include("./runtests_islands_solve.jl") end diff --git a/test/runtests_islands_solve.jl b/test/runtests_islands_solve.jl new file mode 100644 index 00000000..2a56c7cf --- /dev/null +++ b/test/runtests_islands_solve.jl @@ -0,0 +1,227 @@ +# runtests_islands_solve.jl +# +# Islands M2 — Level-0 solve machinery, structural verification gates +# (docs/src/islands/design/05 §A; the M2 milestone contract): +# A5 zero-drive null: g ≡ 0 ⇒ residual = machine zero; Newton converges trivially. +# A1+ assembled solve-MMS: Newton–Krylov recovers the manufactured state at design order. +# A8 y_c matching-block smallest-singular-value conditioning monitor. +# A4 (L0) exact discrete particle conservation + entropy sign of the mimetic +# pitch-angle collision operator. +# A3 Δ_cos even / Δ_sin odd parity under ξ-reflection (manufactured J̄_∥). +# A7 the coefficient-free closure identity ⟨∂²h/∂x²⟩_Ω = 0. +# + preconditioner quality (GMRES iteration reduction), far-field BC rows, +# pseudo-arclength fold detection, species/frames plumbing. +# +# STRUCTURAL (pre-physics) checks: manufactured order-unity coefficients only. +# Every physics coefficient in src/ is a [VERIFY]-gated supplied parameter +# (QUESTIONS Q2–Q4); the York-gate physics benchmarks stay skipped in +# benchmarks/islands/ until human clearance. + +using LinearAlgebra + +const IslM2 = GeneralizedPerturbedEquilibrium.Islands +const PS2 = IslM2.PhaseSpace +const Op2 = IslM2.Operators +const So2 = IslM2.Solvers +const V2 = IslM2.Verify +const Mo2 = IslM2.Moments +const Fi2 = IslM2.Fields +const Sp2 = IslM2.SpeciesLists +const Fr2 = IslM2.Frames + +_sgrid(n; ny=n) = PS2.IslandGrid(; nx=n, nxi=8, ny=ny, nE=2, halfwidth_x=6.0, clustering_x=1.0, + y_max=4.0, y_c=1.0, clustering_y=0.8, order=4) + +@testset "Islands L0 solve machinery (M2)" begin + + @testset "species plumbing (02 §1, D3)" begin + bg = Sp2.Maxwellian(; n=1.0, T=1.0, dlnn_dr=-1.0) + ion = Sp2.Species(; name=:D, Z=1.0, m=1.0, background=bg, role=Sp2.Bulk) + trc = Sp2.Species(; name=:Dtrace, Z=1.0, m=1.0, background=Sp2.Maxwellian(; n=1e-4, T=1.0), role=Sp2.Trace) + @test Sp2.validate_species([ion, trc]) == [ion, trc] + @test Sp2.is_bulk(ion) && Sp2.is_trace(trc) + @test Sp2.bulk_species([ion, trc]) == [ion] + # the L0 test pair: trace deuterium copy of the bulk passes the criteria + @test isempty(Sp2.check_trace_criteria([ion, trc])) + # a heavy trace at bulk-like density violates and is flagged (warn, never degrade) + w = Sp2.Species(; name=:W, Z=40.0, m=92.0, background=Sp2.Maxwellian(; n=0.01, T=1.0), role=Sp2.Trace) + @test Sp2.check_trace_criteria([ion, w]) == [:W] + # structural validation + @test_throws ArgumentError Sp2.validate_species([trc]) # no Bulk + @test_throws ArgumentError Sp2.validate_species([ion, ion]) # duplicate names + end + + @testset "frames: gated conventions poison, mechanics work (01 §5)" begin + conv = Fr2.FrameConvention() # all NaN until human-cleared (Q3) + @test !Fr2.is_cleared(conv) + @test isnan(Fr2.omega_dia_form(2, 1.0, -1.0, 2.0, conv)) + @test isnan(Fr2.effective_dlnn_form(-1.0, 1.0, 0.5, conv)) + # the frame-invariant combination is pure bookkeeping + @test Fr2.frame_shift(1.7, 0.4) ≈ 1.3 + # parameter-vector validation + p = Fr2.Level0Parameters(; w_hat=1.0, omega_E_hat=0.0, epsilon=0.1, inv_Lq_hat=1.0, q_s=2.0) + @test p.tau == 1.0 + @test_throws ArgumentError Fr2.Level0Parameters(-1.0, 0.0, 0.1, 1.0, 2.0, 1.0, Dict{Symbol,Float64}()) + end + + @testset "A4 — mimetic pitch operator: exact conservation + entropy sign" begin + g = _sgrid(9) + P = @. g.y.nodes * (4.0 - g.y.nodes) + 0.1 # arbitrary positive test profile + wm = @. 1.0 + 0.1 * g.y.nodes + K, Wq = Op2.conservative_pitch_operator(g.y, P, wm) + for gv in (sin.(g.y.nodes) .+ 0.3 .* g.y.nodes .^ 2, exp.(-g.y.nodes), g.y.nodes) + @test abs(dot(Wq, K * gv)) < 1e-11 # particle conservation, machine level + @test dot(gv .* Wq, K * gv) <= 1e-13 # entropy sign (0 only for constants) + end + @test abs(dot(ones(g.y.n) .* Wq, K * ones(g.y.n))) < 1e-12 # constants in the null space + @test_throws ArgumentError Op2.conservative_pitch_operator(g.y, -P, wm) + # the term applies c ⋅ (K g) and is allocation-free + nx, nξ, ny, nE, nσ = PS2.nnodes(g) + c4 = fill(2.0, nx, nξ, nE, nσ) + term = Op2.PitchAngleDiffusion(K, c4) + U, = V2.manufactured_state(g) + R = Op2.IslandState(g) + cache = Op2.IslandCache(g) + Op2.apply!(R, term, U, g, cache) + @test R.g[3, 2, :, 1, 1] ≈ 2.0 .* (K * U.g[3, 2, :, 1, 1]) + @test (@allocated Op2.apply!(R, term, U, g, cache)) == 0 + end + + @testset "A5 — zero-drive null test" begin + g = _sgrid(7; ny=7) + setup = V2.zero_drive_setup(g) + r = ones(setup.N) + setup.f(r, zeros(setup.N)) + @test maximum(abs, r) == 0.0 # residual is EXACTLY machine zero + # Newton from a small perturbation falls back to the zero state + sol = So2.newton_krylov(setup.f, 1e-3 .* sin.(1:setup.N); rtol=1e-12, atol=1e-12) + @test sol.converged + @test norm(sol.u) < 1e-9 + end + + @testset "A1 — assembled solve-MMS at design order" begin + r17 = V2.solve_mms(17) + r33 = V2.solve_mms(33) + @test r17.converged && r33.converged + # solution error against the manufactured state converges at design order + @test log(r17.err / r33.err) / log(33 / 17) > 3.3 + @test r33.err < 5e-3 + end + + @testset "preconditioner: y-block Jacobi with TSVD cuts GMRES iterations (04 §5)" begin + g = _sgrid(9) + nx, nξ, ny, nE, nσ = PS2.nnodes(g) + P = @. g.y.nodes * (4.0 - g.y.nodes) # degenerate endpoints: zero-flux built in + K, = Op2.conservative_pitch_operator(g.y, P, ones(ny)) + cstiff = fill(30.0, nx, nξ, nE, nσ) # stiff collisional pencil + shift = fill(-1.0, nx, nξ, ny, nE, nσ) + stack = Op2.IslandStack((Op2.PitchAngleDiffusion(K, cstiff), Op2.RadiationSink(shift)), + Op2.Quasineutrality(1.3)) + f0! = So2.flat_residual(stack, g) + N = Op2.statelength(g) + b = sin.((1:N) ./ 7) + f!(out, u) = (f0!(out, u); out .-= b; out) + pc = So2.YBlockJacobi(g, (ix, iξ, iE, iσ) -> I(ny) + cstiff[ix, iξ, iE, iσ] .* K; phi_scale=-1.3) + s0 = So2.newton_krylov(f!, zeros(N); rtol=1e-10, memory=300) + s1 = So2.newton_krylov(f!, zeros(N); rtol=1e-10, memory=300, precond=pc) + @test s0.converged && s1.converged + @test maximum(abs, s0.u .- s1.u) < 1e-8 # same solution + @test s1.gmres_iters < s0.gmres_iters ÷ 2 # preconditioner earns its keep + end + + @testset "far-field BCs replace the boundary residual rows (01 §3)" begin + g = _sgrid(7; ny=7) + nx, nξ, ny, nE, nσ = PS2.nnodes(g) + U, = V2.manufactured_state(g) + bc = Op2.FarFieldConditions(0.1 .+ zeros(nξ, ny, nE, nσ), 0.2 .+ zeros(nξ, ny, nE, nσ), + fill(0.3, nξ), fill(0.4, nξ)) + stack = V2.build_stack(g) + R = Op2.IslandState(g) + cache = Op2.IslandCache(g) + Op2.residual!(R, U, stack, g, cache, bc) + @test R.g[1, 2, 3, 1, 1] ≈ U.g[1, 2, 3, 1, 1] - 0.1 + @test R.g[nx, 2, 3, 1, 2] ≈ U.g[nx, 2, 3, 1, 2] - 0.2 + @test R.Φ[1, 5] ≈ U.Φ[1, 5] - 0.3 + @test R.Φ[nx, 5] ≈ U.Φ[nx, 5] - 0.4 + # interior rows are untouched by the BC application + R2 = Op2.IslandState(g) + Op2.residual!(R2, U, stack, g, cache) + @test R.g[2:(nx-1), :, :, :, :] == R2.g[2:(nx-1), :, :, :, :] + end + + @testset "A8 — y_c matching-block conditioning monitor" begin + g = _sgrid(7; ny=7) + setup = V2.zero_drive_setup(g) + J = So2.dense_jacobian(setup.f, zeros(setup.N)) + mon = V2.yc_block_sigma_min(J, g) + @test isfinite(mon.sigma_min) && mon.sigma_min > 0 + @test all(mon.pencil .>= 1) + # the monitor detects an artificially singularized pencil block (the + # silent-noise regression of L23 §4.2 must be *tested for*) + idx = [Op2.g_flat_index(g, 1, 1, iy, 1, 1) for iy in 3:5] + Jsing = copy(J) + Jsing[idx, :] .= 0.0 + @test V2.yc_block_sigma_min(Jsing, g).sigma_min < 1e-14 + end + + @testset "A3 — Δ_cos even / Δ_sin odd under ξ-reflection" begin + g = _sgrid(9) + J = [exp(-x^2) * (2.0 + 1.5 * cos(ξ) + 0.7 * sin(ξ)) for x in g.x.nodes, ξ in g.ξ.nodes] + Jr = [exp(-x^2) * (2.0 + 1.5 * cos(-ξ) + 0.7 * sin(-ξ)) for x in g.x.nodes, ξ in g.ξ.nodes] + d = Mo2.delta_moments(J, g; prefactor_cos=1.0, prefactor_sin=1.0) + dr = Mo2.delta_moments(Jr, g; prefactor_cos=1.0, prefactor_sin=1.0) + @test d.Δcos ≈ dr.Δcos atol = 1e-12 # even + @test d.Δsin ≈ -dr.Δsin atol = 1e-12 # odd + @test abs(d.Δsin) > 0.1 # the projection actually sees the sin part + # ξ-projection is spectrally exact: a pure cos ξ current has zero sin moment + Jc = [cos(ξ) for x in g.x.nodes, ξ in g.ξ.nodes] + @test abs(Mo2.delta_moments(Jc, g; prefactor_cos=1.0, prefactor_sin=1.0).Δsin) < 1e-13 + end + + @testset "moments: J̄_∥ assembly and Ω-average diagnostics" begin + g = _sgrid(9) + nx, nξ, ny, nE, nσ = PS2.nnodes(g) + U, = V2.manufactured_state(g) + W = ones(ny, nE, nσ) + ion = Sp2.Species(; name=:D, Z=1.0, m=1.0, background=Sp2.Maxwellian(; n=1.0, T=1.0), role=Sp2.Bulk) + anti = Sp2.Species(; name=:A, Z=-1.0, m=1.0, background=Sp2.Maxwellian(; n=1.0, T=1.0), role=Sp2.Bulk) + Jp = zeros(nx, nξ) + # equal & opposite charges with identical g cancel exactly + Mo2.parallel_current!(Jp, [U.g, U.g], [ion, anti], [W, W], g) + @test maximum(abs, Jp) < 1e-14 + # single species with W ≡ 1 reduces to the plain velocity moment + Mo2.parallel_current!(Jp, [U.g], [ion], [W], g) + Mref = zeros(nx, nξ) + Op2.velocity_moment!(Mref, U.g, g) + @test Jp ≈ Mref + # ⟨1⟩_Ω = 1 outside and inside the separatrix; constant J has no polarization part + @test Mo2.omega_average((x, ξ) -> 1.0, 1.5, 1.0) ≈ 1.0 rtol = 1e-8 + @test Mo2.omega_average((x, ξ) -> 1.0, 0.2, 1.0) ≈ 1.0 rtol = 1e-6 + cs = Mo2.channel_split((x, ξ) -> 3.0, 2.0, 1.0) + @test cs.bs ≈ 3.0 rtol = 1e-8 + @test abs(cs.pol(0.5, 1.0)) < 1e-8 + @test Mo2.omega_label(0.0, 0.0, 1.0) ≈ -1.0 # O-point + @test Mo2.omega_label(1.0, float(π), 1.0) ≈ 3.0 + end + + @testset "A7 — flattened-electron identity ⟨∂²h/∂x²⟩_Ω = 0 (coefficient-free)" begin + for Ω in (1.2, 2.0, 5.0), pref in (1.0, 3.7) # prefactor-independent + @test abs(Fi2.flat_average_d2h_dx2(Ω, 1.0; prefactor=pref)) < 1e-10 + end + @test Fi2.h_profile(0.5; prefactor=1.0) == 0.0 # exactly flat inside the separatrix + @test Fi2.h_profile(2.0; prefactor=1.0) > 0.0 + @test Fi2.Q_omega(3.0) > Fi2.Q_omega(1.5) # Q grows with Ω + @test !Fi2.is_cleared(Fi2.ElectronClosure()) # closure constants stay NaN-gated (Q3) + @test isnan(Fi2.ElectronClosure().k_HS) + end + + @testset "pseudo-arclength continuation detects the toy fold (03 §3)" begin + ftoy!(out, u, p) = (out[1] = u[1]^2 + p; out) + pa = So2.pseudo_arclength(ftoy!, [1.0], -1.0; ds=0.3, nsteps=15, rtol=1e-12, atol=1e-12) + @test length(pa.ps) > 10 # stepped through, no stall at the fold + @test !isempty(pa.folds) # the fold at p = 0 is detected + @test maximum(pa.ps) < 0.05 # never steps past the fold parameter + # both branches visited: u > 0 before the fold, u < 0 after + @test any(z -> z[1] > 0.5, pa.us) && any(z -> z[1] < -0.5, pa.us) + end +end