ISLANDS - M2: Level-0 solve machinery (Newton-Krylov, Δ moments, species/frames/fields) — structure, gated physics

.gitignore

@@ -17,6 +17,7 @@
 !docs/resources/**/*.pdf
 *.png
 !docs/src/assets/**/*.png
+!docs/src/islands/figures/*.png
 *.jld2
 .gitattributes
 Manifest.toml

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"

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 |

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).")

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).")

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.")

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

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)

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",
     ],