From f943f7a182105001f6ae997df958d1aa877e3b2c Mon Sep 17 00:00:00 2001 From: logan-nc Date: Wed, 8 Jul 2026 09:21:38 -0400 Subject: [PATCH 1/9] ISLANDS - NEW FEATURE - Add ForwardDiff dependency for operator-stack AD MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit M1 needs forward-mode AD for the JVP/AD-compatibility checks (design docs/src/islands/design/04-numerics.md §9). Adds ForwardDiff to Project.toml (already present transitively in the manifest); removes nothing. Co-Authored-By: Claude Opus 4.8 (1M context) --- Project.toml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Project.toml b/Project.toml index fd2c512d..31cd4deb 100644 --- a/Project.toml +++ b/Project.toml @@ -14,6 +14,7 @@ DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78" FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838" FastInterpolations = "9ea80cae-fc13-4c00-8066-6eaedb12f34b" +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" IMASdd = "c5a45a97-b3f9-491c-b9a7-aa88c3bc0067" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" @@ -43,6 +44,7 @@ DoubleFloats = "1.6.2" FFTW = "1.9.0" FastGaussQuadrature = "1.1.0" FastInterpolations = "0.4.10" +ForwardDiff = "1.4.1" HDF5 = "0.17.2" IMASdd = "8" JLD2 = "0.6.3" From de51b28c179a83189f3f20631f9fa3930e2d25d5 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Wed, 8 Jul 2026 09:21:39 -0400 Subject: [PATCH 2/9] ISLANDS - NEW FEATURE - M1 skeleton: phase-space grids + operator stack + verify harness MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Land the milestone-M1 numerical skeleton (design 03 §1-2, 04; ladder A1/A2), structure only — no [VERIFY] physics coefficients assigned in src/. - phasespace/PhaseSpace.jl: the (x, ξ, y, E, σ) grids with layer-clustered maps. Fourier spectral ∂ξ; Fornberg high-order finite differences for ∂x/∂y on sinh-stretched grids (per-derivative window widths so D1 and D2 are both 4th-order including boundaries); composite-Simpson quadrature weights; Gauss-Laguerre energy nodes. Pure numerics. - operators/Operators.jl: AbstractTerm + apply! + residual! and the term structs of 03 §2 (ParallelStreaming, MagneticDrift with the :original/:improved toggle, ExBDrift as the (x,ξ) Poisson bracket, Collisions, GradientDrive, PerpTransport and RadiationSink L4 stubs, Quasineutrality field residual). Every physics coefficient is a supplied data field, never a literal; allocation-free and generic over eltype so ForwardDiff duals flow through. - verify/Verify.jl: manufactured-solution (MMS) + AD-vs-finite-difference JVP harness, plus allocation probes. Manufactured coefficients are arbitrary order-unity test values (not physics), exercising the discretization only. - Islands.jl: wire the three submodules. physics-verifier: PASS (no [VERIFY]-policy violation). Co-Authored-By: Claude Opus 4.8 (1M context) --- src/Islands/Islands.jl | 26 +- src/Islands/operators/Operators.jl | 468 +++++++++++++++++++++++++++ src/Islands/phasespace/PhaseSpace.jl | 339 +++++++++++++++++++ src/Islands/verify/Verify.jl | 350 ++++++++++++++++++++ 4 files changed, 1173 insertions(+), 10 deletions(-) create mode 100644 src/Islands/operators/Operators.jl create mode 100644 src/Islands/phasespace/PhaseSpace.jl create mode 100644 src/Islands/verify/Verify.jl diff --git a/src/Islands/Islands.jl b/src/Islands/Islands.jl index f783e515..f6693c97 100644 --- a/src/Islands/Islands.jl +++ b/src/Islands/Islands.jl @@ -18,15 +18,21 @@ milestone M1 proceeds. """ module Islands -# Submodules land here as M1 proceeds, following the layout in -# docs/src/islands/design/03-architecture.md §1: -# include("phasespace/PhaseSpace.jl") # grids (x, ξ, λ, E, σ), layer-clustered maps -# include("species/Species.jl") # Species, backgrounds, roles -# include("frames/Frames.jl") # THE frequency/frame conversion module -# include("operators/Operators.jl") # the AbstractTerm stack -# include("fields/Fields.jl") # Φ̃ quasineutrality (A_∥ Ampère at L3) -# include("moments/Moments.jl") # Δ_cos, Δ_sin, profiles, channel decompositions -# include("solvers/Solvers.jl") # Newton–Krylov, continuation, trace pass -# include("verify/Verify.jl") # MMS + analytic-limit hooks, named configs +# 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. +include("phasespace/PhaseSpace.jl") # grids (x, ξ, λ→y, E, σ), layer-clustered maps +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) + +import .PhaseSpace +import .Operators +import .Verify end # module Islands diff --git a/src/Islands/operators/Operators.jl b/src/Islands/operators/Operators.jl new file mode 100644 index 00000000..230dad73 --- /dev/null +++ b/src/Islands/operators/Operators.jl @@ -0,0 +1,468 @@ +""" + Operators + +The Islands operator stack (`03 §2`): the state vector, the residual assembly, +and the `AbstractTerm` family whose sum is the drift-kinetic residual. + +**Milestone-M1 status: allocation-free, AD-compatible *structural* stubs.** Each +term implements its discretized differential/algebraic *structure* — which +derivatives act, in which coordinate, with which sign pattern — but takes every +physics coefficient as *supplied data* (`term.a_xi`, `term.c_D`, …). No literal +physics coefficient, sign, or normalization is written here: those carry +`[VERIFY]`/`[CHECKED]` tags (module `CLAUDE.md`) and are populated only after +human clearance. The verification harness (`Verify`) exercises the discretization +with arbitrary manufactured test coefficients, which is legitimate — it tests +the numerics, not the physics. + +Operator-stack rules enforced here (`03 §2`, `CLAUDE.md`): + + - terms are independent — no term inspects which others are active; + - `apply!` is generic over `eltype(U)` (ForwardDiff duals flow through) and + allocation-free on the hot path (regression-tested); + - orderings that *remove* structure are phase-space configurations, not terms. + +Term ↔ physics map (`03 §2`, `01 §2`): `ParallelStreaming` (island-induced +streaming, `∂ξ` + `∂x`), `MagneticDrift` (precession, `∂ξ`, `:original`/`:improved` +toggle), `ExBDrift` (the `E×B` Poisson bracket in `(x, ξ)`), `Collisions` +(pitch-angle diffusion in `y`), `GradientDrive` (`(v·∇F₀)` source), +`PerpTransport`/`RadiationSink` (Level-4 stubs), and `Quasineutrality` (the +Level-0 field residual, `01 §3`). +""" +module Operators + +using LinearAlgebra +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! + +# --------------------------------------------------------------------------- +# State, cache +# --------------------------------------------------------------------------- +""" + IslandState(g, Φ) + +The Level-0 unknowns (`03 §2`): the orbit-averaged distribution +`g[ix, iξ, iy, iE, iσ]` per grid point and the electrostatic potential +`Φ[ix, iξ]`. Parametric in the element type so ForwardDiff duals flow through. + +## Fields + + - `g` — 5D array over `(x, ξ, y, E, σ)`. + - `Φ` — 2D array over `(x, ξ)`. +""" +struct IslandState{T,A5<:AbstractArray{T,5},A2<:AbstractArray{T,2}} + g::A5 + Φ::A2 +end + +function IslandState{T}(grid::IslandGrid) where {T} + nx, nξ, ny, nE, nσ = nnodes(grid) + return IslandState(zeros(T, nx, nξ, ny, nE, nσ), zeros(T, nx, nξ)) +end +IslandState(grid::IslandGrid) = IslandState{Float64}(grid) + +Base.eltype(::IslandState{T}) where {T} = T + +function Base.similar(U::IslandState{T}) where {T} + return IslandState(similar(U.g), similar(U.Φ)) +end + +function fill_state!(U::IslandState, v) + fill!(U.g, v) + fill!(U.Φ, v) + return U +end + +""" + IslandCache{T}(grid) + +Per-solve scratch buffers, typed to the state element type `T` so the hot path +allocates nothing (and so AD duals get dual-typed scratch). Currently holds the +two potential-gradient fields consumed by the `E×B` bracket. +""" +struct IslandCache{T} + dΦdx::Matrix{T} + dΦdξ::Matrix{T} +end + +function IslandCache{T}(grid::IslandGrid) where {T} + nx, nξ, = nnodes(grid) + return IslandCache{T}(zeros(T, nx, nξ), zeros(T, nx, nξ)) +end +IslandCache(grid::IslandGrid) = IslandCache{Float64}(grid) + +# --------------------------------------------------------------------------- +# Term family +# --------------------------------------------------------------------------- +""" + AbstractTerm + +Supertype of every operator-stack term (`03 §2`). Each concrete term implements +`apply!(R, term, U, grid, cache)` to accumulate its contribution to the residual. +Terms are independent (no term inspects which others are active), generic over +`eltype(U)`, and allocation-free on the hot path. +""" +abstract type AbstractTerm end + +""" + ParallelStreaming(a_xi, a_x) + +Island-induced parallel streaming (`01 §2`, the `∂ξ`/`∂x` channel): adds +`a_xi ∂g/∂ξ + a_x ∂g/∂x` to the residual. `a_xi`, `a_x` are supplied coefficient +arrays shaped like `g` (physics values `[VERIFY]`-gated; never literals here). +""" +struct ParallelStreaming{A} <: AbstractTerm + a_xi::A + a_x::A +end + +""" + MagneticDrift(c_D; variant=:original) + +Orbit-averaged magnetic (precession) drift (`01 §2.1`): adds `c_D ∂g/∂ξ`. The +`variant` selects the `:original` (finite `L̂_B⁻¹`, I19) vs `:improved` +(`L̂_B⁻¹ → 0` proxy, D21) drift-frequency structure — the 8.73→1.46 ρ_bi toggle +(`docs/05 E1`). Structure only: `c_D` (the `ω̂_D` coefficient over `(y, E, σ)`) +is supplied data. +""" +struct MagneticDrift{A} <: AbstractTerm + c_D::A + variant::Symbol +end +MagneticDrift(c_D; variant::Symbol=:original) = MagneticDrift(c_D, variant) + +""" + ExBDrift(c_E) + +`E×B` advection as the Poisson bracket of `Φ` and `g` in `(x, ξ)` +(`01 §2`): adds `c_E (∂Φ/∂ξ · ∂g/∂x − ∂Φ/∂x · ∂g/∂ξ)`. This is the one Level-0 +kinetic term nonlinear in the state (couples `g` and `Φ`); `c_E` is a supplied +scalar coupling. +""" +struct ExBDrift{S} <: AbstractTerm + c_E::S +end + +""" + Collisions(a_y, b_y; model=:pitch_angle) + +Pitch-angle (Lorentz) collision operator (`01 §2.3`) as a second-order operator +in `y`: adds `a_y ∂²g/∂y² + b_y ∂g/∂y`. `model = :pitch_angle` is the Level-0 +form; `:fokker_planck` (Level 1) reuses the same slot. `a_y`, `b_y` are supplied +coefficient arrays (the `ν̂`-weighted diffusion structure, `[VERIFY]`-gated). +""" +struct Collisions{A} <: AbstractTerm + a_y::A + b_y::A + model::Symbol +end +Collisions(a_y, b_y; model::Symbol=:pitch_angle) = Collisions(a_y, b_y, model) + +""" + GradientDrive(drive) + +The `(v_E + v_D + v_ψ̃)·∇F₀` gradient drive (`03 §2`, `01 §2`): a state-independent +source, adds `drive` to the residual. `drive` is a supplied array shaped like +`g`; at Level 0 it is built from the background Maxwellian gradients (`[VERIFY]`). +""" +struct GradientDrive{A} <: AbstractTerm + drive::A +end + +""" + PerpTransport(χ) + +Perpendicular transport (Level-4 closure stub, `03 §2`): adds `χ ∂²g/∂x²`. +Present as structure only; the closure value `χ` is a supplied knob. +""" +struct PerpTransport{S} <: AbstractTerm + χ::S +end + +""" + RadiationSink(κ) + +Radiative energy sink (Level-4 closure stub, `03 §2`): adds `-κ g`. Structure +only; `κ` is a supplied coefficient array. +""" +struct RadiationSink{A} <: AbstractTerm + κ::A +end + +""" + Quasineutrality(α) + +The Level-0 field residual (`01 §3`): `R_Φ = M[g] − α Φ`, where `M[g]` is the +velocity moment `∫dy ∫dE Σ_σ g` (Gauss in `E`, Simpson in `y`) and `α` encodes +the flattened-electron closure scaling (structure of `e_iΦ̂/T_i/(2 L̂_{n0})`, +supplied). This closes `Φ(x, ξ)` inside the global Newton system rather than the +sources' fragile nested Picard loop (`01 §3`, `03 §3`). +""" +struct Quasineutrality{S} <: AbstractTerm + α::S +end + +# --------------------------------------------------------------------------- +# Discrete kernels (allocation-free, generic over eltype) +# +# Directional-derivative kernels accumulate `coef · ∂g/∂dir` into `Rg`. The +# dense matrices `D` are Float64; `D·g` with `g::Dual` promotes correctly, so +# the whole stack is AD-transparent. +# --------------------------------------------------------------------------- +# adds coef .* (∂g/∂ξ) along dim 2 +@inline function _add_dxi!(Rg, coef, g, D) + nx, nξ, ny, nE, nσ = size(g) + @inbounds for iσ in 1:nσ, iE in 1:nE, iy in 1:ny, ix in 1:nx + for a in 1:nξ + acc = zero(eltype(Rg)) + for b in 1:nξ + acc += D[a, b] * g[ix, b, iy, iE, iσ] + end + Rg[ix, a, iy, iE, iσ] += coef[ix, a, iy, iE, iσ] * acc + end + end + return Rg +end + +# adds coef .* (∂g/∂x) along dim 1 +@inline function _add_dx!(Rg, coef, g, D) + nx, nξ, ny, nE, nσ = size(g) + @inbounds for iσ in 1:nσ, iE in 1:nE, iy in 1:ny, iξ in 1:nξ + for a in 1:nx + acc = zero(eltype(Rg)) + for b in 1:nx + acc += D[a, b] * g[b, iξ, iy, iE, iσ] + end + Rg[a, iξ, iy, iE, iσ] += coef[a, iξ, iy, iE, iσ] * acc + end + end + return Rg +end + +# adds coef .* (Dⁿ g along dim 3, y). D is D1 or D2. +@inline function _add_dy!(Rg, coef, g, D) + 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 + for a in 1:ny + acc = zero(eltype(Rg)) + for b in 1:ny + acc += D[a, b] * g[ix, iξ, b, iE, iσ] + end + Rg[ix, iξ, a, iE, iσ] += coef[ix, iξ, a, iE, iσ] * acc + end + end + return Rg +end + +# scalar-coefficient ∂²/∂x² (PerpTransport): adds χ .* D2x g along dim 1 +@inline function _add_dx2_scalar!(Rg, χ, g, D) + nx, nξ, ny, nE, nσ = size(g) + @inbounds for iσ in 1:nσ, iE in 1:nE, iy in 1:ny, iξ in 1:nξ + for a in 1:nx + acc = zero(eltype(Rg)) + for b in 1:nx + acc += D[a, b] * g[b, iξ, iy, iE, iσ] + end + Rg[a, iξ, iy, iE, iσ] += χ * acc + end + end + return Rg +end + +# ∂Φ/∂x and ∂Φ/∂ξ into cache buffers (allocation-free) +@inline function _potential_gradients!(cache::IslandCache, Φ, Dx, Dξ) + nx, nξ = size(Φ) + dΦdx, dΦdξ = cache.dΦdx, cache.dΦdξ + @inbounds for iξ in 1:nξ, ix in 1:nx + ax = zero(eltype(Φ)) + for b in 1:nx + ax += Dx[ix, b] * Φ[b, iξ] + end + dΦdx[ix, iξ] = ax + end + @inbounds for ix in 1:nx, iξ in 1:nξ + aξ = zero(eltype(Φ)) + for b in 1:nξ + aξ += Dξ[iξ, b] * Φ[ix, b] + end + dΦdξ[ix, iξ] = aξ + end + return cache +end + +# --------------------------------------------------------------------------- +# apply!(R, term, U, grid, cache) — accumulate the term's contribution. +# --------------------------------------------------------------------------- +""" + apply!(R, term, U, grid, cache) + +Accumulate `term`'s contribution to the residual `R` at state `U` on `grid`, +using `cache` for scratch. Each `AbstractTerm` defines a method; the assembly in +[`residual!`](@ref) walks the stack. Allocation-free and generic over `eltype(U)` +so ForwardDiff duals flow through (verification ladder A2). +""" +function apply!(R::IslandState, t::ParallelStreaming, U::IslandState, grid::IslandGrid, ::IslandCache) + _add_dxi!(R.g, t.a_xi, U.g, grid.ξ.D1) + _add_dx!(R.g, t.a_x, U.g, grid.x.D1) + return R +end + +function apply!(R::IslandState, t::MagneticDrift, U::IslandState, grid::IslandGrid, ::IslandCache) + _add_dxi!(R.g, t.c_D, U.g, grid.ξ.D1) + return R +end + +function apply!(R::IslandState, t::Collisions, U::IslandState, grid::IslandGrid, ::IslandCache) + _add_dy!(R.g, t.a_y, U.g, grid.y.D2) + _add_dy!(R.g, t.b_y, U.g, grid.y.D1) + return R +end + +function apply!(R::IslandState, t::GradientDrive, U::IslandState, ::IslandGrid, ::IslandCache) + @inbounds @. R.g += t.drive + return R +end + +function apply!(R::IslandState, t::PerpTransport, U::IslandState, grid::IslandGrid, ::IslandCache) + _add_dx2_scalar!(R.g, t.χ, U.g, grid.x.D2) + return R +end + +function apply!(R::IslandState, t::RadiationSink, U::IslandState, ::IslandGrid, ::IslandCache) + @inbounds @. R.g += -t.κ * U.g + return R +end + +function apply!(R::IslandState, t::ExBDrift, U::IslandState, grid::IslandGrid, cache::IslandCache) + _potential_gradients!(cache, U.Φ, grid.x.D1, grid.ξ.D1) + dΦdx, dΦdξ = cache.dΦdx, cache.dΦdξ + g = U.g + Dx, Dξ = grid.x.D1, grid.ξ.D1 + cE = t.c_E + nx, nξ, ny, nE, nσ = size(g) + @inbounds for iσ in 1:nσ, iE in 1:nE, iy in 1:ny, iξ in 1:nξ, ix in 1:nx + dgdx = zero(eltype(R.g)) + for b in 1:nx + dgdx += Dx[ix, b] * g[b, iξ, iy, iE, iσ] + end + dgdξ = zero(eltype(R.g)) + for b in 1:nξ + dgdξ += Dξ[iξ, b] * g[ix, b, iy, iE, iσ] + end + R.g[ix, iξ, iy, iE, iσ] += cE * (dΦdξ[ix, iξ] * dgdx - dΦdx[ix, iξ] * dgdξ) + end + return R +end + +function apply!(R::IslandState, t::Quasineutrality, U::IslandState, grid::IslandGrid, ::IslandCache) + velocity_moment!(R.Φ, U.g, grid; accumulate=true) + @inbounds @. R.Φ += -t.α * U.Φ + return R +end + +# --------------------------------------------------------------------------- +# Velocity moment M[g](x,ξ) = ∫dy ∫dE Σ_σ g (Simpson in y, Gauss in E). +# --------------------------------------------------------------------------- +""" + velocity_moment!(M, g, grid; accumulate=false) + +Accumulate the phase-space velocity moment `∫dy ∫dE Σ_σ g` into `M[ix, iξ]` +using the Simpson `y`-weights and Gauss `E`-weights of `grid` (`03 §2`). Set +`accumulate=true` to add into `M` (residual assembly); otherwise `M` is zeroed +first. +""" +function velocity_moment!(M, g, grid::IslandGrid; accumulate::Bool=false) + accumulate || fill!(M, zero(eltype(M))) + wy = grid.y.wq + wE = grid.E.weights + nx, nξ, ny, nE, nσ = size(g) + @inbounds for iσ in 1:nσ, iE in 1:nE, iy in 1:ny + w = wy[iy] * wE[iE] + for iξ in 1:nξ, ix in 1:nx + M[ix, iξ] += w * g[ix, iξ, iy, iE, iσ] + end + end + return M +end + +# --------------------------------------------------------------------------- +# Stack + residual assembly +# --------------------------------------------------------------------------- +""" + IslandStack(kinetic, field) + +A named configuration (`03 §2`): the tuple of kinetic terms acting on `R.g` and +the `field` term (`Quasineutrality`) closing `R.Φ`. Stored as a tuple for type +stability so `residual!` stays allocation-free. + +## Fields + + - `kinetic` — tuple of `AbstractTerm`s contributing to the kinetic residual. + - `field` — the field-equation term. +""" +struct IslandStack{K<:Tuple,F<:AbstractTerm} + kinetic::K + field::F +end + +""" + residual!(R, U, stack, grid, cache) + +Assemble the full residual `R = (R_g, R_Φ)` in place: zero `R`, apply every +kinetic term, then the field term. Allocation-free and AD-transparent. +""" +function residual!(R::IslandState, U::IslandState, stack::IslandStack, grid::IslandGrid, cache::IslandCache) + fill_state!(R, zero(eltype(R))) + _apply_kinetic!(R, stack.kinetic, U, grid, cache) + apply!(R, stack.field, U, grid, cache) + 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) + apply!(R, terms[1], U, grid, cache) + return _apply_kinetic!(R, Base.tail(terms), U, grid, cache) +end + +# --------------------------------------------------------------------------- +# Flatten / unflatten for the Newton–Krylov / JVP interface +# --------------------------------------------------------------------------- +""" + statelength(grid) + +Number of scalar unknowns in the flattened state (`g` plus `Φ`). +""" +function statelength(grid::IslandGrid) + nx, nξ, ny, nE, nσ = nnodes(grid) + return nx * nξ * ny * nE * nσ + nx * nξ +end + +""" + flatten!(v, U) + +Copy state `U` into the flat vector `v` (`g` block then `Φ` block). +""" +function flatten!(v, U::IslandState) + ng = length(U.g) + copyto!(view(v, 1:ng), vec(U.g)) + copyto!(view(v, (ng+1):(ng+length(U.Φ))), vec(U.Φ)) + return v +end + +""" + unflatten!(U, v) + +Copy the flat vector `v` back into state `U`. +""" +function unflatten!(U::IslandState, v) + ng = length(U.g) + copyto!(vec(U.g), view(v, 1:ng)) + copyto!(vec(U.Φ), view(v, (ng+1):(ng+length(U.Φ)))) + return U +end + +end # module Operators diff --git a/src/Islands/phasespace/PhaseSpace.jl b/src/Islands/phasespace/PhaseSpace.jl new file mode 100644 index 00000000..9a9c241e --- /dev/null +++ b/src/Islands/phasespace/PhaseSpace.jl @@ -0,0 +1,339 @@ +""" + PhaseSpace + +Phase-space grids and discretization operators for the Islands drift-kinetic +solver: the `(x, ξ, λ→y, E, σ)` coordinates of design doc `03 §1` with the +layer-clustered mappings of `04 §1`. + +This module is **pure numerics** — grid node placement, spectral/finite-difference +differentiation matrices, and quadrature weights. It contains no physics +coefficients: nothing here carries a `[VERIFY]`/`[CHECKED]` tag, and the +milestone-M1 discipline (build the discretization, not the physics numbers) +lives one layer up in `Operators`. + +Design-order summary (verified by the MMS ladder A1, `Verify`): + + - `ξ` (helical angle): Fourier pseudo-spectral on the periodic domain `[0, L)` + — exponential convergence for smooth periodic data (`04 §1`). + - `x` (radial) and `y` (pitch): high-order finite differences on a stretched, + layer-clustered grid — algebraic convergence at the requested `order` + (`04 §1`; the layer packing targets `x = 0` and `y = y_c = 1`). + - `E` (energy): Gauss quadrature on the `F₀`-weighted semi-infinite domain + (`04 §1`, Maxwellian weight at Level 0 via Gauss–Laguerre). + - `σ = ±1`: the two `sgn(v_∥)` sheets. +""" +module PhaseSpace + +using LinearAlgebra +import FastGaussQuadrature + +export FourierGrid, MappedFDGrid, GaussGrid, IslandGrid +export nnodes, differentiate_fourier, fd_weights + +# --------------------------------------------------------------------------- +# Finite-difference weights (Fornberg 1988, Math. Comp. 51, 699) — weights for +# derivatives 0..m of a function sampled at arbitrary `nodes`, evaluated at x0. +# Returns `w[j+1, d+1]` = weight of node j for the d-th derivative. +# --------------------------------------------------------------------------- +""" + fd_weights(m, nodes, x0) + +Finite-difference weights for derivatives `0:m` of a function sampled at +`nodes`, evaluated at `x0`, via Fornberg's recurrence. `w[j, d+1]` multiplies +`f(nodes[j])` in the approximation of the `d`-th derivative. Exact for +polynomials up to degree `length(nodes) - 1`. +""" +function fd_weights(m::Int, nodes::AbstractVector{<:Real}, x0::Real) + n = length(nodes) + @assert m < n "need at least m+1 nodes for the m-th derivative" + w = zeros(Float64, n, m + 1) + c1 = 1.0 + c4 = nodes[1] - x0 + w[1, 1] = 1.0 + for i in 2:n + mn = min(i, m + 1) + c2 = 1.0 + c5 = c4 + c4 = nodes[i] - x0 + for j in 1:(i-1) + c3 = nodes[i] - nodes[j] + c2 *= c3 + if j == i - 1 + for k in mn:-1:2 + w[i, k] = c1 * ((k - 1) * w[i-1, k-1] - c5 * w[i-1, k]) / c2 + end + w[i, 1] = -c1 * c5 * w[i-1, 1] / c2 + end + for k in mn:-1:2 + w[j, k] = (c4 * w[j, k] - (k - 1) * w[j, k-1]) / c3 + end + w[j, 1] = c4 * w[j, 1] / c3 + end + c1 = c2 + end + return w +end + +# --------------------------------------------------------------------------- +# ξ: Fourier pseudo-spectral, periodic on [0, L). +# --------------------------------------------------------------------------- +""" + FourierGrid(n; L=2π) + +Uniform periodic grid of `n` nodes on `[0, L)` with the Fourier spectral +first-derivative matrix `D1` (`04 §1`, `ξ` coordinate). `n` must be even. + +## Fields + + - `n` — number of nodes. + - `L` — period length. + - `nodes` — node positions `j·L/n`, `j = 0:n-1`. + - `D1` — `n×n` spectral first-derivative matrix. +""" +struct FourierGrid + n::Int + L::Float64 + nodes::Vector{Float64} + D1::Matrix{Float64} +end + +function FourierGrid(n::Int; L::Real=2π) + iseven(n) || throw(ArgumentError("FourierGrid needs an even node count (got $n)")) + h = 2π / n # computational grid spacing on [0, 2π) + nodes = collect((0:(n-1)) .* (L / n)) + D1 = zeros(Float64, n, n) + @inbounds for i in 1:n, j in 1:n + if i != j + k = i - j + D1[i, j] = 0.5 * (-1.0)^k / tan(k * h / 2) + end + end + D1 .*= (2π / L) # rescale d/dξ from [0,2π) to [0,L) + return FourierGrid(n, Float64(L), nodes, D1) +end + +""" + differentiate_fourier(g, grid) + +Spectral first derivative of a periodic sample vector `g` on `grid` (allocating +convenience wrapper around `grid.D1 * g`; hot paths use the matrix directly). +""" +differentiate_fourier(g::AbstractVector, grid::FourierGrid) = grid.D1 * g + +# --------------------------------------------------------------------------- +# x, y: high-order finite differences on a layer-clustered stretched grid. +# +# A uniform computational coordinate s ∈ [-1, 1] is mapped to the physical +# coordinate by a smooth, monotonic, layer-clustering map; physical derivative +# matrices follow by the chain rule so the FD order is preserved (`04 §1`). +# --------------------------------------------------------------------------- +""" + MappedFDGrid(n; halfwidth, clustering=0.0, center=0.0, domain=:symmetric, order=4) + +High-order finite-difference grid of `n` nodes with layer clustering. + +The uniform computational coordinate `s ∈ [-1, 1]` is mapped to the physical +coordinate by `sinh` stretching (`clustering = β`): points cluster toward the +map center as `β` grows, and the map degenerates to uniform as `β → 0`. This is +the design-doc packing that targets the internal layers — `x = 0` (rational +surface) and `y = y_c = 1` (trapped–passing boundary), `04 §1–2`. + +`domain = :symmetric` gives `[center - halfwidth, center + halfwidth]` with +clustering at `center`; `domain = :half` gives `[0, halfwidth]` with clustering +at `center` (used for `y ∈ [0, y_max]` packed at `y_c`). + +Physical first/second-derivative matrices `D1`, `D2` are built with Fornberg +weights on windows sized per derivative (`order + d` points for the `d`-th +derivative, shifted one-sided near the boundaries) so both reach `order`-th +order accuracy *uniformly*, including at the boundary rows; the chain rule uses +the analytic map Jacobian `x'(s)` and `x''(s)`. + +`wq` holds composite-Simpson quadrature weights on the same nodes (built on the +uniform computational grid and pushed through the map Jacobian), so +`∫ f dx ≈ Σ wq[j] f(nodes[j])` at fourth order — matching the FD order for the +velocity-moment integrals (`03 §2`, moments). Requires an odd `n`. + +## Fields + + - `n`, `order` — node count and nominal FD order. + - `nodes` — physical node positions. + - `D1`, `D2` — `n×n` physical first/second-derivative matrices. + - `wq` — composite-Simpson quadrature weights on `nodes`. +""" +struct MappedFDGrid + n::Int + order::Int + nodes::Vector{Float64} + D1::Matrix{Float64} + D2::Matrix{Float64} + wq::Vector{Float64} +end + +function MappedFDGrid(n::Int; halfwidth::Real, clustering::Real=0.0, center::Real=0.0, + domain::Symbol=:symmetric, order::Int=4) + isodd(n) || throw(ArgumentError("MappedFDGrid needs odd n for composite Simpson (got $n)")) + # widest window is for D2 (order+2 points); need enough nodes for it. + n > order + 2 || throw(ArgumentError("need n > order+2 ($n ≤ $(order + 2))")) + + s = collect(range(-1.0, 1.0; length=n)) # uniform computational grid + hw = Float64(halfwidth) + β = Float64(clustering) + + # Map s → physical, with analytic first/second derivatives of the map. + if domain === :symmetric + # x(s) = center + hw·sinh(β s)/sinh(β): monotone, clusters at s=0 ↔ x=center. + if abs(β) < 1e-12 + xs = center .+ hw .* s + dxds = fill(hw, n) + d2xds2 = zeros(n) + else + sb = sinh(β) + xs = center .+ hw .* sinh.(β .* s) ./ sb + dxds = hw .* β .* cosh.(β .* s) ./ sb + d2xds2 = hw .* β^2 .* sinh.(β .* s) ./ sb + end + elseif domain === :half + # Map [-1,1] → [0, hw] with clustering at s* ↔ x=center via sinh about s*. + # u(s) = sinh(β(s - s*)); x = hw·(u - u(-1))/(u(1) - u(-1)). + sstar = clamp(2 * (center / hw) - 1, -1.0, 1.0) + if abs(β) < 1e-12 + xs = hw .* (s .+ 1) ./ 2 + dxds = fill(hw / 2, n) + d2xds2 = zeros(n) + else + u(z) = sinh(β * (z - sstar)) + um1 = u(-1.0) + span = u(1.0) - um1 + xs = hw .* (u.(s) .- um1) ./ span + dxds = hw .* β .* cosh.(β .* (s .- sstar)) ./ span + d2xds2 = hw .* β^2 .* sinh.(β .* (s .- sstar)) ./ span + end + else + throw(ArgumentError("unknown domain $domain (use :symmetric or :half)")) + end + + D1s = _fd_matrix(s, 1, order) + D2s = _fd_matrix(s, 2, order) + + # Chain rule: d/dx = (1/x')·d/ds ; d²/dx² = (1/x'²)·d²/ds² − (x''/x'³)·d/ds. + invp = 1.0 ./ dxds + D1 = Diagonal(invp) * D1s + D2 = Diagonal(invp .^ 2) * D2s - Diagonal(d2xds2 .* invp .^ 3) * D1s + + # Composite Simpson on uniform s (ds = 2/(n-1)), times the map Jacobian dx/ds. + ds = 2.0 / (n - 1) + wq = similar(dxds) + @inbounds for j in 1:n + c = (j == 1 || j == n) ? 1.0 : (iseven(j) ? 4.0 : 2.0) + wq[j] = c * ds / 3 * dxds[j] + end + + return MappedFDGrid(n, order, xs, Matrix(D1), Matrix(D2), wq) +end + +# Dense derivative matrix for the `deriv`-th derivative on an arbitrary 1D node +# set at the requested `order`, using windows of `order + deriv` points (the +# minimum for `order`-th accuracy of the `deriv`-th derivative), rounded up to +# an odd width so the interior window is centered; near the ends the window is +# shifted one-sided. Weights are Fornberg's, exact for polynomials up to the +# window degree. +function _fd_matrix(nodes::AbstractVector{<:Real}, deriv::Int, order::Int) + n = length(nodes) + stencil = order + deriv + isodd(stencil) || (stencil += 1) # centered interior window needs odd width + stencil = min(stencil, n) + D = zeros(Float64, n, n) + half = stencil ÷ 2 + @inbounds for i in 1:n + lo = clamp(i - half, 1, n - stencil + 1) + hi = lo + stencil - 1 + w = fd_weights(deriv, @view(nodes[lo:hi]), nodes[i]) + for (col, j) in enumerate(lo:hi) + D[i, j] = w[col, deriv+1] + end + end + return D +end + +# --------------------------------------------------------------------------- +# E: Gauss quadrature on the F₀-weighted semi-infinite energy domain. +# --------------------------------------------------------------------------- +""" + GaussGrid(n; kind=:laguerre) + +Gauss quadrature nodes/weights for velocity-space (energy) integrals over the +`F₀`-weighted semi-infinite domain (`04 §1`). `kind = :laguerre` gives the +Level-0 Maxwellian weight `∫₀^∞ f(E) e^{-E} dE ≈ Σ wᵢ f(Eᵢ)` (Gauss–Laguerre); +the slowing-down weight (Level 2) only changes `kind`, not the machinery. + +## Fields + + - `n` — number of quadrature nodes. + - `nodes` — abscissae `Eᵢ`. + - `weights` — quadrature weights `wᵢ` (the `F₀` weight is folded in). +""" +struct GaussGrid + n::Int + nodes::Vector{Float64} + weights::Vector{Float64} +end + +function GaussGrid(n::Int; kind::Symbol=:laguerre) + if kind === :laguerre + nodes, weights = FastGaussQuadrature.gausslaguerre(n) + elseif kind === :legendre + nodes, weights = FastGaussQuadrature.gausslegendre(n) + else + throw(ArgumentError("unknown quadrature kind $kind")) + end + return GaussGrid(n, collect(nodes), collect(weights)) +end + +# --------------------------------------------------------------------------- +# Bundle: the full Level-0 orbit-averaged 4D phase space (x, ξ, y, E) plus σ. +# --------------------------------------------------------------------------- +""" + IslandGrid(; nx, nxi, ny, nE, halfwidth_x, clustering_x, y_max, y_c, + clustering_y, xi_period=2π, energy_kind=:laguerre, order=4) + +The assembled Level-0 phase-space grid: radial `x`, helical `ξ`, pitch `y`, +energy `E`, and the two `σ = ±1` sheets (`03 §1`). Grids are packed at the +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]`. +""" +struct IslandGrid + x::MappedFDGrid + ξ::FourierGrid + y::MappedFDGrid + E::GaussGrid + σ::Vector{Float64} +end + +function IslandGrid(; nx::Int, nxi::Int, ny::Int, nE::Int, + halfwidth_x::Real, clustering_x::Real=0.0, + y_max::Real, y_c::Real=1.0, clustering_y::Real=0.0, + xi_period::Real=2π, energy_kind::Symbol=:laguerre, order::Int=4) + x = MappedFDGrid(nx; halfwidth=halfwidth_x, clustering=clustering_x, center=0.0, + domain=:symmetric, order=order) + ξ = FourierGrid(nxi; L=xi_period) + 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]) +end + +""" + nnodes(grid::IslandGrid) + +Tuple `(nx, nξ, ny, nE, nσ)` of the phase-space grid dimensions. +""" +nnodes(g::IslandGrid) = (g.x.n, g.ξ.n, g.y.n, g.E.n, length(g.σ)) + +end # module PhaseSpace diff --git a/src/Islands/verify/Verify.jl b/src/Islands/verify/Verify.jl new file mode 100644 index 00000000..4c585b04 --- /dev/null +++ b/src/Islands/verify/Verify.jl @@ -0,0 +1,350 @@ +""" + Verify + +The Islands verification harness (`04 §4`, `05 §A`): manufactured-solution (MMS) +convergence checks and AD-vs-finite-difference JVP checks for the operator +stack. Callable from the test suite (`test/runtests_islands_*.jl`) and from +scripts. + +**These are structural (pre-physics) checks — ladder A1/A2.** They exercise the +*discretization order* and the *AD plumbing* using arbitrary, order-unity +manufactured coefficients. That is deliberate and policy-clean: nothing here is +a physics coefficient, so nothing here carries a `[VERIFY]` tag. Physics +benchmarks (ladder B+) live in `benchmarks/islands/` and stay `[VERIFY]`-gated +until human-cleared. + +Design orders checked: + + - `ξ` derivatives — Fourier spectral (a bandlimited manufactured `ξ`-profile is + differentiated to machine precision). + - `x`, `y` derivatives — high-order finite differences at the grid `order` + (default 4): halving the mesh cuts the error by `2^order`. + - assembled kinetic residual — the min of the above (algebraic, set by `x`/`y`). +""" +module Verify + +using LinearAlgebra +import ForwardDiff +import ..PhaseSpace: IslandGrid, MappedFDGrid, GaussGrid, nnodes +import ..Operators: IslandState, IslandCache, IslandStack, residual!, velocity_moment!, + apply!, statelength, flatten!, unflatten!, + ParallelStreaming, MagneticDrift, ExBDrift, Collisions, GradientDrive, + PerpTransport, Quasineutrality + +export manufactured_state, + test_coefficients, build_stack, + mms_operator_error, mms_assembled_error, estimate_order, + jvp_fd_maxerror, moment_selfconvergence, + term_allocations, residual_allocations + +# --------------------------------------------------------------------------- +# Manufactured solution: smooth, separable, ξ-periodic and bandlimited. +# Each factor's analytic first/second derivatives are known in closed form, +# so the continuous value of every operator is exact at the nodes. +# --------------------------------------------------------------------------- +_Gx(x) = exp(-x^2 / 2) +_Gx′(x) = -x * _Gx(x) +_Gx″(x) = (x^2 - 1) * _Gx(x) +_Gξ(ξ) = sin(ξ) + 0.5 * cos(2ξ) +_Gξ′(ξ) = cos(ξ) - sin(2ξ) +_Gy(y) = exp(-(y - 1)^2 / 2) +_Gy′(y) = -(y - 1) * _Gy(y) +_Gy″(y) = ((y - 1)^2 - 1) * _Gy(y) +_GE(E) = exp(-0.3 * E) +_Gσ(σ) = 1 + 0.25 * σ + +_Px(x) = exp(-x^2 / 4) +_Px′(x) = -(x / 2) * _Px(x) +_Pξ(ξ) = cos(ξ) +_Pξ′(ξ) = -sin(ξ) + +# Fill a 5D array over (x, ξ, y, E, σ) from a scalar function of the node coords. +function _fill5(grid::IslandGrid, f) + nx, nξ, ny, nE, nσ = nnodes(grid) + A = Array{Float64}(undef, nx, nξ, ny, nE, nσ) + @inbounds for iσ in 1:nσ, iE in 1:nE, iy in 1:ny, iξ in 1:nξ, ix in 1:nx + A[ix, iξ, iy, iE, iσ] = f(grid.x.nodes[ix], grid.ξ.nodes[iξ], grid.y.nodes[iy], + grid.E.nodes[iE], grid.σ[iσ]) + end + return A +end + +function _fill2(grid::IslandGrid, f) + nx, nξ, = nnodes(grid) + A = Array{Float64}(undef, nx, nξ) + @inbounds for iξ in 1:nξ, ix in 1:nx + A[ix, iξ] = f(grid.x.nodes[ix], grid.ξ.nodes[iξ]) + end + return A +end + +""" + manufactured_state(grid) + +Return `(U, deriv)` where `U::IslandState` holds the manufactured `g*`, `Φ*` on +`grid` and `deriv` is a NamedTuple of the analytic partial-derivative arrays +(`dgdx`, `dgdξ`, `dgdy`, `d2gdy2`, `d2gdx2`, `dΦdx`, `dΦdξ`) used to form exact +continuous operator values. +""" +function manufactured_state(grid::IslandGrid) + g = _fill5(grid, (x, ξ, y, E, σ) -> _Gx(x) * _Gξ(ξ) * _Gy(y) * _GE(E) * _Gσ(σ)) + Φ = _fill2(grid, (x, ξ) -> _Px(x) * _Pξ(ξ)) + U = IslandState(g, Φ) + deriv = ( + dgdx=_fill5(grid, (x, ξ, y, E, σ) -> _Gx′(x) * _Gξ(ξ) * _Gy(y) * _GE(E) * _Gσ(σ)), + dgdξ=_fill5(grid, (x, ξ, y, E, σ) -> _Gx(x) * _Gξ′(ξ) * _Gy(y) * _GE(E) * _Gσ(σ)), + dgdy=_fill5(grid, (x, ξ, y, E, σ) -> _Gx(x) * _Gξ(ξ) * _Gy′(y) * _GE(E) * _Gσ(σ)), + d2gdy2=_fill5(grid, (x, ξ, y, E, σ) -> _Gx(x) * _Gξ(ξ) * _Gy″(y) * _GE(E) * _Gσ(σ)), + d2gdx2=_fill5(grid, (x, ξ, y, E, σ) -> _Gx″(x) * _Gξ(ξ) * _Gy(y) * _GE(E) * _Gσ(σ)), + dΦdx=_fill2(grid, (x, ξ) -> _Px′(x) * _Pξ(ξ)), + dΦdξ=_fill2(grid, (x, ξ) -> _Px(x) * _Pξ′(ξ)) + ) + return U, deriv +end + +# --------------------------------------------------------------------------- +# Test coefficients: arbitrary, smooth, order-unity. NOT physics — they exercise +# the discretization. `a_y > 0` keeps the pitch operator a sensible diffusion. +# --------------------------------------------------------------------------- +""" + test_coefficients(grid) + +Build a NamedTuple of arbitrary smooth manufactured coefficient arrays/scalars +for every term (`a_xi`, `a_x`, `c_D`, `a_y`, `b_y`, `drive`, `c_E`, `χ`, `α`). +These are MMS test values, not physics coefficients. +""" +function test_coefficients(grid::IslandGrid) + return ( + a_xi=_fill5(grid, (x, ξ, y, E, σ) -> 1.0 + 0.2 * cos(ξ) + 0.1 * x), + a_x=_fill5(grid, (x, ξ, y, E, σ) -> 0.8 + 0.1 * sin(2ξ) + 0.05 * y), + c_D=_fill5(grid, (x, ξ, y, E, σ) -> 0.5 * σ * (1 + 0.1 * E)), + a_y=_fill5(grid, (x, ξ, y, E, σ) -> 0.3 + 0.05 * x^2), + b_y=_fill5(grid, (x, ξ, y, E, σ) -> 0.1 * (y - 1)), + drive=_fill5(grid, (x, ξ, y, E, σ) -> 0.2 * exp(-x^2) * cos(ξ)), + c_E=0.7, + χ=0.15, + α=1.3 + ) +end + +# --------------------------------------------------------------------------- +# Continuous (exact) operator values from the analytic manufactured derivatives. +# --------------------------------------------------------------------------- +function _continuous(term::Symbol, deriv, c) + if term === :streaming + return c.a_xi .* deriv.dgdξ .+ c.a_x .* deriv.dgdx + elseif term === :drift + return c.c_D .* deriv.dgdξ + elseif term === :collisions + return c.a_y .* deriv.d2gdy2 .+ c.b_y .* deriv.dgdy + elseif term === :perp + return c.χ .* deriv.d2gdx2 + elseif term === :drive + return c.drive + elseif term === :exb + # broadcast the 2D potential gradients over (y, E, σ) + dΦdx = reshape(deriv.dΦdx, size(deriv.dΦdx, 1), size(deriv.dΦdx, 2), 1, 1, 1) + dΦdξ = reshape(deriv.dΦdξ, size(deriv.dΦdξ, 1), size(deriv.dΦdξ, 2), 1, 1, 1) + return c.c_E .* (dΦdξ .* deriv.dgdx .- dΦdx .* deriv.dgdξ) + else + throw(ArgumentError("unknown term $term")) + end +end + +function _term_object(term::Symbol, c) + if term === :streaming + return ParallelStreaming(c.a_xi, c.a_x) + elseif term === :drift + return MagneticDrift(c.c_D; variant=:original) + elseif term === :collisions + return Collisions(c.a_y, c.b_y; model=:pitch_angle) + elseif term === :perp + return PerpTransport(c.χ) + elseif term === :drive + return GradientDrive(c.drive) + elseif term === :exb + return ExBDrift(c.c_E) + else + throw(ArgumentError("unknown term $term")) + end +end + +# --------------------------------------------------------------------------- +# MMS error drivers +# --------------------------------------------------------------------------- +""" + mms_operator_error(grid, term) + +Max-norm error between the discrete `apply!` of a single kinetic `term` +(`:streaming`, `:drift`, `:exb`, `:collisions`, `:perp`, `:drive`) on the +manufactured state and its exact continuous value on `grid`. +""" +function mms_operator_error(grid::IslandGrid, term::Symbol) + U, deriv = manufactured_state(grid) + c = test_coefficients(grid) + R = IslandState(grid) + cache = IslandCache(grid) + apply!(R, _term_object(term, c), U, grid, cache) + exact = _continuous(term, deriv, c) + return maximum(abs, R.g .- exact) +end + +""" + mms_assembled_error(grid; terms=(:streaming,:drift,:exb,:collisions,:perp,:drive)) + +Max-norm error between the assembled discrete kinetic residual (sum of `terms`) +on the manufactured state and the summed continuous values. +""" +function mms_assembled_error(grid::IslandGrid; + terms=(:streaming, :drift, :exb, :collisions, :perp, :drive)) + U, deriv = manufactured_state(grid) + c = test_coefficients(grid) + R = IslandState(grid) + cache = IslandCache(grid) + exact = zeros(size(R.g)) + for t in terms + apply!(R, _term_object(t, c), U, grid, cache) + exact .+= _continuous(t, deriv, c) + end + return maximum(abs, R.g .- exact) +end + +""" + estimate_order(errors, refine) + +Observed convergence order from a sequence of `errors` at mesh-refinement +factors `refine` (ratio of successive node counts): `log(eₖ/eₖ₊₁)/log(rₖ)`. +Returns the vector of per-step orders. +""" +function estimate_order(errors::AbstractVector, refine::AbstractVector) + return [log(errors[k] / errors[k+1]) / log(refine[k]) for k in 1:(length(errors)-1)] +end + +# --------------------------------------------------------------------------- +# Assembled stack for the JVP check (includes the E×B nonlinearity + field). +# --------------------------------------------------------------------------- +""" + build_stack(grid; coeffs=test_coefficients(grid)) + +Assemble an `IslandStack` with the full Level-0-shaped term set (streaming, +drift, E×B, collisions, drive) plus the quasineutrality field term, using the +supplied (manufactured) coefficients. Used by the JVP check. +""" +function build_stack(grid::IslandGrid; coeffs=test_coefficients(grid)) + kinetic = ( + ParallelStreaming(coeffs.a_xi, coeffs.a_x), + MagneticDrift(coeffs.c_D; variant=:original), + ExBDrift(coeffs.c_E), + Collisions(coeffs.a_y, coeffs.b_y; model=:pitch_angle), + GradientDrive(coeffs.drive) + ) + return IslandStack(kinetic, Quasineutrality(coeffs.α)) +end + +""" + jvp_fd_maxerror(grid; h=1e-6, seed=1) + +Compare the forward-mode-AD Jacobian–vector product of the assembled residual +against a central finite-difference directional derivative, at the manufactured +state in a deterministic direction. Returns the max-norm difference (ladder A2). +The residual is nonlinear (E×B bracket + g/Φ coupling), so this is a genuine +check of both the AD plumbing and its correctness. +""" +function jvp_fd_maxerror(grid::IslandGrid; h::Float64=1e-6, seed::Int=1) + stack = build_stack(grid) + cache_f = IslandCache{Float64}(grid) + U, = manufactured_state(grid) + + N = statelength(grid) + u0 = Vector{Float64}(undef, N) + flatten!(u0, U) + # deterministic direction (no RNG dependence) + v = Float64[0.5 + 0.5 * sin(seed * k) for k in 1:N] + v ./= norm(v) + + # residual as a flat-vector function, generic over eltype + function resid_vec!(out, uvec) + T = eltype(uvec) + Ut = IslandState(reshape(view(uvec, 1:length(U.g)), size(U.g)), + reshape(view(uvec, (length(U.g)+1):length(uvec)), size(U.Φ))) + Rt = IslandState(reshape(view(out, 1:length(U.g)), size(U.g)), + reshape(view(out, (length(U.g)+1):length(out)), size(U.Φ))) + cache = IslandCache{T}(grid) + residual!(Rt, Ut, stack, grid, cache) + return out + end + + # AD JVP: d/dε residual(u0 + ε v) |_{ε=0} + jvp_ad = ForwardDiff.derivative(0.0) do ε + out = Vector{typeof(ε)}(undef, N) + resid_vec!(out, u0 .+ ε .* v) + end + + # central-difference JVP + rp = similar(u0) + rm = similar(u0) + resid_vec!(rp, u0 .+ h .* v) + resid_vec!(rm, u0 .- h .* v) + jvp_fd = (rp .- rm) ./ (2h) + + return maximum(abs, jvp_ad .- jvp_fd) +end + +# --------------------------------------------------------------------------- +# Velocity-moment quadrature self-convergence (Simpson in y, at grid order). +# --------------------------------------------------------------------------- +""" + moment_selfconvergence(grid_coarse, grid_fine) + +Max-norm difference between the velocity moment of the manufactured `g*` on two +`y`-resolutions. The two grids must share `(nx, nξ)` so the moments live on the +same `(x, ξ)` grid and compare pointwise; refine `ny` (and/or `nE`) between them +to probe the `y`-quadrature order. +""" +function moment_selfconvergence(grid_coarse::IslandGrid, grid_fine::IslandGrid) + nnodes(grid_coarse)[1:2] == nnodes(grid_fine)[1:2] || + throw(ArgumentError("moment_selfconvergence needs grids sharing (nx, nξ)")) + Uc, = manufactured_state(grid_coarse) + Uf, = manufactured_state(grid_fine) + Mc = zeros(nnodes(grid_coarse)[1], nnodes(grid_coarse)[2]) + Mf = zeros(nnodes(grid_fine)[1], nnodes(grid_fine)[2]) + velocity_moment!(Mc, Uc.g, grid_coarse) + velocity_moment!(Mf, Uf.g, grid_fine) + return maximum(abs, Mc .- Mf) +end + +# --------------------------------------------------------------------------- +# Allocation probes (hot-path allocation regression, `04 §9`). +# --------------------------------------------------------------------------- +""" + term_allocations(grid, term) + +Bytes allocated by a single warmed `apply!` of kinetic `term` on `grid`. Should +be zero for the allocation-free hot path. +""" +function term_allocations(grid::IslandGrid, term::Symbol) + c = test_coefficients(grid) + U, = manufactured_state(grid) + R = IslandState(grid) + cache = IslandCache(grid) + t = _term_object(term, c) + apply!(R, t, U, grid, cache) # warm up compilation + return @allocated apply!(R, t, U, grid, cache) +end + +""" + residual_allocations(grid; coeffs=test_coefficients(grid)) + +Bytes allocated by a single warmed full-`residual!` assembly on `grid`. Should +be zero for the allocation-free hot path. +""" +function residual_allocations(grid::IslandGrid; coeffs=test_coefficients(grid)) + stack = build_stack(grid; coeffs=coeffs) + U, = manufactured_state(grid) + R = IslandState(grid) + cache = IslandCache(grid) + residual!(R, U, stack, grid, cache) # warm up compilation + return @allocated residual!(R, U, stack, grid, cache) +end + +end # module Verify From 5e103e19c998fc6bc6a7e84239518417acafbab1 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Wed, 8 Jul 2026 09:21:56 -0400 Subject: [PATCH 3/9] ISLANDS - NEW FEATURE - M1 verification tests (ladder A1/A2) wired into runtests MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit test/runtests_islands_grids.jl and test/runtests_islands_operators.jl, included in test/runtests.jl: - A1 per-operator MMS: 4th-order convergence for the ∂x/∂y differential terms, machine precision for the Fourier ∂ξ (drift) term; assembled kinetic residual converges at 4th order. - A2: forward-mode-AD JVP of the (nonlinear) assembled residual matches the central finite-difference directional derivative to ~1e-9. - Allocation regression: every apply! and residual! hot path allocates 0 bytes. - Grid unit tests: Fourier exactness, mapped-FD order, Simpson/Gauss quadrature, layer-clustered packing, IslandGrid assembly and input validation. All 53 Islands tests pass. Co-Authored-By: Claude Opus 4.8 (1M context) --- test/runtests.jl | 2 + test/runtests_islands_grids.jl | 93 ++++++++++++++++++++++++++++ test/runtests_islands_operators.jl | 97 ++++++++++++++++++++++++++++++ 3 files changed, 192 insertions(+) create mode 100644 test/runtests_islands_grids.jl create mode 100644 test/runtests_islands_operators.jl diff --git a/test/runtests.jl b/test/runtests.jl index 2ae4125e..7679f942 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -36,4 +36,6 @@ else include("./runtests_coils.jl") include("./runtests_imas.jl") include("./runtests_rerun_from_h5.jl") + include("./runtests_islands_grids.jl") + include("./runtests_islands_operators.jl") end diff --git a/test/runtests_islands_grids.jl b/test/runtests_islands_grids.jl new file mode 100644 index 00000000..10dd3833 --- /dev/null +++ b/test/runtests_islands_grids.jl @@ -0,0 +1,93 @@ +# runtests_islands_grids.jl +# +# Islands module — phase-space grid / discretization unit tests (src/Islands/phasespace). +# Pure numerics: spectral and finite-difference differentiation and quadrature. +# No physics coefficients here (nothing [VERIFY]-tagged); these back the MMS +# ladder A1 checks in runtests_islands_operators.jl. + +const PS = GeneralizedPerturbedEquilibrium.Islands.PhaseSpace + +@testset "Islands phase-space grids" begin + + @testset "Fourier spectral ∂ξ is exact for bandlimited data" begin + fg = PS.FourierGrid(16; L=2π) + x = fg.nodes + g = @. sin(3x) + cos(2x) - 0.5 * sin(x) + dg_exact = @. 3cos(3x) - 2sin(2x) - 0.5cos(x) + @test maximum(abs, fg.D1 * g .- dg_exact) < 1e-12 + # odd node count is rejected + @test_throws ArgumentError PS.FourierGrid(15) + end + + @testset "Mapped FD converges at design order (uniform grid)" begin + # smooth decaying test function on [-6, 6] + f(x) = exp(-x^2 / 2) + f′(x) = -x * exp(-x^2 / 2) + f″(x) = (x^2 - 1) * exp(-x^2 / 2) + err1 = Float64[] + err2 = Float64[] + ns = [17, 33, 65] + for n in ns + g = PS.MappedFDGrid(n; halfwidth=6.0, order=4) + xn = g.nodes + push!(err1, maximum(abs, g.D1 * f.(xn) .- f′.(xn))) + push!(err2, maximum(abs, g.D2 * f.(xn) .- f″.(xn))) + end + # 4th-order: halving h cuts error by ≳ 2^4 (allow margin for pre-asymptotics) + @test log(err1[2] / err1[3]) / log(ns[3] / ns[2]) > 3.7 + @test log(err2[2] / err2[3]) / log(ns[3] / ns[2]) > 3.7 + @test err1[end] < 1e-3 + @test err2[end] < 1e-3 + end + + @testset "Layer-clustered map preserves order and packs the center" begin + # a strongly clustered grid still differentiates a smooth function accurately + g = PS.MappedFDGrid(65; halfwidth=6.0, clustering=2.0, order=4) + # node spacing is smallest near the clustering center (x = 0) + Δ = diff(g.nodes) + icenter = argmin(abs.(g.nodes[1:(end-1)] .+ g.nodes[2:end]) ./ 2) + @test Δ[icenter] < Δ[1] + @test Δ[icenter] < Δ[end] + f(x) = sin(x) * exp(-x^2 / 8) + f′(x) = (cos(x) - x / 4 * sin(x)) * exp(-x^2 / 8) + @test maximum(abs, g.D1 * f.(g.nodes) .- f′.(g.nodes)) < 1e-3 + end + + @testset "Half-domain grid packs at y_c and spans [0, y_max]" begin + g = PS.MappedFDGrid(17; halfwidth=4.0, clustering=1.0, center=1.0, domain=:half, order=4) + @test g.nodes[1] ≈ 0.0 atol = 1e-12 + @test g.nodes[end] ≈ 4.0 atol = 1e-12 + @test issorted(g.nodes) + end + + @testset "Simpson quadrature weights integrate at design order" begin + # ∫_0^4 exp(-(y-1)^2/2) dy — self-convergence (no closed form needed): + # successive-refinement differences must shrink at ≳ 4th order. + quad(n) = + let g = PS.MappedFDGrid(n; halfwidth=4.0, clustering=1.0, center=1.0, domain=:half, order=4) + sum(g.wq .* exp.(-(g.nodes .- 1) .^ 2 ./ 2)) + end + q = quad.([17, 33, 65, 129]) + d1, d2, d3 = abs(q[2] - q[1]), abs(q[3] - q[2]), abs(q[4] - q[3]) + @test log(d1 / d2) / log(2) > 3.3 + @test log(d2 / d3) / log(2) > 3.3 + end + + @testset "Gauss–Laguerre quadrature is exact on polynomials × e^{-E}" begin + gg = PS.GaussGrid(6; kind=:laguerre) + # ∫_0^∞ E^k e^{-E} dE = k! + for (k, want) in ((0, 1.0), (1, 1.0), (2, 2.0), (3, 6.0), (4, 24.0)) + @test sum(gg.weights .* gg.nodes .^ k) ≈ want rtol = 1e-10 + end + end + + @testset "IslandGrid assembles all five coordinates" begin + ig = PS.IslandGrid(; nx=17, nxi=16, ny=9, nE=4, halfwidth_x=6.0, clustering_x=1.5, + y_max=4.0, y_c=1.0, clustering_y=1.0) + @test PS.nnodes(ig) == (17, 16, 9, 4, 2) + @test ig.σ == [1.0, -1.0] + @test length(ig.x.wq) == 17 + # even nx is rejected (Simpson needs odd) + @test_throws ArgumentError PS.MappedFDGrid(16; halfwidth=6.0) + end +end diff --git a/test/runtests_islands_operators.jl b/test/runtests_islands_operators.jl new file mode 100644 index 00000000..007a3ee7 --- /dev/null +++ b/test/runtests_islands_operators.jl @@ -0,0 +1,97 @@ +# runtests_islands_operators.jl +# +# Islands operator-stack skeleton — verification ladder A1/A2 (docs/src/islands/design/05). +# A1 MMS: per-operator and assembled-system convergence at design order. +# A2 JVP (forward-mode AD) vs. finite-difference residual directional derivative. +# + allocation regression for the `apply!` / `residual!` hot paths (docs/04 §9). +# +# These are STRUCTURAL (pre-physics) checks. The manufactured coefficients they +# use are arbitrary order-unity test values — not physics — so nothing here is +# [VERIFY]-gated. Physics benchmarks (ladder B+) stay skipped until human-cleared. + +const Isl = GeneralizedPerturbedEquilibrium.Islands +const V = Isl.Verify +const PSg = Isl.PhaseSpace +const Op = Isl.Operators + +# nx = ny refined together; ξ bandlimited (nxi small is exact); nE small. +_grid(n) = PSg.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) +_order(e_coarse, e_fine, n_coarse, n_fine) = log(e_coarse / e_fine) / log(n_fine / n_coarse) + +@testset "Islands operator stack (A1/A2)" begin + + @testset "A1 — per-operator MMS convergence at design order" begin + # x/y differential operators: fourth-order finite differences. + for term in (:streaming, :exb, :collisions, :perp) + e17 = V.mms_operator_error(_grid(17), term) + e33 = V.mms_operator_error(_grid(33), term) + @test _order(e17, e33, 17, 33) > 3.3 + @test e33 < e17 + end + # ξ-only operator (magnetic drift): Fourier spectral → machine precision + # on the bandlimited manufactured ξ-profile. + @test V.mms_operator_error(_grid(17), :drift) < 1e-10 + # state-independent source is reproduced exactly (no discretization). + @test V.mms_operator_error(_grid(17), :drive) < 1e-12 + end + + @testset "A1 — assembled kinetic-residual MMS convergence" begin + e17 = V.mms_assembled_error(_grid(17)) + e33 = V.mms_assembled_error(_grid(33)) + @test _order(e17, e33, 17, 33) > 3.3 + @test e33 < 5e-3 + end + + @testset "A1 — velocity-moment quadrature convergence (refine y)" begin + mky(ny) = PSg.IslandGrid(; nx=9, nxi=8, ny=ny, nE=3, halfwidth_x=6.0, clustering_x=1.0, + y_max=4.0, y_c=1.0, clustering_y=0.8, order=4) + d1 = V.moment_selfconvergence(mky(17), mky(33)) + d2 = V.moment_selfconvergence(mky(33), mky(65)) + @test _order(d1, d2, 33, 65) > 3.3 + # grids must share (nx, nξ) + @test_throws ArgumentError V.moment_selfconvergence(mky(17), _grid(17)) + end + + @testset "A2 — AD JVP matches finite-difference directional derivative" begin + # residual is nonlinear (E×B bracket couples g and Φ), so this exercises + # the AD plumbing and its correctness together. + @test V.jvp_fd_maxerror(_grid(9)) < 1e-6 + @test V.jvp_fd_maxerror(_grid(9); seed=7) < 1e-6 + end + + @testset "allocation regression — apply!/residual! are allocation-free" begin + g = _grid(9) + for term in (:streaming, :drift, :exb, :collisions, :perp, :drive) + @test V.term_allocations(g, term) == 0 + end + @test V.residual_allocations(g) == 0 + end + + @testset "structural — state, stack, flatten/unflatten" begin + g = _grid(9) + U, = V.manufactured_state(g) + @test eltype(U) == Float64 + @test size(U.g) == PSg.nnodes(g) + @test size(U.Φ) == (9, 8) + + # flatten/unflatten round-trips + n = Op.statelength(g) + @test n == prod(PSg.nnodes(g)) + 9 * 8 + v = zeros(n) + Op.flatten!(v, U) + U2 = Op.IslandState(g) + Op.unflatten!(U2, v) + @test U2.g == U.g && U2.Φ == U.Φ + + # residual! runs and produces a finite, correctly-shaped result + stack = V.build_stack(g) + R = Op.IslandState(g) + cache = Op.IslandCache(g) + Op.residual!(R, U, stack, g, cache) + @test all(isfinite, R.g) && all(isfinite, R.Φ) + + # the magnetic-drift variant toggle is carried on the term + @test Op.MagneticDrift(zeros(1, 1, 1, 1, 1); variant=:improved).variant == :improved + end +end From d63aa74ba7b4fcff0c308f8943cf560844fe3733 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Wed, 8 Jul 2026 09:21:56 -0400 Subject: [PATCH 4/9] ISLANDS - MINOR - Docs: M1 status page, LOG entry, resolve Q1 (julia path) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - docs/src/islands.md: document the M1 skeleton (PhaseSpace/Operators/Verify), restating the structure-only [VERIFY] discipline. - LOG.md: M1 session entry (what moved / blocked / next). - QUESTIONS.md: Q1 RESOLVED — julia is at the ncl2128 software dir and must be invoked with a clean LD_LIBRARY_PATH (OMFIT contamination); used here to run the suite locally. Co-Authored-By: Claude Opus 4.8 (1M context) --- docs/src/islands.md | 26 +++++++++++++++++++++++--- docs/src/islands/LOG.md | 30 ++++++++++++++++++++++++++++++ docs/src/islands/QUESTIONS.md | 26 ++++++++++++++------------ 3 files changed, 67 insertions(+), 15 deletions(-) diff --git a/docs/src/islands.md b/docs/src/islands.md index 99c8ada6..52d78479 100644 --- a/docs/src/islands.md +++ b/docs/src/islands.md @@ -3,9 +3,29 @@ 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 (milestone M1). The design documents live under `islands/design/` in -the documentation source tree, and the module conventions in -`src/Islands/CLAUDE.md`. +development. The design documents live under `islands/design/` in the +documentation source tree, and the module conventions in `src/Islands/CLAUDE.md`. + +## Status — milestone M1 (skeleton) + +M1 lands the numerical skeleton, not the 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`. + +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. ## API Reference diff --git a/docs/src/islands/LOG.md b/docs/src/islands/LOG.md index 5c5ef305..8fb07e59 100644 --- a/docs/src/islands/LOG.md +++ b/docs/src/islands/LOG.md @@ -8,6 +8,36 @@ relevant. --- +## 2026-07-08 — M1 skeleton: phase-space grids + operator stack + MMS/AD harness + +- **Moved**: Landed the M1 core (design `03 §1–2`, `04`, ladder `A1/A2`). Three + `src/Islands/` submodules, all structure-only (no `[VERIFY]` physics numbers): + - `phasespace/PhaseSpace.jl` — the `(x, ξ, y, E, σ)` grids with layer-clustered + maps: Fourier spectral `∂ξ`, Fornberg high-order FD `∂x`/`∂y` on `sinh`-stretched + grids (window sized per-derivative so `D1`/`D2` are both 4th-order incl. + boundaries), composite-Simpson quadrature weights, Gauss–Laguerre energy nodes. + - `operators/Operators.jl` — `AbstractTerm` + `apply!` + `residual!`; the term + structs of `03 §2` (`ParallelStreaming`, `MagneticDrift` with the + `:original/:improved` toggle, `ExBDrift` as the `(x,ξ)` Poisson bracket, + `Collisions`, `GradientDrive`, `PerpTransport`/`RadiationSink` L4 stubs, + `Quasineutrality` field residual). Every physics coefficient is a **supplied + data field** — no literal in `src/`. Allocation-free, AD-generic. + - `verify/Verify.jl` — manufactured-solution + AD-vs-FD JVP harness. + - Tests `test/runtests_islands_{grids,operators}.jl` (wired into `runtests.jl`): + A1 per-operator MMS → 4th order for `∂x/∂y` terms, machine-precision for the + `∂ξ` term; assembled kinetic residual → 4th order; A2 JVP-vs-FD agree to ~6e-9; + **allocation regression = 0 bytes** for every `apply!` and `residual!`. All + 53 Islands tests green. Added `ForwardDiff` to `Project.toml` (design `04 §9`). +- **physics-verifier**: PASS — audited all six new/changed files, no + `[VERIFY]`-policy violation; the flagged literature numbers (8.73/1.46 ρ_bi, + k=−1.173, …) appear only in docstring prose, never assigned to a coefficient. +- **Blocked**: nothing. **Q1 RESOLVED**: julia is at + `/mnt/homes_global/ncl2128/software/julia-1.11.7/bin/julia`; must be run with + `env -u LD_LIBRARY_PATH` (OMFIT contamination). Used it to run the suite here. +- **Next**: M2 — wire moments (`Δ_cos`, `Δ_sin`), `frames/`, `species/`, and the + Newton–Krylov solver toward the L0 single-species solve; every physics + coefficient stays `[VERIFY]`-gated with a skipped benchmark until cleared. + ## 2026-07-08 — Harden Stop hook against OMFIT LD_LIBRARY_PATH contamination - **Moved**: Diagnosed why the Stop hook's package-load check fails on this box. diff --git a/docs/src/islands/QUESTIONS.md b/docs/src/islands/QUESTIONS.md index 77ecadf4..7592aa28 100644 --- a/docs/src/islands/QUESTIONS.md +++ b/docs/src/islands/QUESTIONS.md @@ -21,21 +21,23 @@ supervising sessions. --- -## Q1 — Julia not on the automation shell PATH — OPEN +## Q1 — Julia not on the automation shell PATH — RESOLVED (by Claude, 2026-07-08) - **Context**: Phase A bootstrap. Verifying the `Islands` module skeleton loads (`using GeneralizedPerturbedEquilibrium`) and running the test suite requires - `julia`, but it is not on the non-interactive shell's PATH (no `julia` module, + `julia`, but it was not on the non-interactive shell's PATH (no `julia` module, none in `$HOME`; other users have installs under `/mnt/homes*/…/julia*`). - **Question**: What is the canonical `julia` invocation for automation on this cluster, and will the overnight loop's scratch-clone environment expose it? - The unattended loop **cannot run `julia --project` tests without it** — this is - a hard prerequisite for Phase B, not just a local convenience. -- **Options**: (a) a `module load ` line prepended in the loop's shell - rc / launch script; (b) an absolute path to a `juliaup`/`julia` binary the - user owns; (c) install juliaup under `ncl2128` and pin 1.11. -- **Recommendation**: user provides the module/path; bake it into the tmux - launch script (and the Stop hook, which runs the fast test subset). Until - then, CI (`test.yaml`) is the only validation of Julia changes. -- **Gated work**: local verification of every Julia change; the Phase-B overnight - loop's ability to run tests / meet its definition-of-done. +- **Resolution**: the `ncl2128`-owned install is at + `/mnt/homes_global/ncl2128/software/julia-1.11.7/bin/julia` (option (b): an + absolute path to a user-owned binary), and it is on this session's PATH. The + M1 run used it to build and run `test/runtests.jl` locally. The **only caveat** + is the OMFIT `LD_LIBRARY_PATH` contamination already documented in LOG + (2026-07-08): the binary must be invoked with a clean loader path — + `env -u LD_LIBRARY_PATH /mnt/homes_global/ncl2128/software/julia-1.11.7/bin/julia + --project=. …` — or the conda libs shadow Julia's bundled artifacts. The Stop + hook already applies `env -u LD_LIBRARY_PATH`; the overnight loop's launch + script must do the same for its own gpec runs. +- **Gated work (now unblocked)**: local verification of every Julia change; the + overnight loop's ability to run tests / meet its definition-of-done. From 525155f95d4e27cbe5d9b469eaa02dca278a445c Mon Sep 17 00:00:00 2001 From: logan-nc Date: Wed, 8 Jul 2026 09:28:09 -0400 Subject: [PATCH 5/9] ISLANDS - MINOR - LOG: reference M1 PR #320 Co-Authored-By: Claude Opus 4.8 (1M context) --- docs/src/islands/LOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/src/islands/LOG.md b/docs/src/islands/LOG.md index 8fb07e59..3c7dfba9 100644 --- a/docs/src/islands/LOG.md +++ b/docs/src/islands/LOG.md @@ -10,6 +10,7 @@ relevant. ## 2026-07-08 — M1 skeleton: phase-space grids + operator stack + MMS/AD harness +- **PR**: #320 (`feature/islands-m1` → `feature/islands`); full suite green. - **Moved**: Landed the M1 core (design `03 §1–2`, `04`, ladder `A1/A2`). Three `src/Islands/` submodules, all structure-only (no `[VERIFY]` physics numbers): - `phasespace/PhaseSpace.jl` — the `(x, ξ, y, E, σ)` grids with layer-clustered From db9a8f11ef93dcc9fa84bc108f8f8d2cf393dfa8 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Wed, 8 Jul 2026 12:57:39 -0400 Subject: [PATCH 6/9] ISLANDS - MINOR - Author M2 launch prompt + seed clearance queue MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Set up the next milestone (M2 — L0 solve machinery) for the autonomous loop, mirroring the M1 launch-prompt bootstrap. - docs/src/islands/design/M2-launch-prompt.md: the M2 milestone contract. Scope = full L0 solve machinery (solvers/, moments/, frames/, fields/, species/, neoclassical-matching BCs) as AD-compatible allocation-free structure with every physics coefficient a [VERIFY]-gated parameter. DoD = the physics-free structural gates (A5 null, assembled solve-MMS, A8 y_c conditioning, A4 conservation, A3 parity, A7 ⟨∂²h/∂x²⟩=0) + suite + PR + Paper-I OUTLINE + the clearance queue. The York gates (B5a/b/c, B2, B4) are explicitly OUT — they need human-cleared physics; reaching one by hardcoding is a policy violation, not completion. - QUESTIONS.md: seed the parallel-human clearance queue — Q2 (ratify D7/D8), Q3 (clear the L0 [CHECKED] coefficient set with exact cites), Q4 (open [VERIFY]s: the psi-tilde q_s'/q_s typo, B5a collisionality, acquire WCHH96 + Park 2022). Co-Authored-By: Claude Opus 4.8 (1M context) --- docs/src/islands/QUESTIONS.md | 68 ++++++++ docs/src/islands/design/M2-launch-prompt.md | 174 ++++++++++++++++++++ 2 files changed, 242 insertions(+) create mode 100644 docs/src/islands/design/M2-launch-prompt.md diff --git a/docs/src/islands/QUESTIONS.md b/docs/src/islands/QUESTIONS.md index 7592aa28..9c6c6df7 100644 --- a/docs/src/islands/QUESTIONS.md +++ b/docs/src/islands/QUESTIONS.md @@ -41,3 +41,71 @@ supervising sessions. script must do the same for its own gpec runs. - **Gated work (now unblocked)**: local verification of every Julia change; the overnight loop's ability to run tests / meet its definition-of-done. + +## Q2 — Ratify Decisions D7 and D8 — OPEN + +- **Context**: M2 setup. The L0 equation set and benchmark targets rest on two + Decision-Log proposals (`docs/00`) that are dated 2026-07-07 and still marked + "needs human ratification". Until ratified, M2 cannot fix the L0 physics form or + pin its benchmark tolerances. +- **Question**: Ratify (or amend) **D7** — implement Level-0 physics from an + *independent re-derivation* cross-checked against the L23-amended equation set, + treating I19 Eq. (A.1) as printed as known-errata (L23 §2.6), with ω_E a scanned + input from day one — and **D8** — pin the benchmark grid as the three-code + triangle (DK-NTM / RDK-NTM / kokuchou) with B5a/b/c configs, superseding the + single "York thresholds" gate item. +- **Options**: (a) ratify both as written; (b) ratify with amendments; (c) direct a + different L0 derivation strategy. +- **Recommendation**: ratify both — they encode the project's core anti-guessing + thesis (published O(1) coefficients in this lineage are demonstrably wrong) and + the benchmark structure the ladder already assumes. +- **Gated work**: pinning any L0 physics coefficient; the York gates B5a/b/c; the + Paper-I physics claims. + +## Q3 — Clear the Level-0 coefficient set ([CHECKED] → human sign-off) — OPEN + +- **Context**: M2 builds the L0 solve machinery with every physics coefficient a + parameterized `[VERIFY]` stub. Reaching the York gates needs these `[CHECKED]` + (AI-transcribed, cited, but not human-signed-off) items cleared. `[CHECKED]` is + not permission to hardcode (`docs/01` header). Each is implemented structurally; + none has a value in `src/`. +- **Question**: Check each against its source PDF and sign off (record paper + Eq./p. + in `docs/01`), or flag a discrepancy: + - `ω̂_D` magnetic-drift frequency + the `:original`/`:improved` `L̂_B⁻¹` toggle + `[CHECKED: I19 Eq. (32) def.; D21 Eqs. 15, B1; D21 Eq. A2, p. 16]` — drives the + 8.73 → 1.46 ρ_bi threshold shift. + - Pitch-angle collision kernel `[CHECKED: I19 Eqs. (9)–(12); Diss19 Eqs. 2.25–2.30; + WCHH96 Eq. (62)]` + the analytic velocity average + `⟨ν̂_ii⟩_u = (4ε^{3/2}ν_★/3√π)(√2 − ln(1+√2))` `[CHECKED: L23 Eq. 4.1.6, p. 88]` + + normalization `ν_★ = ν_jj Rq/(ε^{3/2}v_th)` `[CHECKED: L23 Eq. (2.3.40)]`. + - Electron-closure constants `k ≃ −1.173` (Hirshman–Sigmar) and `f_p ≃ 1 − 1.46√ε` + `[CHECKED: I19 Eq. (22); L23 Eqs. 2.5.5–2.5.8]`. + - Quasineutrality closure `e_iΦ̂/T_i = [δn̄_i/n₀ + x − ĥ(Ω)]/(2 L̂_{n0})` + `[CHECKED: I19 Eq. (A.11); L23 Eq. (2.4.14)]` and the Picard form + `δΦ̂ = (δn̂_i − δn̂_e)/2` `[CHECKED: Diss19 Eq. 2.45]`. +- **Options**: (a) clear item-by-item after PDF check; (b) require an independent + re-derivation first (couples to Q2/D7) before clearing. +- **Recommendation**: clear via re-derivation (per D7) rather than transcription + alone — L23 §2.6 documents concrete errors in the published set. +- **Gated work**: populating any of these coefficients; the York/large-w/polarization + gates (B5a/b/c, B2, B4); the A7 number-bearing identities (`k`, `f_p`, `⟨ν̂_ii⟩`). + +## Q4 — Resolve open [VERIFY]s and acquire two missing sources — OPEN + +- **Context**: Independent of the `[CHECKED]` clearances above, three genuinely + open `[VERIFY]` items block pinning M2's moment normalization and B5a tolerance. +- **Question**: + - `ψ̃` amplitude: is it `(w_ψ²/4)(q_s′/q_s)` (Diss19/D21/L23, dimensional analysis) + or `(w_ψ²/4)(q_s/q_s′)` (one I19 extraction)? `[VERIFY: check I19 as printed — + possible typo in the paper itself]` — sets the `Δ_cos/Δ_sin` prefactor. + - B5a run collisionality: I19 §4.2 states `ν_★ = 0.01`; L23 p. 82 quotes DK-NTM at + `ν_★ = 10⁻³`. Resolve before pinning the B5a tolerance. + - Acquire **WCHH96** (analytic electron closure / large-w limits, B2) and **Park + PoP 29 (2022)** (SLAYER Q-convention map, D1) into the `docs/08` reference + library — both are currently cited only via transcription. +- **Options**: (a) user resolves from the source PDFs; (b) an independent + re-derivation pins `ψ̃` and the collisionality is read from the reproduced run. +- **Recommendation**: resolve `ψ̃` by re-derivation (it is a clean dimensional check); + read the B5a collisionality from whichever run the D8 triangle pins as canonical. +- **Gated work**: the `Δ_cos/Δ_sin` normalization; the B5a threshold tolerance; the + B2 large-w and D1 SLAYER comparisons. diff --git a/docs/src/islands/design/M2-launch-prompt.md b/docs/src/islands/design/M2-launch-prompt.md new file mode 100644 index 00000000..530358d2 --- /dev/null +++ b/docs/src/islands/design/M2-launch-prompt.md @@ -0,0 +1,174 @@ +# M2 launch prompt (Islands overnight autonomous run) + +> This file is fed verbatim to the overnight loop: +> `claude --permission-mode dontAsk --continue -p "$(cat docs/src/islands/design/M2-launch-prompt.md)"`. +> It is the milestone contract (design doc `06 §2.1`). Keep it stable across +> relaunches so `--continue` resumes the same objective. + +You are working milestone **M2** of the Islands module (`src/Islands/`), a +steady-state drift-kinetic island/layer solver, autonomously and unattended. M1 +(phase-space grids + operator-stack skeleton + MMS/AD harness) is complete +(PR #320); M2 builds the **Level-0 solve machinery** on top of it. + +## Read first (every session) + +`src/Islands/CLAUDE.md` (module conventions + the [VERIFY] policy), +`docs/src/islands/LOG.md` and `docs/src/islands/QUESTIONS.md` (session memory + +open blockers), then the design docs +`docs/src/islands/design/{00-roadmap,01-physics-level0,02-species-and-eps,03-architecture,04-numerics,05-verification,06-autonomy,07-documentation-and-papers}.md` +(read by number `NN`; the file globs match). The repo-root `CLAUDE.md` governs +GPEC-wide conventions. + +## The gating reality (read this before you plan) + +M2's roadmap headline is "L0 single-species solve, Δ moments, **York gates** → +Paper I." **The York gates are NOT reachable in this run** and reaching one is a +policy violation, not completion. Every L0 physics coefficient the gates need +(the `ω̂_D`/`L̂_B` toggle, the collision kernel, `k=−1.173`, `f_p=1−1.46√ε`, +`⟨ν̂_ii⟩_u`, the quasineutrality closure, the York numbers `8.73`/`1.46` ρ_bi) is +at most **`[CHECKED]`** — none is human-cleared — and Decisions **D7** (implement +L0 from an independent re-derivation vs the L23-amended set) and **D8** (pin the +B5a/b/c benchmark triangle) are **unratified** (`docs/00` Decision Log). `[CHECKED]` +is *not* permission to hardcode a number (`docs/01` header; `CLAUDE.md` policy). + +So M2 here = **build the full L0 solve *structure* with every physics coefficient +a `[VERIFY]`-gated parameter, close the physics-free structural gates, and escalate +a prioritized clearance queue** for the human. A thin follow-up run fills the +numbers and hits the York gates *after* a human clears the tags. + +## Goal (set this as your `/goal` completion condition) + +**M2 is done when all of the following hold:** + +1. **The L0 solve machinery exists** as AD-compatible, allocation-free structure, + with **every physics coefficient a supplied `[VERIFY]`-gated parameter (no + literal in `src/`)** — new submodules per design `03 §1`: + - `solvers/` — matrix-free **Newton–Krylov** (GMRES on the M1 ForwardDiff JVP; + inexact-Newton / Eisenstat–Walker forcing; line search; a pseudo-arclength + continuation scaffold that detects folds from day one; a physics-block + preconditioner that treats the `y_c` matching block explicitly), `04 §3, §5`. + - `moments/` — `Δ_cos`/`Δ_sin` assembly: species-charge-weighted velocity moment + of `g` → `J̄_∥(x, ξ)` → the `∮ J̄_∥ {cos, sin} ξ dξ ∫ dψ` projections (`01 §4`), + plus the bootstrap/polarization channel decomposition *structure*. The `μ₀R/2ψ̃` + prefactor and `ψ̃` stay **gated** (open `[VERIFY]`); the projection + quadrature + is pure numerics. + - `frames/` — THE `ω`/normalization conversion module (`01 §5`): the conversion + *forms* with every sign and normalization flagged and gated. This module owns + all `ω`-sign conventions (the polarization sign disputes are frame disputes — + do not reproduce them elsewhere). Unit-test only the mechanical, sign-free + identities. + - `fields/` — formalize the quasineutrality residual (already stubbed in + `operators/`); the flattened-electron closure *structure* gated. + - `species/` — `Species`, `AbstractBackground`, `SpeciesRole {Bulk, Trace}` + plumbing (`02 §1`, Decision D3) — pure data structures, fully unblocked. The L0 + test is a single bulk ion + a trace-deuterium copy. + - **neoclassical-matching far-field BCs** (`01 §3`, `04 §1`): structure only; the + no-island neoclassical far-field solution is gated. **Never bare Neumann** (the + L23 "winged" spurious-branch failure). + +2. **Physics-free structural gates green** (`05 §A`), in `test/runtests_islands_*.jl`, + passing in the full suite: + - **A5** zero-drive null: gradients and `Φ̃` off ⇒ `g ≡ 0` exactly, residual = + machine zero, Newton converges trivially. + - **Assembled solve-MMS**: a manufactured `g*` with `GradientDrive` set so `g*` is + the exact Newton solution; the solver recovers `g*` at design order (extends + M1's residual-MMS to a full converged solve). + - **A8** `y_c` matching-block smallest-singular-value conditioning monitor + (regression = the silent-noise failure mode of L23 §4.2). + - **A4** (L0): particle conservation + discrete entropy sign `∫ g C[g]/F_M ≤ 0` + of the discretized collision operator (structural — holds for any valid + discretization, independent of the physical `ν` value). + - **A3** parity: `Δ_cos` even / `Δ_sin` odd under the appropriate `(ξ, σ, ω_E)` + reflection, tested on a manufactured `J̄_∥` (coefficient-free). + - **A7** the coefficient-free identity `⟨∂²h/∂x²⟩_Ω = 0` (defer the number-bearing + `k`, `f_p`, `⟨ν̂_ii⟩` identities to a post-clearance run). + +3. The full suite passes: `julia --project=. test/runtests.jl`. + +4. **`docs/src/islands/papers/paper-1/OUTLINE.md`** written as the Paper-I figure + contract (`07 §3`): claims → figures → ladder IDs (B5a/b/c, B2, B4). This is the + level-*start* deliverable, not the level-*gate*. + +5. **The clearance queue** (the parallel-human deliverable): a consolidated, + prioritized set of `QUESTIONS.md` entries covering every coefficient the York + gates need + D7/D8 ratification + the open `[VERIFY]`s, **each paired with a + skipped B-series benchmark** in `benchmarks/islands/` that references its tag — so + the human can clear in parallel and the follow-up run fills numbers and hits gates. + +6. A PR is open onto `feature/islands` (a sub-branch → `feature/islands`, as M1's + PR #320 did). + +7. **No `[VERIFY]`/uncleared-`[CHECKED]` coefficient has been assigned a value** + without a `docs/src/islands/QUESTIONS.md` entry. + +**Explicitly NOT in the M2 DoD (gated on human clearance — do not attempt):** the +York gates B5a/b/c, B2, B4; B1 (needs an external NEO/NCLASS run); any physics +number. A "green York gate" reached by hardcoding is the exact failure this project +exists to prevent. + +## Scope + +- **Reuse existing GPEC machinery; do not reimplement** (repo-root CLAUDE.md + minimal-change discipline): `FastInterpolations.integrate` / + `cumulative_integrate` for the ψ/flux integrals (`src/Equilibrium/Equilibrium.jl`, + `src/ForceFreeStates/Resist.jl`); the `src/KineticForces/BounceAveraging.jl` + velocity-space λ-averaging pattern for the `(y, E, σ)` moments; the allocation-free + `QuadGK.quadgk!` pattern (`src/ForceFreeStates/Galerkin/GalerkinAssembly.jl`) for + hot quadrature; the `EquilibriumConfig` / `build_inputs_from_toml` TOML-config + pattern (`src/Equilibrium/EquilibriumTypes.jl`, + `src/GeneralizedPerturbedEquilibrium.jl`) for an `[Islands]` `gpec.toml` section in + a new `io/`. +- **Add `Krylov.jl`** to `Project.toml` `[deps]` + `[compat]` (matrix-free GMRES; + `04 §9` names Krylov.jl / LinearSolve.jl — prefer Krylov.jl: lighter and + purpose-built). The JVP is the M1 ForwardDiff operator; form **no** global sparse + Jacobian except in a tiny-grid debug mode. +- Create `benchmarks/islands/` (+ `figures/`) and at least one `islands_*` regression + case integrated with `regression-harness/` — every physics benchmark ships + **skipped**, referencing its `[VERIFY]`/`QUESTIONS` id. + +## The hard rule (non-negotiable) + +**Never assign a value to a `[VERIFY]` or uncleared-`[CHECKED]` physics coefficient, +sign, or normalization.** The moment you would need a specific number/sign from the +literature that isn't human-cleared: (a) implement the *structure* with the +coefficient as a named parameter, (b) add a skipped benchmark referencing the tag, +(c) write a `QUESTIONS.md` entry (context / question / options / recommendation / +gated work), and (d) **switch to the next unblocked task**. Guessing a coefficient +is the exact failure this project exists to prevent. Manufactured, order-unity test +coefficients in the MMS/verification harness are legitimate (they test numerics, not +physics) and carry no tag — do not confuse the two. + +## Working discipline + +- Before committing any physics-adjacent change, run the **`physics-verifier`** + subagent; if it returns BLOCK, fix or escalate — do not commit. +- **Run every new kernel once under `--check-bounds=yes`** before trusting it. (M1 + lesson: an `@inbounds` index-swap corrupted memory and passed silently until a + forced bounds-checked run caught it. Assume new nested-index kernels are guilty + until a bounds-checked run clears them.) +- Invoke julia with a clean loader path: + `env -u LD_LIBRARY_PATH /mnt/homes_global/ncl2128/software/julia-1.11.7/bin/julia + --project=. …` (Q1 resolution: the OMFIT `LD_LIBRARY_PATH` shadows Julia's + artifacts otherwise). +- Commit granularly to a milestone sub-branch with messages in the repo format + (`ISLANDS - - `); reference `QUESTIONS.md`/ladder IDs where relevant. + Push only Islands branches (never `develop`/`main`; the hooks enforce this). +- Never weaken a tolerance or re-baseline a target to reach "done" — that is a + blocker, not a fix. +- Append a `LOG.md` entry (what moved / blocked / next) before ending each session. + The `Stop` hook will keep you from ending with a dirty tree or a broken build. A + session that ends without a pushed branch and a status note has failed its exit + criteria (`06 §2.5`). +- If you exhaust the milestone's unblocked work (everything remaining is gated on + `QUESTIONS.md`/human clearance), commit, log, and let the session end — the outer + loop and the human pick it up. + +## Definition of NOT done (do not stop early) + +Do not declare M2 done if: any structural gate (A5, solve-MMS, A8, A4, A3, A7's +`⟨∂²h/∂x²⟩=0`) is failing or skipped-as-a-shortcut; the L0 machinery submodules are +absent; the suite is red; the tree is dirty; the Paper-I OUTLINE or the clearance +queue is missing; or any coefficient was guessed. Those are blockers, not +completion. Conversely, do **not** keep working past a green structural ladder in an +attempt to reach the York gates — those are gated on human clearance and out of +scope for this run. From 42fffb62fe8857ff1be30ad2770458543a27668c Mon Sep 17 00:00:00 2001 From: logan-nc Date: Wed, 8 Jul 2026 12:58:18 -0400 Subject: [PATCH 7/9] ISLANDS - MINOR - Fix design-doc filename in M2 prompt read-first list Co-Authored-By: Claude Opus 4.8 (1M context) --- docs/src/islands/design/M2-launch-prompt.md | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/docs/src/islands/design/M2-launch-prompt.md b/docs/src/islands/design/M2-launch-prompt.md index 530358d2..58a50184 100644 --- a/docs/src/islands/design/M2-launch-prompt.md +++ b/docs/src/islands/design/M2-launch-prompt.md @@ -15,9 +15,8 @@ steady-state drift-kinetic island/layer solver, autonomously and unattended. M1 `src/Islands/CLAUDE.md` (module conventions + the [VERIFY] policy), `docs/src/islands/LOG.md` and `docs/src/islands/QUESTIONS.md` (session memory + open blockers), then the design docs -`docs/src/islands/design/{00-roadmap,01-physics-level0,02-species-and-eps,03-architecture,04-numerics,05-verification,06-autonomy,07-documentation-and-papers}.md` -(read by number `NN`; the file globs match). The repo-root `CLAUDE.md` governs -GPEC-wide conventions. +`docs/src/islands/design/{00-roadmap,01-physics-level0,02-species-and-eps,03-architecture,04-numerics,05-verification,06-autonomy-and-tooling,07-documentation-and-papers}.md`. +The repo-root `CLAUDE.md` governs GPEC-wide conventions. ## The gating reality (read this before you plan) From 7e6bcf6dd5591bdc8ebf20f864feb4b8a92ad80d Mon Sep 17 00:00:00 2001 From: logan-nc Date: Wed, 8 Jul 2026 13:37:24 -0400 Subject: [PATCH 8/9] ISLANDS - MINOR - Reframe M2 prompt as an interactive /goal contract MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Drop the overnight-loop framing (dontAsk/--continue header, unattended language, Stop-hook exit-criteria) — M2 is driven interactively via /goal in auto mode. Blocked coefficients are surfaced to the user to clear live rather than only parked in QUESTIONS.md. The contract (scope, DoD, [VERIFY] hard rule, gated York gates) is unchanged. Co-Authored-By: Claude Opus 4.8 (1M context) --- docs/src/islands/design/M2-launch-prompt.md | 30 +++++++++------------ 1 file changed, 13 insertions(+), 17 deletions(-) diff --git a/docs/src/islands/design/M2-launch-prompt.md b/docs/src/islands/design/M2-launch-prompt.md index 58a50184..2b406ce1 100644 --- a/docs/src/islands/design/M2-launch-prompt.md +++ b/docs/src/islands/design/M2-launch-prompt.md @@ -1,14 +1,12 @@ -# M2 launch prompt (Islands overnight autonomous run) +# M2 milestone contract (Islands) -> This file is fed verbatim to the overnight loop: -> `claude --permission-mode dontAsk --continue -p "$(cat docs/src/islands/design/M2-launch-prompt.md)"`. -> It is the milestone contract (design doc `06 §2.1`). Keep it stable across -> relaunches so `--continue` resumes the same objective. +> The milestone contract for Islands M2 (design doc `06 §2.1`). Set it as the +> `/goal` completion condition when working M2. You are working milestone **M2** of the Islands module (`src/Islands/`), a -steady-state drift-kinetic island/layer solver, autonomously and unattended. M1 -(phase-space grids + operator-stack skeleton + MMS/AD harness) is complete -(PR #320); M2 builds the **Level-0 solve machinery** on top of it. +steady-state drift-kinetic island/layer solver. M1 (phase-space grids + +operator-stack skeleton + MMS/AD harness) is complete (PR #320); M2 builds the +**Level-0 solve machinery** on top of it. ## Read first (every session) @@ -132,8 +130,9 @@ sign, or normalization.** The moment you would need a specific number/sign from literature that isn't human-cleared: (a) implement the *structure* with the coefficient as a named parameter, (b) add a skipped benchmark referencing the tag, (c) write a `QUESTIONS.md` entry (context / question / options / recommendation / -gated work), and (d) **switch to the next unblocked task**. Guessing a coefficient -is the exact failure this project exists to prevent. Manufactured, order-unity test +gated work), and (d) **ask the user to clear it if they are available, otherwise +move to the next unblocked task**. Guessing a coefficient is the exact failure this +project exists to prevent. Manufactured, order-unity test coefficients in the MMS/verification harness are legitimate (they test numerics, not physics) and carry no tag — do not confuse the two. @@ -154,13 +153,10 @@ physics) and carry no tag — do not confuse the two. Push only Islands branches (never `develop`/`main`; the hooks enforce this). - Never weaken a tolerance or re-baseline a target to reach "done" — that is a blocker, not a fix. -- Append a `LOG.md` entry (what moved / blocked / next) before ending each session. - The `Stop` hook will keep you from ending with a dirty tree or a broken build. A - session that ends without a pushed branch and a status note has failed its exit - criteria (`06 §2.5`). -- If you exhaust the milestone's unblocked work (everything remaining is gated on - `QUESTIONS.md`/human clearance), commit, log, and let the session end — the outer - loop and the human pick it up. +- Append a `LOG.md` entry (what moved / blocked / next) before wrapping up a working + session, and keep the tree clean and the build green (`06 §2.5`). +- If the remaining work is all gated on human clearance, surface the blockers to the + user (they can clear a tag or ratify D7/D8 live) rather than stalling. ## Definition of NOT done (do not stop early) From 0780eb87f286cea6a07859021f9b654504dcc531 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Wed, 8 Jul 2026 15:31:35 -0400 Subject: [PATCH 9/9] ISLANDS - BUG FIX - Cover Islands submodule docstrings in docs (checkdocs missing_docs) Documenter's checkdocs=:exports requires every exported docstring in the package (including submodules) to appear in a page; the islands.md @autodocs block only listed the top-level Islands module, so the PhaseSpace/Operators/ Verify exports failed the docs build. Add per-submodule @autodocs sections. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_016pUwoFRo9a5AjF8HuwcZMr --- docs/src/islands.md | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/docs/src/islands.md b/docs/src/islands.md index 52d78479..3317325a 100644 --- a/docs/src/islands.md +++ b/docs/src/islands.md @@ -32,3 +32,21 @@ their `[VERIFY]` tags are cleared. ```@autodocs Modules = [GeneralizedPerturbedEquilibrium.Islands] ``` + +### Phase-space grids (`Islands.PhaseSpace`) + +```@autodocs +Modules = [GeneralizedPerturbedEquilibrium.Islands.PhaseSpace] +``` + +### Operator stack (`Islands.Operators`) + +```@autodocs +Modules = [GeneralizedPerturbedEquilibrium.Islands.Operators] +``` + +### Verification harness (`Islands.Verify`) + +```@autodocs +Modules = [GeneralizedPerturbedEquilibrium.Islands.Verify] +```