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 case ggj_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]

@matt-pharr what's the plan for this one? can we merge it soon so we have fewer active PRs to juggle?

[matt-pharr]

@logan-nc I have kept it up to date with develop. It can be merged soon but needs some cleanup. The GGJ portion is validated but we want to see what Via's last couple of weeks of work gets us