InnerLayer (GGJ) - NEW - rotated-ray collocation backend (:ray), robust to |Q| ~ 500

[matt-pharr]

Summary

Adds a third interchangeable GGJ inner-layer solver — GGJModel{:ray}, a rotated-contour spectral-element collocation method — and makes it the default. It returns the same parity-projected matching data (Δ_odd, Δ_even) in the same deltac convention as :galerkin and :shooting, but stays accurate far past where both existing backends fail: validated to Q = 500i, the regime resistivity scans (|Q| ∝ η^(−1/3)) and rotation physics actually visit. The branch then integrates the layer into the RPEC outer↔inner matching pipeline: a backend-agnostic solve_inner_profile interface, the composite (cut + layer) inner-region solution, a corrected layer-based penetrated field, and new HDF5 outputs.

The :ray backend

Method. The layer equations are continued analytically onto the ray x = e^{iθ}s with θ = arg(Q)/4, which makes the parabolic-cylinder WKB exponent exactly real and lifts the coefficient pseudo-resonance x² ≈ −Q²(G+KF) — which sits on the real axis for imaginary Q — off the contour. On [0, s_m] a global Chebyshev spectral-element collocation BVP is solved with parity conditions at the origin and (Δ, c₁, c₂) as bordered unknowns: one sparse LU per mesh plus a rank-3 Woodbury correction for the second parity. The far-field condition reuses the existing validated inps Wasow kernel (evaluated at complex x by the new RayAsymptotics.jl; construction untouched), applied at the series radius S and transported inward by an L-stable 2-stage Radau IIA march in the quotient modulo the decaying exponential pair. The damped-zone march runs in Complex{Double64}: at large S, Float64 linear algebra floors Δ at ~1e-4 mixing — a limit generic to any inps-based method, removed by extended precision. Since Δ is an analytic invariant of the contour angle, a θ-perturbed re-solve doubles as a built-in error bar (delta_convergence).

Validation.

• Manufactured bordered problem with prescribed Δ*: 3e-14.
• Fortran rmatch pins and :galerkin cross-checks at moderate Q: 1e-5–1e-9.
• Campaign scale (Slurm, 102 equilibria × 65-point Q-grids incl. DIII-D and ITER PFPO, ~14k solves): zero θ-consistency flags above 1e-3; worst Δ₁ drift 5.8e-4, worst Δ₂ 3.7e-6. The legacy Galerkin self-disagrees at O(1) for |Q| ≳ 1 on the same points.
• In-repo: test/runtests_innerlayer.jl gains :ray testsets (GW2020 Eq. 55 paper point vs :galerkin; q=4 physical pins at Q = 500i; θ-invariance; internal machinery units; solve_inner_profile contract for both backends).

Default change.GGJModel() now constructs GGJModel{:ray}(). There are no bare GGJModel() call sites in the repo — all benchmark and harness uses pass solver= explicitly — and the backends take disjoint numerical-knob keywords, so stale galerkin-knob calls fail loudly rather than returning a wrong Δ.

RPEC matching integration

• solve_inner_profile interface (InnerLayerInterface.jl): matching data + reconstructed layer profiles + unit conversions, so matching drivers never touch model internals. The :ray implementation takes certified Δ from the optimal-θ contour and profiles from a θ = 0 re-solve (the on-axis pseudo-resonance is a regular point for the BVP); their Δ agreement is logged as a runtime certificate (certΔ).
• gal_inner_solver control ("ray" default | "galerkin"), with the galerkin numerical knobs exposed as gal_inner_* TOML keys.
• Composite inner-region solution: the Galerkin cut solution (sing_get_ua_res_cut, xi_cut/cut_range) supplies the smooth outer background the layer solution grafts onto — matching Fortran match.f intotsol, and required for the inner b^ψ profile to overlap the outer eigenfunction.
• Penetrated field from the layer (bug fix + new): the physical penetrated (reconnected) resonant field is the layer-center Ψ(0). The ported match.f-style scale omitted a −2π/v₁ factor (3.5–7× surface-dependent error on the DIII-D benchmark); the corrected scale is derived from the code's own outer b-convention chained with GWP2016 Eqs. 16/A17 and validated by an inner↔outer overlap test on all four DIII-D surfaces. B_pen now flows from ForceFreeStates into PerturbedEquilibrium via inner_bpen and contracts through the standard C_coeffs machinery. Ideal runs produce exactly zero (perfect shielding). Kinetic runs currently also report zero — the pointwise fallback was retired and a kinetic-appropriate replacement is deferred to a follow-up.
• Δ_mn derivative accuracy: Δ_mn is evaluated from the finite-difference jump across the rational surface as before; gal-matched runs now use the analytic Hermite-cubic derivative (ud_store carries the exact gal-basis ξ′), while shooting runs keep the chord-slope evaluation.
• delta_coil: edge coil-response block from the Riccati S-axis BVP (unit source per edge poloidal mode; one factorization, all modes as RHS columns).

Outputs, docs, infrastructure

• New HDF5: galerkin/match/inner_params/{E,F,G,H,K,M,taua,taur,v1}, galerkin/match/inner/b_$i, galerkin/solution/{xi_cut, cut_range}, singular/delta_coil, singular_coupling/{forcing_solution_weights, rational_area}, response/psi_n.
• Docs: JCP-style documentation standard (docs/DOC_STANDARD.md) with provenance-stamped, manifest-tracked figures; full InnerLayer methods page with four generated figures; citation fixes.
• Regression: new ggj_ray_q500i computed case (Δ at Q = 500i, 1e-4 absolute noise floor reflecting BLAS-dependent damped-zone roundoff at ~1e-5 relative).
• Examples: LAR_resistive_match_test / LAR_ideal_match_test — a second geometry for the matched solve, with its ideal-limit companion.
• Benchmarks: hardcoded local Fortran paths replaced with required CLI args / env vars (no on-disk defaults). Library code reads no environment variables.

Dependencies

• Random (stdlib) promoted from test-extras to [deps] — fixed-seed decaying-pair extraction (provably Δ-neutral; seeded for bit-reproducibility).
• DoubleFloats was already a dependency; only the using is new.

Physics motivation

On a DIII-D k_η scan the legacy galerkin layer produces a non-monotonic-in-η jump of the matched resonant field at the lowest resistivity (|Q| ≳ 1); the ray backend restores a smooth, monotonic shielding family — the inner-layer solver choice materially changes a physical prediction exactly where |Q| leaves the legacy regime. A replication-grade methods write-up exists and will accompany the paper in preparation.

[matt-pharr]

@logan-nc@d-burg this is something I have been working on with claude for a few weeks. It has a number of advantages over the galerkin GGJ solver. It uses collocation which is also a finite element method, but I learned that the reason that the galerkin solver sometimes freaks out is because there are poles to the GGJ solutions in Q space, and when you get near those poles the solutions become super oscillatory on the real axis. However, since we are solving this in a complex domain, you can show that the GGJ solutions should be entire in the complex plane, so the mercier terms in the large |X| limit are the same polynomial + coefficient even off of the real-axis (which is where the physically-relevant part of the solution lies). So you can get a much cleaner/quicker exponential decay of the non-limiting terms if you integrate on a ray that is rotated off of the real axis. This also gives way to an entirely different numerical stability check, which is to re-integrate for multiple rotation angles and see that you get the same Delta_even and Delta_odd. Overall this method is far superior to the galerkin one in all of my testing so far, and can safely go up to Q ~ 500i with a maximum error in the final Delta of ~1e-4 which is pretty good for this problem. Error stays ~1e-7 for more reasonable Q. This enables scanning to lower resistivities (Q ~ 1/eta^1/3) and higher rotation rates (or to absurd growth rates). It also makes the solve more robust for lower |Q|. I will do some more testing on this branch before I merge but I am opening it for review now.

[logan-nc]

Documentation for the :ray backend + a repeatable docs standard

This adds full user-facing documentation for the new rotated-ray inner-layer backend, and uses it as the occasion to establish a reusable Journal of Computational Physics-style documentation standard for the repo.

📖 Rendered preview (live now for this PR)

• Inner Layer page (the main deliverable, with all four figures inline): https://openfusiontoolkit.github.io/GPEC/previews/PR323/inner_layer/
• Ballooning page (migrated to the new figure layout): https://openfusiontoolkit.github.io/GPEC/previews/PR323/ballooning/
• Preview home: https://openfusiontoolkit.github.io/GPEC/previews/PR323/

What's here

Documentation standard — docs/DOC_STANDARD.md defines a module-page template (governing equations → numerical method → validation figures → API) and the figure policy. Linked from Developer Notes; ballooning.md and the new inner_layer.md are the reference exemplars.

Figure infrastructure — docs/figure_tools.jl provides save_doc_figure, which stamps every PNG with its generating git commit + date and records provenance (script, commit, date, depends) in docs/src/figures/manifest.toml.

• Content figures now live in docs/src/figures/<module>/, each PNG committed next to the make_*.jl that produced it. docs/src/assets/ is reserved for Documenter chrome.
• Figures are pre-generated and committed, never rebuilt at docs-build time — regenerate only when a file in a figure's depends list changes the numbers it shows (see DOC_STANDARD.md).

Inner Layer page (docs/src/inner_layer.md, stub → full page) with four computed figures on the q=4 benchmark surface: backend accuracy vs :ray along the imaginary axis (shooting reaches |Q|≈1, galerkin ≈4, ray to 500), the rotated-ray/pseudo-resonance geometry, the seamless collocation↔asymptotic field match, and the Q=500i convergence error bar.

Ballooning migration — the two existing figures moved to docs/src/figures/ballooning/ as the second exemplar, registered legacy (their generators predate the provenance system; one compared the since-removed Mercier.jl).

Review-related fixes folded in

• rescale_delta docstring correction (GGJParameters.jl): the physical-rescale exponent was mislabeled X₀^{+2√(−D_I)}; the code correctly applies X₀^{−2√(−D_I)}·v₁^{2√(−D_I)} (verified against the Fortran rmatch pins). Docstring now matches the code. No behavior change.
• New regression caseggj_ray_q500i pins the :ray matching data at Q=500i on the q=4 surface (Δ_odd = 2.472008 + 13.35402i, Δ_even = 0.1374966 + 0.7427542i).
• Unit tests for the internal machinery (cheblobatto, ode_matrix, parity_rows, decaying_pair, profile diagnostics, delta_convergence).

Verification

• Docs build clean (no missing-docs / invalid-link / cross-ref errors); all six figures embed — confirmed live in the preview above.
• test/runtests_innerlayer.jl: 37/37 pass.
• Regression harness: both GGJ cases compute; the :ray case matches the test pins.

_{Note: committed figure stamps read 3a0837e3-dirty (generated before this commit existed). Content is final; the stamps can optionally be refreshed to the clean hash by re-running the make_*.jl scripts.}

🤖 Generated with Claude Code

[logan-nc]

@matt-pharr I was having trouble understanding/trusting all your claims from the pure text descriptions above. As such, I am taking this opportunity to initiate a new, more-complete docs standard that includes figures. Figures were already creeping in (see legacy ballooning figures) and a standard was needed anyways - I figured I might as well lump it in here so i could actually see what was happening 😝

[logan-nc]

@matt-pharr can you check the "Review-related fixes folded in".
Also read through https://openfusiontoolkit.github.io/GPEC/previews/PR323/inner_layer for (a) accuracy and (b) any improvements you want to make.

[logan-nc]

@claude provide a concise review of this, identifying major issues only

[claude]

Claude finished @logan-nc's task in 3m 39s —— View job

---

Review: rotated-ray GGJ inner-layer backend + RPEC matching integration

• Gather context (diff vs origin/develop, 69 files)
• Review new :ray backend (Ray.jl, RayAsymptotics.jl) — physics/numerics pass
• Review RPEC matching integration (solve_inner_profile, penetrated-field scale fix, composite inner solution)
• Review defaults/interface changes and regression coverage
• Post review

Scope note: this is a large PR (4139/187, 44 commits). I focused on the physics-critical new code (Ray.jl, RayAsymptotics.jl, the matching/penetrated-field wiring) rather than docs/benchmark infrastructure.

Major issues

1. Production solve_inner(::GGJModel{:ray}, ...) drops the solve-quality diagnostics it computes.Ray.jl:936-939 builds RaySolveResult with series_ok, bc_cond, and resid fields, but the solve_inner wrapper — the entry point actually called by gal_match_rpec/resonant_match_rpec in production — discards everything except the two Δ values (InnerLayerResponse(res.Δ[2], res.Δ[1])). A poorly-converged or ill-conditioned surface in a large batch run (the PR cites a 14k-solve campaign) has no signal beyond a @warn line, easy to miss at scale and invisible to any downstream automated check. Worth threading at least bc_cond/series_ok through to the matching driver (or asserting on them) before this becomes the default solver.

2. The flagship |Q|~500i validation is self-referential in-repo. The q=4 benchmark test at Q=500i (test/runtests_innerlayer.jl) pins against values from "the pre-port validation suite" (i.e. the branch's own prior computed output), not an independent Fortran/analytic ground truth — unlike the Q=0.1234 real-axis case, which does reuse a genuine Fortran-rmatch-cross-checked pin. The only in-repo cross-check at the high-|Q| imaginary regime is θ-invariance (delta_convergence). That check is blind to any θ-independent systematic error (asymptotic-basis column order, rescale_delta, the parity swap at Ray.jl:915) — such a bug would reproduce identically at every θ and pass. The 102-equilibrium/14k-solve campaign referenced in the PR body isn't reproducible from the repo. Consider adding at least one independent ground-truth point at large |Q| (even a slow high-precision shooting cross-check) before this ships as the default.

3. Extended precision claim is incomplete in the damped-zone march. The PR states the march runs in Complex{Double64} specifically to remove the ~1e-4 Δ floor from Float64 mixing at large S. But the ray parameter s driving the Radau IIA stages (Ray.jl:391-392, 442, stage abscissae and the s += h accumulation) is plain Float64, accumulated over potentially many steps — this reintroduces the same error channel the extended-precision LU was meant to close. Worth carrying s as Double64 too, or explaining why it doesn't matter.

4. Likely DimensionMismatch for multi-toroidal-mode (npert>1) runs when the inner layer isn't gal-matched.GeneralizedPerturbedEquilibrium.jl:595 sizes the inner_bpen fallback as (msing, numpert_total), but SingularCoupling.jl:422 consumes it as transpose(C_coeffs) * inner_bpen[s,:] where C_coeffs is mpert × numpert_total — i.e. it needs columns of length mpert, matching gal_match_rpec's convention (mcoil = mpert, GalerkinMatch.jl:31). This only happens to work when npert == 1 (single-n, the common case); for nlow ≠ nhigh runs without RPEC gal-matching this will throw at runtime rather than silently zero-filling as intended.

Minor

• benchmarks/deltacoil_sensitivity/scripts/plot_deltacoil_heatmap_vary.jl:11 still has a hardcoded personal path (/Users/viaweber/Desktop/...) as the ARGS-empty fallback, despite the PR description's claim that benchmark scripts no longer have on-disk defaults.
• src/ForceFreeStates/Resist.jl is missing a trailing newline (CLAUDE.md file-hygiene rule; pre-commit should catch this but flagging since it's visible in the diff).
• Ray.jlsolve_inner_profile(::GGJModel{:ray}, ...) calls solve_ray(params, Q; θ=0.0, kwargs...) (Ray.jl:961) — if a caller's kwargs ever included θ, Julia's keyword-splat ordering means it silently overrides the explicit θ=0.0, defeating the "certified vs θ=0" cross-check the certΔ certificate relies on. No current call site passes θ this way, so this is latent, not active.

Not major, but worth flagging for the record

GGJModel() now defaults to :ray instead of :galerkin — confirmed there are no bare GGJModel() call sites elsewhere in the repo, so this is safe internally, but it's a real behavior change for any external/downstream user of the package.

Physics fundamentals of the rotated-ray method itself (θ=arg(Q)/4 driving the parabolic-cylinder pair off the real-axis pseudo-resonance, entire-function analytic continuation, Chebyshev collocation with bordered Δ/c₁/c₂ unknowns + rank-3 Woodbury, Radau IIA far-field matching) check out on manual inspection — no sign/index errors or placeholder terms found in Ray.jl/RayAsymptotics.jl.

---

⚠️ Per repo policy (CLAUDE.md): this PR must not be merged without a human reviewer's approval — that requirement is non-negotiable regardless of any AI review outcome. Also per CLAUDE.md, please run the regression harness (regress --cases <relevant-cases> --refs develop,local) and confirm it's clean before merge, even though a ggj_ray_q500i regression case has already been added in this PR.
· branch GGJ_colocation_backend