Stride deltacoil

.claude/agent-memory/fortran-physics-reviewer/MEMORY.md

@@ -4,3 +4,4 @@
 - [Known Physics Issues in PerturbedEquilibrium](known_issues_perteq.md) — xss_mn not implemented, xi_modes placeholders
 - [KineticForces (NTV) Audit Checklist](kinetic_ntv_map.md) — pentrc->KineticForces map, Logan 2015 matrices, what to verify
 - [InnerLayer (Resistive) Audit Checklist](resistive_layer_map.md) — rmatch->InnerLayer map, GGJ Wasow basis / Δ′, what to verify
+- [Galerkin Δ′ Assembly Map](galerkin_assembly_map.md) — gal.f<->GalerkinAssembly.jl; resonant sign chain verified; PASS

.claude/agent-memory/fortran-physics-reviewer/galerkin_assembly_map.md

@@ -0,0 +1,37 @@
+---
+name: galerkin-assembly-map
+description: RDCON singular Galerkin Δ′ solver — gal.f ↔ GalerkinAssembly.jl element-level weak-form correspondence and audit notes
+metadata:
+  type: project
+---
+
+# RDCON Galerkin Δ′ assembly: gal.f ↔ GalerkinAssembly.jl
+
+File: `src/ForceFreeStates/Galerkin/GalerkinAssembly.jl` ↔ `GPEC_fortran/rdcon/gal.f`.
+On branch `feature/rdcon-galerkin`. Reviewed 2026-06-16.
+
+## Correspondence (verified, not just docstring-trusted)
+- gal_hermite (222-262) ↔ gal_hermite: pb/qb EXACT.
+- gal_get_fkg (1013-1071) ↔ gal_get_fkg: F=Fbar·sf[i]·sf[j], K=K·sf[i], G=G. singfac=mlow-nn*q+(0..mpert-1) = DIRECT (m-nq), not reciprocal. Correct.
+- gal_gauss_quad (1079-1138) ↔ gal_gauss_quad!: F·qq + K·qp + conj(K[jpert,ipert])·pq + G·pp. swap pb(2)/pb(3) (0-based) = Julia pb[3]/pb[4]. Correct.
+- gal_extension (1146-1322) ↔ gal_extension!: emat/ediag/rhs/erhs all match incl surface terms 1277-1300 and w1/w2.
+- gal_resonant! ↔ gal_lsode_int(780-900)+gal_lsode_der(908-954)+gal_make_arrays signs(1366-1369). QuadGK replaces LSODE.
+- gal_assemble_mat (694-772) ↔ gal_assemble_mat!: LU offset=kl+ku+1, chol offset=1, conj(e) Hermitian placement. Correct.
+- gal_assemble_rhs (642-686) ↔ gal_assemble_rhs!: isol++ at ext2-left / res-right. Correct.
+- gal_set_boundary (1554-1620) ↔ gal_set_boundary!: axis idx, edge 3-way (rpec/vac/fixed). DOF index 3(0-based)=Julia 4.
+
+## Resonant sign chain (the trap — verified correct)
+Fortran: LSODE integrates x0→x_lsode giving raw u; then line 882 `IF right: u=-u`; then make_arrays:
+erhs=-u_res(1), ediag=+u_res(2), rhs=-u_hermite1, emat=+u_hermite2.
+Net with s=+1(left)/-1(right): erhs=-s·res1, ediag=s·res2, rhs=-s·hbig, emat=s·hsmall.
+Julia sgn_u = s (=-1 right/+1 left): erhs=-sgn_u·raw1, ediag=sgn_u·raw2, rhs=-sgn_u·hbig, emat=sgn_u·hsmall. MATCHES.
+QuadGK limits x0→x_lsode preserve LSODE orientation so raw == pre-negation u. Correct.
+du_res conj-weight uses the SMALL solution (isol(2)) in both. Correct.
+
+## Only open item (MINOR, out-of-file)
+sing_get_ua_gal/dua_gal two-sided z<0 → sig=-1 series at -z with derivative ×(-1) is a REIMPLEMENTATION
+of sing.f's internal singular-surface handling, not a line port. The ×(-1) chain is internally consistent
+but worth a dedicated check against sing.f sing_get_ua/dua if Δ′ ever disagrees with Fortran. Lives outside
+GalerkinAssembly.jl.
+
+## Verdict: PASS WITH REQUIRED ANNOTATIONS (citation comments only; no physics defects found).

benchmarks/compare_gal_vs_el.jl

@@ -0,0 +1,78 @@
+# Overlay the IDEAL gal matched ξ(ψ) against the EL total-energy eigenmode ξ(ψ), for the same edge
+# eigenvector w. The ideal gal solution should reproduce the EL (DCON) ideal solution.
+#
+# EL :  ξ_EL(ψ)  = U_EL(ψ) · (U_EL_edge \ w)          (fundamental matrix, integration/xi_psi)
+# gal:  ξ_gal(ψ) = U_gal(ψ) · w                       (identity-at-edge ⇒ coefficient is w itself)
+# w = eigenvector of the total energy operator W = W_plasma + W_vacuum (FreeBoundaryStability/XiNorm).
+#
+# Needs ONE gpec.h5 from a run with populate_dense_xi=true (EL dense u_store), gal_match_flag=true,
+# gal_ideal_flag=true (gal ideal matched), vac_flag=true (energy operator).
+# Usage: julia --project=. benchmarks/compare_gal_vs_el.jl [gpec.h5] [out.png] [mode]
+#   mode = "highest" (default, most stable), "lowest" (most unstable), or an integer eigenmode index.
+
+using HDF5, LinearAlgebra, Plots, Printf
+
+h5path = length(ARGS) >= 1 ? ARGS[1] : "/tmp/gal_ideal_test/gpec.h5"
+outpng = length(ARGS) >= 2 ? ARGS[2] : joinpath(@__DIR__, "compare_gal_vs_el.png")
+ksel = length(ARGS) >= 3 ? ARGS[3] : "highest"
+
+to_c(a) = eltype(a) <: Complex ? ComplexF64.(a) : map(x -> ComplexF64(x.re, x.im), a)
+
+et, wt, u1, psiE, gxi, psiG, issing, mlow, sing_psi = h5open(h5path) do f
+    (to_c(read(f["FreeBoundaryStability/XiNorm/eigenmode_energies"])),
+        to_c(read(f["FreeBoundaryStability/XiNorm/W_freeboundary_eigenmodes"])),
+        to_c(read(f["integration/xi_psi"])), read(f["integration/psi"]),
+        to_c(read(f["galerkin/match/xi"])), read(f["galerkin/solution/psi"]),
+        Bool.(read(f["galerkin/solution/issing"])), read(f["info/mlow"]),
+        read(f["galerkin/sing_psi"]))
+end
+
+mpert = size(u1, 1)
+k = ksel == "highest" ? argmax(real.(et)) :
+    ksel == "lowest" ? argmin(real.(et)) :
+    parse(Int, ksel)
+w = wt[:, k]
+@printf("comparing eigenmode k=%d   δW = %.4e %+.4ei\n", k, real(et[k]), imag(et[k]))
+
+# EL profile
+cEL = u1[:, :, end] \ w
+xiE = reduce(hcat, (u1[:, :, ip] * cEL for ip in 1:size(u1, 3)))   # (mpert, nE)
+
+# gal ideal profile (identity-at-edge ⇒ coefficient = w); drop on-surface points
+keep = .!issing
+psiGk = psiG[keep]
+xiG = reduce(hcat, (gxi[:, ip, :] * w for ip in findall(keep)))    # (mpert, nGk)
+
+ms = mlow .+ (0:mpert-1)
+peak = [maximum(abs, @view xiE[i, :]) for i in 1:mpert]
+order = sortperm(peak; rev=true)
+ndom = min(5, mpert)
+
+# quantitative bulk agreement on the dominant harmonic (linear-interp EL onto gal grid, off rationals/edge)
+lininterp(xq, x, y) = (j = clamp(searchsortedlast(x, xq), 1, length(x) - 1);
+t = (xq - x[j]) / (x[j+1] - x[j]); y[j] * (1 - t) + y[j+1] * t)
+idom = order[1]
+inbulk(p) = 0.1 <= p <= 0.92 && all(abs(p - ps) > 8e-3 for ps in sing_psi)
+sel = [ip for ip in eachindex(psiGk) if inbulk(psiGk[ip])]
+num = sum(abs2(xiG[idom, ip] - lininterp(psiGk[ip], psiE, xiE[idom, :])) for ip in sel)
+den = sum(abs2(xiG[idom, ip]) for ip in sel)
+@printf("dominant m=%d:  bulk relative L2 ‖gal−EL‖/‖gal‖ = %.3e  (n=%d, off-layer/off-edge)\n",
+    ms[idom], sqrt(num / den), length(sel))
+
+plt = plot(; size=(1000, 650), xlabel="ψ_N", ylabel="|ξ^ψ_m(ψ)|",
+    title=@sprintf("Ideal gal vs EL total-energy eigenmode  (δW=%.2e)  — solid: EL, dashed: gal", real(et[k])),
+    legend=:topleft, left_margin=12Plots.mm, bottom_margin=5Plots.mm, right_margin=4Plots.mm)
+palette = theme_palette(:auto)
+for (c, i) in enumerate(order[1:ndom])
+    col = palette[mod1(c, length(palette))]
+    plot!(plt, psiE, abs.(@view xiE[i, :]); color=col, lw=2.2, label="EL m=$(ms[i])")
+    plot!(plt, psiGk, abs.(@view xiG[i, :]); color=col, lw=1.6, ls=:dash, label="gal m=$(ms[i])")
+end
+for (j, ps) in enumerate(sing_psi)
+    vline!(plt, [ps]; color=:black, ls=:dot, lw=1, label=(j == 1 ? "rational q" : ""))
+end
+
+display(plt)
+savefig(plt, outpng)
+println("dominant harmonics (m): ", [ms[i] for i in order[1:ndom]])
+println("saved: ", abspath(outpng))

benchmarks/compare_jbgradpsi_m2.jl

@@ -0,0 +1,61 @@
+# Compare the area-normalized b^ψ (perturbed_equilibrium/response/psi_area = b^ψ/⟨J·|∇ψ|⟩_θ) for one
+# poloidal harmonic between two GPEC runs that are identical except for which ξ feeds PerturbedEquilibrium:
+#   (1) IDEAL galerkin matched ξ   (gal_match_flag=true, gal_ideal_flag=true)
+#   (2) SHOOTING ξ                  (gal_match_flag=false)
+#
+# PE writes no ψ grid, so it's reconstructed: gal-ideal → galerkin/solution/psi minus issing points;
+# shooting → integration/psi.
+# Usage: julia --project=. benchmarks/compare_jbgradpsi_m2.jl [gal_h5] [shoot_h5] [out.png] [m]
+
+using HDF5, Plots, Printf
+
+gal_h5 = length(ARGS) >= 1 ? ARGS[1] : "/tmp/gal_ideal_test/gpec.h5"
+sh_h5 = length(ARGS) >= 2 ? ARGS[2] : "/tmp/shooting_test/gpec.h5"
+outpng = length(ARGS) >= 3 ? ARGS[3] : joinpath(@__DIR__, "cmp_jbgradpsi_m2.png")
+mtarget = length(ARGS) >= 4 ? parse(Int, ARGS[4]) : 2
+
+to_c(a) = eltype(a) <: Complex ? ComplexF64.(a) : map(x -> ComplexF64(x.re, x.im), a)
+
+# gal-ideal run: PE grid = gal solution grid with the on-surface (issing) points dropped
+pa_g, psi_g, mlow, sing_psi, sing_m = h5open(gal_h5) do f
+    pa = to_c(read(f["perturbed_equilibrium/response/psi_area"]))   # [npsi, mpert]
+    iss = Bool.(read(f["galerkin/solution/issing"]))
+    (pa, read(f["galerkin/solution/psi"])[.!iss], read(f["info/mlow"]),
+        read(f["galerkin/sing_psi"]), read(f["galerkin/sing_m"]))
+end
+# shooting run: PE grid = integration/psi
+pa_s, psi_s = h5open(sh_h5) do f
+    (to_c(read(f["perturbed_equilibrium/response/psi_area"])), read(f["integration/psi"]))
+end
+
+size(pa_g, 1) == length(psi_g) || error("gal grid mismatch: npsi=$(size(pa_g,1)) vs grid=$(length(psi_g))")
+size(pa_s, 1) == length(psi_s) || error("shoot grid mismatch: npsi=$(size(pa_s,1)) vs grid=$(length(psi_s))")
+
+col = mtarget - mlow + 1
+@printf("m=%d → column %d (mlow=%d, mpert=%d)\n", mtarget, col, mlow, size(pa_g, 2))
+g = @view pa_g[:, col]
+s = @view pa_s[:, col]
+
+# resonant surface for this m (q = m/n)
+psi_res = mtarget in sing_m ? sing_psi[findfirst(==(mtarget), sing_m)] : NaN
+# value at the grid point nearest the resonant surface, each run
+if !isnan(psi_res)
+    ig = argmin(abs.(psi_g .- psi_res))
+    is = argmin(abs.(psi_s .- psi_res))
+    @printf("m=%d resonant surface ψ=%.4f\n", mtarget, psi_res)
+    @printf("  nearest |b^ψ_area|:  gal-ideal=%.4e (ψ=%.4f)   shooting=%.4e (ψ=%.4f)   ratio=%.3f\n",
+        abs(g[ig]), psi_g[ig], abs(s[is]), psi_s[is], abs(g[ig]) / abs(s[is]))
+end
+
+plt = plot(; size=(1000, 620), xlabel="ψ_N", ylabel="|b^ψ / ⟨J·|∇ψ|⟩|  (area-normalized)",
+    title="m=$mtarget perturbed normal field — gal-ideal vs shooting PE", legend=:topleft,
+    left_margin=13Plots.mm, bottom_margin=5Plots.mm, right_margin=4Plots.mm)
+plot!(plt, psi_g, abs.(g); lw=2.2, color=1, label="gal-ideal ξ → PE")
+plot!(plt, psi_s, abs.(s); lw=1.8, color=2, ls=:dash, label="shooting ξ → PE")
+for (j, ps) in enumerate(sing_psi)
+    ism2 = sing_m[j] == mtarget
+    vline!(plt, [ps]; color=(ism2 ? :red : :black), ls=:dot, lw=(ism2 ? 1.8 : 1),
+        label=(ism2 ? "q=$mtarget surface" : (j == 1 ? "rational q" : "")))
+end
+savefig(plt, outpng)
+println("saved: ", abspath(outpng))

benchmarks/gal_validation/gal_validation/gal_beta_scan.jl

@@ -0,0 +1,111 @@
+#!/usr/bin/env julia
+# Gal-enabled TJ-analytic beta (pc) scan: at each pc, run the full GPEC pipeline with gal_flag=true
+# and capture three Δ′ measures per resonant surface (m=2 q=2, m=3 q=3):
+#   gal      = galerkin/pest3_Delta diagonal      (RDCON singular-Galerkin, this port)
+#   stride   = singular/delta_prime_matrix diag   (STRIDE BVP, pr178)
+#   cajump   = (ca_r-ca_l)/((2pi)^2 psio)          (ca-jump, pr178 primary)
+# Reuses examples/LAR_beta_scan/gpec.toml. Writes /tmp/gal_beta_scan_results.h5 + CSV.
+using Pkg
+Pkg.activate(joinpath(@__DIR__, "..", "Projects", "GPEC_dev", "docs", "GPEC_julia"))
+const REPO = "/Users/pharr/Projects/GPEC_dev/docs/GPEC_julia"
+Pkg.activate(REPO)
+using GeneralizedPerturbedEquilibrium
+using HDF5, TOML, Printf
+
+_warped_grid(x0, x1, N; p=2.0) = [x0 + (x1 - x0) * (1 - (1 - i / (N - 1))^p) for i in 0:N-1]
+const PC_FULL = _warped_grid(0.001, 0.1735, 40; p=2.0)
+const PC_TEST = [0.001, 0.05, 0.10, 0.15, 0.17]
+const BASE = TOML.parsefile(joinpath(REPO, "examples", "LAR_beta_scan", "gpec.toml"))
+
+function extract(h5)
+    h5open(h5, "r") do f
+        out = Dict{String,Any}()
+        out["q0"] = read(f, "equil/q0"); out["qmax"] = read(f, "equil/qmax")
+        out["dW_total"] = real(read(f, "vacuum/et")[1])
+        psio = read(f, "equil/psio"); denom = (2pi)^2 * psio
+        mlow = read(f, "info/mlow"); mpert = read(f, "info/mpert"); nlow = read(f, "info/nlow")
+        # surface m-values (STRIDE side)
+        m_s = vec(read(f, "singular/m")); n_s = vec(read(f, "singular/n")); msing = length(m_s)
+        # STRIDE delta_prime_matrix
+        dpm = haskey(f["singular"], "delta_prime_matrix") ? read(f, "singular/delta_prime_matrix") : nothing
+        cal = read(f, "singular/ca_left"); car = read(f, "singular/ca_right")
+        # gal
+        gm = haskey(f, "galerkin") && haskey(f["galerkin"], "sing_m") ? vec(read(f, "galerkin/sing_m")) : Int[]
+        gp = haskey(f, "galerkin") && haskey(f["galerkin"], "pest3_Delta") ? read(f, "galerkin/pest3_Delta") : nothing
+        for (s, m) in enumerate(m_s)
+            key = "m$(m)"
+            # ca-jump
+            ipr = 1 + m - mlow + (n_s[s] - nlow) * mpert
+            out["cajump_$key"] = (car[ipr, ipr, 2, s] - cal[ipr, ipr, 2, s]) / denom
+            # stride
+            out["stride_$key"] = dpm !== nothing && s <= size(dpm, 1) ? dpm[s, s] : NaN + NaN*im
+            # gal (match by m)
+            gi = findfirst(==(m), gm)
+            out["gal_$key"] = (gp !== nothing && gi !== nothing) ? gp[gi, gi] : NaN + NaN*im
+        end
+        return out
+    end
+end
+
+function run_point(pc)
+    dir = mktempdir(; prefix="gal_beta_")
+    try
+        c = deepcopy(BASE)
+        c["TJ_ANALYTIC_INPUT"]["pc"] = pc
+        c["ForceFreeStates"]["gal_flag"] = true
+        c["ForceFreeStates"]["gal_solver"] = "LU"
+        c["ForceFreeStates"]["HDF5_filename"] = joinpath(dir, "gpec.h5")
+        open(joinpath(dir, "gpec.toml"), "w") do io; TOML.print(io, c); end
+        GeneralizedPerturbedEquilibrium.main([dir])
+        return extract(joinpath(dir, "gpec.h5"))
+    catch e
+        @warn "pc=$pc failed" exception=(e, catch_backtrace())
+        return nothing
+    finally
+        rm(dir; force=true, recursive=true)
+    end
+end
+
+pcs = ("--test" in ARGS) ? PC_TEST : PC_FULL
+@info "Gal beta scan: $(length(pcs)) points (ε=0.2, qc=1.5, qa=3.6)"
+results = Tuple{Float64,Dict{String,Any}}[]
+t0 = time()
+for (i, pc) in enumerate(pcs)
+    ti = time()
+    r = run_point(pc)
+    el = time() - ti; tot = time() - t0
+    if r !== nothing
+        push!(results, (pc, r))
+        g2 = get(r, "gal_m2", NaN); s2 = get(r, "stride_m2", NaN)
+        g3 = get(r, "gal_m3", NaN); s3 = get(r, "stride_m3", NaN)
+        @printf("[%2d/%2d] pc=%.5f  q=2: gal=%+9.3f stride=%+9.3f | q=3: gal=%+9.3f stride=%+9.3f  dW=%+.4f  (%.0fs, tot %.0fs)\n",
+            i, length(pcs), pc, real(g2), real(s2), real(g3), real(s3), get(r, "dW_total", NaN), el, tot)
+        flush(stdout)
+    end
+end
+
+# Save
+out_h5 = "/tmp/gal_beta_scan_results.h5"
+isfile(out_h5) && rm(out_h5)
+h5open(out_h5, "w") do f
+    f["pc"] = [r[1] for r in results]
+    for key in ("gal_m2", "stride_m2", "cajump_m2", "gal_m3", "stride_m3", "cajump_m3", "dW_total")
+        vals = [get(r[2], key, NaN+NaN*im) for r in results]
+        if eltype(vals) <: Complex
+            f["$(key)_re"] = real.(vals); f["$(key)_im"] = imag.(vals)
+        else
+            f[key] = Float64.(vals)
+        end
+    end
+end
+open("/tmp/gal_beta_scan_results.csv", "w") do io
+    println(io, "pc,gal_m2,stride_m2,cajump_m2,gal_m3,stride_m3,cajump_m3,dW_total")
+    for (pc, r) in results
+        @printf(io, "%.6f,%.6f,%.6f,%.6f,%.6f,%.6f,%.6f,%.6f\n", pc,
+            real(get(r,"gal_m2",NaN)), real(get(r,"stride_m2",NaN)), real(get(r,"cajump_m2",NaN)),
+            real(get(r,"gal_m3",NaN)), real(get(r,"stride_m3",NaN)), real(get(r,"cajump_m3",NaN)),
+            get(r,"dW_total",NaN))
+    end
+end
+@info "Saved /tmp/gal_beta_scan_results.h5 and .csv ($(length(results)) points)"
+println("DONE")