diff --git a/.claude/agent-memory/fortran-physics-reviewer/fortran_correspondence_map.md b/.claude/agent-memory/fortran-physics-reviewer/fortran_correspondence_map.md index 95852271..4ba8494e 100644 --- a/.claude/agent-memory/fortran-physics-reviewer/fortran_correspondence_map.md +++ b/.claude/agent-memory/fortran-physics-reviewer/fortran_correspondence_map.md @@ -1,33 +1,36 @@ --- name: Fortran-Julia File Correspondence Map -description: Mapping between Fortran GPEC source files at ~/Code/gpec and Julia JPEC modules for physics review +description: Mapping between Fortran GPEC source files and Julia JPEC modules for physics review type: reference --- +Fortran paths are relative to the root of the Fortran GPEC clone; Julia paths are +relative to this repository. + ## Coil/ForcingTerms -- `~/Code/gpec/coil/coil.F` (coil_read) -> `src/ForcingTerms/CoilGeometry.jl` (apply_transforms) -- `~/Code/gpec/coil/field.F` (field_bs_psi) -> `src/ForcingTerms/BiotSavart.jl` + `src/ForcingTerms/CoilFourier.jl` +- `coil/coil.F` (coil_read) -> `src/ForcingTerms/CoilGeometry.jl` (apply_transforms) +- `coil/field.F` (field_bs_psi) -> `src/ForcingTerms/BiotSavart.jl` + `src/ForcingTerms/CoilFourier.jl` ## PerturbedEquilibrium -- `~/Code/gpec/gpec/gpeq.f` (gpeq_sol, gpeq_contra, gpeq_surface, gpeq_normal) -> `src/PerturbedEquilibrium/FieldReconstruction.jl` + `src/PerturbedEquilibrium/ResponseMatrices.jl` -- `~/Code/gpec/gpec/gpresp.f` (gpresp_pinduct, gpresp_sinduct, gpresp_permeab) -> `src/PerturbedEquilibrium/ResponseMatrices.jl` -- `~/Code/gpec/gpec/gpout.f` (gpout_singcoup, gpout_xbnormal) -> `src/PerturbedEquilibrium/SingularCoupling.jl` + `src/PerturbedEquilibrium/FieldReconstruction.jl` -- `~/Code/gpec/gpec/gpvacuum.f` (gpvacuum_flxsurf) -> `src/PerturbedEquilibrium/SingularCoupling.jl` (compute_surface_inductance_from_greens) +- `gpec/gpeq.f` (gpeq_sol, gpeq_contra, gpeq_surface, gpeq_normal) -> `src/PerturbedEquilibrium/FieldReconstruction.jl` + `src/PerturbedEquilibrium/ResponseMatrices.jl` +- `gpec/gpresp.f` (gpresp_pinduct, gpresp_sinduct, gpresp_permeab) -> `src/PerturbedEquilibrium/ResponseMatrices.jl` +- `gpec/gpout.f` (gpout_singcoup, gpout_xbnormal) -> `src/PerturbedEquilibrium/SingularCoupling.jl` + `src/PerturbedEquilibrium/FieldReconstruction.jl` +- `gpec/gpvacuum.f` (gpvacuum_flxsurf) -> `src/PerturbedEquilibrium/SingularCoupling.jl` (compute_surface_inductance_from_greens) ## KineticForces (NTV, Fortran `pentrc/`) — see kinetic_ntv_map.md for the audit checklist -- `~/Code/gpec/pentrc/torque.F90` -> `src/KineticForces/Torque.jl` (tpsi! single-surface torque) -- `~/Code/gpec/pentrc/pitch.f90` -> `src/KineticForces/PitchIntegration.jl` (pitch-angle lambda integration, bounce-averaged operator assembly) -- `~/Code/gpec/pentrc/energy.f90` -> `src/KineticForces/EnergyIntegration.jl` (energy-space quadrature integrand) -- `~/Code/gpec/pentrc/pentrc.F90` (orchestration) -> `src/KineticForces/{KineticForces.jl,Compute.jl}` -- `~/Code/gpec/pentrc/inputs.f90` + `params.f90` -> `src/KineticForces/KineticForcesStructs.jl` -- `~/Code/gpec/pentrc/dcon_interface.f` -> `src/KineticForces/CalculatedKineticMatrices.jl` (bridge to ForceFreeStates kinetic stability) + `src/ForceFreeStates/{Kinetic.jl,FixedKineticMatrices.jl}` +- `pentrc/torque.F90` -> `src/KineticForces/Torque.jl` (tpsi! single-surface torque) +- `pentrc/pitch.f90` -> `src/KineticForces/PitchIntegration.jl` (pitch-angle lambda integration, bounce-averaged operator assembly) +- `pentrc/energy.f90` -> `src/KineticForces/EnergyIntegration.jl` (energy-space quadrature integrand) +- `pentrc/pentrc.F90` (orchestration) -> `src/KineticForces/{KineticForces.jl,Compute.jl}` +- `pentrc/inputs.f90` + `params.f90` -> `src/KineticForces/KineticForcesStructs.jl` +- `pentrc/dcon_interface.f` -> `src/KineticForces/CalculatedKineticMatrices.jl` (bridge to ForceFreeStates kinetic stability) + `src/ForceFreeStates/{Kinetic.jl,FixedKineticMatrices.jl}` - bounce averaging of the 6 perturbed-action matrices -> `src/KineticForces/BounceAveraging.jl` ## InnerLayer (resistive matched-asymptotics, Fortran `rmatch/`) — see resistive_layer_map.md for the audit checklist -- `~/Code/gpec/rmatch/deltac.f` (Galerkin solver) -> `src/InnerLayer/GGJ/Galerkin.jl` -- `~/Code/gpec/rmatch/deltar.f` (shooting solver) -> `src/InnerLayer/GGJ/Shooting.jl` -- `~/Code/gpec/rmatch/inps.f` + `inpso.f` (Wasow asymptotic basis: T,J,P,B,Q,C,D,Y,Z,U matrices) -> `src/InnerLayer/GGJ/InnerAsymptotics.jl` -- `~/Code/gpec/rmatch/{inner.f,match.f,msing.f,gamma.f}` -> `src/InnerLayer/GGJ/` (parameters, matching data, special functions) +- `rmatch/deltac.f` (Galerkin solver) -> `src/InnerLayer/GGJ/Galerkin.jl` +- `rmatch/deltar.f` (shooting solver) -> `src/InnerLayer/GGJ/Shooting.jl` +- `rmatch/inps.f` + `inpso.f` (Wasow asymptotic basis: T,J,P,B,Q,C,D,Y,Z,U matrices) -> `src/InnerLayer/GGJ/InnerAsymptotics.jl` +- `rmatch/{inner.f,match.f,msing.f,gamma.f}` -> `src/InnerLayer/GGJ/` (parameters, matching data, special functions) - SLAYER drift-MHD two-fluid solver: `src/InnerLayer/SLAYER/` is a placeholder pending implementation ## Key Fortran Conventions diff --git a/.claude/agent-memory/fortran-physics-reviewer/galerkin_assembly_map.md b/.claude/agent-memory/fortran-physics-reviewer/galerkin_assembly_map.md index 3dd7bec8..a24db8f4 100644 --- a/.claude/agent-memory/fortran-physics-reviewer/galerkin_assembly_map.md +++ b/.claude/agent-memory/fortran-physics-reviewer/galerkin_assembly_map.md @@ -49,7 +49,7 @@ lives only in match_delta (Newton eigenvalue search), NOT the RPEC path. Two iss Confirm intended convention vs Fortran before trusting resistive Δ(Q); ideal_flag path (γ unused) is unaffected. RESOLUTION 2026-06-20: author confirms the Hz convention is DELIBERATE (gal_rotation is f[Hz], so γ=2πi·n·f); the overstated comments were corrected. Δ(Q) now guarded by test/runtests_innerlayer.jl (Julia self-pin on -glasser_wang_2020_eq55). Remaining follow-up (mpharr): upgrade that guard to a Fortran deltac_run / paper +glasser_wang_2020_eq55). Remaining follow-up: upgrade that guard to a Fortran deltac_run / paper cross-check to fully close the convention question. ## Verdict: PASS WITH REQUIRED ANNOTATIONS for outer solve; RPEC γ convention is SHOULD-FIX (verify vs match.f). diff --git a/.claude/agent-memory/fortran-physics-reviewer/kinetic_ntv_map.md b/.claude/agent-memory/fortran-physics-reviewer/kinetic_ntv_map.md index 414be7a9..5f5dfb56 100644 --- a/.claude/agent-memory/fortran-physics-reviewer/kinetic_ntv_map.md +++ b/.claude/agent-memory/fortran-physics-reviewer/kinetic_ntv_map.md @@ -6,7 +6,7 @@ type: reference The KineticForces module computes neoclassical toroidal viscosity (NTV) torque and the kinetic energy/matrices from trapped-particle nonambipolar transport. Fortran counterpart is -`~/Code/gpec/pentrc/` (file mapping in fortran_correspondence_map.md). +the Fortran `pentrc/` sources (file mapping in fortran_correspondence_map.md). ## Governing theory - Logan & Park (2013) PoP 20, 122507: diamagnetic frequency (Eq. 7), energy integrand (Eq. 8), torque normalization Im(T) = 2n·δW_k (Eq. 19). diff --git a/.claude/agent-memory/fortran-physics-reviewer/resistive_layer_map.md b/.claude/agent-memory/fortran-physics-reviewer/resistive_layer_map.md index 1a794e7f..57e3b1d7 100644 --- a/.claude/agent-memory/fortran-physics-reviewer/resistive_layer_map.md +++ b/.claude/agent-memory/fortran-physics-reviewer/resistive_layer_map.md @@ -7,7 +7,7 @@ type: reference The InnerLayer module performs matched-asymptotic analysis of resistive MHD stability at rational surfaces — solving the resistive inner-layer equations and computing the tearing stability parameter Δ via parity-projected splitting (Δ_odd, Δ_even). Fortran counterpart is -`~/Code/gpec/rmatch/` (file mapping in fortran_correspondence_map.md). +the Fortran `rmatch/` sources (file mapping in fortran_correspondence_map.md). ## Governing theory - Glasser (2016) PoP 23, 072505: computation of resistive instabilities by matched asymptotic expansions; shooting method, Frobenius exponent splitting. diff --git a/.claude/agents/fortran-physics-reviewer.md b/.claude/agents/fortran-physics-reviewer.md index 986f8eb3..8def0bc2 100644 --- a/.claude/agents/fortran-physics-reviewer.md +++ b/.claude/agents/fortran-physics-reviewer.md @@ -35,7 +35,7 @@ You operate under a hard budget to protect the user's token quota: - **KineticForces** (NTV, formerly PENTRC): Logan & Park (2013) PoP 20, 122507; Logan (2015) PhD Thesis (Ch. 7 Eqs 7.30–7.35, App. C/D); Park et al. (2009) PRL 102, 065002 - **InnerLayer** (resistive matched-asymptotics, GGJ/SLAYER): Glasser (2016) PoP 23, 072505; Glasser (2018) PoP 25, 032501; Glasser (2020) "Asymptotic solutions and convergence studies of the resistive inner region equations"; Wang et al. (2020) PoP 27, 122509 - **Analysis**: visualization/post-processing only — no physics kernels; verify it reads the correct HDF5 groups and applies correct conventions rather than auditing equations. -- Locate the corresponding Fortran GPEC source at `~/Code/gpec` (ask the user for the path if not found there). KineticForces ↔ `~/Code/gpec/pentrc/`, InnerLayer ↔ `~/Code/gpec/rmatch/`. Read the relevant Fortran routines carefully. Your project memory holds a maintained correspondence map and per-domain audit checklists — consult it first. +- Locate the corresponding Fortran GPEC source in the user's local Fortran GPEC clone (ask the user for its path). KineticForces ↔ `pentrc/`, InnerLayer ↔ `rmatch/`. Read the relevant Fortran routines carefully. Your project memory holds a maintained correspondence map and per-domain audit checklists — consult it first. ### Step 2: Read the Fortran implementation - Open and carefully read the relevant Fortran source files. @@ -132,7 +132,7 @@ Examples of what to record: # Persistent Agent Memory -You have a persistent, file-based memory system at `~/.claude/agent-memory/fortran-physics-reviewer/` (shared across all clones of this repo on this machine). When reading or writing files, expand `~` to the absolute home directory path (e.g. use `/Users/yourname/.claude/agent-memory/fortran-physics-reviewer/`). Create the directory if it does not exist before writing. +You have a persistent, file-based memory system in this repository at `.claude/agent-memory/fortran-physics-reviewer/` (versioned with the repo, so findings are shared with everyone working on it). Resolve it relative to the repository root, and create the directory if it does not exist before writing. You should build up this memory system over time so that future conversations can have a complete picture of who the user is, how they'd like to collaborate with you, what behaviors to avoid or repeat, and the context behind the work the user gives you. diff --git a/.gitignore b/.gitignore index 7e1ff4f6..0729e730 100644 --- a/.gitignore +++ b/.gitignore @@ -16,9 +16,13 @@ !docs/resources/*.pdf *.png !docs/src/assets/**/*.png +# committed documentation figures (generated by docs/src/figures/**/make_*.jl) +!docs/src/figures/**/*.png *.jld2 .gitattributes Manifest.toml +# the docs figure-provenance manifest is not a Julia Manifest.toml (case-insensitive FS collides) +!docs/src/figures/manifest.toml *.log docs/build/ anaconda_projects/ diff --git a/CLAUDE.md b/CLAUDE.md index 476b1ba0..ac675430 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -8,7 +8,7 @@ GPEC (Generalized Perturbed Equilibrium Code, Julia implementation) is a compreh **Relationship to Fortran GPEC**: This Julia GPEC is an evolution of the Fortran GPEC code suite, available at https://github.com/PrincetonUniversity/GPEC. When users reference "the Fortran code", "the original GPEC", "Fortran GPEC", or "legacy VACUUM", they are referring to that Fortran codebase — not a runtime dependency of this Julia implementation. This Julia implementation reimplements and extends GPEC's functionality with improved performance and maintainability. -**Local GPEC Repository**: For code conversion or comparison with the original Fortran implementation, check for a local GPEC repository at `~/Code/gpec`. If not found at this location, ask the user for the correct path. +**Local GPEC Repository**: For code conversion or comparison with the original Fortran implementation, a local clone of the Fortran GPEC repository is needed; ask the user for its path. ## Key References diff --git a/Project.toml b/Project.toml index fd2c512d..478b05c1 100644 --- a/Project.toml +++ b/Project.toml @@ -25,6 +25,7 @@ PlotlyJS = "f0f68f2c-4968-5e81-91da-67840de0976a" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" @@ -57,14 +58,11 @@ QuadGK = "2.11.3" Roots = "2.2.13" SparseArrays = "1" SpecialFunctions = "2.5.1" +Random = "1" StaticArrays = "1.9.15" Statistics = "1" TOML = "1" Test = "1" julia = "1.11" -[extras] -Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" -[targets] -test = ["Random"] diff --git a/benchmarks/benchmark_coil_ForcingTerms_against_fortran.jl b/benchmarks/benchmark_coil_ForcingTerms_against_fortran.jl index d50fe72f..22e2b193 100644 --- a/benchmarks/benchmark_coil_ForcingTerms_against_fortran.jl +++ b/benchmarks/benchmark_coil_ForcingTerms_against_fortran.jl @@ -13,7 +13,8 @@ runs the Julia Biot-Savart pipeline, and produces a 4-panel diagnostic figure: Usage: julia --project=. benchmarks/check_coil_pipeline.jl [/path/to/fortran/run] -Default run directory: ~/Code/gpec/docs/examples/DIIID_ideal_example +The Fortran run directory (a DIII-D ideal example run) is taken from the command +line, or from the GPEC_FORTRAN_DIR environment variable if no argument is given. """ using GeneralizedPerturbedEquilibrium @@ -149,7 +150,8 @@ end # Main # --------------------------------------------------------------------------- -run_dir = length(ARGS) > 0 ? ARGS[1] : expanduser("~/Code/gpec/docs/examples/DIIID_ideal_example") +run_dir = length(ARGS) > 0 ? ARGS[1] : get(ENV, "GPEC_FORTRAN_DIR", "") +isempty(run_dir) && error("No Fortran run directory given: pass it as the first argument or set GPEC_FORTRAN_DIR") isdir(run_dir) || error("Run directory not found: $run_dir") mkpath(OUTPUT_DIR) diff --git a/benchmarks/benchmark_diiid_ideal_ntv_torque.jl b/benchmarks/benchmark_diiid_ideal_ntv_torque.jl index 53e7c595..2263e78c 100644 --- a/benchmarks/benchmark_diiid_ideal_ntv_torque.jl +++ b/benchmarks/benchmark_diiid_ideal_ntv_torque.jl @@ -20,8 +20,8 @@ the Fortran example directory. Usage: julia --project=. benchmarks/benchmark_diiid_ideal_ntv_torque.jl [fortran_dir] - fortran_dir defaults to \$GPEC_FORTRAN_DIIID or - ~/Code/gpec/docs/examples/DIIID_ideal_example. + fortran_dir is a Fortran GPEC DIII-D ideal example run directory, taken from + the command line or from \$GPEC_FORTRAN_DIIID. One of the two is required. """ using Printf @@ -36,7 +36,9 @@ const KF = GPE.KineticForces const Eq = GPE.Equilibrium const PE = GPE.PerturbedEquilibrium -const DEFAULT_FORTRAN_DIR = expanduser("~/Code/gpec/docs/examples/DIIID_ideal_example") +"Fortran GPEC DIII-D ideal example run directory, from the environment (no on-disk default)." +default_fortran_dir() = get(() -> error("Set GPEC_FORTRAN_DIIID, or pass the Fortran run directory as the first argument"), + ENV, "GPEC_FORTRAN_DIIID") """ discover_inputs(fortran_dir) → (; eq_file, kin_file, xclebsch, pentrc_nc) @@ -218,7 +220,7 @@ end _p(args...) = (println(stderr, args...); flush(stderr)) _pf(fmt, args...) = (print(stderr, Printf.format(Printf.Format(fmt), args...)); flush(stderr)) -function run_benchmark(fortran_dir::String=DEFAULT_FORTRAN_DIR) +function run_benchmark(fortran_dir::String=default_fortran_dir()) _p("=" ^ 70) _p(" DIIID Kinetic Benchmark: Julia KineticForces vs Fortran PENTRC") _p(" Fortran example: $fortran_dir") @@ -479,6 +481,6 @@ end # Run only when invoked as a script (so other benchmarks can `include()` this # file to reuse `load_fortran_xclebsch` without triggering the full run). if abspath(PROGRAM_FILE) == @__FILE__ - fortran_dir = length(ARGS) >= 1 ? ARGS[1] : get(ENV, "GPEC_FORTRAN_DIIID", DEFAULT_FORTRAN_DIR) + fortran_dir = length(ARGS) >= 1 ? ARGS[1] : default_fortran_dir() results = run_benchmark(fortran_dir) end diff --git a/benchmarks/benchmark_diiid_kinetic_stability.jl b/benchmarks/benchmark_diiid_kinetic_stability.jl index 9d89690c..b5465f1a 100644 --- a/benchmarks/benchmark_diiid_kinetic_stability.jl +++ b/benchmarks/benchmark_diiid_kinetic_stability.jl @@ -24,8 +24,8 @@ the non-Hermitian FKG reduction and the inner kinetic-matrix tolerances. Usage: julia --project=. benchmarks/benchmark_diiid_kinetic_stability.jl [fortran_dir] - fortran_dir defaults to \$GPEC_FORTRAN_DIIID_DCON or - ~/Code/gpec/docs/examples/DIIID_kinetic_example. + fortran_dir is a Fortran GPEC DIII-D kinetic example run directory, taken from + the command line or from \$GPEC_FORTRAN_DIIID_DCON. One of the two is required. """ using Printf @@ -36,7 +36,9 @@ using GeneralizedPerturbedEquilibrium const GPE = GeneralizedPerturbedEquilibrium const AnalysisFFS = GPE.Analysis.ForceFreeStates -const DEFAULT_FORTRAN_DIR = expanduser("~/Code/gpec/docs/examples/DIIID_kinetic_example") +"Fortran GPEC DIII-D kinetic example run directory, from the environment (no on-disk default)." +default_fortran_dir() = get(() -> error("Set GPEC_FORTRAN_DIIID_DCON, or pass the Fortran run directory as the first argument"), + ENV, "GPEC_FORTRAN_DIIID_DCON") """ discover_inputs(fortran_dir) → (; eq_file, kin_file, dcon_nc) @@ -174,7 +176,7 @@ end _p(args...) = (println(stderr, args...); flush(stderr)) _pf(fmt, args...) = (print(stderr, Printf.format(Printf.Format(fmt), args...)); flush(stderr)) -function run_benchmark(fortran_dir::String=DEFAULT_FORTRAN_DIR) +function run_benchmark(fortran_dir::String=default_fortran_dir()) _p("=" ^ 70) _p(" DIIID Kinetic-DCON Benchmark: Julia KF→FFS vs Fortran DCON") _p(" Fortran example: $fortran_dir") @@ -253,7 +255,6 @@ end # Run only when invoked as a script. if abspath(PROGRAM_FILE) == @__FILE__ - fortran_dir = length(ARGS) >= 1 ? ARGS[1] : - get(ENV, "GPEC_FORTRAN_DIIID_DCON", DEFAULT_FORTRAN_DIR) + fortran_dir = length(ARGS) >= 1 ? ARGS[1] : default_fortran_dir() run_benchmark(fortran_dir) end diff --git a/benchmarks/benchmark_solovev_kinetic_stability.jl b/benchmarks/benchmark_solovev_kinetic_stability.jl index 3ced63df..4ab64e5a 100644 --- a/benchmarks/benchmark_solovev_kinetic_stability.jl +++ b/benchmarks/benchmark_solovev_kinetic_stability.jl @@ -31,10 +31,9 @@ the pre-rewrite (`ximag`) value of 34.176 was not. Usage: julia --project=. benchmarks/benchmark_solovev_kinetic_stability.jl -Environment: - GPEC_FORTRAN_DCON path to Fortran `dcon` (default ~/Code/gpec/bin/dcon) +Environment (both required): + GPEC_FORTRAN_DCON path to the Fortran `dcon` executable GPEC_FORTRAN_SOLDIR Fortran solovev_kinetic_example dir used as deck template - (default ~/Code/gpec/docs/examples/solovev_kinetic_example) """ using Printf @@ -46,9 +45,8 @@ const GPE = GeneralizedPerturbedEquilibrium const REPO = normpath(joinpath(@__DIR__, "..")) const JULIA_FIXTURE = joinpath(REPO, "test", "test_data", "regression_solovev_kinetic_calculated") -const FORTRAN_DCON = get(ENV, "GPEC_FORTRAN_DCON", expanduser("~/Code/gpec/bin/dcon")) -const FORTRAN_SOLDIR = get(ENV, "GPEC_FORTRAN_SOLDIR", - expanduser("~/Code/gpec/docs/examples/solovev_kinetic_example")) +const FORTRAN_DCON = get(ENV, "GPEC_FORTRAN_DCON", "") +const FORTRAN_SOLDIR = get(ENV, "GPEC_FORTRAN_SOLDIR", "") _p(args...) = (println(stderr, args...); flush(stderr)) @@ -86,7 +84,7 @@ fields that must match the Julia regression fixture (read from its `gpec.toml` and `sol.toml`). Converts the fixture's `kinetic.dat` into Fortran `kinetic.txt`. """ function build_matched_fortran_deck(workdir::String) - isdir(FORTRAN_SOLDIR) || error("Fortran deck template not found: $FORTRAN_SOLDIR") + isdir(FORTRAN_SOLDIR) || error("Fortran deck template not found: '$FORTRAN_SOLDIR' (set GPEC_FORTRAN_SOLDIR)") for f in ("sol.in", "equil.in", "dcon.in", "pentrc.in", "vac.in") cp(joinpath(FORTRAN_SOLDIR, f), joinpath(workdir, f); force=true) end diff --git a/docs/DOC_STANDARD.md b/docs/DOC_STANDARD.md new file mode 100644 index 00000000..a1af2043 --- /dev/null +++ b/docs/DOC_STANDARD.md @@ -0,0 +1,126 @@ +# GPEC documentation standard + +This is the authoring standard for the GPEC user documentation under +`docs/src/`. The goal is that each physics module page reads like a compact +*Journal of Computational Physics* methods paper: the governing equations, the +numerical method, and the validation that the method works — with figures whose +provenance is explicit. Those module pages **are** built by Documenter and +published on the public documentation site. + +This standard file (`docs/DOC_STANDARD.md`) is itself the one exception: it lives +outside `docs/src/` and is not listed in `make.jl`, so Documenter does not +render it to the website. It is a contributor-facing reference, read here in the +repository — which is why `developer_notes.md` links to it by GitHub URL rather +than as a site-internal page. + +The `Ballooning and Mercier Local Stability` and `Inner Layer Module` pages are +the reference exemplars; new or reworked module pages should follow their shape. + +## Module-page template + +A module reference page (`docs/src/.md`) has these sections, in order. +Not every section is large for every module, but the skeleton is the same. + +1. **Overview** — what the module computes and where it sits in the + equilibrium → stability → perturbed-equilibrium pipeline. +2. **Governing equations** — the physics the module discretizes, written as + numbered display equations. Each equation or block cites the paper it comes + from (the PDFs in `docs/resources/`), e.g. *(GWP2016 Eq. 11)*. Citing + equation numbers is required, not optional — it is what makes the + implementation auditable against theory. +3. **Numerical method** — the discretization and algorithm: what is expanded in + what basis, how the linear/eigenvalue/BVP problem is assembled and solved, + and the *design choices* that matter (why this contour, this ordering, this + precision). This is the part that distinguishes a methods page from an API + dump. +4. **Validation & benchmarks** — evidence the method is correct: a table of + pinned values cross-checked against an independent code or analytic result, + and figures (see below) showing convergence, invariance, or agreement. +5. **Practical usage** — the configuration keys and the public entry points a + user calls, with a short worked example. +6. **API Reference** — an `@autodocs` block for the module so every exported + docstring is rendered and `checkdocs=:exports` stays satisfied. + +## Figures + +### Where figures live + +- **Content figures** (anything shown on a page) live in + `docs/src/figures//`. Each committed `.png` sits next to the + script `make_.jl` that produced it — script and figure are always + committed together, so it is always clear how a figure was made. +- `docs/src/assets/` is reserved for **Documenter chrome only** (the site logo, + themes, custom CSS). Content figures never go there. +- The shared helper `docs/figure_tools.jl` (outside `src/`, so it is not + published) provides `save_doc_figure`, provenance stamping, the manifest, and + the `step_series` spectrum helper. Every `make_.jl` `include`s it. + +### Where a figure's data comes from + +Documentation figures get their **own** generator scripts under +`docs/src/figures/`. They are not forced to route through `benchmarks/`. When a +figure genuinely *is* a benchmark or Fortran cross-check result, its +`make_.jl` may call the relevant `benchmarks/` script, and its `depends` +list records that coupling — but the committed PNG and a runnable script still +live together under `docs/src/figures//`. + +`benchmarks/` remains the home of performance comparisons and Fortran +cross-checks (its outputs are git-ignored and may be heavy). `docs/src/figures/` +is the home of the small, purpose-built, committed figures that illustrate a +page. + +### Figures do not build with the docs + +Figure scripts run **manually**, in the root project: + +```bash +julia --project=. docs/src/figures//make_.jl +``` + +They are **never** run at Documenter build time. The build only embeds the +committed PNGs, so `docs/Project.toml` stays minimal (no `Plots`) and the build +stays fast. This is also why we do not regenerate every figure on every docs +change — see the regeneration policy below. + +### Provenance + +`save_doc_figure` stamps a small `GPEC · ` mark in the corner +of every figure (so a figure's age is visible when browsing the rendered site) +and writes a machine-readable entry to `docs/src/figures/manifest.toml`: + +```toml +[inner_layer.rotated_ray_contour] +script = "make_rotated_ray_contour.jl" +commit = "3a0837e3" +date = "2026-07-08" +depends = ["src/InnerLayer/GGJ/Ray.jl"] +``` + +`depends` lists the repo-relative source files whose *numbers* the figure +visualizes. Generate figures as the last step before committing a change; the +stamp then reflects the state they were made against (a `-dirty` suffix marks a +figure generated with an uncommitted working tree, which is honest, not an +error). + +### When to regenerate a figure + +Regenerate a figure (and its manifest entry) only when: + +1. a file in its `depends` list changed the numbers it shows — the regression + harness is the trigger: if a tracked quantity for that module moved, its + figures are suspect; +2. the physics or method the figure illustrates changed; or +3. the figure script itself changed. + +Do **not** regenerate for prose edits, formatting, or changes to unrelated +modules. Because each figure carries its generating commit and a `depends` +list, a reviewer (or a future `docs/check_figures.jl`) can tell when a figure +predates the last change to a file it depends on and is therefore stale. + +### Legacy figures + +A figure migrated from before this system, whose original generator was not +preserved and which cannot be faithfully reproduced from current code, is +registered with `register_legacy_figure` (`commit = "legacy"`, a `note` +explaining why). This is honest bookkeeping for genuinely irreproducible +figures — it is never a substitute for writing a generator for a new figure. diff --git a/docs/figure_tools.jl b/docs/figure_tools.jl new file mode 100644 index 00000000..ca8001be --- /dev/null +++ b/docs/figure_tools.jl @@ -0,0 +1,166 @@ +# figure_tools.jl +# +# Shared helpers for GPEC documentation figures. Every figure generator +# (`docs/src/figures//make_.jl`) `include`s this file and calls +# `save_doc_figure` to stamp, save, and register its output. See +# `docs/DOC_STANDARD.md` for the figure-organization and provenance policy. +# +# This file lives OUTSIDE `docs/src/` so Documenter does not publish it. The +# `make_.jl` scripts do live under `docs/src/figures/` (co-located with +# their PNGs) and are published alongside the figures for reproducibility. +# +# Figure scripts run MANUALLY in the root project, never at docs-build time: +# julia --project=. docs/src/figures//make_.jl +# so the Documenter build stays fast and Plots-free. + +using Plots +using Dates +using TOML + +const _DOCS_DIR = @__DIR__ # .../docs +const _REPO_ROOT = normpath(joinpath(_DOCS_DIR, "..")) # repo root +const FIGURES_ROOT = joinpath(_DOCS_DIR, "src", "figures") # docs/src/figures +const MANIFEST_PATH = joinpath(FIGURES_ROOT, "manifest.toml") + +""" + figure_provenance() -> (; commit, date, dirty) + +Short git commit and ISO date stamped onto every figure. `commit` is the +`git rev-parse --short HEAD` hash (reusing the idiom from +`benchmarks/benchmark_git_branches.jl`), with a `-dirty` suffix when the +working tree has uncommitted changes so a figure never claims a cleaner +provenance than it has. Falls back to `"unknown"` outside a git checkout. +""" +function figure_provenance() + date = string(Dates.today()) + try + h = strip(read(`git -C $(_REPO_ROOT) rev-parse --short HEAD`, String)) + dirty = !isempty(strip(read(`git -C $(_REPO_ROOT) status --porcelain`, String))) + return (; commit=(dirty ? "$h-dirty" : h), date=date, dirty=dirty) + catch + return (; commit="unknown", date=date, dirty=false) + end +end + +""" + stamp!(p; commit, date, subplot=length(p.subplots)) + +Annotate a small `GPEC · ` provenance mark in the bottom-right +corner of subplot `subplot` (default: the last panel), so the figure's age is +visible when browsing the rendered docs. Works on linear or log axes +(positions in data coordinates). Never throws: a plotting hiccup warns and +leaves the figure unstamped rather than losing the plot. +""" +function stamp!(p::Plots.Plot; commit::AbstractString, date::AbstractString, + subplot::Int=length(p.subplots)) + try + sp = p[subplot] + xl = Plots.xlims(sp) + yl = Plots.ylims(sp) + xpos = xl[2] # right edge (data coords) + ypos = yl[1] # bottom edge + txt = Plots.text("GPEC $(commit) · $(date)", 6, RGBA(0.5, 0.5, 0.5, 0.9), + :right, :bottom) + annotate!(sp, xpos, ypos, txt) + catch err + @warn "stamp! failed; saving figure without provenance mark" exception = err + end + return p +end + +""" + save_doc_figure(p, mod, name; script, depends, npx...) -> String + +Stamp `p` with the current git provenance, save it to +`docs/src/figures//.png`, and upsert its `manifest.toml` entry +(`script`, `commit`, `date`, `depends`). `depends` lists the repo-relative +source files whose numbers the figure visualizes — the regeneration policy in +`docs/DOC_STANDARD.md` keys off them. Returns the absolute PNG path (printed +so it can be opened directly). Extra keyword args are forwarded to `savefig`. +""" +function save_doc_figure(p::Plots.Plot, mod::AbstractString, name::AbstractString; + script::AbstractString, depends::AbstractVector{<:AbstractString}=String[]) + prov = figure_provenance() + stamp!(p; commit=prov.commit, date=prov.date) + + outdir = joinpath(FIGURES_ROOT, mod) + mkpath(outdir) + outpath = joinpath(outdir, name * ".png") + savefig(p, outpath) + + _update_manifest!(mod, name; script=String(script), commit=prov.commit, + date=prov.date, depends=String.(depends)) + + println("Saved doc figure: $(abspath(outpath))") + return abspath(outpath) +end + +# Load / merge-write the figures manifest, keyed [.]. Kept sorted and +# human-diffable; one entry per committed figure. +function _load_manifest() + return isfile(MANIFEST_PATH) ? TOML.parsefile(MANIFEST_PATH) : Dict{String,Any}() +end + +function _write_manifest(manifest) + open(MANIFEST_PATH, "w") do io + println(io, "# GPEC documentation-figure provenance manifest.") + println(io, "# One [.
] entry per committed PNG under docs/src/figures/.") + println(io, "# Regenerate a figure (and this entry) only when a file in its `depends`") + println(io, "# list changes the numbers it shows — see docs/DOC_STANDARD.md.") + println(io) + TOML.print(io, manifest; sorted=true) + end + return manifest +end + +function _update_manifest!(mod::AbstractString, name::AbstractString; + script::String, commit::String, date::String, depends::Vector{String}) + manifest = _load_manifest() + modtbl = get!(manifest, mod, Dict{String,Any}()) + modtbl[name] = Dict{String,Any}( + "script" => script, + "commit" => commit, + "date" => date, + "depends" => depends + ) + return _write_manifest(manifest) +end + +""" + register_legacy_figure(mod, name; note, depends=String[]) + +Record a manifest entry for a figure that predates this provenance system and +whose original generator script was not preserved (`commit = "legacy"`). Use +this only for migrated figures that cannot be faithfully reproduced from the +current code — never to skip writing a generator for a new figure. `note` +states why the figure is legacy (e.g. "compared the since-removed Mercier.jl"). +""" +function register_legacy_figure(mod::AbstractString, name::AbstractString; + note::AbstractString, depends::AbstractVector{<:AbstractString}=String[]) + manifest = _load_manifest() + modtbl = get!(manifest, mod, Dict{String,Any}()) + modtbl[name] = Dict{String,Any}( + "script" => "(not preserved)", + "commit" => "legacy", + "date" => "unknown", + "depends" => String.(depends), + "note" => String(note) + ) + _write_manifest(manifest) + println("Registered legacy figure: $(mod).$(name)") + return manifest +end + +""" + step_series(m_vals, amps) -> (m_ext, amp_ext) + +Canonical spectrum-plot helper (CLAUDE.md plotting convention): pad a zero on +each end of a discrete-mode series so a `seriestype=:steppre` plot draws clean +boxes that fall to zero at the edges. Centralized here so doc figures and +benchmarks share one definition. +""" +function step_series(m_vals, amps) + m_ext = [m_vals[1] - 1; m_vals; m_vals[end] + 1] + amp_ext = [0.0; amps; 0.0] + return m_ext, amp_ext +end diff --git a/docs/src/ballooning.md b/docs/src/ballooning.md index 6207b225..3313e535 100644 --- a/docs/src/ballooning.md +++ b/docs/src/ballooning.md @@ -854,7 +854,7 @@ numerically fragile. The comparison below shows that the `Baloo.f` style, which sets the displacement to vanish at ``\pm\theta_{\max}`` and computes ``\Delta'`` from the resulting solution, gives the most stable continuous profile. -![Comparison of numerical stability among methods for extracting ``\Delta'`` and ``c_{a1}``.](assets/ballooning/delta_prime_matching.png) +![Comparison of numerical stability among methods for extracting ``\Delta'`` and ``c_{a1}``.](figures/ballooning/delta_prime_matching.png) `Bal.jl` therefore adopts the `Baloo.f`-style boundary replacement for the ballooning ``\Delta'`` calculation. @@ -905,7 +905,7 @@ The previous standalone `Mercier.jl` path can therefore be removed from the loca stability calculation. The comparison below shows agreement between the previous Mercier calculation and the new `Bal.jl` calculation. -![Comparison of Mercier criterion calculations from `Mercier.jl` and `Bal.jl`.](assets/ballooning/mercier_comparison.png) +![Comparison of Mercier criterion calculations from `Mercier.jl` and `Bal.jl`.](figures/ballooning/mercier_comparison.png) The local-stability output now stores ballooning ``\Delta'`` in the fourth `locstab_fs` entry. In the HDF5 output this is written as diff --git a/docs/src/citations.md b/docs/src/citations.md index dfb65fb6..143b9a49 100644 --- a/docs/src/citations.md +++ b/docs/src/citations.md @@ -104,7 +104,7 @@ Provides an efficient algorithm for computing the full Δ' matrix (coupling betw > *Physics of Plasmas* **27**, 012506 (2020). > DOI: [10.1063/1.5134999](https://doi.org/10.1063/1.5134999) -Showcases a different basis that significantly aids convergence of the Galerkin solver used in the GGJ InnerLayer module. +Constructs the Wasow large-``x`` asymptotic basis (the `inps` kernel) that supplies the far-field boundary condition for every GGJ InnerLayer backend — the Galerkin solver and the rotated-contour collocation (`:ray`) solver alike. --- diff --git a/docs/src/developer_notes.md b/docs/src/developer_notes.md index 9a77b751..95f586da 100644 --- a/docs/src/developer_notes.md +++ b/docs/src/developer_notes.md @@ -16,6 +16,18 @@ CODE - TAG - Detailed message Where CODE is the module name (EQUIL, ForceFreeStates, VAC, PERTURBED EQUILIBRIUM, etc.) and TAG describes the type of change (WIP, MINOR, IMPROVEMENT, BUG FIX, NEW FEATURE, REFACTOR, CLEANUP, etc.). This format is used for compiling release notes — tags should be human-readable but are not enforced to a fixed set. +## Documentation standard + +Module reference pages follow a shared *Journal of Computational Physics*-style +structure (governing equations → numerical method → validation figures → API), +and every documentation figure is committed together with the script that made +it and a provenance stamp. The full policy — page template, figure +organization under `docs/src/figures//`, provenance, and the +regenerate-only-when-`depends`-change rule — is in +[`docs/DOC_STANDARD.md`](https://github.com/OpenFUSIONToolkit/GPEC/blob/develop/docs/DOC_STANDARD.md). +The `Ballooning and Mercier Local Stability` and `Inner Layer Module` pages are +the reference exemplars. + ## Regression Testing The regression harness **must be run on every pull request before merging into `develop`**. It is the project's primary safeguard for tracking how numerical results evolve across changes, so it is only useful if every PR exercises it. When you open a PR, paste the regression report into the PR thread so reviewers can see what moved (and what did not). If your change touches a quantity that is not yet tracked, add a new regression case — or extend an existing one — in the same PR. diff --git a/docs/src/assets/ballooning/delta_prime_matching.png b/docs/src/figures/ballooning/delta_prime_matching.png similarity index 100% rename from docs/src/assets/ballooning/delta_prime_matching.png rename to docs/src/figures/ballooning/delta_prime_matching.png diff --git a/docs/src/assets/ballooning/mercier_comparison.png b/docs/src/figures/ballooning/mercier_comparison.png similarity index 100% rename from docs/src/assets/ballooning/mercier_comparison.png rename to docs/src/figures/ballooning/mercier_comparison.png diff --git a/docs/src/figures/inner_layer/backend_regime_map.png b/docs/src/figures/inner_layer/backend_regime_map.png new file mode 100644 index 00000000..d86221bc Binary files /dev/null and b/docs/src/figures/inner_layer/backend_regime_map.png differ diff --git a/docs/src/figures/inner_layer/convergence_Sinvariance.png b/docs/src/figures/inner_layer/convergence_Sinvariance.png new file mode 100644 index 00000000..d8b55082 Binary files /dev/null and b/docs/src/figures/inner_layer/convergence_Sinvariance.png differ diff --git a/docs/src/figures/inner_layer/make_backend_regime_map.jl b/docs/src/figures/inner_layer/make_backend_regime_map.jl new file mode 100644 index 00000000..163b5049 --- /dev/null +++ b/docs/src/figures/inner_layer/make_backend_regime_map.jl @@ -0,0 +1,65 @@ +# make_backend_regime_map.jl +# +# Figure: accuracy of the three GGJ backends along the imaginary-Q axis. The +# :ray backend is the reference (its own error stays below ~1e-5 out to +# |Q| = 500 — see the convergence figure). The relative error of :shooting and +# :galerkin against it is plotted versus |Q|: :shooting holds to |Q| ~ 1, +# :galerkin to |Q| ~ 4, and both then lose all accuracy, while :ray continues. +# The crossovers of each curve with the 1% line are the practical reach of each +# method — this is what motivates the :ray backend. A few tens of solves; ~1 min. +# +# Run manually: julia --project=. docs/src/figures/inner_layer/make_backend_regime_map.jl + +include(joinpath(@__DIR__, "..", "..", "..", "figure_tools.jl")) + +using GeneralizedPerturbedEquilibrium +using LinearAlgebra +using Printf + +const IL = GeneralizedPerturbedEquilibrium.InnerLayer +const GGJ = IL.GGJ + +p = IL.q4_surface_benchmark() +γ(aq) = im * aq * GGJ.q0(p) # imaginary axis: Q = i|Q| + +# Worst-case relative error over the two parity components vs the :ray reference. +relerr(Δ, Δref) = maximum(abs.(Δ .- Δref) ./ abs.(Δref)) + +absQ = 10 .^ range(-1, log10(500); length=22) +ref = [IL.solve_inner(IL.GGJModel(; solver=:ray), p, γ(aq)) for aq in absQ] + +gal_x, gal_e = Float64[], Float64[] +shoot_x, shoot_e = Float64[], Float64[] +for (k, aq) in enumerate(absQ) + try + Δ = IL.solve_inner(IL.GGJModel(; solver=:galerkin), p, γ(aq)) + all(isfinite, Δ) && (push!(gal_x, aq); push!(gal_e, max(relerr(Δ, ref[k]), 1e-16))) + catch + end + aq > 10 && continue # :shooting is meaningless past its regime + try + Δ = IL.solve_inner(IL.GGJModel(; solver=:shooting), p, γ(aq)) + all(isfinite, Δ) && (push!(shoot_x, aq); push!(shoot_e, max(relerr(Δ, ref[k]), 1e-16))) + catch + end +end + +@printf("computed %d :ray refs, %d :galerkin, %d :shooting points\n", + length(absQ), length(gal_x), length(shoot_x)) + +plt = plot(; size=(820, 560), legend=:bottomright, xscale=:log10, yscale=:log10, + xlabel="|Q| (imaginary axis, Q = i|Q|)", + ylabel="relative error vs :ray reference", + title="Backend accuracy along the imaginary-Q axis (q=4 surface)", + ylims=(1e-7, 3e0), left_margin=9Plots.mm, bottom_margin=5Plots.mm) + +plot!(plt, shoot_x, shoot_e; lw=2, marker=:diamond, ms=5, color=:seagreen, label=":shooting") +plot!(plt, gal_x, gal_e; lw=2, marker=:circle, ms=5, color=:orange, label=":galerkin") +hline!(plt, [1e-2]; color=:red, ls=:dash, lw=1.5, label="1% accuracy") +annotate!(plt, 8.0, 3e-6, + Plots.text(":ray reference\n(error < 1e-5 to |Q| = 500)", 8, :center, RGB(0.2, 0.3, 0.6))) + +save_doc_figure(plt, "inner_layer", "backend_regime_map"; + script=basename(@__FILE__), + depends=["src/InnerLayer/GGJ/Ray.jl", "src/InnerLayer/GGJ/Galerkin.jl", + "src/InnerLayer/GGJ/Shooting.jl"]) diff --git a/docs/src/figures/inner_layer/make_convergence_Sinvariance.jl b/docs/src/figures/inner_layer/make_convergence_Sinvariance.jl new file mode 100644 index 00000000..3e2a71b2 --- /dev/null +++ b/docs/src/figures/inner_layer/make_convergence_Sinvariance.jl @@ -0,0 +1,45 @@ +# make_convergence_Sinvariance.jl +# +# Figure: honest numerical error bar on the matching data Δ at Q = 500i. The +# `delta_convergence` battery re-solves with each numerical knob perturbed on +# an independent axis (contour angle θ, spectral order p, series order/radius, +# refinement depth, march tolerance, handoff radius, purification e-folds). Δ +# is an analytic invariant of the contour, so the worst-case spread across +# these orthogonal perturbations bounds the true error. Eight solves; ~30 s. +# +# Run manually: julia --project=. docs/src/figures/inner_layer/make_convergence_Sinvariance.jl + +include(joinpath(@__DIR__, "..", "..", "..", "figure_tools.jl")) + +using GeneralizedPerturbedEquilibrium +using Printf + +const IL = GeneralizedPerturbedEquilibrium.InnerLayer + +p = IL.q4_surface_benchmark() +Q = 500.0im + +conv = IL.delta_convergence(p, Q; verbose=false) +names = [r.name for r in conv.table] +d1 = [max(r.d1, 1e-16) for r in conv.table] # relative change of Δ_odd +d2 = [max(r.d2, 1e-16) for r in conv.table] # relative change of Δ_even +n = length(names) + +plt = plot(; size=(820, 540), legend=:topright, yscale=:log10, + xticks=(1:n, names), xrotation=30, ylims=(1e-12, 1e-2), + ylabel="relative change of Δ per knob", xlabel="perturbed numerical knob", + title=@sprintf("Convergence & contour-invariance of Δ, Q = %gi (q=4)", imag(Q)), + left_margin=10Plots.mm, bottom_margin=14Plots.mm) + +scatter!(plt, 1:n, d1; marker=:circle, ms=7, color=1, label="δΔ_odd") +scatter!(plt, 1:n, d2; marker=:diamond, ms=7, color=2, label="δΔ_even") + +# Worst-case spread — the reported error bar on each parity. +hline!(plt, [conv.spread[1]]; color=1, ls=:dash, lw=1.5, + label=@sprintf("worst-case Δ_odd %.1e", conv.spread[1])) +hline!(plt, [conv.spread[2]]; color=2, ls=:dash, lw=1.5, + label=@sprintf("worst-case Δ_even %.1e", conv.spread[2])) + +save_doc_figure(plt, "inner_layer", "convergence_Sinvariance"; + script=basename(@__FILE__), + depends=["src/InnerLayer/GGJ/Ray.jl"]) diff --git a/docs/src/figures/inner_layer/make_rotated_ray_contour.jl b/docs/src/figures/inner_layer/make_rotated_ray_contour.jl new file mode 100644 index 00000000..15abe167 --- /dev/null +++ b/docs/src/figures/inner_layer/make_rotated_ray_contour.jl @@ -0,0 +1,65 @@ +# make_rotated_ray_contour.jl +# +# Figure: the rotated integration ray in the complex layer-coordinate (x) plane. +# Shows why the :ray backend rotates the contour by θ = arg(Q)/4: the real-axis +# contour (θ = 0) runs straight through the pseudo-resonance at +# x² = −Q²(G + K F) (real and large on the imaginary-Q axis), while the rotated +# ray clears it. Geometry only — no BVP solve — so it is cheap to regenerate. +# +# Run manually: julia --project=. docs/src/figures/inner_layer/make_rotated_ray_contour.jl + +include(joinpath(@__DIR__, "..", "..", "..", "figure_tools.jl")) + +using GeneralizedPerturbedEquilibrium +using Printf + +const IL = GeneralizedPerturbedEquilibrium.InnerLayer +const GGJ = IL.GGJ + +# q = 4 rational-surface benchmark on the imaginary-Q axis (the regime that +# defeats :galerkin). Q is the scaled growth rate; θ makes the parabolic +# exponent real. +p = IL.q4_surface_benchmark() +Q = 500.0im +θ = angle(Q) / 4 # = π/8 = 22.5° + +# Pseudo-resonance location: x² = −Q²(G + K F). On the imaginary-Q axis this is +# real and large, so it sits ON the real-x contour. +x_pr = sqrt(-Q^2 * (p.G + p.K * p.F)) +xpr = real(x_pr) # imaginary part ≈ 0 here + +# Matching radius from the series-residual criterion along the ray (annotation). +S, _, ok = IL.pick_smax(p, Q; θ=θ) + +# Draw both contours out to ~1.8× the pseudo-resonance radius so it is in frame. +smax = 1.8 * xpr +s = range(0, smax; length=400) +ray = cis(θ) .* s # x = e^{iθ} s + +plt = plot(; size=(720, 560), legend=:topleft, framestyle=:zerolines, + xlabel="Re x", ylabel="Im x", aspect_ratio=:equal, + title="Rotated integration ray, Q = 500i (q=4 surface)", + left_margin=6Plots.mm, bottom_margin=4Plots.mm) + +# Real-axis contour (θ = 0) and the rotated ray. +plot!(plt, [0, smax], [0, 0]; lw=2, ls=:dash, color=:gray, + label="real-axis contour (θ = 0)") +plot!(plt, real.(ray), imag.(ray); lw=3, color=1, + label=@sprintf("rotated ray x = e^{iθ}s, θ = %.1f°", rad2deg(θ))) + +# Pseudo-resonance on the real axis — the point the rotation steps around. +scatter!(plt, [xpr], [0.0]; marker=:xcross, ms=9, color=:red, msw=3, + label="pseudo-resonance x² = −Q²(G+KF)") + +# Clearance note above the ray near the pseudo-resonance radius, plus the true +# matching radius. Positions are tied to xpr so they stay inside the frame. +annotate!(plt, 1.02 * xpr, 0.56 * xpr, + Plots.text("ray stays clear of the\nreal-axis pseudo-resonance", 8, :left, + RGB(0.15, 0.15, 0.15))) +annotate!(plt, 0.45 * xpr, -0.14 * xpr, + Plots.text(@sprintf("matching radius S ≈ %.2g%s", S, ok ? "" : " (series tol not met)"), + 8, :gray, :left)) + +save_doc_figure(plt, "inner_layer", "rotated_ray_contour"; + script=basename(@__FILE__), + depends=["src/InnerLayer/GGJ/Ray.jl"]) diff --git a/docs/src/figures/inner_layer/make_solution_profiles.jl b/docs/src/figures/inner_layer/make_solution_profiles.jl new file mode 100644 index 00000000..2165f358 --- /dev/null +++ b/docs/src/figures/inner_layer/make_solution_profiles.jl @@ -0,0 +1,60 @@ +# make_solution_profiles.jl +# +# Figure: the inner-layer fields (Ψ, Ξ, Υ) reconstructed on the rotated ray for +# the q=4 benchmark surface. The collocation solution on [0, s_m] (solid) and +# the analytic u_small + Δ·u_big asymptotic representation for s ≥ S (dashed) +# share the same power-law tail — the numeric↔asymptotic overlap the outer- +# region matching relies on. One BVP solve; a few minutes at large imaginary Q. +# +# Run manually: julia --project=. docs/src/figures/inner_layer/make_solution_profiles.jl + +include(joinpath(@__DIR__, "..", "..", "..", "figure_tools.jl")) + +using GeneralizedPerturbedEquilibrium +using Printf + +const IL = GeneralizedPerturbedEquilibrium.InnerLayer +const GGJ = IL.GGJ + +p = IL.q4_surface_benchmark() +Q = 2.0im # imaginary axis, beyond the |Q|≪1 shooting regime, +# small enough that the collocation domain reaches +# the series radius directly (no march) — so the +# numeric and asymptotic segments join seamlessly. +isol = 1 # "odd" parity solution (Ψ'(0)=Ξ(0)=Υ(0)=0) + +res = IL.solve_ray(p, Q) +s_m = res.breaks[end] +@printf("solve_ray: Q=%s θ=%.3f S=%.4g s_m=%.4g resid=%.1e\n", + string(Q), res.θ, res.S, s_m, res.resid) + +# Collocation solution over the BVP domain [0, s_m]. +prof = IL.solution_profile(res; npc=8) + +# Analytic tail on [S, few·S] where the inps series is trusted. +srange = 10 .^ range(log10(res.S), log10(4 * res.S); length=80) +asy = IL.asymptotic_profile(p, res, srange) + +fields = ((:Ψ, prof.Ψ, asy.Ψ, 1), (:Ξ, prof.Ξ, asy.Ξ, 2), (:Υ, prof.Υ, asy.Υ, 3)) + +plt = plot(; size=(760, 560), legend=:bottomleft, xscale=:log10, yscale=:log10, + ylims=(1e-3, 1e2), + xlabel="ray parameter s (x = e^{iθ}s)", ylabel="|field| (odd parity)", + title=@sprintf("Inner-layer fields on the rotated ray, Q = %gi (q=4)", imag(Q)), + left_margin=8Plots.mm, bottom_margin=5Plots.mm) + +for (nm, col, acol, ci) in fields + plot!(plt, prof.s[2:end], abs.(col[2:end, isol]); lw=2.5, color=ci, + label="|$(nm)| collocation") + plot!(plt, asy.s, abs.(acol[:, isol]); lw=2, ls=:dash, color=ci, label="|$(nm)| asymptotic") +end + +# Matching radius: with no march the BVP edge s_m coincides with the series +# radius S, so numeric and asymptotic meet at a single point. +vline!(plt, [res.S]; color=:black, ls=:dot, lw=1, + label=@sprintf("S = s_m ≈ %.1f (match point)", res.S)) + +save_doc_figure(plt, "inner_layer", "solution_profiles"; + script=basename(@__FILE__), + depends=["src/InnerLayer/GGJ/Ray.jl", "src/InnerLayer/GGJ/RayAsymptotics.jl", + "src/InnerLayer/GGJ/InnerAsymptotics.jl"]) diff --git a/docs/src/figures/inner_layer/rotated_ray_contour.png b/docs/src/figures/inner_layer/rotated_ray_contour.png new file mode 100644 index 00000000..27884cf1 Binary files /dev/null and b/docs/src/figures/inner_layer/rotated_ray_contour.png differ diff --git a/docs/src/figures/inner_layer/solution_profiles.png b/docs/src/figures/inner_layer/solution_profiles.png new file mode 100644 index 00000000..3a1c2ef5 Binary files /dev/null and b/docs/src/figures/inner_layer/solution_profiles.png differ diff --git a/docs/src/figures/manifest.toml b/docs/src/figures/manifest.toml new file mode 100644 index 00000000..2ad9ef5d --- /dev/null +++ b/docs/src/figures/manifest.toml @@ -0,0 +1,43 @@ +# GPEC documentation-figure provenance manifest. +# One [.
] entry per committed PNG under docs/src/figures/. +# Regenerate a figure (and this entry) only when a file in its `depends` +# list changes the numbers it shows — see docs/DOC_STANDARD.md. + + +[ballooning.delta_prime_matching] +commit = "legacy" +date = "unknown" +depends = ["src/ForceFreeStates/Ballooning.jl"] +note = "Predates the provenance system; compared methods for extracting the ballooning Delta_prime and c_a1. Original generator script not preserved." +script = "(not preserved)" + +[ballooning.mercier_comparison] +commit = "legacy" +date = "unknown" +depends = ["src/ForceFreeStates/Ballooning.jl"] +note = "Predates the provenance system; compared the since-removed standalone Mercier.jl against Bal.jl. Not reproducible from current code (Mercier.jl was removed once Bal.jl subsumed it)." +script = "(not preserved)" + +[inner_layer.backend_regime_map] +commit = "3a0837e3-dirty" +date = "2026-07-08" +depends = ["src/InnerLayer/GGJ/Ray.jl", "src/InnerLayer/GGJ/Galerkin.jl", "src/InnerLayer/GGJ/Shooting.jl"] +script = "make_backend_regime_map.jl" + +[inner_layer.convergence_Sinvariance] +commit = "3a0837e3-dirty" +date = "2026-07-08" +depends = ["src/InnerLayer/GGJ/Ray.jl"] +script = "make_convergence_Sinvariance.jl" + +[inner_layer.rotated_ray_contour] +commit = "3a0837e3-dirty" +date = "2026-07-08" +depends = ["src/InnerLayer/GGJ/Ray.jl"] +script = "make_rotated_ray_contour.jl" + +[inner_layer.solution_profiles] +commit = "3a0837e3-dirty" +date = "2026-07-08" +depends = ["src/InnerLayer/GGJ/Ray.jl", "src/InnerLayer/GGJ/RayAsymptotics.jl", "src/InnerLayer/GGJ/InnerAsymptotics.jl"] +script = "make_solution_profiles.jl" diff --git a/docs/src/inner_layer.md b/docs/src/inner_layer.md index cc80f52e..25c4416d 100644 --- a/docs/src/inner_layer.md +++ b/docs/src/inner_layer.md @@ -1,17 +1,286 @@ # Inner Layer Module -The InnerLayer module provides abstract scaffolding for resistive inner-layer -models used in matched asymptotic expansions for resistive MHD stability -analysis. It currently includes the GGJ (Glasser–Greene–Johnson) shooting -method for computing the inner-layer response. +The `InnerLayer` module solves the **resistive inner-region problem** of +matched-asymptotic resistive-MHD stability: the thin layer around a rational +surface where resistivity, inertia and the ideal singularity balance. Its +output is the pair of parity-projected matching data ``(\Delta_\mathrm{odd}, +\Delta_\mathrm{even})`` that the outer ideal region (DCON/GPEC) matches to in +order to obtain resistive growth rates, tearing ``\Delta'``, and resistive +interchange stability. -## InnerLayer +The module provides an abstract [`InnerLayerModel`](@ref InnerLayer.InnerLayerModel) interface and one +concrete model so far, the **GGJ (Glasser–Greene–Johnson)** layer, with three +interchangeable solver backends. Future inner-layer models (SLAYER, kinetic +layers) plug in through the same interface. + +Two source papers define the equations and the asymptotic construction, and are +cited by their equation numbers throughout the code and below: + +- **GWP2016** — A. H. Glasser, Z. R. Wang and J.-K. Park, *Computation of + resistive instabilities by matched asymptotic expansions*, Phys. Plasmas + **23**, 112506 (2016). +- **GW2020** — A. H. Glasser and Z. R. Wang, *Asymptotic solutions and + convergence studies of the resistive inner region equations*, Phys. Plasmas + **27**, 012506 (2020). + +## Governing equations + +At a rational surface the layer equations couple three fields — the perturbed +flux ``\Psi`` and two auxiliary layer variables ``\Xi`` and ``\Upsilon`` — as +functions of the scaled distance ``x`` from the surface. In the GWP2016 form +(Eq. 11 ≡ GW2020 Eq. 1) they are + +```math +\begin{aligned} +\Psi'' &= H\,\Upsilon' + Q\,(\Psi - x\,\Xi),\\[2pt] +\Xi'' &= \frac{1}{Q^{2}}\Big[\,Q x^{2}\,\Xi - Q x\,\Psi - (E+F)\,\Upsilon - H\,\Psi'\,\Big],\\[2pt] +\Upsilon''&= \frac{1}{Q}\Big[\,x^{2}\,\Upsilon - x\,\Psi - Q^{2}\big(G(\Xi-\Upsilon) - K(E\Xi + F\Upsilon + H\Psi')\big)\Big], +\end{aligned} +``` + +where ``E, F, G, H, K`` are the flux-surface-averaged equilibrium coefficients +of GWP2016 Eq. (A8) and ``Q`` is the dimensionless growth rate defined below. +These coefficients are collected in [`GGJParameters`](@ref InnerLayer.GGJParameters). + +### Dimensionless scales + +The layer is scaled by the local Lundquist number ``S = \tau_R/\tau_A`` (ratio +of resistive to Alfvén time). The displacement and growth-rate scales are + +```math +X_0 = S^{-1/3}, \qquad Q_0 = X_0/\tau_A \qquad \text{(GWP2016 Eqs. A14–A15)}, +``` + +so a physical growth rate ``\gamma`` maps to the scaled eigenvalue +``Q = \gamma/Q_0`` ([`inner_Q`](@ref InnerLayer.inner_Q)) that appears in the equations above, and +the scaled matching data are returned to physical units by the +``X_0^{-2\sqrt{-D_I}} = S^{\,2\sqrt{-D_I}/3}`` rescaling (together with the +``v_1^{\,2\sqrt{-D_I}}`` linear-scale factor; [`rescale_delta`](@ref InnerLayer.rescale_delta)). + +### Mercier index and matching exponents + +The premise of an inner-layer model is local Mercier stability. The Mercier +interchange index and its resistive counterpart are + +```math +D_I = E + F + H - \tfrac14, \qquad +D_R = E + F + H^2 = D_I + \big(H-\tfrac12\big)^2 \qquad \text{(GWP2016 A9–A10)}, +``` + +([`mercier_di`](@ref InnerLayer.mercier_di), [`mercier_dr`](@ref InnerLayer.mercier_dr)). For ``D_I < 0`` the large-``x`` +solutions are the two power laws ``x^{r_\pm}`` with the Frobenius exponents + +```math +r_\pm = \tfrac32 \pm \sqrt{-D_I}, \qquad p_1 \equiv \sqrt{-D_I} +\qquad \text{(GW2020 Eq. 49)}. +``` + +The matching datum ``\Delta`` is the amplitude of the large solution +``x^{r_+}`` when the small solution ``x^{r_-}`` is normalized to unity (i.e. the +ratio of large to small coefficient), in the physical (outer-matching) +normalization. Imposing the two parities at ``x=0`` (odd: +``\Psi'(0)=\Xi(0)=\Upsilon(0)=0``; even: ``\Psi(0)=\Xi'(0)=\Upsilon'(0)=0``) +gives the pair ``(\Delta_\mathrm{odd}, \Delta_\mathrm{even})`` of GWP2016 +Eqs. (34)–(35). + +## Numerical method + +!!! note "Benchmark surface used in the figures" + Every figure on this page is computed on the ``q = 4`` rational-surface + benchmark ([`q4_surface_benchmark`](@ref InnerLayer.q4_surface_benchmark)) — + a fixed set of GGJ layer coefficients (``S = \tau_R/\tau_A \approx 4.6\times10^{6}``, + ``D_I \approx -0.31``) frozen into code from the ``q = 4`` rational surface of + the DIII-D-like (TkMkr H-mode) equilibrium of the `resistive_resmets` + benchmark, at that case's per-surface ``\eta`` and ``\rho``. The well-conditioned + real-``Q`` cross-check in the Validation section instead uses the + Glasser & Wang (2020) Eq. (55) surface. + +### Solver backends + +The `GGJ` model exposes three solvers through the `solver` type parameter of +[`GGJModel`](@ref InnerLayer.GGJModel); all consume the same large-``x`` asymptotic basis and return +``\Delta`` in the same convention: + +| backend | method | robust regime | +|:--------|:-------|:--------------| +| `:shooting` | backward stable shoot from ``X_\mathrm{max}\to 0`` | ``\lvert Q\rvert \ll 1`` | +| `:galerkin` | Hermite-cubic finite-element (real axis) | drift onset ``\lvert Q\rvert \sim 1``; 1% error by ``\lvert Q\rvert \sim 4`` | +| `:ray` (default) | rotated-contour spectral-element collocation | ``\lvert Q\rvert \sim 500`` on/near the imaginary axis | + +The difficulty is the large-``|Q|`` imaginary-``Q`` axis itself, where rotation +and resistivity scans live: there the layer's pseudo-resonance at +``x^2 = -Q^2(G+KF)`` sits directly on the real-``x`` integration path and the +exponential dichotomy weakens, so both older backends degrade — `:shooting` +through the dichotomy of the backward shoot, `:galerkin` through real-axis +oscillation and the on-axis pseudo-resonance. (The poles of ``\Delta(Q)`` +itself lie on and near the **real**-``Q`` axis — they are the layer +eigenvalues; ``\Delta`` is smooth along the imaginary axis.) Measured against +the `:ray` reference along the imaginary axis, `:shooting` holds to +``|Q|\sim 1`` and `:galerkin` to ``|Q|\sim 4`` (the 1% error crossings) +before both lose all accuracy: + +![Relative error of the :shooting and :galerkin backends against the :ray reference along the imaginary-Q axis, on the q=4 benchmark surface. Each curve's crossing of the 1% line marks that method's practical reach.](figures/inner_layer/backend_regime_map.png) + +The `:ray` backend was written to reach the large-``|Q|`` imaginary-axis regime +(rotation and resistivity scans) where these fail. The remainder of this section +describes it. + +### Entire-solution formulation + +`:ray` works with the **plain** state ``v = (\Psi, \Xi, \Upsilon, \Psi', +\Xi', \Upsilon')`` and writes the layer equations as the first-order system + +```math +\frac{dv}{dx} = M(x)\,v . +``` + +The coefficient matrix (the `GGJ` internal `ode_matrix`) is **polynomial** in ``x``, so +``x = 0`` is an ordinary point — the ``x^{-2}/x^{-4}`` singularities of the +GW2020 Eq. (2) scaled form ``(x\Psi,\ \Psi'/x,\ \dots)`` are artifacts of that +scaling, not of the equations. Because the coefficients are entire, the system +continues analytically to complex ``x``, which is what makes the contour +rotation below legitimate. + +### The rotated ray + +The equations are continued onto the ray + +```math +x = e^{i\theta}\, s, \qquad s \in [0, S], \qquad \theta = \tfrac14\arg Q, +``` + +The angle ``\theta = \tfrac14\arg Q`` makes the parabolic-cylinder exponent of +the outer solutions exactly real and clears the pseudo-resonance at +``x^2 = -Q^2(G + K F)``, which on the imaginary-``Q`` axis is real and large and +therefore sits directly on the un-rotated (real-``x``) contour. Rotating the +contour lifts it off that point; ``\theta = 0`` recovers a real-axis solve. + +![The rotated integration ray in the complex layer-coordinate plane at Q = 500i. The real-axis contour runs through the pseudo-resonance x² = −Q²(G+KF); the ray at θ = arg(Q)/4 = 22.5° clears it.](figures/inner_layer/rotated_ray_contour.png) + +### Spectral-element collocation BVP + +On ``[0, s_m]`` the system is discretized by a **global Chebyshev +spectral-element collocation** boundary-value problem: right-biased (Radau-like) +collocation at the Chebyshev–Lobatto nodes of each cell, three parity conditions +at the ordinary-point origin ``s = 0``, and six matching conditions at the outer +edge, + +```math +v(S) - \Delta\, U_b - c_1 E_1 - c_2 E_2 = U_s , +``` + +with the matching datum ``\Delta`` and the two decaying-mode amplitudes +``c_1, c_2`` carried as **bordered unknowns**. Here ``U_s, U_b`` are the small +and large power-like solutions and ``E_{1,2}`` the forward-decaying exponential +pair. No quantity is ever propagated across the layer, so the exponential +dichotomy that limits the shooting backend never enters; the boundary condition +splits it exactly. + +The two parities differ only in their three ``s=0`` rows, i.e. by a rank-3 +update, so both are obtained from a **single sparse LU factorization** plus a +Woodbury correction (the `GGJ` internal `_solve_parities`) rather than two +factorizations. A residual-driven bisection refinement adds cells until the +collocation residual meets tolerance, with a roundoff-plateau guard. + +### Far-field boundary and the inward march + +The large-``x`` boundary data use the **same** `inps` Wasow asymptotic basis as +the other backends (GW2020 Eqs. 3–52), evaluated at complex ``x`` by +`RayAsymptotics.jl` and applied at the series radius ``S`` where the series is +trusted (residual below tolerance; [`pick_smax`](@ref InnerLayer.pick_smax)). For large ``|Q|`` the +trusted radius ``S`` can be far outside the collocation domain ``s_m``, so the +power-pair data are transported inward from ``S`` to ``s_m`` by an **L-stable +2-stage Radau IIA march** in the quotient modulo the decaying exponential pair — +the subspace in which ``\Delta`` is defined (the `GGJ` internal +`march_boundary`). An L-stable implicit method is essential: the decaying pair +grows under backward integration, so any explicit marcher is stability-limited +to ``O(\rho S^2)`` steps, while the Radau march *damps* the unresolvable +backward-growing directions instead of amplifying them. + +The damped-zone march runs in `Complex{Double64}` extended precision. At large +``S`` the near-parallel power-pair geometry amplifies the structured backward +error of the ill-conditioned implicit solves into ``\Delta``-mixing of order +``10^{-4}`` at ``|Q| = 500`` in `Float64`; extended precision removes this +floor, while the well-conditioned resolved band stays in `Float64`. + +The result is a seamless numeric↔asymptotic solution: the collocation solution +on ``[0, s_m]`` and the analytic ``u_\mathrm{small} + \Delta\,u_\mathrm{big}`` +representation for ``s \ge S`` share the same power-law tail — the overlap the +outer-region matching relies on. + +![Inner-layer fields Ψ, Ξ, Υ on the rotated ray for the q=4 surface at Q = 2i. The collocation solution (solid) joins the asymptotic representation (dashed) seamlessly at the match point S = s_m.](figures/inner_layer/solution_profiles.png) + +## Validation and benchmarks + +**Cross-check against the Fortran `rmatch` solver.** At the Glasser & Wang +(2020) Eq. (55) operating point ``Q = 0.1234`` (real) the `:galerkin` backend +reproduces the Fortran `rmatch deltac` solver to ``\sim 10^{-8}``: + +| quantity | `:galerkin` (= Fortran to ``10^{-8}``) | `:ray` | +|:---------|:---------------------------------------|:-------| +| ``\Delta_\mathrm{odd}`` | ``3.698368\times 10^{4}`` | ``3.69789\times 10^{4}`` | +| ``\Delta_\mathrm{even}`` | ``14.759721`` | ``14.759715`` | + +The large ``\lvert\Delta_\mathrm{odd}\rvert \sim 4\times10^{4}`` means this +operating point sits near a pole of ``\Delta(Q)``, where every solver's error +is amplified by the pole geometry; the ``1.3\times10^{-4}`` +`:galerkin`↔`:ray` difference in ``\Delta_\mathrm{odd}`` (versus +``4\times10^{-7}`` in ``\Delta_\mathrm{even}``) is consistent with that +amplification, not a defect of either backend. + +**Physical benchmark on the imaginary axis.** On the ``q=4`` rational-surface +benchmark ([`q4_surface_benchmark`](@ref InnerLayer.q4_surface_benchmark), ``S \approx 4.58\times10^6``, +``D_I \approx -0.312``) the `:ray` backend is pinned at ``Q = 500i``, a regime +entirely beyond `:galerkin`: + +| quantity | value at ``Q = 500i`` | +|:---------|:----------------------| +| ``\Delta_\mathrm{odd}`` | ``2.4720 + 13.3540\,i`` | +| ``\Delta_\mathrm{even}`` | ``0.13750 + 0.74275\,i`` | + +**Convergence and contour invariance.** Because ``\Delta`` is an analytic +invariant of the contour angle, re-solving with each numerical knob perturbed on +an independent axis ([`delta_convergence`](@ref InnerLayer.delta_convergence)) gives an honest error bar: at +``Q = 500i`` the worst-case spread across all knobs is ``\sim 5\times10^{-6}`` +for ``\Delta_\mathrm{odd}`` and ``\sim 6\times10^{-7}`` for +``\Delta_\mathrm{even}``. That is a single-machine error bar: across +machines/BLAS builds the absolute values reproduce to ``\sim 10^{-5}`` +relative (the damped-zone implicit solves carry platform-dependent structured +roundoff), which is why the table above quotes five significant figures. + +![Relative change of Δ at Q = 500i under independent perturbations of each numerical knob (contour angle, spectral order, series order/radius, refinement depth, march tolerance, handoff radius, purification depth). The worst-case spread is the reported error bar.](figures/inner_layer/convergence_Sinvariance.png) + +## Choosing a backend + +`:ray` is the **default** — `GGJModel()` constructs `GGJModel{:ray}()`. It is +the correct choice for ``|Q| \gtrsim 1`` and for any ``Q`` near the imaginary +axis. `:galerkin` remains available and may be faster for very small real +``|Q|``, but note the backends take disjoint numerical-knob keywords: pass +`GGJModel(solver=:galerkin)` explicitly when supplying the Galerkin `nx`/`xfac` +knobs — passing them to the default `:ray` backend throws a `MethodError`. + +```julia +using GeneralizedPerturbedEquilibrium.InnerLayer + +p = q4_surface_benchmark() # GGJParameters for the q=4 benchmark surface +γ = 500im * GGJ.q0(p) # physical rate placing Q on the imaginary axis at 500i + +Δ = solve_inner(GGJModel(), p, γ) # (Δ_odd, Δ_even) with the default :ray backend + +# Full result (raw Δ, contour, mesh, nodal solution, residuals) via solve_ray: +res = solve_ray(p, GGJ.inner_Q(p, γ)) +res.Δ, res.resid, res.bc_cond +``` + +## API Reference + +### InnerLayer ```@autodocs Modules = [GeneralizedPerturbedEquilibrium.InnerLayer] ``` -## GGJ +### GGJ ```@autodocs Modules = [GeneralizedPerturbedEquilibrium.InnerLayer.GGJ] diff --git a/examples/LAR_ideal_match_test/gpec.toml b/examples/LAR_ideal_match_test/gpec.toml new file mode 100644 index 00000000..07f67030 --- /dev/null +++ b/examples/LAR_ideal_match_test/gpec.toml @@ -0,0 +1,80 @@ +# LAR ideal-limit reference: large-aspect-ratio analytic TJ equilibrium, epsilon = a/R0 = 0.2, n = 1. +# +# The ideal-limit companion to LAR_resistive_match_test, identical except that gal_ideal_flag skips +# the resistive inner layer entirely and builds the bare (perfectly shielded) coil columns. Running +# the two side by side isolates what the inner layer contributes on this geometry. The per-surface +# resistivity, density and rotation inputs are omitted here because the ideal limit never reads them. +# This run terminates after the stability stage. + +[Equilibrium] +eq_type = "tj_analytic" # Type of the input 2D equilibrium (analytic large-aspect-ratio model) +jac_type = "hamada" # Coordinate system (hamada, pest, boozer, equal_arc, park, custom) +grid_type = "ldp" # Radial grid packing type (ldp = linear-derivative packing toward rationals) +psilow = 0.01 # Lower limit of normalized flux coordinate +psihigh = 0.995 # Upper limit of normalized flux coordinate +mpsi = 128 # Number of radial grid points (0 = auto-compute from psi_accuracy) +mtheta = 512 # Number of poloidal grid points + +[TJ_ANALYTIC_INPUT] +lar_r0 = 2.0 # Major radius R0 [m] +lar_a = 0.4 # Minor radius a [m] (epsilon = a/R0 = 0.2) +qc = 1.5 # On-axis safety factor +qa = 3.6 # Edge safety factor +pc = 0.01 # Normalized on-axis pressure +mu = 2.0 # Pressure peaking exponent: p2 = pc*(1-r^2)^mu +B0 = 12.0 # On-axis toroidal field [T] +ma = 128 # Radial grid points +mtau = 128 # Poloidal grid points + +[Wall] +shape = "conformal" # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath) +a = 20 # Distance from plasma (conformal) or shape parameter + +[ForceFreeStates] +HDF5_filename = "gpec_lar_ideal.h5" # Output filename for the stability-stage results +local_stability_flag = true # Perform local stability analysis (Mercier and ballooning) across the psi profile +vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes + +qlow = 1.02 # Integration initiated at q determined by min(q0, qlow) +qhigh = 3.6 # Integration terminated at q limit determined by min(qa, qhigh) +sing_start = 0 # Start integration at the sing_start'th rational from the axis (psilow) +nn_low = 1 # Smallest toroidal mode number to include +nn_high = 1 # Largest toroidal mode number to include +delta_mlow = 8 # Expands lower bound of Fourier harmonics +delta_mhigh = 8 # Expands upper bound of Fourier harmonics +mthvac = 960 # Number of points used in splines over poloidal angle at plasma-vacuum interface + +eulerlagrange_tolerance = 1e-12 # Relative tolerance for ODE integration of Euler-Lagrange equations +singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced +ucrit = 1e4 # Maximum fraction of solutions allowed before re-normalized +sing_order = 6 # Power-series order for the singular-surface asymptotics +save_interval = 3 # Save every Nth ODE step (1=all, 10=every 10th). Always saves near rational surfaces. + +use_parallel = true # Run parallel FM-propagator BVP path (unlocks singular/delta_prime_matrix) +parallel_threads = 2 # BVP thread cap (1 = serial/bit-deterministic; 2 is about +20% speedup) +populate_dense_xi = false # Dense axis-basis xi for the FFS HDF5 output. Not needed here: the gal-matched path builds its own dense xi. +set_psilim_via_dmlim = false # Keep psilim at psihigh (do not truncate at last_rational_q + dmlim) +dmlim = 0.2 # Truncate integration at (last_rational_q + dmlim) / n (only used when set_psilim_via_dmlim = true) +force_termination = true # Terminate after the stability stage (no perturbed-equilibrium section here) +write_outputs_to_HDF5 = true # Write stability outputs to HDF5 + +# Outer-region singular Galerkin Delta-prime solver with rpec coil columns. +gal_flag = true # Enable the outer-region singular Galerkin Delta-prime solve +gal_solver = "LU" # Banded solver: "LU" (zgbtrf/zgbtrs) or "cholesky" (zpbtrf/zpbtrs). rpec requires "LU". +gal_nx = 256 # Hermite-cubic elements per interval between singular surfaces +gal_nq = 6 # Gauss-Lobatto quadrature order per element +gal_pfac = 0.03 # Grid packing ratio toward singular surfaces (smaller = tighter packing) +gal_dx0 = 5e-4 # Resonant-element integration truncation distance from the rational (times 1/|n q'|) +gal_dx1 = 1e-3 # Resonant-element size (times 1/|n q'|) +gal_dx2 = 1e-3 # Extension-element size (times 1/|n q'|) +gal_cutoff = 10 # Number of elements carrying the large (driving) solution as a source term +gal_tol = 1e-10 # Resonant-cell adaptive-quadrature (QuadGK) tolerance +gal_gnstep = 5000 # Maximum resonant-integration steps +gal_dx1dx2_flag = true # Enable the special dx1/dx2 element sizing for resonant/extension cells +gal_sing_order = 6 # Base power-series order for the Galerkin singular asymptotics +gal_sing_order_ceiling = true # Auto-raise the series order by ceil(2*Re(alpha)) on high-Mercier-index surfaces +gal_rpec_flag = true # Append mpert coil-response columns (unit edge sources) to the solve, giving the delta_coil block + +# DRIVEN (RPEC) outer-inner matching in the ideal limit: the inner layer is skipped entirely. +gal_match_flag = true # Solve the coil-driven matched xi(psi) from the gal Delta-prime. Requires gal_rpec_flag = true. +gal_ideal_flag = true # true = build the IDEAL solution (skip the inner layer, bare coil columns) diff --git a/examples/LAR_resistive_match_test/gpec.toml b/examples/LAR_resistive_match_test/gpec.toml new file mode 100644 index 00000000..7fb9f565 --- /dev/null +++ b/examples/LAR_resistive_match_test/gpec.toml @@ -0,0 +1,87 @@ +# LAR resistive-match test: large-aspect-ratio analytic TJ equilibrium, epsilon = a/R0 = 0.2, n = 1. +# +# A different equilibrium from the DIII-D-like validation case, used to check that the +# outer-region Galerkin Delta-prime solve and the driven (RPEC) outer-inner match stay self-consistent +# on a second geometry: matching residual near machine precision and convergence under gal_sing_order +# refinement. For n = 1 with q in [1.5, 3.6] the rational surfaces are q = 2 and q = 3, so msing = 2 +# and the per-surface resistive arrays below have length 2 (core to edge). This run terminates after +# the stability stage. The companion LAR_ideal_match_test case is the ideal-limit reference. + +[Equilibrium] +eq_type = "tj_analytic" # Type of the input 2D equilibrium (analytic large-aspect-ratio model) +jac_type = "hamada" # Coordinate system (hamada, pest, boozer, equal_arc, park, custom) +grid_type = "ldp" # Radial grid packing type (ldp = linear-derivative packing toward rationals) +psilow = 0.01 # Lower limit of normalized flux coordinate +psihigh = 0.995 # Upper limit of normalized flux coordinate +mpsi = 128 # Number of radial grid points (0 = auto-compute from psi_accuracy) +mtheta = 512 # Number of poloidal grid points + +[TJ_ANALYTIC_INPUT] +lar_r0 = 2.0 # Major radius R0 [m] +lar_a = 0.4 # Minor radius a [m] (epsilon = a/R0 = 0.2) +qc = 1.5 # On-axis safety factor +qa = 3.6 # Edge safety factor +pc = 0.01 # Normalized on-axis pressure +mu = 2.0 # Pressure peaking exponent: p2 = pc*(1-r^2)^mu +B0 = 12.0 # On-axis toroidal field [T] +ma = 128 # Radial grid points +mtau = 128 # Poloidal grid points + +[Wall] +shape = "conformal" # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath) +a = 20 # Distance from plasma (conformal) or shape parameter + +[ForceFreeStates] +HDF5_filename = "gpec_lar_resistive.h5" # Output filename for the stability-stage results +local_stability_flag = true # Perform local stability analysis (Mercier and ballooning) across the psi profile +vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes + +qlow = 1.02 # Integration initiated at q determined by min(q0, qlow) +qhigh = 3.6 # Integration terminated at q limit determined by min(qa, qhigh) +sing_start = 0 # Start integration at the sing_start'th rational from the axis (psilow) +nn_low = 1 # Smallest toroidal mode number to include +nn_high = 1 # Largest toroidal mode number to include +delta_mlow = 8 # Expands lower bound of Fourier harmonics +delta_mhigh = 8 # Expands upper bound of Fourier harmonics +mthvac = 960 # Number of points used in splines over poloidal angle at plasma-vacuum interface + +eulerlagrange_tolerance = 1e-12 # Relative tolerance for ODE integration of Euler-Lagrange equations +singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced +ucrit = 1e4 # Maximum fraction of solutions allowed before re-normalized +sing_order = 6 # Power-series order for the singular-surface asymptotics +save_interval = 3 # Save every Nth ODE step (1=all, 10=every 10th). Always saves near rational surfaces. + +use_parallel = true # Run parallel FM-propagator BVP path (unlocks singular/delta_prime_matrix) +parallel_threads = 2 # BVP thread cap (1 = serial/bit-deterministic; 2 is about +20% speedup) +populate_dense_xi = false # Dense axis-basis xi for the FFS HDF5 output. Not needed here: the gal-matched path builds its own dense xi. +set_psilim_via_dmlim = false # Keep psilim at psihigh (do not truncate at last_rational_q + dmlim) +dmlim = 0.2 # Truncate integration at (last_rational_q + dmlim) / n (only used when set_psilim_via_dmlim = true) +force_termination = true # Terminate after the stability stage (no perturbed-equilibrium section here) +write_outputs_to_HDF5 = true # Write stability outputs to HDF5 + +# Outer-region singular Galerkin Delta-prime solver with rpec coil columns. +gal_flag = true # Enable the outer-region singular Galerkin Delta-prime solve +gal_solver = "LU" # Banded solver: "LU" (zgbtrf/zgbtrs) or "cholesky" (zpbtrf/zpbtrs). rpec requires "LU". +gal_nx = 256 # Hermite-cubic elements per interval between singular surfaces +gal_nq = 6 # Gauss-Lobatto quadrature order per element +gal_pfac = 0.03 # Grid packing ratio toward singular surfaces (smaller = tighter packing) +gal_dx0 = 5e-4 # Resonant-element integration truncation distance from the rational (times 1/|n q'|) +gal_dx1 = 1e-3 # Resonant-element size (times 1/|n q'|) +gal_dx2 = 1e-3 # Extension-element size (times 1/|n q'|) +gal_cutoff = 10 # Number of elements carrying the large (driving) solution as a source term +gal_tol = 1e-10 # Resonant-cell adaptive-quadrature (QuadGK) tolerance +gal_gnstep = 5000 # Maximum resonant-integration steps +gal_dx1dx2_flag = true # Enable the special dx1/dx2 element sizing for resonant/extension cells +gal_sing_order = 6 # Base power-series order for the Galerkin singular asymptotics +gal_sing_order_ceiling = true # Auto-raise the series order by ceil(2*Re(alpha)) on high-Mercier-index surfaces +gal_rpec_flag = true # Append mpert coil-response columns (unit edge sources) to the solve, giving the delta_coil block + +# DRIVEN (RPEC) outer-inner matching. Per-surface resistive inputs (core to edge) for the two +# rational surfaces q = 2 and q = 3: eta=8e-8, rho=3.3e-7 kg/m^3, rotation=1 Hz, Gamma=5/3. +gal_match_flag = true # Solve the coil-driven matched xi(psi) from the gal Delta-prime + the inner-layer Delta(Q). Requires gal_rpec_flag = true. +gal_ideal_flag = false # true = build the IDEAL solution (skip the inner layer, bare coil columns). false = the resistive matching below. +gal_inner_solver = "ray" # Inner-layer Delta backend: "ray" (rotated-contour collocation, robust at large |Q|) or "galerkin" +gal_eta = [8e-8, 8e-8] # Per-surface resistivity eta (core to edge) +gal_rho = [3.3e-7, 3.3e-7] # Per-surface mass density rho [kg/m^3] (core to edge) +gal_rotation = [1.0, 1.0] # Per-surface rotation frequency f [Hz] (core to edge); forced eigenvalue gamma_s = 2*pi*i*n*f +gal_gamma = 1.6666666666666667 # Ratio of specific heats Gamma for the resistive-layer coefficients (5/3) diff --git a/regression-harness/cases/ggj_ray_q500i.toml b/regression-harness/cases/ggj_ray_q500i.toml new file mode 100644 index 00000000..f5722ee4 --- /dev/null +++ b/regression-harness/cases/ggj_ray_q500i.toml @@ -0,0 +1,56 @@ +# Regression case: GGJ inner-layer rotated-ray collocation backend at Q = 500i. +# A "computed" case (kind = "computed", no example_dir) that calls solve_inner with the +# :ray backend on the physical q=4 benchmark surface at Q = 500i — a regime beyond the +# :galerkin backend — and tracks the Δ_odd / Δ_even matching coefficients. Guards the +# rotated-ray solver (and its extended-precision damped-zone march) against numerical +# drift across changes. Each [quantities.*] block names an HDF5 path in the computed +# output, how to extract it, and the noise floor for diffs. +[case] +name = "ggj_ray_q500i" +description = "GGJ inner-layer :ray backend, q=4 surface, Q = 500i (beyond :galerkin)" +kind = "computed" + +# Δ_odd / Δ_even from solve_inner(GGJModel(:ray), q4_surface_benchmark(), 500im·Q₀). +# Complex-valued; both real and imaginary parts are tracked separately. +# noise_threshold 1e-4 (absolute): the extended-precision damped-zone march carries +# BLAS/platform-dependent structured roundoff at the ~1e-4 absolute level on these +# O(10) components (~1e-5 relative); a tighter floor flags phantom cross-machine drift. +[quantities.delta_odd_real] +h5path = "ggj/delta_odd_real" +type = "real_scalar" +extract = "value" +label = "Re(Δ_odd)" +noise_threshold = 1e-4 +order = 10 + +[quantities.delta_odd_imag] +h5path = "ggj/delta_odd_imag" +type = "real_scalar" +extract = "value" +label = "Im(Δ_odd)" +noise_threshold = 1e-4 +order = 11 + +[quantities.delta_even_real] +h5path = "ggj/delta_even_real" +type = "real_scalar" +extract = "value" +label = "Re(Δ_even)" +noise_threshold = 1e-4 +order = 12 + +[quantities.delta_even_imag] +h5path = "ggj/delta_even_imag" +type = "real_scalar" +extract = "value" +label = "Im(Δ_even)" +noise_threshold = 1e-4 +order = 13 + +[quantities.runtime] +h5path = "" +type = "runtime" +extract = "value" +label = "Runtime (s)" +noise_threshold = 0.0 +order = 90 diff --git a/regression-harness/src/runner.jl b/regression-harness/src/runner.jl index 8705cb4d..84e80838 100644 --- a/regression-harness/src/runner.jl +++ b/regression-harness/src/runner.jl @@ -79,6 +79,31 @@ open(ARGS[2], "w") do f end """ +# GGJ rotated-ray backend at Q = 500i on the q=4 benchmark surface — a regime beyond the +# :galerkin backend. Builds the physical rate γ = 500i·Q₀ so inner_Q lands exactly on the +# imaginary axis at 500i, then writes the parity matching data. Runtime to ARGS[2] as usual. +const COMPUTED_GGJ_RAY_SCRIPT_TEMPLATE = """ +using Pkg +%INSTANTIATE% +using GeneralizedPerturbedEquilibrium +using GeneralizedPerturbedEquilibrium.InnerLayer +using HDF5 +p = q4_surface_benchmark() +γ = 500.0im * InnerLayer.GGJ.q0(p) +t_start = time() +Δ = solve_inner(GGJModel(solver=:ray), p, γ) +elapsed = time() - t_start +h5open(ARGS[1], "w") do fid + fid["ggj/delta_odd_real"] = real(Δ[1]) + fid["ggj/delta_odd_imag"] = imag(Δ[1]) + fid["ggj/delta_even_real"] = real(Δ[2]) + fid["ggj/delta_even_imag"] = imag(Δ[2]) +end +open(ARGS[2], "w") do f + println(f, elapsed) +end +""" + # Self-contained separatrix-finder regression (PR #296). Loads a fixed-boundary EFIT whose # computational box hugs the LCFS (eps=0.05 TokaMaker aspect-scan g-file): outside the prescribed # LCFS the coil-vacuum flux turns back above the boundary value before the grid edge, so the old @@ -139,6 +164,8 @@ Pick the subprocess script template for a `kind="computed"` case. function _computed_script_template(case_spec::CaseSpec) if case_spec.name == "ggj_reference" return COMPUTED_GGJ_SCRIPT_TEMPLATE + elseif case_spec.name == "ggj_ray_q500i" + return COMPUTED_GGJ_RAY_SCRIPT_TEMPLATE elseif case_spec.name == "efit_fixedbdy_separatrix" return COMPUTED_SEPARATRIX_SCRIPT_TEMPLATE end diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl index 1f8fe3bc..8315d63c 100644 --- a/src/ForceFreeStates/EulerLagrange.jl +++ b/src/ForceFreeStates/EulerLagrange.jl @@ -310,8 +310,10 @@ function compute_axis_init(ffit::FourFitVars, profiles::Equilibrium.ProfileSplin m22 = sf * kd * fi # Frobenius matrix A₀_j = ψ_low · M_j [Glasser 2016 Eq. 51] + #! format: off F_eig = eigen([psi_low*m11 psi_low*m12; psi_low*m21 psi_low*m22]) + #! format: on eig_vals = F_eig.values eig_vecs = F_eig.vectors diff --git a/src/ForceFreeStates/ForceFreeStatesStructs.jl b/src/ForceFreeStates/ForceFreeStatesStructs.jl index 7f727ea5..bc0b2f4a 100644 --- a/src/ForceFreeStates/ForceFreeStatesStructs.jl +++ b/src/ForceFreeStates/ForceFreeStatesStructs.jl @@ -194,6 +194,13 @@ A mutable struct holding internal state variables for stability calculations. raw 2msing×2msing BVP solution to produce the PEST3-compatible tearing parameter. """ delta_prime_matrix::Matrix{ComplexF64} = Matrix{ComplexF64}(undef, 0, 0) + """ + Edge coil-response matrix of shape (2msing × numpert_total). Column k is the resonant + small-solution response at each surface side to a unit source on edge poloidal mode k, + built by imposing the Eq. (37) rpec edge boundary condition on the Riccati BVP + (`_solve_bvp_edge_coil`). Empty unless the S-axis vacuum-edge BVP was assembled. + """ + delta_coil::Matrix{ComplexF64} = Matrix{ComplexF64}(undef, 0, 0) end """ @@ -302,6 +309,13 @@ A mutable struct containing control parameters for stability analysis, set by th # --- DRIVEN (RPEC) outer↔inner asymptotic matching (rmatch match_rpec port) --- gal_match_flag::Bool = false # enable the RPEC inner-layer matching: solve the coil-driven matched ξ(ψ) from the gal Δ′ + the inner-layer Δ(Q). Requires gal_rpec_flag=true. gal_ideal_flag::Bool = false # within the match, build the IDEAL solution: skip the inner-layer Δ, use bare coil columns (cout=0). Mirrors Fortran rmatch coil%ideal_flag (the EL reference). eta/rho/rotation ignored. + gal_inner_solver::String = "ray" # inner-layer Δ backend for the match: "ray" (rotated-contour collocation, certified Δ at the optimal θ = arg(Q)/4; robust for |Q| ≳ 1) or "galerkin" (Hermite-cubic inps; drifts for |Q| ≳ 1) + # --- Inner-layer "galerkin" backend knobs (used only when gal_inner_solver = "galerkin") --- + gal_inner_xfac::Float64 = 10.0 # asymptotic-matching radius multiplier (inps_xfac: xmax × 10) + gal_inner_nx::Int = 1280 # inner-layer grid cells (128 · xfac in the reference) + gal_inner_nq::Int = 5 # quadrature order per cell + gal_inner_cutoff::Int = 5 # cells carrying the large solution as driving term + gal_inner_kmax::Int = 8 # large-x asymptotic series order (↔ order_pow) gal_eta::Vector{Float64} = Float64[] # per-surface resistivity η (length msing, core→edge); Fortran rmatch `eta` gal_rho::Vector{Float64} = Float64[] # per-surface mass density ρ [kg/m³] (length msing, core→edge); Fortran rmatch `massden` gal_rotation::Vector{Float64} = Float64[] # per-surface rotation frequency f [Hz] (length msing, core→edge); forced eigenvalue γ_s = 2πi·n·f. Fortran rmatch `rotation` diff --git a/src/ForceFreeStates/Galerkin/GalerkinAssembly.jl b/src/ForceFreeStates/Galerkin/GalerkinAssembly.jl index 6d28ed53..3032e6ef 100644 --- a/src/ForceFreeStates/Galerkin/GalerkinAssembly.jl +++ b/src/ForceFreeStates/Galerkin/GalerkinAssembly.jl @@ -22,6 +22,11 @@ @inline sing_get_dua_gal(g::GalSingAsymp, z::Float64) = z >= 0 ? sing_get_dua_res(g.right, z) : (-1) .* sing_get_dua_res(g.left, -z) +# Leading-order ("cut") resonant basis, used to build the cut outer solution that the inner-layer +# solution is added to. Same two-sided convention as `sing_get_ua_gal`. +@inline sing_get_ua_gal_cut(g::GalSingAsymp, z::Float64) = + z >= 0 ? sing_get_ua_res_cut(g.right, z) : sing_get_ua_res_cut(g.left, -z) + @inline sing_get_ua_gal!(out, g::GalSingAsymp, z::Float64) = z >= 0 ? sing_get_ua_res!(out, g.right, z) : sing_get_ua_res!(out, g.left, -z) @inline function sing_get_dua_gal!(out, g::GalSingAsymp, z::Float64) diff --git a/src/ForceFreeStates/Galerkin/GalerkinMatch.jl b/src/ForceFreeStates/Galerkin/GalerkinMatch.jl index 21fc87ff..e9d91d85 100644 --- a/src/ForceFreeStates/Galerkin/GalerkinMatch.jl +++ b/src/ForceFreeStates/Galerkin/GalerkinMatch.jl @@ -45,6 +45,8 @@ function gal_match_rpec(ctrl::ForceFreeStatesControl, equil, intr::ForceFreeStat bpen = zeros(ComplexF64, msing, mcoil) inner_psi = Vector{Float64}[] # no inner layer in the ideal limit inner_xi = Matrix{ComplexF64}[] + inner_b = Matrix{ComplexF64}[] + inner_params = InnerLayer.GGJParameters[] # inner layer skipped in the ideal limit rpec_eig = zeros(ComplexF64, msing) residual = 0.0 else @@ -57,6 +59,8 @@ function gal_match_rpec(ctrl::ForceFreeStatesControl, equil, intr::ForceFreeStat sings, _, _ = gal_resonant_surfaces(intr, equil) length(sings) == msing || error("gal_match_rpec: re-derived $(length(sings)) surfaces, expected msing=$msing") + ctrl.gal_inner_solver in ("ray", "galerkin") || + error("gal_match_rpec: gal_inner_solver = \"$(ctrl.gal_inner_solver)\" (expected \"ray\" or \"galerkin\")") # --- inner-layer matching data Δ(Q) per surface (deltac_run; match.f) --- # solve_inner_profile returns the same Δ as solve_inner plus the reconstructed inner-layer field, @@ -71,25 +75,47 @@ function gal_match_rpec(ctrl::ForceFreeStatesControl, equil, intr::ForceFreeStat inner_psi = Vector{Vector{Float64}}(undef, msing) inner_odd = Vector{Vector{ComplexF64}}(undef, msing) # Ξ₁, antisymmetric across ψ_s inner_even = Vector{Vector{ComplexF64}}(undef, msing) # Ξ₂, symmetric across ψ_s + inner_bodd = Vector{Vector{ComplexF64}}(undef, msing) # b^ψ₁ = scale·resc·Ψ₁, symmetric (Ψ′₁(0)=0) + inner_beven = Vector{Vector{ComplexF64}}(undef, msing) # b^ψ₂ = scale·resc·Ψ₂, antisymmetric (Ψ₂(0)=0) + inner_params = Vector{InnerLayer.GGJParameters}(undef, msing) for i in 1:msing params = resist_eval(sings[i], equil, intr; eta=ctrl.gal_eta[i], rho=ctrl.gal_rho[i], gamma=ctrl.gal_gamma, ising=i) + inner_params[i] = params γ = 2π * im * nn * ctrl.gal_rotation[i] # forced eigenvalue; gal_rotation is f [Hz], γ = 2πi·n·f rpec_eig[i] = γ - # Inner-solve knobs matched to the Fortran rmatch deltac/inps reference (DELTAC_LIST): inps basis, - # inps_xfac=10 (xmax×10, grid nx = 128·10 = 1280), nq=5, cutoff=5, order_pow=8 (↔ kmax). - Δ, _, prof, _ = InnerLayer.GGJ.solve_inner_profile(params, γ; xfac=10.0, nx=1280, nq=5, cutoff=5, kmax=8) # (Δ₁, Δ₂) in deltac.f convention - deltar[i, 1] = Δ[1] - deltar[i, 2] = Δ[2] - x0i = InnerLayer.GGJ.x0(params) - resc = (params.v1 / x0i)^(0.5 + InnerLayer.GGJ.p1(params)) # deltac.f amplitude rescale - scale = chi1 * im * nn * sings[i].q1 * x0i # match.f b-field scaling - pen[i, 1] = scale * prof.Ψ[1, 1] * resc # layer center X=0, parity 1 (Ψ(0)≠0) - pen[i, 2] = scale * prof.Ψ[1, 2] * resc # parity 2 (Ψ(0)=0 ⇒ ~0) - xvar = prof.x .* (x0i / params.v1) # inner X → ψ-distance (deltac.f:1822) + # gal_inner_solver = "ray" (default) uses the rotated-ray backend (certified optimal-θ Δ + # + θ=0 real-axis profile re-solve with runtime certificate — see the InnerLayer + # solve_inner_profile(::GGJModel{:ray}, ...) docstring); "galerkin" uses the Hermite-FEM + # real-axis solve with the ctrl.gal_inner_* knobs (defaults match the Fortran rmatch + # deltac/inps reference). + inner = if ctrl.gal_inner_solver == "ray" + InnerLayer.solve_inner_profile(InnerLayer.GGJModel(; solver=:ray), params, γ) + else + InnerLayer.solve_inner_profile(InnerLayer.GGJModel(; solver=:galerkin), params, γ; + xfac=ctrl.gal_inner_xfac, nx=ctrl.gal_inner_nx, nq=ctrl.gal_inner_nq, cutoff=ctrl.gal_inner_cutoff, kmax=ctrl.gal_inner_kmax) + end + deltar[i, 1] = inner.Δ[1] + deltar[i, 2] = inner.Δ[2] + profΨ, profΞ, xg = inner.Ψ, inner.Ξ, inner.x + # Amplitude rescale: inner profiles are normalized in X = v1·δψ/x0 (inner_psi below); + # inner.rescale converts a big-branch amplitude to the outer δψ-normalization. + resc = inner.rescale + # b-field scaling, derived from the code's own outer convention (SingularCoupling): + # b_m = 2πi·χ₁·(m−nq)·ξ_m, m−nq = −n·q′·δψ, δψ = dψdx·X, resc·Ξ(X) ↔ ξ_m, + # and the far-field identity Ψ = X·Ξ (GWP2016 Eq. 16; Ψ is the normal-field variable, Eq. A17): + # b_m = −2πi·χ₁·n·q′·dψdx·resc·Ψ. + scale = -2π * chi1 * im * nn * sings[i].q1 * inner.dψdx + pen[i, 1] = scale * profΨ[1, 1] * resc # layer center X=0, parity 1 (Ψ(0)≠0) + pen[i, 2] = scale * profΨ[1, 2] * resc # parity 2 (Ψ(0)=0 ⇒ ~0) + xvar = xg .* inner.dψdx # inner X → ψ-distance (deltac.f:1822) inner_psi[i] = vcat(reverse(sings[i].psifac .- xvar), sings[i].psifac .+ xvar) - inner_odd[i] = resc .* vcat(reverse(.-prof.Ξ[:, 1]), prof.Ξ[:, 1]) # comp 2, parity 1 (odd: −left,+right) - inner_even[i] = resc .* vcat(reverse(prof.Ξ[:, 2]), prof.Ξ[:, 2]) # comp 2, parity 2 (even) + inner_odd[i] = resc .* vcat(reverse(.-profΞ[:, 1]), profΞ[:, 1]) # comp 2, parity 1 (odd: −left,+right) + inner_even[i] = resc .* vcat(reverse(profΞ[:, 2]), profΞ[:, 2]) # comp 2, parity 2 (even) + # b^ψ profiles on the same two-sided grid: Ψ is the normal-field variable (GWP2016 A17), + # b_m(δψ) = scale·resc·Ψ(X) throughout the layer (→ outer frozen-in relation via Ψ = XΞ). + inner_bodd[i] = (scale * resc) .* vcat(reverse(profΨ[:, 1]), profΨ[:, 1]) # parity 1: Ψ even + inner_beven[i] = (scale * resc) .* vcat(reverse(.-profΨ[:, 2]), profΨ[:, 2]) # parity 2: Ψ odd end # --- assemble the 4·msing matching system (match.f) --- @@ -135,6 +161,8 @@ function gal_match_rpec(ctrl::ForceFreeStatesControl, equil, intr::ForceFreeStat # Matched inner-layer ξ_ψ(ψ) per surface, per coil drive (match.f intotsol): odd parity weighted by # cin[2i] (cofin(2·ising)), even parity by cin[2i-1] (cofin(2·ising-1)). inner_xi = [inner_odd[i] * transpose(cin[2i, :]) .+ inner_even[i] * transpose(cin[2i-1, :]) for i in 1:msing] + # Matched inner-layer b^ψ(ψ) per surface, per coil drive. + inner_b = [inner_bodd[i] * transpose(cin[2i, :]) .+ inner_beven[i] * transpose(cin[2i-1, :]) for i in 1:msing] end # --- matched outer ξ/ξ′ per coil drive (match.f); ideal: cout=0 ⇒ bare coil column --- @@ -153,7 +181,57 @@ function gal_match_rpec(ctrl::ForceFreeStatesControl, equil, intr::ForceFreeStat end end - return GalMatchResult(cout, cin, xi, xi_deriv, deltar, bpen, inner_psi, inner_xi, rpec_eig, residual) + # --- composite inner-region solution: cut outer background + layer (match.f intotsol/intotsol_b) --- + # The layer solution alone carries only the resonant content it resolves; the smooth background + # removed by the cut has to be added back for the inner and outer solutions to overlap in the + # matching region. Without this the inner profile does not graft onto the outer eigenfunction. + if !isempty(inner_psi) + sols_cut = gal_result.solution.xi_cut + isempty(sols_cut) && error("gal_match_rpec: solution.xi_cut is empty — the cut solution is " * + "required to build the composite inner-region solution") + cut_range = gal_result.solution.cut_range + keep = .!gal_result.solution.issing + psi_keep = gal_result.solution.psi[keep] + chi1_c = 2π * equil.psio + for i in 1:msing + m_res = round(Int, nn * sings[i].q) + ires = m_res - intr.mlow + 1 + # Clip the layer to the window where the cut is active; outside it the cut removes + # nothing and the composite is undefined (Fortran writes no points there). + lo, hi = cut_range[i, 1], cut_range[i, 2] + inside = findall(p -> lo <= p <= hi, inner_psi[i]) + if isempty(inside) + @warn "gal_match_rpec: inner layer at ψ=$(round(sings[i].psifac, digits=5)) lies outside the " * + "resonant/extension cells (ψ ∈ [$lo, $hi]); raise gal_dx1/gal_dx2 to overlap the layer" surface = i + continue + end + length(inside) < length(inner_psi[i]) && @info "gal_match_rpec: surface $i inner region clipped to the " * + "cut window ($(length(inside)) of $(length(inner_psi[i])) points)" + inner_psi[i] = inner_psi[i][inside] + inner_xi[i] = inner_xi[i][inside, :] + inner_b[i] = inner_b[i][inside, :] + + # cut outer background for this surface's resonant harmonic, per coil drive + cutmn = Matrix{ComplexF64}(undef, length(psi_keep), mcoil) + for j in 1:mcoil + @views cutmn[:, j] .= sols_cut[ires, keep, 2msing+j] + for isol in 1:2msing + @views cutmn[:, j] .+= cout[isol, j] .* sols_cut[ires, keep, isol] + end + end + itp = cubic_interp(psi_keep, Series(cutmn); extrap=ExtendExtrap()) + buf = Vector{ComplexF64}(undef, mcoil) + hint = Ref(1) + for (ip, psi_p) in enumerate(inner_psi[i]) + itp(buf, psi_p; hint=hint) + singfac = m_res - nn * equil.profiles.q_spline(psi_p) + @views inner_xi[i][ip, :] .+= buf + @views inner_b[i][ip, :] .+= (2π * im * chi1_c * singfac) .* buf + end + end + end + + return GalMatchResult(cout, cin, xi, xi_deriv, deltar, bpen, inner_psi, inner_xi, inner_b, inner_params, rpec_eig, residual) end """ diff --git a/src/ForceFreeStates/Galerkin/GalerkinSolution.jl b/src/ForceFreeStates/Galerkin/GalerkinSolution.jl index 596df080..6ad3a5a1 100644 --- a/src/ForceFreeStates/Galerkin/GalerkinSolution.jl +++ b/src/ForceFreeStates/Galerkin/GalerkinSolution.jl @@ -9,7 +9,10 @@ # derivative (sing_get_dua_gal), never by spline-differentiating ξ at the packed edge. # # Fortran module flags restore_uh/us/ul default .TRUE. (gal.f); the full reconstruction is always -# performed. The cut path (gal.f, used only for the matching diagnostic) is NOT ported here. +# performed. The `cut=true` path swaps the resonant Frobenius series for its leading-order term +# (`sing_get_ua_gal_cut`); the resulting cut solution supplies the regular background that the +# resistive inner-layer solution is added to when forming the composite solution at a rational +# surface (see `gal_match_rpec`). # Sampling points per cell, matching Fortran interp_np_res / interp_np (gal.f): coarse in regular cells, # dense in resonant/extension cells to resolve the near-singular asymptotic series. @@ -37,16 +40,24 @@ const GAL_INTERP_NP_RES = 60 end """ - gal_get_solution(ws, asymps, sings, intr, x, iintvl, icell, isol; want_d=true) + gal_get_solution(ws, asymps, sings, intr, x, iintvl, icell, isol; want_d=true, cut_delta=nothing) -> (sol, dsol, iintvl, icell) Reconstruct the displacement `sol` (length `mpert`) — and, when `want_d`, its analytic radial derivative `dsol = d(sol)/dψ` — at flux `x` for solution column `isol`, from the solved Galerkin -coefficients `ws.sol`. Port of `gal_get_solution` (gal.f, non-cut path). `iintvl`/`icell` +coefficients `ws.sol`. Port of `gal_get_solution` (gal.f). `iintvl`/`icell` are a monotone cell cursor advanced and returned for the next call. `dsol` is `nothing` if `!want_d`. + +Passing `cut_delta` (the `(nsol, 2·msing)` small-solution coefficient matrix, i.e. `GalerkinResult.delta`) +selects the *cut* solution: the leading-order resonant content is subtracted from the reconstruction, +leaving the smooth background that the resistive inner-layer solution is added to when forming the +composite solution at a rational surface. `want_d` is forced off in that case. """ function gal_get_solution(ws::GalWorkspace, asymps::Vector{GalSingAsymp}, sings::Vector{SingType}, - intr::ForceFreeStatesInternal, x::Float64, iintvl::Int, icell::Int, isol::Int; want_d::Bool=true) + intr::ForceFreeStatesInternal, x::Float64, iintvl::Int, icell::Int, isol::Int; want_d::Bool=true, + cut_delta::Union{Nothing,AbstractMatrix{ComplexF64}}=nothing) + + cut_delta === nothing || (want_d = false) msing = length(ws.intvl) - 1 mpert = intr.numpert_total @@ -107,8 +118,9 @@ function gal_get_solution(ws::GalWorkspace, asymps::Vector{GalSingAsymp}, sings: end # large solution (coeff 1) only for the column driving this surface/side (gal.f) - if (isol == 2jsing - 1 && cell.extra == GAL_SIDE_LEFT) || - (isol == 2jsing && cell.extra == GAL_SIDE_RIGHT) + driving = (isol == 2jsing - 1 && cell.extra == GAL_SIDE_LEFT) || + (isol == 2jsing && cell.extra == GAL_SIDE_RIGHT) + if driving if cell.etype == GCT_RES || cell.etype == GCT_EXT || cell.etype == GCT_EXT1 @views sol .+= ua[:, 1, 1] want_d && (@views dsol .+= dua[:, 1, 1]) @@ -118,19 +130,34 @@ function gal_get_solution(ws::GalWorkspace, asymps::Vector{GalSingAsymp}, sings: end end + # Cut solution (gal.f cut_flag block): remove the leading-order resonant content, evaluated at x + # for every resonant cell type. The small coefficient comes from the Δ′ matrix rather than + # `cell.emap` because ext1/ext2 cells carry no emap. + if cut_delta !== nothing + jsol = cell.extra == GAL_SIDE_LEFT ? 2jsing - 1 : 2jsing + ua_cut = sing_get_ua_gal_cut(asymp, x - psi_s) + @views sol .-= cut_delta[isol, jsol] .* ua_cut[:, 2, 1] + driving && (@views sol .-= ua_cut[:, 1, 1]) + end + return sol, dsol, iintvl, icell end """ - gal_output_solution(ws, asymps, sings, intr, profiles, psihigh) -> GalerkinSolution + gal_output_solution(ws, asymps, sings, intr, profiles, psihigh; delta=nothing) -> GalerkinSolution Build the gal-native packed radial grid (inner→edge; `GAL_INTERP_NP_RES` points per resonant/extension cell, `GAL_INTERP_NP` elsewhere, plus the edge point ψ=psihigh) and evaluate ξ AND the analytic ξ′ over it for every solution column. Port of `gal_output_solution` (gal.f); the binary/ASCII writes and the `b_flag` ξ→b^ψ conversion (off by default) are not ported — we keep ξ itself. + +When `delta` (the `(nsol, 2·msing)` small-solution coefficient matrix) is supplied, the cut solution +`xi_cut` is evaluated on the same grid; it is the background of the composite inner-region solution +built by `gal_match_rpec`. """ function gal_output_solution(ws::GalWorkspace, asymps::Vector{GalSingAsymp}, sings::Vector{SingType}, - intr::ForceFreeStatesInternal, profiles, psihigh::Float64) + intr::ForceFreeStatesInternal, profiles, psihigh::Float64; + delta::Union{Nothing,AbstractMatrix{ComplexF64}}=nothing) mpert = intr.numpert_total msing = length(ws.intvl) - 1 @@ -160,6 +187,7 @@ function gal_output_solution(ws::GalWorkspace, asymps::Vector{GalSingAsymp}, sin ngrid = length(psi) xi = zeros(ComplexF64, mpert, ngrid, ws.nsol) xi_deriv = zeros(ComplexF64, mpert, ngrid, ws.nsol) + xi_cut = delta === nothing ? zeros(ComplexF64, 0, 0, 0) : zeros(ComplexF64, mpert, ngrid, ws.nsol) for isol in 1:ws.nsol iintvl = 0 icell = 1 @@ -169,7 +197,31 @@ function gal_output_solution(ws::GalWorkspace, asymps::Vector{GalSingAsymp}, sin @views xi[:, ip, isol] .= sol @views xi_deriv[:, ip, isol] .= dsol end + delta === nothing && continue + iintvl = 0 + icell = 1 + for ip in 1:ngrid + issing[ip] && continue + solc, _, iintvl, icell = gal_get_solution(ws, asymps, sings, intr, psi[ip], iintvl, icell, isol; + cut_delta=delta) + @views xi_cut[:, ip, isol] .= solc + end + end + + # ψ span of the resonant + extension cells flanking each surface: outside it the cut subtracts + # nothing, so the composite inner-region solution is only defined within these bounds + # (Fortran `outs%xext`). + cut_range = zeros(Float64, 0, 0) + if delta !== nothing + cut_range = hcat(fill(Inf, msing), fill(-Inf, msing)) + for iintvl in 0:msing, icell in 1:ws.nx + cell = ws.intvl[iintvl+1].cells[icell] + cell.etype == GCT_NONE && continue + js = _cell_jsing(cell, iintvl) + cut_range[js, 1] = min(cut_range[js, 1], cell.x[1]) + cut_range[js, 2] = max(cut_range[js, 2], cell.x[2]) + end end - return GalerkinSolution(psi, q, issing, xi, xi_deriv) + return GalerkinSolution(psi, q, issing, xi, xi_deriv, xi_cut, cut_range) end diff --git a/src/ForceFreeStates/Galerkin/GalerkinSolve.jl b/src/ForceFreeStates/Galerkin/GalerkinSolve.jl index 0b7dd658..aa033a36 100644 --- a/src/ForceFreeStates/Galerkin/GalerkinSolve.jl +++ b/src/ForceFreeStates/Galerkin/GalerkinSolve.jl @@ -179,7 +179,8 @@ function galerkin_solve(ctrl::ForceFreeStatesControl, equil, ffit::FourFitVars, # Reconstruct ξ(ψ) AND analytic ξ′(ψ) on the gal-native grid (gal_output_solution). ctrl.verbose && @info "Reconstructing outer-region ξ and analytic ξ′ on the gal grid" - solution = gal_output_solution(ws, asymps, sings, intr, equil.profiles, psihigh) + solution = gal_output_solution(ws, asymps, sings, intr, equil.profiles, psihigh; + delta=(ctrl.gal_match_flag ? delta : nothing)) sing_psi = [s.psifac for s in sings] sing_q = [s.q for s in sings] @@ -256,6 +257,8 @@ function write_galerkin!(out_h5, result::GalerkinResult) out_h5["galerkin/solution/issing"] = collect(sol.issing) out_h5["galerkin/solution/xi"] = sol.xi out_h5["galerkin/solution/xi_deriv"] = sol.xi_deriv + isempty(sol.xi_cut) || (out_h5["galerkin/solution/xi_cut"] = sol.xi_cut) + isempty(sol.cut_range) || (out_h5["galerkin/solution/cut_range"] = sol.cut_range) end if result.match !== nothing m = result.match @@ -270,8 +273,14 @@ function write_galerkin!(out_h5, result::GalerkinResult) for i in eachindex(m.inner_psi) out_h5["galerkin/match/inner/psi_$i"] = m.inner_psi[i] out_h5["galerkin/match/inner/xi_$i"] = m.inner_xi[i] + out_h5["galerkin/match/inner/b_$i"] = m.inner_b[i] end out_h5["galerkin/match/residual"] = m.residual + if !isempty(m.inner_params) + for f in (:E, :F, :G, :H, :K, :M, :taua, :taur, :v1) + out_h5["galerkin/match/inner_params/$(f)"] = [getfield(pp, f) for pp in m.inner_params] + end + end end return nothing end diff --git a/src/ForceFreeStates/Galerkin/GalerkinStructs.jl b/src/ForceFreeStates/Galerkin/GalerkinStructs.jl index 2fcd9173..cbfd59d1 100644 --- a/src/ForceFreeStates/Galerkin/GalerkinStructs.jl +++ b/src/ForceFreeStates/Galerkin/GalerkinStructs.jl @@ -126,6 +126,14 @@ PerturbedEquilibrium). Populated by `gal_output_solution` (GalerkinSolution.jl). - `psi::Vector{Float64}`, `q::Vector{Float64}` — radial grid (inner→edge) and its safety factor. - `issing::Vector{Bool}` — grid points sitting on a singular surface (skipped; left zero). - `xi::Array{ComplexF64,3}`, `xi_deriv::Array{ComplexF64,3}` — `(mpert, ngrid, nsol)` ξ and dξ/dψ. + - `xi_cut::Array{ComplexF64,3}` — `(mpert, ngrid, nsol)` cut solution: ξ with the leading-order + resonant content removed. This is the smooth background the resistive inner-layer solution is + added to when forming the composite solution at a rational surface, so the inner and outer + solutions overlap in the matching region. Sized `(0,0,0)` when the cut path was not requested. + - `cut_range::Matrix{Float64}` — `(msing, 2)` ψ bounds of the resonant + extension cells flanking + each rational surface. Outside these bounds no resonant content is subtracted, so `xi_cut` + equals `xi` there and the composite solution is undefined; the inner-region solution is only + meaningful inside this window. Sized `(0,0)` when the cut path was not requested. """ struct GalerkinSolution psi::Vector{Float64} @@ -133,6 +141,8 @@ struct GalerkinSolution issing::Vector{Bool} xi::Array{ComplexF64,3} xi_deriv::Array{ComplexF64,3} + xi_cut::Array{ComplexF64,3} + cut_range::Matrix{Float64} end """ @@ -169,6 +179,8 @@ struct GalMatchResult bpen::Matrix{ComplexF64} inner_psi::Vector{Vector{Float64}} inner_xi::Vector{Matrix{ComplexF64}} + inner_b::Vector{Matrix{ComplexF64}} # per surface, matched inner-layer b^ψ(ψ) on inner_psi + inner_params::Vector{InnerLayer.GGJParameters} # per-surface layer coefficients from resist_eval rpec_eig::Vector{ComplexF64} residual::Float64 end diff --git a/src/ForceFreeStates/Riccati.jl b/src/ForceFreeStates/Riccati.jl index 715c6421..9b06ef07 100644 --- a/src/ForceFreeStates/Riccati.jl +++ b/src/ForceFreeStates/Riccati.jl @@ -142,8 +142,10 @@ function assemble_fm_matrix(propagators::Vector{ChunkPropagator}, idx_range; isempty(idx_range) && return Phi for i in idx_range p = propagators[i] + #! format: off Phi_i = [p.block_upper_ic[:,:,1] p.block_lower_ic[:,:,1]; p.block_upper_ic[:,:,2] p.block_lower_ic[:,:,2]] + #! format: on Phi = Phi_i * Phi if condition condition_propagator!(Phi, N) @@ -334,6 +336,12 @@ function compute_delta_prime_matrix!( @info "Δ' BVP: nMat=$nMat, rank(M)=$(rank(M)), cond(M)=$(@sprintf("%.2e", cond(M)))" end + # rpec coil-response block: needs the S-axis row layout that `_solve_bvp_edge_coil` assumes + # for its edge rows, and a vacuum edge in the assembled matrix (col_edge in junc_rows). + if use_S_axis && wv !== nothing + intr.delta_coil = _solve_bvp_edge_coil(M, col_edge, msing, N, ipert_all) + end + intr.delta_prime_matrix = _solve_bvp_and_combine_pest3( M, msing, N, nMat, use_S_axis, ipert_all, col_edge, ctrl, debug) end @@ -628,6 +636,30 @@ function _assemble_bvp_S_axis(uShootR::Vector{Matrix{ComplexF64}}, return M, nMat, col_edge end +# Coil-response block for the Eq. (37) edge [Glasser-Kolemen 2018 PoP 25, 032501]: impose the rpec edge +# boundary condition — identity edge plus a unit source per poloidal mode, matching RDCON's +# `gal_set_boundary` rpec branch — then read the small-solution (+N slot) coefficient at every surface. +# The BC is the same for every edge mode, so the matrix is factorized once and all N modes are solved +# together as columns of one right-hand side. Returns delta_coil (2·msing × N), rows = surface side. +function _solve_bvp_edge_coil(M::Matrix{ComplexF64}, col_edge, msing::Int, N::Int, ipert_all::Vector{Int}) + nMat = size(M, 1) + bot = (nMat-2msing-N+1):(nMat-2msing) # Eq. (38) bottom rows, carrying the W_V block + Mc = copy(M) + Mc[bot, :] .= 0 + Mc[bot, col_edge] .= I(N) # Dirichlet edge: the edge coefficients equal the source + B = zeros(ComplexF64, nMat, N) + B[bot, :] .= I(N) # unit drive, one column per edge poloidal mode + X = lu(Mc) \ B + delta_coil = zeros(ComplexF64, 2msing, N) + for j in 1:msing + row_left = _col_left(j, N)[ipert_all[j]+N] + row_right = _col_right(j, N)[ipert_all[j]+N] + @views delta_coil[2j-1, :] .= X[row_left, :] + @views delta_coil[2j, :] .= X[row_right, :] + end + return delta_coil +end + # Fallback BVP assembly with FM-based axis BC (used when no Riccati S matrices are available). # Uses the conditioned axis propagator Phi_R[1][:,N+1:2N] in place of S-axis matching. function _assemble_bvp_FM_axis(Phi_L_mats::Vector{Matrix{ComplexF64}}, @@ -1497,8 +1529,10 @@ Glasser-Kolemen (2018) Phys. Plasmas 25, 032501 Eq. 33. function apply_propagator_inverse!(odet::OdeState, prop::ChunkPropagator) N = size(odet.u, 1) # Assemble 2N×2N backward FM Φ_bwd - Φ = [prop.block_upper_ic[:,:,1] prop.block_lower_ic[:,:,1]; - prop.block_upper_ic[:,:,2] prop.block_lower_ic[:,:,2]] + #! format: off + Φ = [prop.block_upper_ic[:,:,1] prop.block_lower_ic[:,:,1]; + prop.block_upper_ic[:,:,2] prop.block_lower_ic[:,:,2]] + #! format: on # Φ_bwd maps state at psi_end → psi_start (well-conditioned). # We want Φ_fwd = Φ_bwd⁻¹ to advance state from psi_start → psi_end. # Solving Φ_bwd · x = [U₁_old; U₂_old] gives x = Φ_bwd⁻¹ · [U₁_old; U₂_old]. diff --git a/src/ForceFreeStates/Sing.jl b/src/ForceFreeStates/Sing.jl index 211e6f30..61bd0916 100644 --- a/src/ForceFreeStates/Sing.jl +++ b/src/ForceFreeStates/Sing.jl @@ -944,6 +944,48 @@ Allocating wrapper for [`sing_get_ua_res!`](@ref); returns the big/small columns sing_get_ua_res(sing_asymp::SingAsymptotics, dpsi::Float64) = sing_get_ua_res!(Array{ComplexF64,3}(undef, size(sing_asymp.vmat, 1), 2, 2), sing_asymp, dpsi) +""" + sing_get_ua_res_cut!(out, sing_asymp, dpsi) -> out + +Leading-order ("cut") form of [`sing_get_ua_res!`](@ref): keeps only the `t = 0` term of the +Frobenius series instead of summing all `2·sing_order` terms, then applies the same unshear. + +This is the resonant basis used to build the *cut* outer solution, whose role is to supply the +regular background that the resistive inner-layer solution is added to when forming the composite +solution near a rational surface. Because the full and cut series share the same leading term, the +difference between them is exactly the higher-order content that the layer solution replaces. +""" +function sing_get_ua_res_cut!(out::AbstractArray{ComplexF64,3}, sing_asymp::SingAsymptotics, dpsi::Float64) + vmat = sing_asymp.vmat + N = size(vmat, 1) + sqrtfac = sqrt(dpsi) + ρ = sing_asymp.r1[1] + cbig = sing_asymp.r2[1] + csml = sing_asymp.r2[2] + + @inbounds for k in 1:2, ii in 1:N + out[ii, 1, k] = vmat[ii, cbig, k, 1] + out[ii, 2, k] = vmat[ii, csml, k, 1] + end + + pfac = dpsi^sing_asymp.alpha[1] + @inbounds for k in 1:2, ii in 1:N + out[ii, 1, k] /= pfac + out[ii, 2, k] *= pfac + end + @inbounds out[ρ, 1, 1] /= sqrtfac + @inbounds out[ρ, 2, 1] /= sqrtfac + @inbounds out[ρ, 1, 2] *= sqrtfac + @inbounds out[ρ, 2, 2] *= sqrtfac + return out +end + +""" +Allocating wrapper for [`sing_get_ua_res_cut!`](@ref); returns the big/small columns as `(N, 2, 2)`. +""" +sing_get_ua_res_cut(sing_asymp::SingAsymptotics, dpsi::Float64) = + sing_get_ua_res_cut!(Array{ComplexF64,3}(undef, size(sing_asymp.vmat, 1), 2, 2), sing_asymp, dpsi) + """ Allocating wrapper for [`sing_get_dua_res!`](@ref); returns the big/small columns as `(N, 2, 2)`. """ diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 2e10a1bd..efd5e697 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -167,7 +167,6 @@ function main_from_inputs( @info "\n Equilibrium\n$_SECTION" equil_start = time() - # Build data structures from inputs intr = ForceFreeStatesInternal(; dir_path=path) ffs_table = inputs["ForceFreeStates"] @@ -534,6 +533,10 @@ function main_from_inputs( if ctrl.gal_flag && ctrl.gal_match_flag && gal_data !== nothing && gal_data.match !== nothing @info "PerturbedEquilibrium: using the RPEC-matched gal solution" pe_odet = gal_matched_odestate(gal_data, ffit, intr) + pe_intr.odet_from_gal = true + pe_intr.inner_bpen = gal_data.match.bpen + else + pe_intr.inner_bpen = zeros(ComplexF64, intr.msing, intr.numpert_total) end # Reuse the forcing modes loaded at snapshot time (or injected by @@ -785,6 +788,14 @@ function write_outputs_to_HDF5( out_h5["singular/delta_prime_matrix"] = intr.delta_prime_matrix end + # Edge coil-response matrix, stored (numpert_total × 2msing) = (edge mode, surface-side) to match + # the galerkin/delta_coil layout so H5Web heatmaps share axes (x = edge mode, y = surface-side). + # Internal intr.delta_coil stays (2msing × numpert_total); transpose only at write. + if intr.msing > 0 && !isempty(intr.delta_coil) + dc = permutedims(intr.delta_coil) + out_h5["singular/delta_coil"] = dc + end + # Write kinetic singular surface data (det(F̄) near-zeros) and the cond(F̄) scan # used to find them. Populated only when kinetic crossings were searched for. out_h5["singular/kinetic/kmsing"] = intr.kmsing diff --git a/src/InnerLayer/GGJ/GGJ.jl b/src/InnerLayer/GGJ/GGJ.jl index e7f0488a..3a69ce79 100644 --- a/src/InnerLayer/GGJ/GGJ.jl +++ b/src/InnerLayer/GGJ/GGJ.jl @@ -1,15 +1,18 @@ # GGJ.jl # -# Glasser–Greene–Johnson resistive inner-layer model. Provides two +# Glasser–Greene–Johnson resistive inner-layer model. Provides three # interchangeable solvers selected via the `solver` type-parameter of # `GGJModel`: # -# - `:shooting` – stable backward shoot from X_max → 0 (Phase 3) -# - `:galerkin` – Hermite-cubic finite element method (Phase 4) +# - `:shooting` – stable backward shoot from X_max → 0; +# numerically stable only for |Q| ≪ 1, should not be used +# - `:galerkin` – Hermite-cubic finite element method; +# real-axis method, degrades for |Q| ≳ 1 off the real axis +# - `:ray` – rotated-contour spectral-element collocation (Ray.jl); +# robust to |Q| ~ 500 on/near the imaginary axis # -# Both solvers share the same `inps` Wasow asymptotic-basis kernel -# (`InnerAsymptotics.jl`) for the large-x boundary condition. They return -# the parity-projected matching data `(Δ_odd, Δ_even)` of GWP2016 Eqs. (34)–(35). +# They return the parity-projected matching data +# `(Δ_odd, Δ_even)` of GWP2016 Eqs. (34)–(35) in the same (deltac) convention. # # Equation references throughout this module use two source papers: # @@ -31,32 +34,39 @@ module GGJ using LinearAlgebra using StaticArrays +using SparseArrays +using Random +using Printf +using DoubleFloats: Double64 -import ..InnerLayerModel, ..solve_inner +import ..InnerLayerModel, ..solve_inner, ..solve_inner_profile """ GGJModel{S} <: InnerLayerModel -Glasser–Greene–Johnson resistive inner-layer model. The type parameter `S` -selects the solver: `:galerkin` (default) for the Hermite-cubic finite element -solver and `:shooting` for the backward stable-shoot solver. Both -implementations consume the same `inps` asymptotic-basis kernel and return -the parity-projected matching data. +Glasser–Greene–Johnson resistive inner-layer model. `S` selects the solver +backend: `:ray` (default; robust at large |Q| on/near the imaginary axis), +`:galerkin` (real-axis Hermite FEM; degrades for |Q| ≳ 1), or `:shooting` +(|Q| ≪ 1 only). The backends take different numerical-knob keywords. """ struct GGJModel{S} <: InnerLayerModel end -GGJModel(; solver::Symbol=:galerkin) = GGJModel{solver}() +GGJModel(; solver::Symbol=:ray) = GGJModel{solver}() include("GGJParameters.jl") include("InnerAsymptotics.jl") +include("RayAsymptotics.jl") include("Reference.jl") include("Shooting.jl") include("Galerkin.jl") +include("Ray.jl") export GGJModel, GGJParameters export mercier_di, mercier_dr, inner_Q, rescale_delta export build_asymptotics, evaluate_asymptotics, pick_xmax export InnerAsymptoticsCache export glasser_wang_2020_eq55 +export solve_ray, RaySolveResult, pick_smax, physical_ua_dua +export delta_convergence, solution_profile, asymptotic_profile, q4_surface_benchmark end # module GGJ diff --git a/src/InnerLayer/GGJ/GGJParameters.jl b/src/InnerLayer/GGJ/GGJParameters.jl index 732ab781..579d2e04 100644 --- a/src/InnerLayer/GGJ/GGJParameters.jl +++ b/src/InnerLayer/GGJ/GGJParameters.jl @@ -8,13 +8,10 @@ """ GGJParameters -Dimensionless parameters of the Glasser–Greene–Johnson inner-layer model -at a single rational surface, plus the local Alfvén/resistive timescales -needed to scale the matching data back to physical Δ. The equilibrium -coefficients `E, F, G, H, K, M` are the flux-surface averages defined in -GWP2016 Eq. (A8); they enter the inner-region equations (Eq. 11). - -Fields are the same as the Fortran `resist_type`: +Glasser–Greene–Johnson inner-layer parameters at one rational surface: the +flux-surface-averaged equilibrium coefficients of GWP2016 Eq. (A8) plus the +local timescales that scale the matching data back to physical Δ. Same +fields as the Fortran `resist_type`: | field | meaning | |:------- |:-------------------------------------------------------------- | @@ -29,8 +26,7 @@ Fields are the same as the Fortran `resist_type`: | `v1` | Linear scale factor used in the V₁ rescaling | | `ising` | Index of the singular surface (traceability only) | -The complex growth rate `γ` is **not** stored here; it is passed as a -separate argument to `solve_inner`. +The growth rate `γ` is not stored here; it is a separate argument to `solve_inner`. """ Base.@kwdef struct GGJParameters E::Float64 @@ -48,26 +44,22 @@ end """ mercier_di(p::GGJParameters) -> Float64 -Mercier interchange index `D_I = E + F + H − 1/4` (GWP2016 Eq. A9; also the -quantity under the root in the GW2020 Eq. 49 power exponents). `D_I > 0` -indicates local ideal interchange instability. +Mercier interchange index `D_I = E + F + H − 1/4` (GWP2016 Eq. A9); `D_I > 0` means local ideal interchange instability. """ mercier_di(p::GGJParameters) = p.E + p.F + p.H - 0.25 """ mercier_dr(p::GGJParameters) -> Float64 -Resistive interchange index `D_R = E + F + H² = D_I + (H − 1/2)²` -(GWP2016 Eq. A10). `D_R > 0` indicates local resistive interchange instability. +Resistive interchange index `D_R = E + F + H²` (GWP2016 Eq. A10); `D_R > 0` means local resistive interchange instability. """ mercier_dr(p::GGJParameters) = p.E + p.F + p.H * p.H """ p1(p::GGJParameters) -> Float64 -`p₁ = √(−D_I)`, the Mercier power that sets the large-x Frobenius exponents -`r± = 3/2 ± √(−D_I)` (GW2020 Eq. 49; μ = −1/2 ± √(−D_I) in GWP2016 Eq. 26). -The Mercier-stable branch requires `D_I < 0`; this function asserts it. +`p₁ = √(−D_I)`, setting the large-x Frobenius exponents `r± = 3/2 ± p₁` +(GW2020 Eq. 49). Throws unless `D_I < 0` (Mercier-stable required). """ function p1(p::GGJParameters) di = mercier_di(p) @@ -78,43 +70,36 @@ end """ sfac(p::GGJParameters) -> Float64 -Lundquist number `S = τ_R / τ_A` (GWP2016 Sec. I) that sets the inner-layer -length and time scales `X₀ ∝ S^(−1/3)`, `Q₀ ∝ S^(−1/3)/τ_A`. +Lundquist number `S = τ_R / τ_A`. """ sfac(p::GGJParameters) = p.taur / p.taua """ x0(p::GGJParameters) -> Float64 -Inner-layer displacement scale factor `X₀ = S^(−1/3)` (GWP2016 Eq. A14), -relating scaled `X` to physical `x = ψ − ψ₀` via `x = X₀ X`. +Inner-layer length scale `X₀ = S^(−1/3)` (GWP2016 Eq. A14); physical `x = X₀ X`. """ x0(p::GGJParameters) = sfac(p)^(-1.0 / 3.0) """ q0(p::GGJParameters) -> Float64 -Inner-layer growth-rate scale factor `Q₀ = X₀ / τ_A` (GWP2016 Eq. A15), -relating the scaled growth rate `Q` to the physical rate `s` via `s = Q₀ Q`. +Growth-rate scale `Q₀ = X₀ / τ_A` (GWP2016 Eq. A15); physical rate `= Q₀ Q`. """ q0(p::GGJParameters) = x0(p) / p.taua """ inner_Q(p::GGJParameters, γ::Number) -> ComplexF64 -Dimensionless scaled inner-layer growth rate `Q = γ / Q₀` (GWP2016 Eq. A15), -the eigenvalue appearing in the inner-region equations (Eq. 11) and the inps -Wasow basis. The argument `γ` may be real or complex; the result is always complex. +Scaled inner-layer growth rate `Q = γ / Q₀` (GWP2016 Eq. A15). """ inner_Q(p::GGJParameters, γ::Number) = ComplexF64(γ) / q0(p) """ rescale_delta(Δ, p::GGJParameters) -> SVector{2,ComplexF64} -Map the scaled inner-layer matching data back to physical Δ at the rational -surface by the `X₀^(2√(−D_I))` rescaling implied by the power-like matching -(GWP2016 Sec. IV; the inner solution ~ X^{r±} converts to physical x = X₀ X). -Operates element-wise on a 2-vector of `(Δ_odd, Δ_even)`. +Rescale the matching data to physical Δ by `S^(2√(−D_I)/3) · v₁^(2√(−D_I))` +(GWP2016 Sec. IV), element-wise on `(Δ_odd, Δ_even)`. """ function rescale_delta(Δ::AbstractVector, p::GGJParameters) s = sfac(p) @@ -122,3 +107,9 @@ function rescale_delta(Δ::AbstractVector, p::GGJParameters) fac = s^(2.0 * pp / 3.0) * p.v1^(2.0 * pp) return SVector{2,ComplexF64}(Δ[1] * fac, Δ[2] * fac) end + +# Profile conversions shared by the solve_inner_profile backends: δψ per unit inner +# coordinate X = v₁·δψ/X₀, and the big-branch (μ₋ = −1/2−p₁) amplitude rescale to the +# outer δψ-normalization, resc = (v₁/X₀)^(−μ₋) — the companion of rescale_delta's +# (v₁/X₀)^(2p₁) = (v₁/X₀)^(μ₊−μ₋). +_profile_conversions(p::GGJParameters) = (; dψdx=x0(p) / p.v1, rescale=(p.v1 / x0(p))^(0.5 + p1(p))) diff --git a/src/InnerLayer/GGJ/Galerkin.jl b/src/InnerLayer/GGJ/Galerkin.jl index 76b7578f..dc4ec5a8 100644 --- a/src/InnerLayer/GGJ/Galerkin.jl +++ b/src/InnerLayer/GGJ/Galerkin.jl @@ -38,12 +38,9 @@ using QuadGK: quadgk # ----------------------------------------------------------------------- """ -Convert the inps 6×2 output U_inps (the Wasow asymptotic basis U = TPQSY of -GW2020 Eq. 53) at coordinate `x` to the physical (Ψ, Ξ, Υ) and (Ψ', Ξ', Υ') -representation used by deltac/inpso. The 6-component first-order state packs -(Ψ, Ξ, Υ, Ψ', Ξ', Υ') with the GW2020 Eq. (2) scaling 𝚿 ≡ (xΨ, Ξ, Υ), hence -the /x and ·x factors below. Returns `(ua, dua)` each 3×2 complex, where columns -are the two power-like (Mercier) solutions and rows are the components (Ψ, Ξ, Υ). +Convert the inps 6×2 basis (GW2020 Eq. 53) at `x` to physical `(Ψ, Ξ, Υ)` +values and derivatives, each 3×2 (columns = the two power-like solutions). +The /x and ·x factors undo the GW2020 Eq. (2) scaling 𝚿 ≡ (xΨ, Ξ, Υ). """ function _physical_ua_dua(cache::InnerAsymptoticsCache, x::Real) U, dU = evaluate_asymptotics(cache, x; derivative=true, apply_T=true) @@ -61,30 +58,18 @@ function _physical_ua_dua(cache::InnerAsymptoticsCache, x::Real) end """ -Build the (I, U, V) coefficient matrices of the second-order system -`I·u'' − V·u' − U·u = 0` for u = (Ψ, Ξ, Υ) at coordinate `x`. Port of -inpso_get_uv. All matrices are 3×3 complex. - -These are the matrices A, B, C of GWP2016 Eq. (12), `A Ψ'' + B Ψ' + C Ψ = 0`, -with A given in Eq. (14), B in Eq. (14), and C in Eq. (15). The code's weak-form -layout flips the off-diagonal signs: (I, V, U) = (A, −B, −C). - -Each matrix row is one GGJ inner-region equation (GWP2016 Eq. 11 ≡ GW2020 Eq. 1); -reading `I u'' − V u' − U u` across a row reproduces the corresponding equation: - -row 1 (Ψ): Ψ_xx − H Υ_x − Q(Ψ − x Ξ) = 0 -⇒ I=(1,0,0) V=(0,0,H) U=(Q, −Qx, 0) -row 2 (Ξ): Q² Ξ_xx − Q x² Ξ + Q x Ψ + (E+F) Υ + H Ψ_x = 0 -⇒ I=(0,Q²,0) V=(−H,0,0) U=(−Qx, Q x², −(E+F)) -row 3 (Υ): Q Υ_xx − x² Υ + x Ψ + Q²[G(Ξ−Υ) − K(E Ξ + F Υ + H Ψ_x)] = 0 -⇒ I=(0,0,Q) V=(K Q² H,0,0) U=(−x, −Q²(G−KE), x²+Q²(G+KF)) +Coefficient matrices of the second-order system `I·u'' − V·u' − U·u = 0` for +u = (Ψ, Ξ, Υ) at `x`; port of inpso_get_uv, 3×3 complex. These are A, B, C of +GWP2016 Eqs. (12)–(15) with `(I, V, U) = (A, −B, −C)` and rows scaled by +(1, Q², Q). The inline comments below give the row-by-row map to the inner +equations. """ function _physical_uv(params::GGJParameters, Q::ComplexF64, x::Real) - e = ComplexF64(params.E); + e = ComplexF64(params.E) f = ComplexF64(params.F) - h = ComplexF64(params.H); + h = ComplexF64(params.H) g = ComplexF64(params.G) - k = ComplexF64(params.K); + k = ComplexF64(params.K) q = Q q2 = q * q x2 = x * x @@ -98,10 +83,11 @@ function _physical_uv(params::GGJParameters, Q::ComplexF64, x::Real) # row 1 (Ψ): (Q, −Q x, 0) # row 2 (Ξ): (−Q x, Q x², −(E+F)) # row 3 (Υ): (−x, −Q²(G−KE), x²+Q²(G+KF)) + #! format: off U = @SMatrix ComplexF64[ - q (-q * x) 0 - (-x / q) * q2 (x2 / q) * q2 (-(e + f) / q2) * q2 - (-x / q) * q (-(g - k * e) * q) * q (x2 / q + (g + k * f) * q) * q + q (-q * x) 0 + (-x / q) * q2 (x2 / q) * q2 (-(e + f) / q2) * q2 + (-x / q) * q (-(g - k * e) * q) * q (x2 / q + (g + k * f) * q) * q ] # V = −B of GWP2016 Eq. (14): coefficients of −u' in each equation (same row scaling as U): @@ -109,10 +95,11 @@ function _physical_uv(params::GGJParameters, Q::ComplexF64, x::Real) # row 2 (Ξ): (−H, 0, 0) ⇒ +H Ψ_x # row 3 (Υ): (K Q² H, 0, 0) ⇒ −K Q² H Ψ_x V = @SMatrix ComplexF64[ - 0 0 h - (-h / q2) * q2 0 0 - (h * k * q) * q 0 0 + 0 0 h + (-h / q2) * q2 0 0 + (h * k * q) * q 0 0 ] + #! format: on return Imat, U, V end @@ -127,7 +114,7 @@ function _hermite(x::Real, x0::Real, x1::Real) dx = x1 - x0 t0 = (x - x0) / dx t1 = 1 - t0 - t02 = t0 * t0; + t02 = t0 * t0 t12 = t1 * t1 pb = SVector{4,Float64}( t12 * (1 + 2t0), @@ -226,7 +213,7 @@ function _xmax_3level(params::GGJParameters, Q::ComplexF64; if !any(set) break end - x_prev = x; + x_prev = x delta_prev = dmax x *= dxfac end @@ -310,9 +297,9 @@ function _build_grid_and_workspace(nx::Int, xmax::Float64, dx1::Float64, dx2::Fl x_nodes[nx-1] = xmax - (dx1 + dx2) ixmax = nx - 2 # packed region is 0..ixmax - x0 = x_nodes[1]; + x0 = x_nodes[1] x1_packed = x_nodes[ixmax+1] - xm = (x1_packed + x0) / 2; + xm = (x1_packed + x0) / 2 dxp = (x1_packed - x0) / 2 mx = ixmax ÷ 2 packed = xm .+ dxp .* _pack(mx, pfac, side) @@ -327,7 +314,7 @@ function _build_grid_and_workspace(nx::Int, xmax::Float64, dx1::Float64, dx2::Fl imap = 1 for ix in 1:nx et = etypes[ix] - xl = x_nodes[ix]; + xl = x_nodes[ix] xr = x_nodes[ix+1] cell_np = if et == CT_NONE || et == CT_EXT1 || et == CT_EXT2 @@ -552,7 +539,7 @@ function _resonant_integral(cell::GalerkinCell, params::GGJParameters, function integrand_11(x) ua, dua = _physical_ua_dua(cache, x) Imat, Umat, Vmat = _physical_uv(params, Q, x) - ua1 = ua[:, 1]; + ua1 = ua[:, 1] dua1 = dua[:, 1] return transpose(dua1) * Imat * dua1 + transpose(ua1) * Vmat * dua1 + transpose(ua1) * Umat * ua1 end @@ -561,8 +548,8 @@ function _resonant_integral(cell::GalerkinCell, params::GGJParameters, ua, dua = _physical_ua_dua(cache, x) Imat, Umat, Vmat = _physical_uv(params, Q, x) dua1 = dua[:, 1] - ua1 = ua[:, 1]; - ua2 = ua[:, 2]; + ua1 = ua[:, 1] + ua2 = ua[:, 2] dua2 = dua[:, 2] return transpose(dua1) * Imat * dua2 + transpose(ua1) * Vmat * dua2 + transpose(ua1) * Umat * ua2 end @@ -581,7 +568,7 @@ function _assemble_and_solve!(ws::GalerkinWorkspace, params::GGJParameters, Q::ComplexF64, cache::InnerAsymptoticsCache; nq::Int=4, tol_res::Float64=1e-5) - mpert = 3; + mpert = 3 np = 3 quad_nodes, quad_weights = gausslobatto(nq + 1) offset = ws.kl + ws.kl + 1 # kl + ku + 1 since ku = kl @@ -603,14 +590,14 @@ function _assemble_and_solve!(ws::GalerkinWorkspace, for ip in 0:np_eff, ipert in 1:mpert i = cell.map[ipert, ip+1] if i > ws.ndim - ; - continue; + + continue end for jp in 0:np_eff, jpert in 1:mpert j = cell.map[jpert, jp+1] if j > ws.ndim - ; - continue; + + continue end ws.mat[offset+i-j, j, 1] += cell_mat[ipert, jpert, ip+1, jp+1] end @@ -641,8 +628,8 @@ function _assemble_and_solve!(ws::GalerkinWorkspace, continue end if i > ws.ndim - ; - continue; + + continue end for jp in 0:npp, jpert in 1:mpert j = jp < size(cell.map, 2) ? cell.map[jpert, jp+1] : cell.emap[1] @@ -650,8 +637,8 @@ function _assemble_and_solve!(ws::GalerkinWorkspace, continue end if j > ws.ndim - ; - continue; + + continue end ws.mat[offset+i-j, j, 1] += cell_mat_ext[ipert, jpert, ip+1, jp+1] end @@ -718,8 +705,8 @@ function _assemble_and_solve!(ws::GalerkinWorkspace, for ipert in 1:mpert i = cell1.map[ipert, 1] # ip=0 DOFs if i > ws.ndim - ; - continue; + + continue end for jj in max(1, i-ws.kl):min(ws.ndim, i+ws.kl) ws.mat[offset+i-jj, jj, isol] = 0 @@ -734,20 +721,20 @@ function _assemble_and_solve!(ws::GalerkinWorkspace, # → row=Ξ(ip=0), col=Ξ(ip=1): A[map[2,1], map[2,2]] = 1 # → row=Υ(ip=0), col=Υ(ip=1): A[map[3,1], map[3,2]] = 1 if isol == 1 - i = cell1.map[1, 1]; + i = cell1.map[1, 1] j = cell1.map[1, 2] ws.mat[offset+i-j, j, isol] = 1 for ipert in 2:3 - i = cell1.map[ipert, 1]; + i = cell1.map[ipert, 1] j = cell1.map[ipert, 1] ws.mat[offset+i-j, j, isol] = 1 end else - i = cell1.map[1, 1]; + i = cell1.map[1, 1] j = cell1.map[1, 1] ws.mat[offset+i-j, j, isol] = 1 for ipert in 2:3 - i = cell1.map[ipert, 1]; + i = cell1.map[ipert, 1] j = cell1.map[ipert, 2] ws.mat[offset+i-j, j, isol] = 1 end @@ -761,8 +748,8 @@ function _assemble_and_solve!(ws::GalerkinWorkspace, end # Solve for each parity using LAPACK banded LU (gbtrf! + gbtrs!) - n = ws.ndim; - kl = ws.kl; + n = ws.ndim + kl = ws.kl ku = kl for isol in 1:2 ab = copy(ws.mat[:, :, isol]) @@ -815,9 +802,12 @@ function _solution_profile(ws::GalerkinWorkspace; npc::Int=10) end k += 1 xs[k] = x - Ψ[k, 1] = vals[1, 1]; Ψ[k, 2] = vals[1, 2] - Ξ[k, 1] = vals[2, 1]; Ξ[k, 2] = vals[2, 2] - Υ[k, 1] = vals[3, 1]; Υ[k, 2] = vals[3, 2] + Ψ[k, 1] = vals[1, 1] + Ψ[k, 2] = vals[1, 2] + Ξ[k, 1] = vals[2, 1] + Ξ[k, 2] = vals[2, 2] + Υ[k, 1] = vals[3, 1] + Υ[k, 2] = vals[3, 2] end end p = sortperm(xs) @@ -845,6 +835,19 @@ function solve_inner_profile(params::GGJParameters, γ::Number; return Δ, Q, _solution_profile(ws), (; xmax=xmax_info.xmax, x0=x0(params), sfac=sfac(params)) end +""" + solve_inner_profile(::GGJModel{:galerkin}, params::GGJParameters, γ::Number; kwargs...) + -> (; Δ, x, Ψ, Ξ, dψdx, rescale) + +Hermite-FEM implementation of the [`solve_inner_profile`](@ref) interface: +real-axis solve, so `Δ` and the profiles come from the same solution. Same +numerics/kwargs as `solve_inner(GGJModel(; solver=:galerkin), ...)`. +""" +function solve_inner_profile(::GGJModel{:galerkin}, params::GGJParameters, γ::Number; kwargs...) + Δ, _, prof, _ = solve_inner_profile(params, γ; kwargs...) + return (; Δ=Δ, x=prof.x, Ψ=prof.Ψ, Ξ=prof.Ξ, _profile_conversions(params)...) +end + """ solve_inner(::GGJModel{:galerkin}, params::GGJParameters, γ::Number; kmax::Int=8, nx::Int=512, nq::Int=4, pfac::Float64=1.0, @@ -893,20 +896,15 @@ end # ISSUES AND IS NOT RELIABLE IN BROAD SCANS. IT IS LEFT HERE FOR REFERENCE AND MAY BE REWORKED LATER. """ solve_inner_converged(::GGJModel{:galerkin}, params::GGJParameters, γ::Number; - rtol=1e-2, kmax0=12, xfac0=1.5, cells_per_unit=3.0, - nx_min=1024, nx_max=8192, kmax_step=2, kmax_max=28, - xfac_growth=1.5, max_levels=6, nq=6, pfac=1.0, cutoff=8, tol_res=1e-4) + rtol=1e-2, max_levels=6, kwargs...) -> (; delta, converged, err, kmax, xfac, nx, nlevels) -Convergence-guarded GGJ inner-layer solve: only returns a Δ once it is stable under joint refinement of -the three coupled accuracy knobs — series order `kmax`, asymptotic reach `xfac`, and grid resolution `nx`. -Successively refines all three (raising `kmax` clears the high-|Q| series floor; raising `xfac` clears the -reach floor; `nx` is scaled with `xmax` to hold cells-per-unit-x ≈ `cells_per_unit`, which prevents the -grid-starvation breakup that otherwise corrupts Δ at large `xfac`/high |Q|) until the per-component -relative change of (Δ₁, Δ₂) drops below `rtol`, or `max_levels` is hit (then `converged=false`). - -The metric is per real/imag component with a significance floor (5% of |Δᵢ|), so a converged norm cannot -mask a wrong reactive part Re(Δ₁) — the failure mode under-reach produces (Im dominates |Δ₁|). +Convergence-guarded Galerkin solve: jointly refines the three coupled accuracy +knobs (series order `kmax`, asymptotic reach `xfac`, grid `nx` scaled to hold +cells-per-unit-x ≈ `cells_per_unit`) until the per-component relative change +of (Δ₁, Δ₂) drops below `rtol`, or `max_levels` is hit (`converged=false`). +The metric is per real/imag component with a 5% significance floor, so a +converged norm cannot mask a wrong small component. """ function solve_inner_converged(model::GGJModel{:galerkin}, params::GGJParameters, γ::Number; rtol::Float64=1e-2, kmax0::Int=12, xfac0::Float64=1.5, diff --git a/src/InnerLayer/GGJ/InnerAsymptotics.jl b/src/InnerLayer/GGJ/InnerAsymptotics.jl index 7b6b838c..4dfd641a 100644 --- a/src/InnerLayer/GGJ/InnerAsymptotics.jl +++ b/src/InnerLayer/GGJ/InnerAsymptotics.jl @@ -33,13 +33,9 @@ const _R2 = SVector(3, 4, 5, 6) """ InnerAsymptoticsCache -Frozen state of the `inps` Wasow asymptotic-basis construction for a single -`(GGJParameters, Q)` pair. All matrices are stored as `SMatrix`/`SVector` -so the evaluator can run allocation-free on a hot path. - -Index convention: `P[k+1]` holds the k-th-order matrix `P_k`, `B[k+1]` -holds `B_k`, etc., for k = 0, 1, …, the upper bound documented in each -field. +Frozen `inps` Wasow asymptotic-basis construction for one `(GGJParameters, Q)` +pair, stored as `SMatrix`/`SVector` for allocation-free evaluation. Index +convention: `P[k+1]` holds the k-th-order matrix `P_k`, etc. Fields: @@ -97,6 +93,7 @@ function _build_tjmat(p::GGJParameters, Q::ComplexF64) # T (Eq. 7) — Fortran constructs this with column-major RESHAPE; the # listing below is row-major Julia order matching that layout. + #! format: off T = @SMatrix ComplexF64[ 1 0 h*q q2/λ h*q -q2/λ 0 0 0 -1/λ 0 1/λ @@ -114,6 +111,7 @@ function _build_tjmat(p::GGJParameters, Q::ComplexF64) 0 0 λ/2 0 0 1/2 0 λ/2 0 0 1/2 0 ] + #! format: on # A_0, A_1, A_2 — GW2020 Eqs. (4), (5), and the A₂ matrix of Eq. (3). Build mutable then freeze. A0 = zeros(ComplexF64, 6, 6) @@ -176,9 +174,9 @@ function _lyap_solve(K::SMatrix{6,6,ComplexF64}, λ::ComplexF64) Pm = zeros(ComplexF64, 6, 6) # B is block-diagonal in the (r1, r2) split. - Bm[1, 1] = K[1, 1]; + Bm[1, 1] = K[1, 1] Bm[1, 2] = K[1, 2] - Bm[2, 1] = K[2, 1]; + Bm[2, 1] = K[2, 1] Bm[2, 2] = K[2, 2] for i in 3:6, j in 3:6 Bm[i, j] = K[i, j] @@ -314,7 +312,7 @@ function _coefs(B::Vector{SMatrix{6,6,ComplexF64,36}}, end # Lowest-order Y solution and inverse — GW2020 Eq. (48), with exponents r± from Eq. (49). - r1 = ComplexF64(R[1]); + r1 = ComplexF64(R[1]) r2 = ComplexF64(R[2]) Y0 = @SMatrix ComplexF64[ 1 1 @@ -357,13 +355,9 @@ end """ build_asymptotics(params::GGJParameters, Q::ComplexF64; kmax::Int=8) -> InnerAsymptoticsCache -Construct the `inps` Wasow asymptotic basis for the given GGJ parameters -and dimensionless growth rate `Q`. Truncates each power series at order -`kmax` (default `8`). The returned cache can be evaluated at any `x > 0` -via [`evaluate_asymptotics`](@ref) and queried for an adaptive `X_max` -via [`pick_xmax`](@ref). - -Reference: Glasser & Wang, Phys. Plasmas **27**, 012506 (2020), Eqs. 7–53. +Construct the `inps` Wasow asymptotic basis (GW2020 Eqs. 7–53), truncating +each power series at order `kmax`. Evaluate with [`evaluate_asymptotics`](@ref); +pick a cutoff with [`pick_xmax`](@ref). """ function build_asymptotics(params::GGJParameters, Q::ComplexF64; kmax::Int=8) p1v = p1(params) @@ -498,19 +492,12 @@ end # ----------------------------------------------------------------------- """ - evaluate_asymptotics(cache::InnerAsymptoticsCache, x::Real; - derivative::Bool=true, apply_T::Bool=true) - -> (U, dU) - -Evaluate the inps asymptotic basis at `x > 0`. Returns the 6×2 complex -matrix `U` whose two columns are the algebraically-decaying ("small") -asymptotic solutions of the GGJ system, and (if `derivative=true`) the -6×2 matrix `dU` of their derivatives `dU/dx`. - -If `apply_T=false`, the result is left in the J-rotated coordinate basis -(used by `inps_delta` for residual checks). The default `apply_T=true` -returns the solutions in the original 6-component first-order-system -basis used by `inpso_get_uv` and the shooting / Galerkin solvers. + evaluate_asymptotics(cache::InnerAsymptoticsCache, x::Real; derivative=true, apply_T=true) -> (U, dU) + +Evaluate the inps basis at `x > 0`: `U` is 6×2 with columns the two power-like +solutions (GW2020 Eq. 53), `dU` their x-derivatives (`nothing` unless +`derivative=true`). `apply_T=false` leaves the result in the J-rotated basis +used by the residual check. """ function evaluate_asymptotics(cache::InnerAsymptoticsCache, x::Real; derivative::Bool=true, apply_T::Bool=true) @@ -532,7 +519,7 @@ function evaluate_asymptotics(cache::InnerAsymptoticsCache, x::Real; # Splitting matrix pp (6×2): top 2×2 = I, bottom 4×2 = p21. pp_m = zeros(ComplexF64, 6, 2) - pp_m[1, 1] = 1; + pp_m[1, 1] = 1 pp_m[2, 2] = 1 @inbounds for i in 1:4, j in 1:2 pp_m[i+2, j] = p21[i, j] @@ -590,11 +577,8 @@ end """ asymptotic_residual(cache::InnerAsymptoticsCache, x::Real) -> SVector{2,Float64} -Compute the convergence measure `Δ±` of the asymptotic basis at `x` for each -of the two algebraic columns (GW2020 Eq. 54). Mirrors `inps_delta`: returns -`‖dU − x·matrix·U‖∞ / max(‖dU‖∞, ‖x·matrix·U‖∞)` per column, where -`matrix = J₀ + xfac·J₁ + xfac²·J₂` is the J-rotated coefficient matrix (the -residual of `v' = xJv`, GW2020 Eq. 6). +Convergence measure Δ± of the two series columns at `x` (GW2020 Eq. 54): +`‖dU − x·J(x)·U‖∞ / max(‖dU‖∞, ‖x·J(x)·U‖∞)` per column. Mirrors `inps_delta`. """ function asymptotic_residual(cache::InnerAsymptoticsCache, x::Real) U, dU = evaluate_asymptotics(cache, x; derivative=true, apply_T=false) @@ -616,8 +600,8 @@ function asymptotic_residual(cache::InnerAsymptoticsCache, x::Real) delta = MVector{2,Float64}(0.0, 0.0) @inbounds for j in 1:2 - n0 = 0.0; - n1 = 0.0; + n0 = 0.0 + n1 = 0.0 n2 = 0.0 for i in 1:6 n0 = max(n0, abs(matvec0[i, j])) @@ -631,19 +615,11 @@ function asymptotic_residual(cache::InnerAsymptoticsCache, x::Real) end """ - pick_xmax(params::GGJParameters, Q::ComplexF64; - eps::Float64=1e-7, kmax::Int=8, - xlogmin::Float64=-1.0, xlogmax::Float64=4.0, - dxlog::Float64=0.01) -> (Float64, InnerAsymptoticsCache) - -Sweep `x` log-uniformly upward from `10^xlogmin` and return the smallest -`x` at which `max(asymptotic_residual(cache, x)) < eps` — the cutoff `x_max` -where the GW2020 Eq. (54) convergence measure drops below tolerance -(GW2020 Sec. III, Fig. 3). Also returns the `InnerAsymptoticsCache` it built -so callers can reuse it. Mirrors `inps_xmax`. - -Throws an `ErrorException` if no `x` in the sweep range achieves the -target tolerance. + pick_xmax(params, Q; eps=1e-7, kmax=8, xlogmin=-1.0, xlogmax=4.0, dxlog=0.01) -> (x_max, cache) + +Sweep `x` log-uniformly and return the smallest `x` where the GW2020 Eq. (54) +residual drops below `eps`, plus the cache built along the way. Throws if the +tolerance is never reached in the sweep range. Mirrors `inps_xmax`. """ function pick_xmax(params::GGJParameters, Q::ComplexF64; eps::Float64=1e-7, kmax::Int=8, diff --git a/src/InnerLayer/GGJ/Ray.jl b/src/InnerLayer/GGJ/Ray.jl new file mode 100644 index 00000000..ea0e62ee --- /dev/null +++ b/src/InnerLayer/GGJ/Ray.jl @@ -0,0 +1,1115 @@ +# Ray.jl +# +# Rotated-contour spectral-element collocation backend (`GGJModel{:ray}`) +# for the GGJ inner-layer matching data — robust at large |Q| on and near +# the imaginary axis, where both `:shooting` (exponential dichotomy) and +# `:galerkin` (real-axis oscillation + on-axis pseudo-resonance) degrade. +# +# Validated against the Fortran rmatch pins, the `:galerkin` backend at +# moderate Q, and physical benchmarks to Q = 500i. Summarized method write-up: +# docs/src/inner_layer.md, "Numerical method". In one paragraph: +# +# The equations are continued analytically onto the ray x = e^{iθ}s with +# θ = arg(Q)/4 (parabolic-cylinder exponent exactly real; pseudo-resonance +# cleared). On [0, s_m] a global Chebyshev spectral-element collocation BVP +# in the plain variables v = (Ψ, Ξ, Υ, Ψ', Ξ', Υ') is solved with parity +# conditions at the (ordinary-point) origin and Δ, c₁, c₂ as bordered +# unknowns of one sparse LU per parity (rank-3 Woodbury for the second). +# The far-field condition uses the SAME inps Wasow basis as the other +# backends (RayAsymptotics.jl evaluates it at complex x; the rdcon +# normalization of Δ is retained), applied at the series radius S and +# transported S → s_m by an L-stable 2-stage Radau IIA march in the +# quotient modulo the decaying exponential pair. The damped-zone march +# arithmetic runs in Complex{Double64}: at large S the near-parallel +# power-pair geometry amplifies the structured backward error of the +# ill-conditioned implicit solves (Σ eps·z ≈ eps·ρS²/2) into Δ-mixing +# ~1e-4 at |Q| = 500 in Float64 (docs/src/inner_layer.md, "Far-field +# boundary and the inward march"). +# +# File layout: +# 1. system — plain-variable first-order ODE matrix, parity sets +# 2. mesh — Chebyshev cells, barycentric interpolation, WKB grading +# 3. boundary — decaying pair, Radau IIA quotient march +# 4. solve — assembly, bordered solve, refinement, driver, diagnostics + +# system.jl +# +# The GGJ inner-region equations (GWP2016 Eq. 11 ≡ GW2020 Eq. 1) as a plain +# first-order system in v = (Ψ, Ξ, Υ, Ψ', Ξ', Υ'), dv/dx = M(x)·v. +# All coefficients are POLYNOMIAL in x, so x = 0 is an ordinary point (the +# x⁻²/x⁻⁴ singularities of the GW2020 Eq. 2 form are artifacts of the +# (xΨ, Ψ'/x) scaling) and the system continues analytically to complex x. +# +# Ψ'' = H Υ' + Q(Ψ − xΞ) +# Ξ'' = [Q x² Ξ − Q x Ψ − (E+F) Υ − H Ψ'] / Q² +# Υ'' = [x² Υ − x Ψ − Q²(G(Ξ−Υ) − K(EΞ + FΥ + HΨ'))] / Q +# +# Row-by-row this matches the deltac `inpso_get_uv` matrices (I, V, U) after +# dividing each equation by its diagonal second-derivative coefficient +# (1, Q², Q). + +""" + ode_matrix([CT,] p::GGJParameters, Q, x) -> SMatrix{6,6,CT} + +Coefficient matrix `M(x)` of `dv/dx = M v`, `v = (Ψ, Ξ, Υ, Ψ', Ξ', Υ')`, for +real or complex `x`. `CT` selects the element type (default `ComplexF64`; +`Complex{Double64}` for the extended-precision march). +""" +ode_matrix(p::GGJParameters, Q::ComplexF64, x::Number) = + ode_matrix(ComplexF64, p, Q, x) + +function ode_matrix(::Type{CT}, p::GGJParameters, Q::Number, x::Number) where {CT<:Complex} + e = CT(p.E) + f = CT(p.F) + g = CT(p.G) + h = CT(p.H) + k = CT(p.K) + q = CT(Q) + q2 = q * q + xc = CT(x) + x2 = xc * xc + + m = zeros(MMatrix{6,6,CT}) + # rows 1–3: (Ψ, Ξ, Υ)' = (Ψ', Ξ', Υ') + m[1, 4] = 1 + m[2, 5] = 1 + m[3, 6] = 1 + # row 4: Ψ'' = Q Ψ − Q x Ξ + H Υ' + m[4, 1] = q + m[4, 2] = -q * xc + m[4, 6] = h + # row 5: Ξ'' = −(x/Q) Ψ + (x²/Q) Ξ − ((E+F)/Q²) Υ − (H/Q²) Ψ' + m[5, 1] = -xc / q + m[5, 2] = x2 / q + m[5, 3] = -(e + f) / q2 + m[5, 4] = -h / q2 + # row 6: Υ'' = −(x/Q) Ψ − Q(G−KE) Ξ + (x²/Q + Q(G+KF)) Υ + HKQ Ψ' + m[6, 1] = -xc / q + m[6, 2] = -q * (g - k * e) + m[6, 3] = x2 / q + q * (g + k * f) + m[6, 4] = h * k * q + + return SMatrix{6,6,CT}(m) +end + +""" + parity_rows(isol) -> SVector{3,Int} + +Component indices of `v = (Ψ, Ξ, Υ, Ψ', Ξ', Υ')` that vanish at s = 0 for +each parity solution, mirroring deltac_set_boundary: + + - `isol = 1` ("odd"): Ψ'(0) = 0, Ξ(0) = 0, Υ(0) = 0 → components (4, 2, 3) + - `isol = 2` ("even"): Ψ(0) = 0, Ξ'(0) = 0, Υ'(0) = 0 → components (1, 5, 6) +""" +parity_rows(isol::Int) = isol == 1 ? SVector(4, 2, 3) : SVector(1, 5, 6) +# mesh.jl +# +# Graded spectral-element mesh on the ray parameter s ∈ [0, S], plus the +# Chebyshev–Lobatto reference-cell machinery (nodes, differentiation matrix, +# barycentric interpolation). +# +# The initial grading tracks the local stiffness scales of the rotated +# system: +# - the fast Ohm's-law pair e^{±√Q x}: active near the origin, decaying +# at rate Re(√Q e^{iθ}) — resolved on a fine zone [0, s₁]; +# - the parabolic-cylinder pair e^{±ρ s²/2}, ρ = Re(e^{2iθ}/√Q): decaying +# content is above relative machine floor only for s ≲ s_d = √(2D/ρ) — +# resolved on a rate-adapted zone [s₁, s_d]; +# - beyond s_d the physical solution is a smooth power law: geometric +# (log-spaced) cells suffice out to S. +# Residual-driven bisection refinement (solve.jl) then repairs whatever the +# initial guess missed. + +""" + cheblobatto(p) -> (t, D) + +Chebyshev–Lobatto nodes on [-1, 1] in ASCENDING order and the spectral +differentiation matrix `D` for those nodes (Trefethen's `cheb`, reflected). +""" +function cheblobatto(p::Int) + p >= 1 || throw(ArgumentError("polynomial order p must be ≥ 1")) + # Standard descending Trefethen nodes x_j = cos(jπ/p), j = 0..p. + x = [cos(pi * j / p) for j in 0:p] + c = [j == 0 || j == p ? 2.0 : 1.0 for j in 0:p] .* [(-1.0)^j for j in 0:p] + D = zeros(p + 1, p + 1) + for i in 1:(p+1), j in 1:(p+1) + if i != j + D[i, j] = (c[i] / c[j]) / (x[i] - x[j]) + end + end + for i in 1:(p+1) + D[i, i] = -sum(D[i, j] for j in 1:(p+1) if j != i) + end + # Reflect to ascending order: y = -x is ascending with the SAME index + # order (x is descending), and for g(y) := f(-y) sampled on y the chain + # rule gives D_y = -D. No index permutation. + t = -x + Dt = -D + return t, Dt +end + +""" + barycentric_weights(t) -> w + +Barycentric interpolation weights for nodes `t` (any distribution). +""" +function barycentric_weights(t::AbstractVector{Float64}) + n = length(t) + w = ones(Float64, n) + for j in 1:n + for k in 1:n + k == j && continue + w[j] /= (t[j] - t[k]) + end + end + return w +end + +""" + barycentric_eval(t, w, vals, τ) -> value + +Evaluate the interpolant of `vals` (matrix: nodes × components) at `τ`. +Returns a vector over components. +""" +function barycentric_eval(t::AbstractVector{Float64}, w::AbstractVector{Float64}, + vals::AbstractMatrix{ComplexF64}, τ::Float64) + n = length(t) + for j in 1:n + if τ == t[j] + return vals[j, :] + end + end + num = zeros(ComplexF64, size(vals, 2)) + den = 0.0 + for j in 1:n + fac = w[j] / (τ - t[j]) + den += fac + @views num .+= fac .* vals[j, :] + end + return num ./ den +end + +""" + initial_breaks(params, Q, θ, S, p; cfl=0.4, drop=40.0, ratio=1.3, + hmax_frac=0.08) -> Vector{Float64} + +Build the initial cell-boundary vector `0 = s₀ < s₁ < … = S` using the +three-zone grading described in the header. `cfl` ≈ resolved-exponential +rate per node; `drop` = e-folds after which decaying exponential content is +considered machine-dead; `ratio` = geometric growth in the outer zone. +""" +function initial_breaks(params::GGJParameters, Q::ComplexF64, θ::Float64, + S::Float64, p::Int; cfl::Float64=0.4, drop::Float64=40.0, + ratio::Float64=1.3, hmax_frac::Float64=0.08) + sq = sqrt(Q) + ρ = real(cis(2θ) / sq) # parabolic pair: exponent ρ s²/2 + ρ > 0 || throw(ArgumentError("contour angle θ=$θ gives non-decaying parabolic pair (ρ=$ρ ≤ 0); use θ ≈ arg(Q)/4")) + # Fast local WKB scales: the Ohm pair |√Q| and the Υ-family + # √|Q(G+KF)| (dominant on extreme-coefficient surfaces, K ~ 10⁶ — + # dormant almost everywhere but must be resolved; see solve.jl header). + μ = abs(sq) + sqrt(abs(Q) * abs(params.G + params.K * params.F)) + μdec = max(real(sq * cis(θ)), 0.1 * μ) # its decay rate along the ray + + s1 = min(drop / μdec, 0.7 * S) # fast pair dead beyond s1 + sd = min(sqrt(2 * drop / ρ), 0.8 * S) # parabolic pair dead beyond sd + sd = max(sd, min(1.05 * s1, S)) + hmax = hmax_frac * S + + # Single grading rule: every cell resolves the full local WKB rate + # μ + ρs (fast pair + parabolic pair), whether the corresponding mode + # amplitudes are alive or dormant — unresolved dormant modes alias into + # spurious bounded discrete modes and poison the global solve. The + # outer power-law region is not discretized at all (handled by the + # boundary march), so this never gets expensive: the BVP domain ends at + # s_m ~ max(s_d, s₁) where μ dominates the rate budget. + _ = (s1, sd, ratio) # zone markers retained in signature for tuning + breaks = [0.0] + while breaks[end] < S + s = breaks[end] + h = cfl * p / (μ + ρ * s) + h = min(h, hmax, S / 8) + snew = s + h + if snew >= S - 0.25 * h + snew = S + end + push!(breaks, min(snew, S)) + end + # Guard against a degenerate final sliver. + if length(breaks) >= 3 && (breaks[end] - breaks[end-1]) < 0.2 * (breaks[end-1] - breaks[end-2]) + deleteat!(breaks, length(breaks) - 1) + end + return breaks +end +# boundary.jl +# +# Large-s boundary data at the matching point on the ray x = e^{iθ}s. +# +# The admissible (bounded-as-s→∞) solution space is spanned by +# { u_small, u_big } — the two power-like Mercier solutions, taken from +# the inps series in the standard normalization, and +# { E₁, E₂ } — the two forward-DECAYING exponential solutions. +# +# The inps series is only trusted at radius S (residual ≤ smax_tol), which +# for extreme-coefficient surfaces (K ~ 10⁶) can be S ~ 10⁴–10⁵, while the +# collocation BVP ends at s_m ~ 30–150. The gap is bridged by MARCHING the +# power-pair data W = [u_small, u_big] backward S → s_m. +# +# Marching design (the hard-won part): +# - Backward, the decaying pair GROWS at rate ρs (ρ = Re(e^{2iθ}/√Q)), up +# to ~10³/unit at large s. Any explicit integrator — including +# projected/continuous-QR variants, whose unprojected stage points +# reintroduce the transverse stiffness — is stability-limited to +# h ≲ 3/(ρs), i.e. O(ρS²) steps: hopeless at S ~ 10⁵. +# - The system is LINEAR, so we use a stiffly-accurate L-STABLE implicit +# method (2-stage Radau IIA, order 3; one 12×12 solve per step). For +# resolved slow content it is 3rd-order accurate; for the unresolvable +# backward-growing exponential directions |R(z)| < 1 at large |z| — +# the method DAMPS what it cannot resolve instead of amplifying it. +# That is normally an accuracy vice; here it is exactly the required +# behavior, because Δ is defined modulo the decaying pair. +# - In the partially-resolved band (|z| = ρs·h ≲ 10) genuine backward +# growth of rounding-injected decaying-pair content still occurs, so W +# is PURGED against a locally-extracted decaying frame every +# Φ = ∫ρs ds ≈ 8 e-folds. Purging only ever removes span{E} content — +# legitimate by the definition of Δ — so the transported coefficients +# of u_small/u_big are exact up to integrator tolerance. + +""" + decaying_pair(params, Q, θ, S; growth=30.0, seed=20260707, dp_res=20.0) -> E (6×2) + +Orthonormal basis of the forward-decaying pair at s = S on the ray: random +seeds integrated backward over [S, S+ΔS] with fully-resolved RK4, ΔS sized so +the pair separates from everything else by ≈ `growth` e-folds. +""" +function decaying_pair(params::GGJParameters, Q::ComplexF64, θ::Float64, S::Float64; + growth::Float64=30.0, seed::Int=20260707, dp_res::Float64=20.0) + sq = sqrt(Q) + ρ = real(cis(2θ) / sq) + ρ > 0 || throw(ArgumentError("θ=$θ gives ρ=$ρ ≤ 0; decaying pair undefined")) + ph = cis(θ) + # Purification window: BOTH backward-growing directions must separate + # from the power pair by ≥ `growth` e-folds over [S, S+ΔS], so the + # budget is set by the SLOWEST grower's rate margin over the algebraic + # power rates (from the frozen spectrum) — not by ρS. At moderate |Q| + # the fast-Ohm grower can be ~0.6/unit while ρS ~ 10²: budgeting on ρS + # leaves the second frame column O(1)-contaminated with power content, + # which then poisons every march purge that uses the frame. + g2 = Inf + gmax_loc = 0.0 + for λ in eigvals(Matrix(ph * ode_matrix(params, Q, ph * S))) + g = -real(λ) + g > 20.0 / S && (g2 = min(g2, g)) + gmax_loc = max(gmax_loc, abs(real(λ))) + end + gsep = isfinite(g2) ? max(g2 - 20.0 / S, 0.05 * g2) : ρ * S + ΔS = min(growth / gsep, 3.0 * S) + + # Step resolution must use the TRUE local spectrum, not the analytic + # ρs + |√Q| estimate: on extreme-coefficient surfaces the Υ-family rate + # √|Q(G+KF)| dominates at small s and the analytic estimate under-counts + # it 10×, leaving h·λ ~ 0.2 — stable but phase-inaccurate, which corrupts + # the extracted frame directions (observed as O(1) θ-drift downstream). + for λ in eigvals(Matrix(ph * ode_matrix(params, Q, ph * (S + ΔS)))) + gmax_loc = max(gmax_loc, abs(real(λ))) + end + rate = max(ρ * (S + ΔS) + abs(sq), gmax_loc) + nsteps = clamp(ceil(Int, dp_res * rate * ΔS), 200, 2_000_000) + h = -ΔS / nsteps + + rng = MersenneTwister(seed) + v = randn(rng, ComplexF64, 6, 2) + v = Matrix(qr(v).Q) + + f(s, y) = (ph * ode_matrix(params, Q, ph * s)) * y + + s = S + ΔS + for istep in 1:nsteps + k1 = f(s, v) + k2 = f(s + h / 2, v + (h / 2) * k1) + k3 = f(s + h / 2, v + (h / 2) * k2) + k4 = f(s + h, v + h * k3) + v = v + (h / 6) * (k1 + 2k2 + 2k3 + k4) + s += h + if istep % 50 == 0 + v = Matrix(qr(v).Q) + end + end + return Matrix(qr(v).Q) +end + +""" + march_boundary(params, Q, θ, S, s_m, Us, Ub; rtol=1e-9, growth=30.0, seed=20260707) + -> (Us_m, Ub_m, F) + +Transport the power-pair boundary data from the series radius `S` inward to +`s_m` along the ray, modulo the decaying exponential pair (the quotient in +which Δ is defined). 2-stage Radau IIA (L-stable) with Φ-budgeted purges +against a carried decaying frame; explicit marchers cannot do this at large S +(see file header). +""" +function march_boundary(params::GGJParameters, Q::ComplexF64, θ::Float64, + S::Float64, s_m::Float64, Us::Vector{ComplexF64}, Ub::Vector{ComplexF64}; + rtol::Float64=1e-9, growth::Float64=30.0, seed::Int=20260707, + ztrk::Float64=0.3, purge_budget::Float64=2.5, dp_res::Float64=20.0) + ph = cis(θ) + 𝔐(s) = ph * ode_matrix(params, Q, ph * s) + + W = hcat(copy(Us), copy(Ub)) + function purge!(W, s) + Floc = decaying_pair(params, Q, θ, s; growth=growth, seed=seed, dp_res=dp_res) + W .-= Floc * (Floc' * W) + return Floc + end + + I6 = Matrix{ComplexF64}(I, 6, 6) + # One Radau IIA step: A = [5/12 -1/12; 3/4 1/4], c = (1/3, 1), b = A[2,:]. + function radau_step(s, W, h) + M1 = 𝔐(s + h / 3) + M2 = 𝔐(s + h) + Ablk = [I6-(5h/12)*M1 (h/12)*M1; -(3h/4)*M2 I6-(h/4)*M2] + rhs = [M1 * W; M2 * W] + K = Ablk \ rhs + return W + h * ((3 / 4) * @view(K[1:6, :]) + (1 / 4) * @view(K[7:12, :])) + end + + # Extended-precision (Double64) Radau step for the damped-zone W march. + # At z = ρ·c·s² ≫ 1 the step matrices have cond ~ z and the LU backward + # error has a SYSTEMATIC direction set by the elimination structure — it + # does not dither away, accumulating a relative error Σ eps·z ≈ eps·ρS²/2 + # on the power columns that is INDEPENDENT of the step law (c cancels). + # The near-parallel power-pair geometry at large S turns that into + # u_small/u_big coefficient mixing β ~ eps·ρS²/2 · (‖u_big‖/‖u_small‖)(S) + # ~ 1e-4 at |Q| = 500 — saturating Δ of the near-pole parity at ~1/β and + # flooring the other at |Δ|·β — the Δ-mixing defect extended precision removes. + # Double64 arithmetic (eps ~ 5e-33) removes it outright; the resolved + # band stays Float64 (z ≤ ztrk there, solves well-conditioned). + CTd = Complex{Double64} + phd = CTd(cos(Double64(θ)), sin(Double64(θ))) + I6d = Matrix{CTd}(I, 6, 6) + 𝔐d(s) = phd * Matrix(ode_matrix(CTd, params, Q, phd * Double64(s))) + function radau_step_d(s, Wd, h) + hd = Double64(h) + M1 = 𝔐d(s + h / 3) + M2 = 𝔐d(s + h) + Ablk = [I6d-(5hd/12)*M1 (hd/12)*M1; -(3hd/4)*M2 I6d-(hd/4)*M2] + rhs = [M1 * Wd; M2 * Wd] + K = Ablk \ rhs + return Wd + hd * ((CTd(3) / 4) * K[1:6, :] + (CTd(1) / 4) * K[7:12, :]) + end + + # --- Schedule: geometric where damped, growth-capped where resolved. --- + # No adaptive controller (three generations died on estimator artifacts); + # the step law is derived from the local mode structure instead. At each + # step the eigenvalues of 𝔐(s) give the two backward-growing exponential + # rates g (Re λ < 0, excluding the algebraic power-pair rates ~ r/s): + # + # damped zone (g_min·c·s > 8 for all growers): h = -c·s. All growing + # content — including the O(smax_tol) fast-direction contamination of + # the series data — sits at z ≫ 1 where L-stable Radau DAMPS it + # (|R| ~ 6/z²). Slow-mode global error ≈ 10²c³ → c = (rtol/100)^{1/3}; + # step count ~ ln(S/s_m)/c independent of S. + # + # resolved band (some grower at z ≲ 8): genuine backward growth must be + # tracked and removed. h = -min(c·s, 0.3/g_max) keeps z ≤ 0.3 so the + # DISCRETE growing eigendirection matches the continuous one to + # O(z²) ~ 1% (an orthogonal purge against a frame misaligned by δ + # leaks δ·content — with e-fold budget Φ ≤ 2.5 the cycle gain + # e^Φ·δ < 0.2 and contamination stays at the rounding floor). The + # frame F is carried CONTINUOUSLY with the same Radau steps + QR + # (backward-attracting ⇒ self-correcting), so no re-extraction cost. + c = clamp(cbrt(rtol / 100.0), 2e-5, 2e-3) + + # --- Phase 1: damped zone, W in Double64 (see radau_step_d comment). --- + # First purge (at S) and all damped-zone steps run in extended precision; + # on resolved-band entry W drops back to Float64 for phase 2. + Wd = CTd.(W) + let Floc = decaying_pair(params, Q, θ, S; growth=growth, seed=seed, dp_res=dp_res) + Fd = CTd.(Floc) + Wd .-= Fd * (Fd' * Wd) + end + s = S + nstep = 0 + while s > s_m + 1e-12 * S + λs = eigvals(Matrix(𝔐(s))) + gmin = Inf + for λ in λs + g = -real(λ) + g > 10.0 / s && (gmin = min(gmin, g)) + end + (isfinite(gmin) && gmin * c * s < 8.0) && break # resolved-band entry + h = -min(c * s, s - s_m) + Wd = radau_step_d(s, Wd, h) + s += h + nstep += 1 + nstep > 5_000_000 && error("march_boundary: step limit exceeded (s=$s)") + end + W = ComplexF64.(Wd) + + # --- Phase 2: resolved band (Float64, carried frame + purges). --- + Φ = 0.0 + F = Matrix{ComplexF64}(undef, 6, 2) + have_F = false + while s > s_m + 1e-12 * S + # Backward-growing exponential rates from the frozen-coefficient + # spectrum (Re λ < 0 ⇒ grows as s decreases); rates ≲ 10/s are the + # algebraic power pair, not exponential growers. + λs = eigvals(Matrix(𝔐(s))) + gmin, gmax = Inf, 0.0 + for λ in λs + g = -real(λ) + if g > 10.0 / s + gmin = min(gmin, g) + gmax = max(gmax, g) + end + end + resolved = isfinite(gmin) && gmin * c * s < 8.0 + h = resolved ? -min(c * s, ztrk / gmax, s - s_m) : -min(c * s, s - s_m) + if resolved && !have_F + F = decaying_pair(params, Q, θ, s; growth=growth, seed=seed, dp_res=dp_res) + W .-= F * (F' * W) + have_F = true + Φ = 0.0 + end + W = radau_step(s, W, h) + if have_F + F = radau_step(s, F, h) + F = Matrix(qr(F).Q) + Φ += gmax * abs(h) + if Φ > purge_budget + W .-= F * (F' * W) + Φ = 0.0 + end + end + s += h + nstep += 1 + nstep > 5_000_000 && error("march_boundary: step limit exceeded (s=$s)") + end + # The boundary basis must span the traces AT s_m of the solutions that + # decay AT INFINITY — which, after the mode families morph around + # s ~ |Q|, is NOT the same subspace as the locally-decaying pair at s_m. + # The carried frame F is exactly the marched ∞-decaying pair (it enters + # the resolved band before the morphing region and tracks the subspace + # through it), so it IS the correct E-basis; re-extracting locally at + # s_m substitutes the wrong subspace and corrupts Δ at O(1) on + # extreme-coefficient surfaces. Local extraction remains only as the + # fallback when the march never entered the resolved band (no morphing + # crossed under resolution — the subspaces then coincide). + if have_F + W .-= F * (F' * W) + else + F = purge!(W, s_m) + end + return W[:, 1], W[:, 2], F +end + +""" + boundary_basis(params, Q, cache, θ, S; growth=30.0) + -> (Us, Ub, E, cond_est) + +Boundary data at s = S directly from the inps series (no march; used when +the series is already valid at the BVP edge): `Us`, `Ub` are the small/large +power-like solutions as 6-component states in the inps normalization, `E` +the numerically-generated decaying pair, `cond_est` the condition number of +the column-normalized basis. +""" +function boundary_basis(params::GGJParameters, Q::ComplexF64, + cache::InnerAsymptoticsCache, θ::Float64, S::Float64; + growth::Float64=30.0, seed::Int=20260707, dp_res::Float64=20.0) + xS = cis(θ) * S + ua, dua = physical_ua_dua(cache, xS) + Ub = ComplexF64[ua[1, 1], ua[2, 1], ua[3, 1], dua[1, 1], dua[2, 1], dua[3, 1]] + Us = ComplexF64[ua[1, 2], ua[2, 2], ua[3, 2], dua[1, 2], dua[2, 2], dua[3, 2]] + E = decaying_pair(params, Q, θ, S; growth=growth, seed=seed, dp_res=dp_res) + + Bmat = hcat(Us ./ norm(Us), Ub ./ norm(Ub), E) + cond_est = cond(Bmat) + return Us, Ub, E, cond_est +end + + +# solve.jl +# +# Global spectral-element collocation solve of the GGJ inner-layer BVP on +# the rotated ray, with the matching data Δ as a bordered unknown. +# +# For each parity (deltac isol = 1 "odd", 2 "even"): +# +# unknowns : v at all global nodes (6 components each), plus (Δ, c₁, c₂) +# equations: dv/ds = e^{iθ} M(e^{iθ}s) v collocated at the p non-left +# Lobatto nodes of every cell (right-biased collocation — +# Radau-IIA-like, stable for the stiff modes); +# 3 parity conditions at s = 0; +# 6 matching conditions at s = S: +# v(S) − Δ·Ub − c₁E₁ − c₂E₂ = Us . +# +# One sparse LU per parity per mesh; residual-driven bisection refinement. +# The exponential dichotomy never enters: no quantity is ever propagated +# across the layer, and the boundary condition splits the dichotomy exactly. + +""" + RaySolveResult + +Output of [`solve_ray`](@ref). `Δ` is the rescaled matching data in the +deltac output ordering (index 1 ↔ even-BC solution, 2 ↔ odd-BC solution, +i.e. the same swap deltac_solve applies); `Δraw` is (isol=1, isol=2) in raw +inps normalization before `rescale_delta`. +""" +struct RaySolveResult + Δ::SVector{2,ComplexF64} + Δraw::SVector{2,ComplexF64} + cexp::SMatrix{2,2,ComplexF64,4} # decaying-pair coefficients (col = isol) + Q::ComplexF64 + θ::Float64 + S::Float64 + series_ok::Bool # pick_smax reached its tolerance + bc_cond::Float64 # boundary-basis condition number + breaks::Vector{Float64} + p::Int + sols::Matrix{ComplexF64} # ndof × 2 nodal solution (isol = 1, 2) + resid::Float64 # final max relative collocation residual + δΔ_refine::Float64 # |Δ change| in last refinement round (rel) + nrounds::Int +end + +# ----------------------------------------------------------------------- +# Assembly. +# ----------------------------------------------------------------------- + +function _assemble_base(params::GGJParameters, Q::ComplexF64, θ::Float64, + breaks::Vector{Float64}, p::Int, + t::Vector{Float64}, D::Matrix{Float64}, + Us::Vector{ComplexF64}, Ub::Vector{ComplexF64}, E::Matrix{ComplexF64}) + N = length(breaks) - 1 + Ng = N * p + 1 + ndof = 6 * Ng + 3 + ph = cis(θ) + βb = norm(Ub) + + nnz_est = N * p * 6 * (p + 7) + 24 + 3 + Ic = Vector{Int}() + Jc = Vector{Int}() + Vc = Vector{ComplexF64}() + sizehint!(Ic, nnz_est) + sizehint!(Jc, nnz_est) + sizehint!(Vc, nnz_est) + + @inbounds for c in 1:N + a, b = breaks[c], breaks[c+1] + Dscale = 2 / (b - a) + for j in 1:p + s_j = a + (b - a) * (t[j+1] + 1) / 2 + M = ph * ode_matrix(params, Q, ph * s_j) + row0 = 6 * ((c - 1) * p + (j - 1)) + gj = (c - 1) * p + j + 1 + for i in 1:6 + row = row0 + i + for k in 0:p + gk = (c - 1) * p + k + 1 + push!(Ic, row) + push!(Jc, 6 * (gk - 1) + i) + push!(Vc, Dscale * D[j+1, k+1]) + end + for i2 in 1:6 + Mv = M[i, i2] + if Mv != 0 + push!(Ic, row) + push!(Jc, 6 * (gj - 1) + i2) + push!(Vc, -Mv) + end + end + end + end + end + + # Matching rows at s = S. + rowm0 = 6 * N * p + @inbounds for i in 1:6 + row = rowm0 + i + push!(Ic, row) + push!(Jc, 6 * (Ng - 1) + i) + push!(Vc, one(ComplexF64)) + push!(Ic, row) + push!(Jc, 6 * Ng + 1) + push!(Vc, -Ub[i] / βb) + push!(Ic, row) + push!(Jc, 6 * Ng + 2) + push!(Vc, -E[i, 1]) + push!(Ic, row) + push!(Jc, 6 * Ng + 3) + push!(Vc, -E[i, 2]) + end + + rhs = zeros(ComplexF64, ndof) + rhs[rowm0 .+ (1:6)] .= Us + + return Ic, Jc, Vc, rhs, βb, ndof, Ng +end + +""" + _solve_parities(Ic, Jc, Vc, rhs, ndof, N, p) -> (sol_odd, sol_even) + +Solve both parity systems from one factorization: the two matrices differ only +in the 3 parity rows at s = 0 (a rank-3 update), so one sparse LU of the +row-equilibrated odd system plus a Woodbury correction replaces the second +factorization, which dominates the linear-algebra cost. +""" +function _solve_parities(Ic::Vector{Int}, Jc::Vector{Int}, Vc::Vector{ComplexF64}, + rhs::Vector{ComplexF64}, ndof::Int, N::Int, p::Int) + rowp0 = 6 * N * p + 6 + c1 = parity_rows(1) + c2 = parity_rows(2) + I2 = copy(Ic) + J2 = copy(Jc) + V2 = copy(Vc) + for k in 1:3 + push!(I2, rowp0 + k) + push!(J2, c1[k]) + push!(V2, one(ComplexF64)) + end + A = sparse(I2, J2, V2, ndof, ndof) + + # Row equilibration. + rmax = zeros(Float64, ndof) + rows = rowvals(A) + vals = nonzeros(A) + for col in 1:ndof + for idx in nzrange(A, col) + r = rows[idx] + a = abs(vals[idx]) + a > rmax[r] && (rmax[r] = a) + end + end + @inbounds for i in 1:ndof + rmax[i] = rmax[i] > 0 ? 1 / rmax[i] : 1.0 + end + Dr = Diagonal(rmax) + F = lu(Dr * A) + x1 = F \ (Dr * rhs) + + # Woodbury: (DrA + U Vᵀ)⁻¹ b = x1 − A⁻¹U (I₃ + Vᵀ A⁻¹U)⁻¹ Vᵀ x1, + # with U[:,k] = Dr[rₖ,rₖ]·e_{rₖ} and Vᵀ[k,:] = (e_{c2ₖ} − e_{c1ₖ})ᵀ. + U = zeros(ComplexF64, ndof, 3) + for k in 1:3 + U[rowp0+k, k] = rmax[rowp0+k] + end + AinvU = F \ U + Vx1 = ComplexF64[x1[c2[k]] - x1[c1[k]] for k in 1:3] + Smat = Matrix{ComplexF64}(I, 3, 3) + for k in 1:3, j in 1:3 + Smat[k, j] += AinvU[c2[k], j] - AinvU[c1[k], j] + end + x2 = x1 - AinvU * (Smat \ Vx1) + # UMFPACK numeric factorizations live in C malloc, invisible to the GC's + # heap accounting — under the refinement loop, lazily-finalized multi-GB + # factorizations pile up (observed: tens of GB on 24k-cell meshes). Free + # eagerly now that both solutions are extracted. + finalize(F) + return x1, x2 +end + +# ----------------------------------------------------------------------- +# Residual estimation for adaptive refinement. +# ----------------------------------------------------------------------- + +function _cell_errors(params::GGJParameters, Q::ComplexF64, θ::Float64, + breaks::Vector{Float64}, p::Int, + t::Vector{Float64}, D::Matrix{Float64}, w::Vector{Float64}, + sol::Vector{ComplexF64}) + N = length(breaks) - 1 + ph = cis(θ) + errs = zeros(Float64, N) + + # Reference-cell test points: left endpoint + inter-node midpoints. + τs = Float64[-1.0] + for j in 1:p + push!(τs, (t[j] + t[j+1]) / 2) + end + + vals = Matrix{ComplexF64}(undef, p + 1, 6) + for c in 1:N + a, b = breaks[c], breaks[c+1] + Dscale = 2 / (b - a) + for j in 0:p + g = (c - 1) * p + j + 1 + for i in 1:6 + vals[j+1, i] = sol[6*(g-1)+i] + end + end + dvals = (Dscale .* D) * vals + + emax = 0.0 + for τ in τs + vτ = barycentric_eval(t, w, vals, τ) + dvτ = barycentric_eval(t, w, dvals, τ) + s_τ = a + (b - a) * (τ + 1) / 2 + Mv = (ph * ode_matrix(params, Q, ph * s_τ)) * vτ + rnorm = 0.0 + scale = 1e-300 + for i in 1:6 + rnorm = max(rnorm, abs(dvτ[i] - Mv[i])) + scale = max(scale, abs(dvτ[i]), abs(Mv[i])) + end + emax = max(emax, rnorm / scale) + end + errs[c] = emax + end + return errs +end + +# ----------------------------------------------------------------------- +# Driver. +# ----------------------------------------------------------------------- + +""" + solve_ray(params::GGJParameters, Q; θ=angle(Q)/4, kmax=12, p=12, S=nothing, + smax_tol=1e-9, refine_tol=1e-8, max_rounds=8, max_cells=6000, + verbose=false, kwargs...) -> RaySolveResult + +Solve the GGJ inner-layer problem on the rotated ray x = e^{iθ}s. +`Δ` follows the deltac output convention (swap + `rescale_delta`), directly +comparable to the `:galerkin` backend. θ = 0 reproduces a real-axis solve; +`S` defaults to the inps series-residual radius (`pick_smax`); `p` is the +polynomial order per spectral element; remaining keywords are march/mesh +knobs (see `march_boundary`, `initial_breaks`, `decaying_pair`). +""" +function solve_ray(params::GGJParameters, Q::ComplexF64; + θ::Float64=angle(Q) / 4, kmax::Int=12, p::Int=12, + S::Union{Nothing,Float64}=nothing, smax_tol::Float64=1e-9, + growth::Float64=30.0, refine_tol::Float64=1e-8, + march_rtol::Float64=1e-9, sm_fac::Float64=1.0, cfl::Float64=0.4, + max_rounds::Int=8, max_cells::Int=6000, seed::Int=20260707, + ztrk::Float64=0.3, purge_budget::Float64=2.5, dp_res::Float64=20.0, + handoff::Union{Nothing,NTuple{2,Vector{ComplexF64}}}=nothing, + boundary::Union{Nothing,Tuple{Vector{ComplexF64},Vector{ComplexF64},Matrix{ComplexF64}}}=nothing, + verbose::Bool=false) + + cache = build_asymptotics(params, Q; kmax=kmax) + series_ok = true + if S === nothing + Sv, cache, series_ok = pick_smax(params, Q; θ=θ, eps=smax_tol, kmax=kmax, cache=cache) + series_ok || @warn "pick_smax: series residual target $smax_tol not reached; using best point S=$Sv" + else + Sv = S + end + + # Decide the BVP outer radius: the layer physics (exponential content + # above machine floor) ends at s_d; beyond that the solution is pure + # power law and is handled by the projected quotient march, not by + # collocation. For small |Q| (s_d ≳ S) no march is needed. + sq = sqrt(Q) + ρ = real(cis(2θ) / sq) + # BVP outer radius s_m: where every exponential solution family has + # decayed ≥ 45 e-folds from the origin, computed from the actual + # frozen-coefficient spectrum (the asymptotic ρs rate is wrong near the + # origin for extreme-coefficient surfaces, where the Υ-family rate + # √|Q(G+KF)| dominates and kills exponential content within s ~ O(1)). + s_m = let acc = 0.0, sprev = 0.0, out = Sv + for sk in exp10.(range(-2, log10(Sv); length=400)) + gdec = Inf + for λ in eigvals(Matrix(cis(θ) * ode_matrix(params, Q, cis(θ) * sk))) + g = -real(λ) # backward-growth = forward-decay rate + g > 10.0 / sk && (gdec = min(gdec, g)) + end + isfinite(gdec) && (acc += gdec * (sk - sprev)) + sprev = sk + if acc >= 45.0 + out = sk + break + end + end + sm_fac * max(out, 2.0) + end + + local Us, Ub, E, bc_cond, Sbvp + if boundary !== nothing + # diagnostic: externally supplied boundary data (Us, Ub, E) at s_m, + # e.g. from an extended-precision march + Us, Ub, E = boundary + Sbvp = s_m + bc_cond = cond(hcat(Us ./ norm(Us), Ub ./ norm(Ub), E)) + elseif s_m < 0.98 * Sv + local UsS, UbS + if handoff === nothing + ua, dua = physical_ua_dua(cache, cis(θ) * Sv) + UbS = ComplexF64[ua[1, 1], ua[2, 1], ua[3, 1], dua[1, 1], dua[2, 1], dua[3, 1]] + UsS = ComplexF64[ua[1, 2], ua[2, 2], ua[3, 2], dua[1, 2], dua[2, 2], dua[3, 2]] + else + # diagnostic: externally supplied series hand-off (Us, Ub) at + # x = e^{iθ}S, e.g. evaluated in extended precision + UsS, UbS = handoff + end + Us, Ub, E = march_boundary(params, Q, θ, Sv, s_m, UsS, UbS; + rtol=march_rtol, growth=growth, seed=seed, + ztrk=ztrk, purge_budget=purge_budget, dp_res=dp_res) + Sbvp = s_m + bc_cond = cond(hcat(Us ./ norm(Us), Ub ./ norm(Ub), E)) + verbose && @printf(" march: S=%.4g → s_m=%.4g, bc_cond=%.1e\n", Sv, s_m, bc_cond) + else + Us, Ub, E, bc_cond = boundary_basis(params, Q, cache, θ, Sv; + growth=growth, seed=seed, dp_res=dp_res) + Sbvp = Sv + end + bc_cond < 1e12 || @warn "boundary basis poorly conditioned (cond ≈ $bc_cond)" + + t, D = cheblobatto(p) + w = barycentric_weights(t) + breaks = initial_breaks(params, Q, θ, Sbvp, p; cfl=cfl) + + Δraw = MVector{2,ComplexF64}(0, 0) + Δprev = MVector{2,ComplexF64}(NaN, NaN) + cexp = zeros(MMatrix{2,2,ComplexF64}) + local sols + resid = Inf + resid_prev = Inf + nstall = 0 + δΔ = Inf + round_used = 0 + + for round in 0:max_rounds + round_used = round + N = length(breaks) - 1 + Ic, Jc, Vc, rhs, βb, ndof, Ng = _assemble_base(params, Q, θ, breaks, p, t, D, Us, Ub, E) + sols = Matrix{ComplexF64}(undef, ndof, 2) + errs = zeros(Float64, N) + sol_odd, sol_even = _solve_parities(Ic, Jc, Vc, rhs, ndof, N, p) + for (isol, sol) in ((1, sol_odd), (2, sol_even)) + sols[:, isol] = sol + Δraw[isol] = sol[6*Ng+1] / βb + cexp[1, isol] = sol[6*Ng+2] + cexp[2, isol] = sol[6*Ng+3] + errs = max.(errs, _cell_errors(params, Q, θ, breaks, p, t, D, w, sol)) + end + resid = maximum(errs) + δΔ = maximum(abs.(Δraw .- Δprev) ./ (abs.(Δraw) .+ 1e-300)) + Δprev .= Δraw + verbose && @printf(" round %d: %d cells, resid %.2e, Δraw[1] = %.6e %+.6eim, δΔ = %.1e\n", + round, N, resid, real(Δraw[1]), imag(Δraw[1]), δΔ) + + (resid < refine_tol || round == max_rounds) && break + # Stagnation guard: on extreme-coefficient surfaces the residual + # bottoms out at the roundoff floor ε·‖M‖ (~10⁻⁸ ≫ refine_tol); + # refining past that plateau burns cells without informing Δ. Two + # consecutive rounds with < 2× improvement ⇒ converged-to-plateau. + nstall = resid > 0.5 * resid_prev ? nstall + 1 : 0 + resid_prev = resid + if nstall >= 2 + verbose && @printf(" refinement stagnated at resid %.2e (roundoff plateau) — stopping\n", resid) + break + end + marked = findall(>(refine_tol), errs) + isempty(marked) && break + if N + length(marked) > max_cells + @warn "solve_ray: cell budget ($max_cells) reached with resid=$resid" + break + end + newbreaks = Float64[breaks[1]] + for c in 1:N + if c in marked + push!(newbreaks, (breaks[c] + breaks[c+1]) / 2) + end + push!(newbreaks, breaks[c+1]) + end + breaks = newbreaks + end + + Δswapped = SVector{2,ComplexF64}(Δraw[2], Δraw[1]) # deltac_solve convention + Δ = rescale_delta(Δswapped, params) + + return RaySolveResult(Δ, SVector(Δraw), SMatrix(cexp), Q, θ, Sv, series_ok, + bc_cond, breaks, p, sols, resid, δΔ, round_used) +end + +solve_ray(params::GGJParameters, Q::Number; kwargs...) = + solve_ray(params, ComplexF64(Q); kwargs...) + +""" + solve_inner(::GGJModel{:ray}, params::GGJParameters, γ::Number; kwargs...) + -> SVector{2,ComplexF64} + +Rotated-ray backend: converts γ via `inner_Q`, forwards keywords to +[`solve_ray`](@ref), and returns `(Δ₁, Δ₂)` in the shared backend convention +(deltac swap + physical rescale). Preferred for |Q| ≳ 1 and near the +imaginary axis; use `solve_ray` directly when the full +[`RaySolveResult`](@ref) is wanted. +""" +function solve_inner(::GGJModel{:ray}, params::GGJParameters, γ::Number; kwargs...) + res = solve_ray(params, inner_Q(params, γ); kwargs...) + return res.Δ +end + +""" + solve_inner_profile(::GGJModel{:ray}, params::GGJParameters, γ::Number; + npc=8, certify_rtol=1e-3, kwargs...) + -> (; Δ, x, Ψ, Ξ, dψdx, rescale, certΔ) + +Rotated-ray implementation of the [`solve_inner_profile`](@ref) interface. The +certified `Δ` comes from the optimal-contour solve at θ = arg(Q)/4 (robust for +|Q| ≳ 1, where real-axis methods drift); the profiles come from a θ = 0 +re-solve on the real axis, valid at physical (RPEC) |Q| since the on-axis +pseudo-resonance is a regular point resolved by the BVP refinement. The +θ=0-vs-optimal-θ Δ agreement doubles as a runtime certificate (two maximally +different contours through the same entire solutions): the relative mismatch +is returned as `certΔ` and warns above `certify_rtol`. `npc` sets the profile +points per mesh cell; extra keywords forward to both [`solve_ray`](@ref) calls +(θ is fixed by the method — do not pass it). +""" +function solve_inner_profile(::GGJModel{:ray}, params::GGJParameters, γ::Number; + npc::Int=8, certify_rtol::Float64=1e-3, kwargs...) + Q = inner_Q(params, γ) + rr = solve_ray(params, Q; kwargs...) # certified Δ (optimal θ) + r0 = solve_ray(params, Q; θ=0.0, kwargs...) # real-axis profile solve + certΔ = abs(r0.Δ[1] - rr.Δ[1]) / max(abs(rr.Δ[1]), 1e-300) + certΔ > certify_rtol && @warn "GGJ ray: θ=0 profile solve disagrees with certified Δ" ising = params.ising Q certΔ + prof = solution_profile(r0; npc=npc) + return (; Δ=rr.Δ, x=prof.s, Ψ=prof.Ψ, Ξ=prof.Ξ, _profile_conversions(params)..., certΔ=certΔ) +end + +""" + evaluate_solution(res::RaySolveResult, s::Real; isol=1) -> SVector{6,ComplexF64} + +Evaluate the solved state v = (Ψ, Ξ, Υ, Ψ', Ξ', Υ') at ray parameter `s` +(diagnostic; barycentric interpolation on the containing cell). +""" +function evaluate_solution(res::RaySolveResult, s::Real; isol::Int=1) + breaks = res.breaks + p = res.p + (breaks[1] <= s <= breaks[end]) || throw(ArgumentError("s=$s outside [0, $(breaks[end])]")) + c = clamp(searchsortedlast(breaks, s), 1, length(breaks) - 1) + a, b = breaks[c], breaks[c+1] + t, _ = cheblobatto(p) + w = barycentric_weights(t) + vals = Matrix{ComplexF64}(undef, p + 1, 6) + for j in 0:p + g = (c - 1) * p + j + 1 + for i in 1:6 + vals[j+1, i] = res.sols[6*(g-1)+i, isol] + end + end + τ = 2 * (s - a) / (b - a) - 1 + return SVector{6,ComplexF64}(barycentric_eval(t, w, vals, τ)...) +end + +""" + solution_profile(res::RaySolveResult; npc=8) -> (; s, x, Ψ, Ξ, Υ) + +Reconstruct (Ψ, Ξ, Υ) on the solution contour, `npc` points per cell; each +field is `npts × 2` with columns (odd, even) and `x = e^{iθ}s`. For rotated +solves the fields live on the ray — continuation back to real x is a +separate (unstable) problem, by design. +""" +function solution_profile(res::RaySolveResult; npc::Int=8) + breaks = res.breaks + p = res.p + t, _ = cheblobatto(p) + w = barycentric_weights(t) + N = length(breaks) - 1 + npts = N * npc + 1 + s = Vector{Float64}(undef, npts) + Ψ = Matrix{ComplexF64}(undef, npts, 2) + Ξ = Matrix{ComplexF64}(undef, npts, 2) + Υ = Matrix{ComplexF64}(undef, npts, 2) + vals = Matrix{ComplexF64}(undef, p + 1, 6) + k = 0 + for c in 1:N + a, b = breaks[c], breaks[c+1] + for m in 0:(npc-1) + k += 1 + τ = -1 + 2 * m / npc + s[k] = a + (b - a) * (τ + 1) / 2 + for isol in 1:2 + for j in 0:p, i in 1:6 + vals[j+1, i] = res.sols[6*((c-1)*p+j)+i, isol] + end + v = barycentric_eval(t, w, vals, τ) + Ψ[k, isol] = v[1] + Ξ[k, isol] = v[2] + Υ[k, isol] = v[3] + end + end + end + # right endpoint + k += 1 + s[k] = breaks[end] + for isol in 1:2 + v = evaluate_solution(res, breaks[end]; isol=isol) + Ψ[k, isol] = v[1] + Ξ[k, isol] = v[2] + Υ[k, isol] = v[3] + end + return (; s=s, x=cis(res.θ) .* s, Ψ=Ψ, Ξ=Ξ, Υ=Υ) +end + +""" + asymptotic_profile(params, res::RaySolveResult, srange) -> (; s, x, Ψ, Ξ, Υ) + +Evaluate the analytic representation `u_small + Δ·u_big` on the ray for +`s ≥ res.S`; concatenated with [`solution_profile`](@ref) it should join +the numeric solution seamlessly across the march zone. +""" +function asymptotic_profile(params::GGJParameters, res::RaySolveResult, + srange::AbstractVector{<:Real}) + cache = build_asymptotics(params, res.Q; kmax=12) + n = length(srange) + Ψ = Matrix{ComplexF64}(undef, n, 2) + Ξ = Matrix{ComplexF64}(undef, n, 2) + Υ = Matrix{ComplexF64}(undef, n, 2) + for (k, s) in enumerate(srange) + ua, _ = physical_ua_dua(cache, cis(res.θ) * s) + for isol in 1:2 + Δ = res.Δraw[isol] + Ψ[k, isol] = ua[1, 2] + Δ * ua[1, 1] + Ξ[k, isol] = ua[2, 2] + Δ * ua[2, 1] + Υ[k, isol] = ua[3, 2] + Δ * ua[3, 1] + end + end + return (; s=collect(float.(srange)), x=cis(res.θ) .* srange, Ψ=Ψ, Ξ=Ξ, Υ=Υ) +end + +""" + delta_convergence(params, Q; verbose=true, refine_tol=1e-9, + max_rounds=16, max_cells=12000, kwargs...) + -> (; Δ, spread, base, table) + +Error bar for the matching data: solve at a strict baseline, then re-solve +with each numerical knob perturbed on an independent axis (contour angle θ, +element order p, series order/radius, refinement depth, march tolerance, +handoff radius, purification e-folds — see `variations` below) and report the +relative change of (Δ₁, Δ₂) per knob. θ is perturbed outward only: inward +moves toward the x² ≈ −Q²(G+KF) pseudo-resonance and measures the perturbed +solve, not the baseline. `spread` is the per-parity worst case over all knobs. +""" +function delta_convergence(params::GGJParameters, Q::ComplexF64; + verbose::Bool=true, refine_tol::Float64=1e-9, + max_rounds::Int=16, max_cells::Int=12000, kwargs...) + base = solve_ray(params, Q; refine_tol=refine_tol, max_rounds=max_rounds, + max_cells=max_cells, kwargs...) + variations = [ + ("θ → 1.2θ", (; θ=(base.θ == 0 ? 0.1 : 1.2 * base.θ))), + ("p → p+4", (; p=base.p + 4)), + ("kmax+4, S↑", (; kmax=16, smax_tol=1e-10)), + ("refine_tol/10", (; refine_tol=refine_tol / 10)), + ("march_rtol/100", (; march_rtol=1e-11)), + ("sm_fac 1→1.4", (; sm_fac=1.4)), + ("growth 30→45", (; growth=45.0)) + ] + table = NamedTuple[] + spread = MVector{2,Float64}(0.0, 0.0) + verbose && @printf("Δ convergence check, Q = %s:\n", string(Q)) + verbose && @printf(" baseline: Δ₁ = %+.8e %+.8eim, Δ₂ = %+.8e %+.8eim\n", + real(base.Δ[1]), imag(base.Δ[1]), real(base.Δ[2]), imag(base.Δ[2])) + for (name, kw) in variations + r = solve_ray(params, Q; refine_tol=refine_tol, max_rounds=max_rounds, + max_cells=max_cells, kwargs..., kw...) + d1 = abs(r.Δ[1] - base.Δ[1]) / abs(base.Δ[1]) + d2 = abs(r.Δ[2] - base.Δ[2]) / abs(base.Δ[2]) + spread[1] = max(spread[1], d1) + spread[2] = max(spread[2], d2) + push!(table, (; name, d1, d2)) + verbose && @printf(" %-16s δΔ₁ = %.2e δΔ₂ = %.2e\n", name, d1, d2) + end + verbose && @printf(" %-16s δΔ₁ = %.2e δΔ₂ = %.2e\n", "WORST CASE", spread[1], spread[2]) + return (; Δ=base.Δ, spread=SVector(spread), base=base, table=table) +end + +delta_convergence(params::GGJParameters, Q::Number; kwargs...) = + delta_convergence(params, ComplexF64(Q); kwargs...) diff --git a/src/InnerLayer/GGJ/RayAsymptotics.jl b/src/InnerLayer/GGJ/RayAsymptotics.jl new file mode 100644 index 00000000..630d5c39 --- /dev/null +++ b/src/InnerLayer/GGJ/RayAsymptotics.jl @@ -0,0 +1,237 @@ +# RayAsymptotics.jl +# +# Complex-x extensions of the inps Wasow asymptotic kernel +# (InnerAsymptotics.jl) for the rotated-ray solver (Ray.jl). +# +# The inps construction (GW2020 Eqs. 3–52) depends only on Q — the contour +# enters only through the EVALUATION point. These methods therefore reuse the +# existing, validated `InnerAsymptoticsCache` and coefficient packing verbatim +# and add evaluation at complex x = e^{iθ}s (principal branch of x^r, valid +# for |arg x| < π), the plain-variable conversion used by the ray boundary +# data, the GW2020 Eq. (54) residual along the ray, and the series-radius +# selector `pick_smax` (the ray analog of `pick_xmax`). +# +# The real-x methods of InnerAsymptotics.jl are untouched. + +# Complex-argument Horner evaluator: same contract as _horner(x::Real, ...) +# in InnerAsymptotics.jl, with x^rvec on the principal branch. +function _horner(x::Complex, c::AbstractMatrix{ComplexF64}; + rvec::Union{Nothing,AbstractVector{<:Real}}=nothing, + derivative::Bool=false) + nrows, ncols = size(c) + n = ncols - 1 + + y = Vector{ComplexF64}(undef, nrows) + @inbounds for i in 1:nrows + y[i] = c[i, ncols] + end + @inbounds for k in (n-1):-1:0 + for i in 1:nrows + y[i] = y[i] * x + c[i, k+1] + end + end + + if rvec !== nothing + @inbounds for i in 1:nrows + y[i] *= ComplexF64(x)^rvec[i] + end + end + + if !derivative + return y, nothing + end + + dy = similar(y) + if rvec !== nothing + @inbounds for i in 1:nrows + dy[i] = c[i, ncols] * (rvec[i] + n) + end + @inbounds for k in (n-1):-1:0 + for i in 1:nrows + dy[i] = dy[i] * x + c[i, k+1] * (rvec[i] + k) + end + end + @inbounds for i in 1:nrows + dy[i] = dy[i] * ComplexF64(x)^rvec[i] / x + end + else + @inbounds for i in 1:nrows + dy[i] = c[i, ncols] * n + end + @inbounds for k in (n-1):-1:0 + for i in 1:nrows + dy[i] = dy[i] * x + c[i, k+1] * k + end + end + @inbounds for i in 1:nrows + dy[i] = dy[i] / x + end + end + return y, dy +end + +""" + evaluate_asymptotics(cache, x::Complex; derivative=true, apply_T=true) -> (U, dU) + +Complex-x version of [`evaluate_asymptotics`](@ref) (GW2020 Eq. 53; principal +branch of `x^r`, valid for |arg x| < π). Columns of the 6×2 `U` are the two +power-like solutions in the GW2020 Eq. (2) state convention +`u = (xΨ, Ξ, Υ, (xΨ)'/x, Ξ'/x, Υ'/x)`. +""" +function evaluate_asymptotics(cache::InnerAsymptoticsCache, x::Complex; + derivative::Bool=true, apply_T::Bool=true) + xfac = 1 / (x * x) + R = cache.R + rvec = SVector(-R[1] / 2, -R[1] / 2, -R[2] / 2, -R[2] / 2) + + cc = _pack_y_coefs(cache) + dd = _pack_qp_coefs(cache) + + yy, dyy = _horner(xfac, cc; rvec=rvec, derivative=derivative) + zz, dzz = _horner(xfac, dd; derivative=derivative) + + y = SMatrix{2,2,ComplexF64}(yy[1], yy[2], yy[3], yy[4]) + q = SMatrix{2,2,ComplexF64}(zz[1], zz[2], zz[3], zz[4]) + p21 = SMatrix{4,2,ComplexF64}(zz[5], zz[6], zz[7], zz[8], + zz[9], zz[10], zz[11], zz[12]) + + pp_m = zeros(ComplexF64, 6, 2) + pp_m[1, 1] = 1 + pp_m[2, 2] = 1 + @inbounds for i in 1:4, j in 1:2 + pp_m[i+2, j] = p21[i, j] + end + pp = SMatrix{6,2,ComplexF64}(pp_m) + + smat = @SMatrix ComplexF64[1 0; 0 xfac] + + qsy = q * smat * y + U = pp * qsy + + if derivative + chain = -2 * xfac / x + dy_v = chain .* dyy + dz_v = chain .* dzz + + dy = SMatrix{2,2,ComplexF64}(dy_v[1], dy_v[2], dy_v[3], dy_v[4]) + dq = SMatrix{2,2,ComplexF64}(dz_v[1], dz_v[2], dz_v[3], dz_v[4]) + dp21 = SMatrix{4,2,ComplexF64}(dz_v[5], dz_v[6], dz_v[7], dz_v[8], + dz_v[9], dz_v[10], dz_v[11], dz_v[12]) + + dpp_m = zeros(ComplexF64, 6, 2) + @inbounds for i in 1:4, j in 1:2 + dpp_m[i+2, j] = dp21[i, j] + end + dpp = SMatrix{6,2,ComplexF64}(dpp_m) + + dsmat = @SMatrix ComplexF64[0 0; 0 (-2*xfac/x)] + dqsy = q * smat * dy + q * dsmat * y + dq * smat * y + dU = pp * dqsy + dpp * qsy + + if apply_T + U = cache.T * U + dU = cache.T * dU + end + return U, dU + else + if apply_T + U = cache.T * U + end + return U, nothing + end +end + +""" + physical_ua_dua(cache, x::Number) -> (ua, dua) + +Convert the inps 6×2 basis at (possibly complex) `x` to physical `(Ψ, Ξ, Υ)` +values and x-derivatives, each 3×2 (column 1 = large power solution, 2 = small). +Same map as the deltac `inpso_get_ua/dua` convention. +""" +function physical_ua_dua(cache::InnerAsymptoticsCache, x::Number) + xc = ComplexF64(x) + U, dU = evaluate_asymptotics(cache, xc; derivative=true, apply_T=true) + ua = zeros(ComplexF64, 3, 2) + dua = zeros(ComplexF64, 3, 2) + @inbounds for j in 1:2 + ua[1, j] = U[1, j] / xc + ua[2, j] = U[2, j] + ua[3, j] = U[3, j] + dua[1, j] = U[4, j] - U[1, j] / (xc * xc) + dua[2, j] = U[5, j] * xc + dua[3, j] = U[6, j] * xc + end + return ua, dua +end + +""" + asymptotic_residual(cache, x::Complex) -> SVector{2,Float64} + +GW2020 Eq. (54) residual of the two series columns at complex `x`. +""" +function asymptotic_residual(cache::InnerAsymptoticsCache, x::Complex) + U, dU = evaluate_asymptotics(cache, x; derivative=true, apply_T=false) + + xfac = 1 / (x * x) + M = cache.J[1] + if cache.kmax > 0 + M = M + xfac * cache.J[2] + end + if cache.kmax > 1 + M = M + xfac * xfac * cache.J[3] + end + + matvec1 = dU + matvec2 = -x * (M * U) + matvec0 = matvec1 + matvec2 + + delta = MVector{2,Float64}(0.0, 0.0) + @inbounds for j in 1:2 + n0 = 0.0 + n1 = 0.0 + n2 = 0.0 + for i in 1:6 + n0 = max(n0, abs(matvec0[i, j])) + n1 = max(n1, abs(matvec1[i, j])) + n2 = max(n2, abs(matvec2[i, j])) + end + denom = max(n1, n2) + delta[j] = denom == 0 ? 0.0 : n0 / denom + end + return SVector{2,Float64}(delta) +end + +""" + pick_smax(params, Q; θ=angle(Q)/4, eps=1e-9, kmax=12, cache=nothing, + slogmin=-1.0, slogmax=6.5, dslog=0.01) -> (S, cache, achieved) + +Ray analog of [`pick_xmax`](@ref): the smallest `s` at which the series +residual along `x = e^{iθ}s` drops below `eps`. If never reached, returns the +residual-minimizing `s` with `achieved = false`. +""" +function pick_smax(params::GGJParameters, Q::ComplexF64; + θ::Float64=angle(Q) / 4, eps::Float64=1e-9, kmax::Int=12, + cache::Union{Nothing,InnerAsymptoticsCache}=nothing, + slogmin::Float64=-1.0, slogmax::Float64=6.5, dslog::Float64=0.01) + c = cache === nothing ? build_asymptotics(params, Q; kmax=kmax) : cache + ph = cis(θ) + best_s = 10.0^slogmin + best_r = Inf + slog = slogmin + while slog <= slogmax + s = 10.0^slog + r = maximum(asymptotic_residual(c, ph * s)) + if r < eps + return s, c, true + end + if r < best_r + best_r = r + best_s = s + end + slog += dslog + end + return best_s, c, false +end + +pick_smax(params::GGJParameters, Q::Number; kwargs...) = + pick_smax(params, ComplexF64(Q); kwargs...) diff --git a/src/InnerLayer/GGJ/Reference.jl b/src/InnerLayer/GGJ/Reference.jl index 6127cbfd..16098a28 100644 --- a/src/InnerLayer/GGJ/Reference.jl +++ b/src/InnerLayer/GGJ/Reference.jl @@ -5,21 +5,13 @@ # corresponding Fortran RMATCH test cases. """ - glasser_wang_2020_eq55() -> GGJParameters + glasser_wang_2020_eq55(; taua=1.0, taur=1e6, v1=1.0) -> GGJParameters -D-shaped aspect-ratio-2, q = 2 surface from Glasser & Wang, Phys. Plasmas -**27**, 012506 (2020), Eq. 55. This is the primary benchmark case for -validating the inps Wasow basis convergence (their Figs. 1–4). This is useful -only for benchmarking the galerkin solver and comparing to published results. - -The five coefficients below are transcribed verbatim from Eq. 55; the paper's -companion operating point is the scaled growth rate `Q = 1.234e-1` (their Fig. 1). -Note Eq. 55 does not tabulate an inner-region matching `Δ(Q)` — its `Δ_±` (Eq. 54) -is a convergence-error norm — so a quantitative `Δ(Q)` cross-check needs a Fortran -rmatch/INPS run, not this paper alone. - -Timescale parameters (taua, taur, v1) are set to canonical normalization; -callers should override them for physical cases. +D-shaped aspect-ratio-2, q = 2 surface from GW2020 Eq. 55 (companion operating +point `Q = 1.234e-1`, their Fig. 1); benchmark for the inps basis convergence +(their Figs. 1–4). Eq. 55 does not tabulate `Δ(Q)` — its `Δ_±` (Eq. 54) is a +convergence norm — so quantitative Δ cross-checks need a Fortran rmatch/inps +run. Timescales default to canonical normalization; override for physical cases. """ function glasser_wang_2020_eq55(; taua::Float64=1.0, taur::Float64=1e6, v1::Float64=1.0) return GGJParameters(; @@ -27,3 +19,18 @@ function glasser_wang_2020_eq55(; taua::Float64=1.0, taur::Float64=1e6, v1::Floa taua=taua, taur=taur, v1=v1 ) end + +""" + q4_surface_benchmark() -> GGJParameters + +Physical q = 4 rational-surface benchmark point (S = τ_R/τ_A ≈ 4.58×10⁶, +D_I ≈ −0.31166, α = √(−D_I) ≈ 0.5583). Primary validation point for the +rotated-ray backend on the imaginary-Q axis (pinned at Q = 100i, 500i in +the test suite). +""" +function q4_surface_benchmark() + return GGJParameters(; + E=-0.13733, F=0.022202, G=7.60633, H=0.053468, K=14.66987, + M=30.26883, taua=2.11226e-7, taur=0.968219, v1=1.55009 + ) +end diff --git a/src/InnerLayer/GGJ/Shooting.jl b/src/InnerLayer/GGJ/Shooting.jl index 7a3162d2..e0f8ce45 100644 --- a/src/InnerLayer/GGJ/Shooting.jl +++ b/src/InnerLayer/GGJ/Shooting.jl @@ -62,11 +62,11 @@ function _build_origin_arrays(p::GGJParameters, Q::ComplexF64; nps::Int=8, rtol: p1v = p1(p) pplus = -0.5 + p1v - e = ComplexF64(p.E); - f = ComplexF64(p.F); + e = ComplexF64(p.E) + f = ComplexF64(p.F) h = ComplexF64(p.H) - g = ComplexF64(p.G); - k = ComplexF64(p.K); + g = ComplexF64(p.G) + k = ComplexF64(p.K) q = Q q3 = q^3 @@ -158,11 +158,11 @@ end function _build_infinity_arrays(p::GGJParameters, Q::ComplexF64, d0inv::SMatrix{4,4,ComplexF64}; rtol::Float64=1e-6, fmax::Float64=1.0) - e = ComplexF64(p.E); - f = ComplexF64(p.F); + e = ComplexF64(p.E) + f = ComplexF64(p.F) h = ComplexF64(p.H) - g = ComplexF64(p.G); - k = ComplexF64(p.K); + g = ComplexF64(p.G) + k = ComplexF64(p.K) q = Q dr = mercier_dr(p) @@ -343,17 +343,10 @@ end rtol_origin::Float64=1e-6, nps::Int=8, fmax::Float64=1.0, solver=Tsit5()) -> SVector{2,ComplexF64} -Solve the GGJ inner-layer matching problem by stable backward shooting in -the origin-diagonalized 4×4 basis. Direct port of the rmatch `deltar.f` -algorithm. - -Returns the parity-projected matching data `(Δ₁, Δ₂)` (already rescaled -back to physical units via `rescale_delta`). Index ordering matches the -Fortran `deltar` output. - -Tolerances `reltol`/`abstol` are the integrator tolerances; `rtol_origin` -controls the truncation error of the origin Frobenius series and the -choice of `tmin`. +Backward shoot in the origin-diagonalized 4×4 basis; direct port of rmatch +`deltar.f`. Returns `(Δ₁, Δ₂)` rescaled to physical units, index ordering as +the Fortran `deltar` output. `rtol_origin` sets the origin-series truncation +and `tmin`. Not for production use — see the file header. """ function solve_inner(::GGJModel{:shooting}, params::GGJParameters, γ::Number; reltol::Float64=1e-6, abstol::Float64=1e-6, diff --git a/src/InnerLayer/InnerLayer.jl b/src/InnerLayer/InnerLayer.jl index 934e8f83..112843b0 100644 --- a/src/InnerLayer/InnerLayer.jl +++ b/src/InnerLayer/InnerLayer.jl @@ -17,13 +17,17 @@ include("GGJ/GGJ.jl") import .GGJ: GGJModel, GGJParameters, build_asymptotics, evaluate_asymptotics, pick_xmax import .GGJ: InnerAsymptoticsCache, mercier_di, mercier_dr, inner_Q, rescale_delta import .GGJ: glasser_wang_2020_eq55, solve_inner_converged # solve_inner_converged: experimental, not exported (reachable as a qualified call only) +import .GGJ: solve_ray, RaySolveResult, pick_smax, physical_ua_dua +import .GGJ: delta_convergence, solution_profile, asymptotic_profile, q4_surface_benchmark # SLAYER imports go here -export InnerLayerModel, solve_inner +export InnerLayerModel, solve_inner, solve_inner_profile export GGJ, GGJModel, GGJParameters export build_asymptotics, evaluate_asymptotics, pick_xmax, InnerAsymptoticsCache export mercier_di, mercier_dr, inner_Q, rescale_delta export glasser_wang_2020_eq55 +export solve_ray, RaySolveResult, pick_smax, physical_ua_dua +export delta_convergence, solution_profile, asymptotic_profile, q4_surface_benchmark # SLAYER exports go here diff --git a/src/InnerLayer/InnerLayerInterface.jl b/src/InnerLayer/InnerLayerInterface.jl index 3c6e9010..a4d6f6d1 100644 --- a/src/InnerLayer/InnerLayerInterface.jl +++ b/src/InnerLayer/InnerLayerInterface.jl @@ -27,3 +27,27 @@ imposed at the rational surface, X = 0). They are the Δ_{j,±}(γ) of Glasser, Wang & Park, Phys. Plasmas **23**, 112506 (2016), Eqs. (34)–(35). """ function solve_inner end + +""" + solve_inner_profile(model::InnerLayerModel, params, γ::Number; kwargs...) + -> (; Δ, x, Ψ, Ξ, dψdx, rescale, ...) + +Compute the inner-layer matching data **and** the reconstructed layer field +profiles for the given `model` — everything an outer↔inner matching driver +needs from the layer, so drivers never touch model internals. Returns a named +tuple with at least: + + - `Δ` — the same `(Δ_odd, Δ_even)` matching data as [`solve_inner`](@ref) + - `x` — real ascending grid in the model's stretched inner coordinate, + `x ≥ 0` with the rational surface at `x = 0` + - `Ψ`, `Ξ` — `length(x) × 2` profiles, columns (odd, even) parity, in the + model's inner normalization: `Ψ` the normal-field + (reconnected-flux) variable, `Ξ` the displacement + - `dψdx` — conversion to poloidal-flux distance, `δψ = dψdx · x` + - `rescale` — amplitude factor converting the inner-normalized profiles to + the outer δψ-normalized convention (companion of the Δ rescale) + +Concrete models may return additional diagnostic fields (e.g. a solve-quality +certificate). Solver-knob keywords are model-specific. +""" +function solve_inner_profile end diff --git a/src/KineticForces/PitchIntegration.jl b/src/KineticForces/PitchIntegration.jl index fa61c595..845cfabc 100644 --- a/src/KineticForces/PitchIntegration.jl +++ b/src/KineticForces/PitchIntegration.jl @@ -166,7 +166,7 @@ in a single pitch integration, sharing one energy integration per (λ, E). Returns a length-`2*nqty` packed buffer: `[wmm | tmm]`. The two halves each reproduce Fortran's independent-pass result at Fortran's element-by-element -convention (verified via matrix-dump comparison vs `~/Code/gpec/dcon/fourfit.F` +convention (verified via matrix-dump comparison vs Fortran `dcon/fourfit.F` `kwmat_l`/`ktmat_l`). Downstream `kwmat ± ktmat` combinations in `ForceFreeStates/Kinetic.jl` then reproduce `sing.f:967-1075` exactly for the non-Hermitian B_k, C_k, E_k diagonals. diff --git a/src/KineticForces/Torque.jl b/src/KineticForces/Torque.jl index af22d165..f9af8648 100644 --- a/src/KineticForces/Torque.jl +++ b/src/KineticForces/Torque.jl @@ -870,7 +870,7 @@ independent integration passes at `torque.F90:842-847`. Per-surface matrix dumps confirm element-by-element agreement with Fortran `fourfit.F:1080-1082` (`kwmat_l`, `ktmat_l`). This is the Fortran convention required by the adjoint combinations `kwmat ± ktmat` in `ForceFreeStates/Kinetic.jl` / -`~/Code/gpec/dcon/sing.f:967-1075` for non-Hermitian B_k, C_k, E_k. +Fortran `dcon/sing.f:967-1075` for non-Hermitian B_k, C_k, E_k. # Arguments - `kwmat::Array{ComplexF64,3}`: Output (mpert×mpert×6), fwmm half, zeroed on entry diff --git a/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl b/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl index 804c6ac8..3a114270 100644 --- a/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl +++ b/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl @@ -74,6 +74,9 @@ Internal state variables for perturbed equilibrium calculations. singular_coupling_metrics::Dict{String,Float64} = Dict{String,Float64}() m_modes::Vector{Int} = Int[] n_modes::Vector{Int} = Int[] + # ForceFreeStates-provided B_pen per (match surface × coil-drive column) from inner layer. + inner_bpen::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0) + odet_from_gal::Bool = false end """ @@ -165,6 +168,8 @@ well-conditioned flux-space inductances L, Λ: island_width_sq::Vector{ComplexF64} = ComplexF64[] penetrated_area_weighted_field::Vector{ComplexF64} = ComplexF64[] delta_prime::Vector{ComplexF64} = ComplexF64[] + forcing_solution_weights::Vector{ComplexF64} = ComplexF64[] + rational_area::Vector{Float64} = Float64[] # Diagnostics [n_rational] island_half_width::Vector{Float64} = Float64[] diff --git a/src/PerturbedEquilibrium/SingularCoupling.jl b/src/PerturbedEquilibrium/SingularCoupling.jl index d9cd1ebc..d0b81c21 100644 --- a/src/PerturbedEquilibrium/SingularCoupling.jl +++ b/src/PerturbedEquilibrium/SingularCoupling.jl @@ -40,6 +40,20 @@ function _hermite_cubic_val(u_a, u_b, du_a, du_b, psi_a, psi_b, psi) return @. h00 * u_a + h * h10 * du_a + h01 * u_b + h * h11 * du_b end +# Cubic Hermite interpolant DERIVATIVE at psi. Used when the endpoint derivatives are analytic +# (gal-matched OdeState: ud_store carries the exact ξ′ including the resonant Frobenius content), +# so the interpolated derivative is consistent with the stored analytic one — unlike a chord slope, +# which differences singular content across nodes. +function _hermite_cubic_deriv(u_a, u_b, du_a, du_b, psi_a, psi_b, psi) + h = psi_b - psi_a + t = (psi - psi_a) / h + d00 = (6t^2 - 6t) / h + d10 = 3t^2 - 4t + 1 + d01 = (-6t^2 + 6t) / h + d11 = 3t^2 - 2t + return @. d00 * u_a + d10 * du_a + d01 * u_b + d11 * du_b +end + # Reflect a periodic theta-space vector θ → -θ (the theta reversal in gpvacuum_flxsurf). _reverse_theta(v::AbstractVector) = circshift(reverse(v), 1) @@ -127,13 +141,20 @@ function compute_singular_coupling_metrics!( state.C_resonant_area_weighted_field = zeros(ComplexF64, n_rational, numpert_total) state.C_resonant_current = zeros(ComplexF64, n_rational, numpert_total) state.C_island_width_sq = zeros(ComplexF64, n_rational, numpert_total) - state.C_penetrated_area_weighted_field = zeros(ComplexF64, n_rational, numpert_total) + have_inner_bpen = !isempty(intr.inner_bpen) + if have_inner_bpen + state.C_penetrated_area_weighted_field = zeros(ComplexF64, n_rational, numpert_total) + else + state.C_penetrated_area_weighted_field = zeros(ComplexF64, 0, 0) + @warn "No inner-layer B_pen supplied; penetrated field not computed." maxlog=1 + end state.C_delta_prime = zeros(ComplexF64, n_rational, numpert_total) state.rational_psi = zeros(Float64, n_rational) state.rational_q = zeros(Float64, n_rational) state.rational_m_res = zeros(Int, n_rational) state.rational_n = zeros(Int, n_rational) state.rational_surface_idx = zeros(Int, n_rational) + state.rational_area = zeros(Float64, n_rational) # Precompute ODE coefficient matrix C_coeffs for all PE forcing modes. # For each forcing mode k: c_k = u_bnd⁻¹ × edge_mn_k @@ -225,12 +246,20 @@ function compute_singular_coupling_metrics!( u_l = _hermite_cubic_val(ua_l, ub_l, dua_l, dub_l, psi_il_l, psi_ir_l, lpsi) u_r = _hermite_cubic_val(ua_r, ub_r, dua_r, dub_r, psi_il_r, psi_ir_r, rpsi) - # Derivative (ud): chord slope from u_store only — ud_store can be systematically off near - # outer surfaces (q=4/5) where the ODE solution varies rapidly; near-cancellation in bwp1 - # then amplifies even a ~10% ud_store error into a large Delta' error. Chord slope avoids - # this by using only u values, which are accurately stored by the ODE integrator. - ud_l = (ub_l .- ua_l) ./ (psi_ir_l - psi_il_l) - ud_r = (ub_r .- ua_r) ./ (psi_ir_r - psi_il_r) + # Derivative (ud): two paths. + # - gal-matched OdeState (intr.odet_from_gal): ud_store is the analytic ξ′ from the gal basis, + # use the analytic Hermite-cubic derivative built from (u, ud). + # - shooting OdeState: chord slope from u_store only — ud_store can be systematically off near + # outer surfaces (q=4/5) where the ODE solution varies rapidly; near-cancellation in bwp1 + # then amplifies even a ~10% ud_store error into a large Delta' error. Chord slope avoids + # this by using only u values, which are accurately stored by the ODE integrator. + if intr.odet_from_gal + ud_l = _hermite_cubic_deriv(ua_l, ub_l, dua_l, dub_l, psi_il_l, psi_ir_l, lpsi) + ud_r = _hermite_cubic_deriv(ua_r, ub_r, dua_r, dub_r, psi_il_r, psi_ir_r, rpsi) + else + ud_l = (ub_l .- ua_l) ./ (psi_ir_l - psi_il_l) + ud_r = (ub_r .- ua_r) ./ (psi_ir_r - psi_il_r) + end q_l = equil.profiles.q_spline(lpsi) q1_l = equil.profiles.q_deriv(lpsi) @@ -249,14 +278,15 @@ function compute_singular_coupling_metrics!( bwp1_l = 2π * im * chi1 * (singfac_l * xsp1_l - nn * q1_l * xsp_l) bwp1_r = 2π * im * chi1 * (singfac_r * xsp1_r - nn * q1_r * xsp_r) jump_vec[k] = bwp1_r - bwp1_l - # C_penetrated_area_weighted_field: midpoint of b^ψ at lpsi/rpsi divided by the scalar surface area. - # Matches Fortran gpout_resp: gpeq_interp_singsurf evaluates bwp_mn at respsi. - # LHS normalization audit (#233): the resonant flux Φ^r divided by the scalar area A^r - # is a genuine field amplitude in tesla and is coordinate-invariant [Pharr 2026; cf. - # the resonant-field definition in the Conventions Reference]. - b_l = chi1 * singfac_l * 2π * im * xsp_l - b_r = chi1 * singfac_r * 2π * im * xsp_r - state.C_penetrated_area_weighted_field[row, k] = (b_l + b_r) / 2 / area + end + + # Inner-layer (cusp-free) penetrated field: bpen[s, j] is linear in the same identity-at-edge + # coil-drive columns as the OdeState solutions, so it contracts with C_coeffs exactly like + # the outer solution values above (xsp = dot(u, ck)); /area matches the area-weighted + # convention of the pointwise row. + if have_inner_bpen && s <= size(intr.inner_bpen, 1) + pen_row = (transpose(C_coeffs) * @view(intr.inner_bpen[s, :])) ./ area + state.C_penetrated_area_weighted_field[row, :] = pen_row end # LHS normalization audit (#233) — output scalar coordinate-invariance per row: @@ -276,6 +306,7 @@ function compute_singular_coupling_metrics!( end state.rational_psi[row] = sing_surf.psifac + state.rational_area[row] = area state.rational_q[row] = sing_surf.q state.rational_m_res[row] = m_res state.rational_n[row] = nn @@ -302,8 +333,9 @@ function compute_singular_coupling_metrics!( state.resonant_area_weighted_field = state.C_resonant_area_weighted_field * forcing_flux state.resonant_current = state.C_resonant_current * forcing_flux state.island_width_sq = state.C_island_width_sq * forcing_flux - state.penetrated_area_weighted_field = state.C_penetrated_area_weighted_field * forcing_flux state.delta_prime = state.C_delta_prime * forcing_flux + have_inner_bpen && (state.penetrated_area_weighted_field = state.C_penetrated_area_weighted_field * forcing_flux) + state.forcing_solution_weights = C_coeffs * forcing_flux # Conform the stored coupling-matrix input basis to the coordinate-invariant root-area-weighted # field (b̃) space (#233 / Pharr 2026): C̃ = C·R, so each stored row acts on the applied field @@ -316,8 +348,8 @@ function compute_singular_coupling_metrics!( state.C_resonant_area_weighted_field = state.C_resonant_area_weighted_field * flux_conform state.C_resonant_current = state.C_resonant_current * flux_conform state.C_island_width_sq = state.C_island_width_sq * flux_conform - state.C_penetrated_area_weighted_field = state.C_penetrated_area_weighted_field * flux_conform state.C_delta_prime = state.C_delta_prime * flux_conform + have_inner_bpen && (state.C_penetrated_area_weighted_field = state.C_penetrated_area_weighted_field * flux_conform) # Phase 5: Island diagnostics from applied resonant vectors compute_island_diagnostics!(state, n_rational) diff --git a/src/PerturbedEquilibrium/Utils.jl b/src/PerturbedEquilibrium/Utils.jl index 9957fe7c..a7f06db8 100644 --- a/src/PerturbedEquilibrium/Utils.jl +++ b/src/PerturbedEquilibrium/Utils.jl @@ -61,6 +61,7 @@ perturbed_equilibrium/ ├── forcing_b / forcing_b_root_area / forcing_b_area # control-surface forcing spectrum (b, b̃, b̄) [numpert_total], tesla ├── response_b / response_b_root_area / response_b_area # control-surface response spectrum (b, b̃, b̄) [numpert_total], tesla ├── response/ +│ ├── psi_n # Radial abscissa ψ_N [npsi] shared by every response profile below │ ├── xi_psi # Radial displacement ξ^ψ = ξ·∇ψ (ComplexF64 [npsi, mpert]) │ ├── xi_psi_J # J·ξ^ψ Jacobian-weighted (from gpeq_contra) │ ├── b_psi_area_weighted # b^ψ / ⟨J·|∇ψ|⟩_θ area-normalized (ComplexF64 [npsi, mpert]) @@ -136,6 +137,7 @@ function write_outputs_to_HDF5( # Response fields (ComplexF64 directly) response_group = haskey(pe_group, "response") ? pe_group["response"] : create_group(pe_group, "response") + !isempty(state.psi_grid) && (response_group["psi_n"] = state.psi_grid) have_xi = !isnothing(state.xi_modes) have_b = have_xi && !isnothing(state.b_modes) response_group["xi_psi"] = have_xi ? state.xi_modes.psi : ComplexF64[] @@ -203,8 +205,10 @@ function write_outputs_to_HDF5( !isempty(state.resonant_area_weighted_field) && (coupling_group["resonant_area_weighted_field"] = state.resonant_area_weighted_field) !isempty(state.resonant_current) && (coupling_group["resonant_current"] = state.resonant_current) !isempty(state.island_width_sq) && (coupling_group["island_width_sq"] = state.island_width_sq) - !isempty(state.penetrated_area_weighted_field) && (coupling_group["penetrated_area_weighted_field"] = state.penetrated_area_weighted_field) + !isempty(state.penetrated_area_weighted_field) && (coupling_group["penetrated_area_weighted_field"] = state.penetrated_area_weighted_field) !isempty(state.delta_prime) && (coupling_group["delta_prime"] = state.delta_prime) + !isempty(state.forcing_solution_weights) && (coupling_group["forcing_solution_weights"] = state.forcing_solution_weights) + !isempty(state.rational_area) && (coupling_group["rational_area"] = state.rational_area) !isempty(state.island_half_width) && (coupling_group["island_half_width"] = state.island_half_width) !isempty(state.chirikov_parameter) && (coupling_group["chirikov_parameter"] = state.chirikov_parameter) diff --git a/test/runtests_eulerlagrange.jl b/test/runtests_eulerlagrange.jl index 6f3ebc54..beb11049 100644 --- a/test/runtests_eulerlagrange.jl +++ b/test/runtests_eulerlagrange.jl @@ -1,6 +1,5 @@ -# TODO: perhaps this isn't the best place for this function? -# Should I do include("../ForceFreeStates/utils.jl") instead? or maybe save these functions in a separate file? -# associated TODO: come up with Gaussian reduction test that doesn't rely on external data +# TODO: this helper may belong in a shared test-utilities file rather than here. +# TODO: come up with a Gaussian reduction test that doesn't rely on external data. function load_u_matrix(filename) lines = readlines(filename) diff --git a/test/runtests_innerlayer.jl b/test/runtests_innerlayer.jl index 7ee36d06..7b087216 100644 --- a/test/runtests_innerlayer.jl +++ b/test/runtests_innerlayer.jl @@ -49,3 +49,118 @@ const GGJ = IL.GGJ @test abs(imag(Δ[2])) < 1e-3 * abs(Δ[2]) end end + +@testset "InnerLayer GGJ :ray backend (rotated-contour collocation)" begin + # The method was validated before landing here: manufactured Δ* to 3e-14, + # Fortran rmatch pins, 96-equilibrium robustness scans to Q = 500i. These + # tests pin the implementation, not the method. + p = IL.glasser_wang_2020_eq55() + + @testset "agrees with :galerkin at the paper point Q = 0.1234" begin + γ = 0.1234 * GGJ.q0(p) + Δ = IL.solve_inner(IL.GGJModel(; solver=:ray), p, γ) + # Same Fortran-cross-checked pins as the Galerkin testset above. + @test real(Δ[1]) ≈ 3.698368e4 rtol = 1e-3 + @test real(Δ[2]) ≈ 14.759721 rtol = 1e-3 + @test abs(imag(Δ[1])) < 1e-3 * abs(Δ[1]) + @test abs(imag(Δ[2])) < 1e-3 * abs(Δ[2]) + # :ray is the default backend. + @test IL.GGJModel() === IL.GGJModel{:ray}() + end + + @testset "q4 physical benchmark at Q = 500i (regime beyond :galerkin)" begin + q4 = IL.q4_surface_benchmark() + γ = 500.0im * GGJ.q0(q4) + Δ = IL.solve_inner(IL.GGJModel(), q4, γ) + # Pins from the pre-port validation suite (post extended-precision + # march fix; S-invariant to 3e-4 / 7e-9 and θ-stable there). + @test Δ[1] ≈ 2.4720608737 + 13.354123514im rtol = 1e-4 + @test Δ[2] ≈ 0.13749694953 + 0.74275468725im rtol = 1e-4 + + # Δ is an analytic invariant of the contour angle: the outward + # θ-check drift is a direct numerical error measurement. + Q = GGJ.inner_Q(q4, γ) + r2 = IL.solve_ray(q4, Q; θ=1.2 * angle(Q) / 4) + @test abs(r2.Δ[2] - Δ[2]) / abs(Δ[2]) < 1e-5 + @test abs(r2.Δ[1] - Δ[1]) / abs(Δ[1]) < 1e-3 + end +end + +@testset "InnerLayer GGJ :ray internal machinery" begin + p = IL.q4_surface_benchmark() + + @testset "cheblobatto nodes and differentiation matrix" begin + t, D = GGJ.cheblobatto(8) + @test length(t) == 9 + @test issorted(t) # ascending, per the reflected convention + @test t[1] ≈ -1 && t[end] ≈ 1 + @test D * ones(9) ≈ zeros(9) atol = 1e-10 # d/dt of a constant is zero + @test D * t ≈ ones(9) atol = 1e-10 # d/dt of the linear function is one + end + + @testset "ode_matrix: ordinary point and type-generic build" begin + Q = 5.0im + # x = 0 is an ordinary point: the coefficient matrix is finite there. + M0 = GGJ.ode_matrix(p, Q, 0.0) + @test all(isfinite, M0) + @test M0[1, 4] == 1 && M0[2, 5] == 1 && M0[3, 6] == 1 # v' = (Ψ',Ξ',Υ') block + # The extended-precision build agrees with the Float64 build. + Md = GGJ.ode_matrix(Complex{GGJ.Double64}, p, Q, 0.3) + @test ComplexF64.(Md) ≈ GGJ.ode_matrix(p, Q, 0.3) rtol = 1e-12 + end + + @testset "parity_rows match the deltac boundary convention" begin + @test GGJ.parity_rows(1) == [4, 2, 3] # odd: Ψ'(0)=Ξ(0)=Υ(0)=0 + @test GGJ.parity_rows(2) == [1, 5, 6] # even: Ψ(0)=Ξ'(0)=Υ'(0)=0 + end + + @testset "decaying_pair is an orthonormal 6×2 frame" begin + Q = 5.0im + θ = angle(Q) / 4 + E = GGJ.decaying_pair(p, Q, θ, 60.0) + @test size(E) == (6, 2) + @test all(isfinite, E) + @test E' * E ≈ I(2) atol = 1e-10 # columns orthonormal + end + + @testset "profile diagnostics reconstruct finite fields" begin + res = IL.solve_ray(p, 5.0im) + prof = IL.solution_profile(res; npc=6) + @test size(prof.Ψ, 2) == 2 && length(prof.s) == size(prof.Ψ, 1) + @test all(isfinite, prof.Ψ) && all(isfinite, prof.Ξ) && all(isfinite, prof.Υ) + # Analytic tail evaluates on the trusted series radius. + asy = IL.asymptotic_profile(p, res, [res.S, 2 * res.S]) + @test all(isfinite, asy.Ψ) + end + + @testset "delta_convergence: small spread, consistent with solve_inner" begin + Q = 5.0im + conv = IL.delta_convergence(p, Q; verbose=false) + Δ = IL.solve_inner(IL.GGJModel(; solver=:ray), p, Q * GGJ.q0(p)) + @test conv.Δ ≈ Δ rtol = 1e-6 # baseline == the plain solve + @test maximum(conv.spread) < 1e-4 # honest error bar is small here + end +end + +@testset "solve_inner_profile interface (matching-driver contract)" begin + p = IL.glasser_wang_2020_eq55() + γ = 0.1234 * GGJ.q0(p) + for model in (IL.GGJModel(; solver=:ray), IL.GGJModel(; solver=:galerkin)) + prof = IL.solve_inner_profile(model, p, γ) + # Δ agrees with the plain matching solve of the same backend (identical solve path). + @test prof.Δ ≈ IL.solve_inner(model, p, γ) rtol = 1e-12 + # Real ascending inner-coordinate grid from the rational surface, profiles npts × 2. + @test issorted(prof.x) + @test prof.x[1] ≈ 0 atol = 1e-12 + @test size(prof.Ψ) == (length(prof.x), 2) && size(prof.Ξ) == (length(prof.x), 2) + @test all(isfinite, prof.Ψ) && all(isfinite, prof.Ξ) + # Parity at the layer center: Ψ(0) ≠ 0 odd-parity column, Ψ(0) = 0 even-parity column. + @test abs(prof.Ψ[1, 2]) < 1e-6 * abs(prof.Ψ[1, 1]) + # Conversion factors match their GGJ definitions. + @test prof.dψdx ≈ GGJ.x0(p) / p.v1 + @test prof.rescale ≈ (p.v1 / GGJ.x0(p))^(0.5 + GGJ.p1(p)) + end + # Ray backend certificate: at real Q the optimal contour is θ = 0, so the two solves coincide. + ray = IL.solve_inner_profile(IL.GGJModel(; solver=:ray), p, γ) + @test ray.certΔ < 1e-12 +end