diff --git a/docs/figure_tools.jl b/docs/figure_tools.jl index ca8001be5..13b91d08f 100644 --- a/docs/figure_tools.jl +++ b/docs/figure_tools.jl @@ -29,13 +29,16 @@ 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. +provenance than it has. Generated figures and the manifest do not count as changes, so +several figures can be regenerated in turn on one clean commit. 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))) + h = String(strip(read(`git -C $(_REPO_ROOT) rev-parse --short HEAD`, String))) + outputs = [":(exclude,glob)docs/src/figures/**/*.png", ":(exclude)docs/src/figures/manifest.toml"] + dirty = !isempty(strip(read(`git -C $(_REPO_ROOT) status --porcelain -- . $outputs`, String))) return (; commit=(dirty ? "$h-dirty" : h), date=date, dirty=dirty) catch return (; commit="unknown", date=date, dirty=false) diff --git a/docs/src/citations.md b/docs/src/citations.md index dd8ef7546..a447cfd34 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) -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. +Constructs the Wasow large-``x`` asymptotic basis (the `inps` kernel) that supplies the far-field boundary condition for every GGJ InnerLayer backend — the traditional Galerkin solver and the rotated-contour (`:ray`) solver alike. --- diff --git a/docs/src/figures/inner_layer/convergence_Sinvariance.png b/docs/src/figures/inner_layer/convergence_Sinvariance.png index d8b550822..cedbc5905 100644 Binary files a/docs/src/figures/inner_layer/convergence_Sinvariance.png 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 index 163b50493..f0978ab30 100644 --- a/docs/src/figures/inner_layer/make_backend_regime_map.jl +++ b/docs/src/figures/inner_layer/make_backend_regime_map.jl @@ -19,7 +19,7 @@ using Printf const IL = GeneralizedPerturbedEquilibrium.InnerLayer const GGJ = IL.GGJ -p = IL.q4_surface_benchmark() +p = IL.GGJ.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. diff --git a/docs/src/figures/inner_layer/make_convergence_Sinvariance.jl b/docs/src/figures/inner_layer/make_convergence_Sinvariance.jl index 3e2a71b2c..2b4129baf 100644 --- a/docs/src/figures/inner_layer/make_convergence_Sinvariance.jl +++ b/docs/src/figures/inner_layer/make_convergence_Sinvariance.jl @@ -2,10 +2,10 @@ # # 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). Δ +# an independent axis (contour angle θ, mesh resolution, series order/radius, +# 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. +# these orthogonal perturbations bounds the true error. Seven solves; ~30 s. # # Run manually: julia --project=. docs/src/figures/inner_layer/make_convergence_Sinvariance.jl @@ -16,13 +16,13 @@ using Printf const IL = GeneralizedPerturbedEquilibrium.InnerLayer -p = IL.q4_surface_benchmark() +p = IL.GGJ.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 +d1 = [max(r.d1, 1e-16) for r in conv.table] # relative change of Δ₊ (interchange) +d2 = [max(r.d2, 1e-16) for r in conv.table] # relative change of Δ₋ (tearing) n = length(names) plt = plot(; size=(820, 540), legend=:topright, yscale=:log10, @@ -31,14 +31,14 @@ plt = plot(; size=(820, 540), legend=:topright, yscale=:log10, 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") +scatter!(plt, 1:n, d1; marker=:circle, ms=7, color=1, label="δΔ₊ (interchange)") +scatter!(plt, 1:n, d2; marker=:diamond, ms=7, color=2, label="δΔ₋ (tearing)") # 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])) + label=@sprintf("worst-case Δ₊ %.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])) + label=@sprintf("worst-case Δ₋ %.1e", conv.spread[2])) save_doc_figure(plt, "inner_layer", "convergence_Sinvariance"; script=basename(@__FILE__), diff --git a/docs/src/figures/inner_layer/make_rotated_ray_contour.jl b/docs/src/figures/inner_layer/make_rotated_ray_contour.jl index 15abe1670..a23c39fba 100644 --- a/docs/src/figures/inner_layer/make_rotated_ray_contour.jl +++ b/docs/src/figures/inner_layer/make_rotated_ray_contour.jl @@ -19,7 +19,7 @@ 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() +p = IL.GGJ.q4_surface_benchmark() Q = 500.0im θ = angle(Q) / 4 # = π/8 = 22.5° diff --git a/docs/src/figures/inner_layer/make_solution_profiles.jl b/docs/src/figures/inner_layer/make_solution_profiles.jl index 2165f358e..6f89b9c66 100644 --- a/docs/src/figures/inner_layer/make_solution_profiles.jl +++ b/docs/src/figures/inner_layer/make_solution_profiles.jl @@ -1,7 +1,7 @@ # 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 q=4 benchmark surface. The Hermite BVP 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. @@ -16,19 +16,19 @@ using Printf const IL = GeneralizedPerturbedEquilibrium.InnerLayer const GGJ = IL.GGJ -p = IL.q4_surface_benchmark() +p = IL.GGJ.q4_surface_benchmark() Q = 2.0im # imaginary axis, beyond the |Q|≪1 shooting regime, -# small enough that the collocation domain reaches +# small enough that the BVP 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) +isol = 2 # odd-parity (tearing, Δ₋) 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) +@printf("solve_ray: Q=%s θ=%.3f S=%.4g s_m=%.4g cells=%d\n", + string(Q), res.θ, res.S, s_m, length(res.breaks) - 1) -# Collocation solution over the BVP domain [0, s_m]. +# BVP solution over [0, s_m]. prof = IL.solution_profile(res; npc=8) # Analytic tail on [S, few·S] where the inps series is trusted. @@ -39,13 +39,13 @@ fields = ((:Ψ, prof.Ψ, asy.Ψ, 1), (:Ξ, prof.Ξ, asy.Ξ, 2), (:Υ, prof.Υ, a 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)", + xlabel="ray parameter s (x = e^{iθ}s)", ylabel="|field| (odd parity, tearing)", 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") + label="|$(nm)| BVP") plot!(plt, asy.s, abs.(acol[:, isol]); lw=2, ls=:dash, color=ci, label="|$(nm)| asymptotic") end @@ -56,5 +56,4 @@ vline!(plt, [res.S]; color=:black, ls=:dot, lw=1, 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"]) + depends=["src/InnerLayer/GGJ/Ray.jl", "src/InnerLayer/GGJ/InnerAsymptotics.jl"]) diff --git a/docs/src/figures/inner_layer/solution_profiles.png b/docs/src/figures/inner_layer/solution_profiles.png index 3a1c2ef57..de394ef0d 100644 Binary files a/docs/src/figures/inner_layer/solution_profiles.png 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 index 2ad9ef5dc..731585bc5 100644 --- a/docs/src/figures/manifest.toml +++ b/docs/src/figures/manifest.toml @@ -25,8 +25,8 @@ depends = ["src/InnerLayer/GGJ/Ray.jl", "src/InnerLayer/GGJ/Galerkin.jl", "src/I script = "make_backend_regime_map.jl" [inner_layer.convergence_Sinvariance] -commit = "3a0837e3-dirty" -date = "2026-07-08" +commit = "5d11aad3e" +date = "2026-10-21" depends = ["src/InnerLayer/GGJ/Ray.jl"] script = "make_convergence_Sinvariance.jl" @@ -37,7 +37,7 @@ 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"] +commit = "5d11aad3e" +date = "2026-10-21" +depends = ["src/InnerLayer/GGJ/Ray.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 7ddfecb95..ff3033c12 100644 --- a/docs/src/inner_layer.md +++ b/docs/src/inner_layer.md @@ -19,8 +19,8 @@ concrete models: the **GGJ (Glasser–Greene–Johnson)** layer documented below interchangeable solver backends, and the pressureless **SLAYER** Fitzpatrick Riccati layer. Further inner-layer models (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: +Two papers define the equations and the asymptotic construction; the code and this page cite +them by equation number rather than restating them: - **GWP2016** — A. H. Glasser, Z. R. Wang and J.-K. Park, *Computation of resistive instabilities by matched asymptotic expansions*, Phys. Plasmas @@ -29,192 +29,108 @@ cited by their equation numbers throughout the code and below: 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). +## Equations and matching data + +The GGJ layer equations (GWP2016 Eq. 11 ≡ GW2020 Eq. 1) couple the perturbed flux ``\Psi`` and +the layer variables ``\Xi``, ``\Upsilon`` in the scaled distance ``x`` from the rational surface. +Their equilibrium coefficients ``E, F, G, H, K`` (GWP2016 Eq. A8) are held in +[`GGJParameters`](@ref InnerLayer.GGJParameters). A physical growth rate ``\gamma`` maps to the +scaled ``Q`` by [`inner_Q`](@ref InnerLayer.inner_Q) (GWP2016 Eqs. A11, A14–A16), and +[`rescale_delta`](@ref InnerLayer.rescale_delta) returns the matching data to physical units. The +model requires Mercier stability, ``D_I < 0`` ([`mercier_di`](@ref InnerLayer.mercier_di), +GWP2016 Eq. A9), so that at large ``x`` the solutions include two power laws (GW2020 Eq. 49). + +The matching datum ``\Delta`` is the ratio of the large to the small power-law coefficient. +Imposing the two parities at ``x = 0`` gives the pair ``(\Delta_+, \Delta_-)`` of GWP2016 +Eqs. (34)–(35). Parity is that of the displacement ``\Xi`` (GWP2016 Eq. 34; Pletzer, Bondeson & +Dewar, J. Comput. Phys. **115**, 530 (1994), Eqs. 32–33): ``\Delta_+`` is even, +``\Psi(0)=\Xi'(0)=\Upsilon'(0)=0``, the **interchange** channel; ``\Delta_-`` is odd, +``\Psi'(0)=\Xi(0)=\Upsilon(0)=0``, the **tearing** channel. Positional pairs, such as `deltar` in +the matching and `RaySolveResult.Δ`, are always ordered ``(\Delta_+, \Delta_-)``, like the +unknowns ``d_\pm`` of GWP2016 Eq. (38) and Pletzer, Bondeson & Dewar Eqs. 83, 85; +`InnerLayerResponse` names them `interchange` and `tearing`. ## 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. + The figures use the ``q = 4`` DIII-D-like surface of + [`q4_surface_benchmark`](@ref InnerLayer.GGJ.q4_surface_benchmark) + (``S \approx 4.6\times10^{6}``, ``D_I \approx -0.31``). ### 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 +[`GGJModel`](@ref InnerLayer.GGJModel); all use 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) +| `:galerkin` | Hermite-cubic finite-element (real axis) | surface- and ``\arg Q``-dependent; see [Choosing a backend](@ref) | +| `:ray` (default) | rotated contour: inward-marched far-field condition + Hermite-cubic BVP | ``\lvert Q\rvert \sim 500`` on/near the imaginary axis | -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. +Rotation and resistivity scans live on the imaginary-``Q`` axis, where the pseudo-resonance at +``x^2 = -Q^2(G+KF)`` sits on the real-``x`` path. On the ``q = 4`` surface there `:shooting` and +`:galerkin` lose accuracy beyond ``|Q| \sim 1`` and ``|Q| \sim 4`` respectively (the poles of ``\Delta(Q)`` +themselves lie near the real-``Q`` axis): -### 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. +![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 rotated ray -The equations are continued onto the ray +`:ray` writes the layer equations as ``dv/dx = M(x)\,v`` in the plain state +``v = (\Psi, \Xi, \Upsilon, \Psi', \Xi', \Upsilon')``. ``M`` is polynomial in ``x``, so +``x = 0`` is an ordinary point and the system continues analytically onto the ray ```math -x = e^{i\theta}\, s, \qquad s \in [0, S], \qquad \theta = \tfrac14\arg Q, +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. +This angle makes the parabolic-cylinder exponent of the outer solutions real and lifts the +pseudo-resonance off the contour; ``\theta = 0`` is 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 +### Hermite-Galerkin 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, +On ``[0, s_m]`` the system is solved with the `:galerkin` backend's Hermite-cubic assembly +(GWP2016 Eq. 32) written on the ray, ``A u_{ss} + e^{i\theta}B u_s + e^{2i\theta}C u = 0``, +with the parity conditions at ``s = 0``. The far end carries the point condition ```math -v(S) - \Delta\, U_b - c_1 E_1 - c_2 E_2 = U_s , +v(s_m) = U_s + \Delta\, U_b + c_1 E_1 + c_2 E_2 , ``` -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. +where ``U_s, U_b`` are the small and large power-like solutions and ``E_{1,2}`` the +forward-decaying exponential pair; ``\Delta, c_1, c_2`` replace the last node's six unknowns. +Each cell resolves the local WKB rate, with a minimum cell count for small ``|Q|``. ### 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) +The boundary data come from the same `inps` Wasow basis as the other backends (GW2020 +Eqs. 3–53), evaluated at the radius ``S`` where its residual is below tolerance +([`pick_smax`](@ref InnerLayer.pick_smax)). At large ``|Q|``, ``S`` lies far beyond ``s_m``, and +the power pair is marched inward by 2-stage Radau IIA modulo the decaying pair, the subspace in +which ``\Delta`` is defined. The decaying pair grows under backward integration at a rate +``\sim\rho s``. Radau's stability function decays as ``2/z``, so steps with every growing mode at +``z > 8`` damp those directions while the slowly varying power pair stays resolved; the step count +grows only like ``\ln(S/s_m)``. Nearer ``s_m`` the march resolves the growing directions, carries +the decaying frame, and projects it out periodically. The damped zone runs in `Complex{Double64}` +because its ill-conditioned implicit solves mix the power pair at a level ``\sim\epsilon\rho S^2``. + +On the real axis with imaginary ``Q`` the switch to resolved steps lands on the +pseudo-resonance, a turning point where the decaying directions cannot be identified; the +default ``\theta = \arg(Q)/4`` keeps it off the contour. + +For matching, `solve_inner_profile` returns ``\Delta`` and the layer-center fields +``\Psi(0), \Xi(0)`` from the rotated solve; ``x = 0`` lies on every contour, so these are the +real-axis values. Fields across the layer on the real ``\psi`` axis are an opt-in diagnostic +(`GGJModel(; real_profiles=true)`, not exposed in the deck): they come from an extra ``\theta = 0`` +solve, while the ``x = 0`` values, and so the penetrated field, stay those of the rotated solve. +The ``\theta = 0`` solve's ``\Delta_\pm`` are checked against the rotated ones (`certΔ`, which +warns on disagreement, as can happen at large ``|Q|`` near the imaginary axis). + +![Inner-layer fields Ψ, Ξ, Υ on the rotated ray for the q=4 surface at Q = 2i. The BVP solution (solid) joins the asymptotic representation (dashed) seamlessly at the match point S = s_m.](figures/inner_layer/solution_profiles.png) ## Validation and benchmarks @@ -224,58 +140,65 @@ 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`` | +| ``\Delta_+`` (interchange) | ``3.698368\times 10^{4}`` | ``3.69789\times 10^{4}`` | +| ``\Delta_-`` (tearing) | ``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. +This point sits near a pole of ``\Delta_+(Q)``, which amplifies every solver's error; hence the +``1.3\times10^{-4}`` `:galerkin`↔`:ray` difference in ``\Delta_+`` against ``4\times10^{-7}`` in ``\Delta_-``. **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``, +benchmark ([`q4_surface_benchmark`](@ref InnerLayer.GGJ.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`` | +| ``\Delta_+`` (interchange) | ``2.4720 + 13.3540\,i`` | +| ``\Delta_-`` (tearing) | ``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. +**Convergence and contour invariance.** ``\Delta`` does not depend on the contour angle, so +re-solving with each numerical knob perturbed in turn +([`delta_convergence`](@ref InnerLayer.delta_convergence)) gives an error bar: at ``Q = 500i`` the +worst case is ``\sim 5\times10^{-6}`` for ``\Delta_+`` and ``\sim 6\times10^{-7}`` for ``\Delta_-``. +Across machines and BLAS builds the values reproduce to ``\sim 10^{-5}`` relative, hence the five +significant figures above. -![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) +![Relative change of Δ at Q = 500i under independent perturbations of each numerical knob (contour angle, mesh resolution, series order/radius, march tolerance, BVP edge s_m, decaying-pair separation). 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`. +`:ray` is the **default** — `GGJModel()` constructs `GGJModel{:ray}()` — and the recommended +backend for every ``Q``. It uses the same Hermite-cubic weak form as `:galerkin`; they differ in +where and how the far-field condition is imposed. `:galerkin` imposes it on the real axis through +resonant and extension cells at the radius where the series residual reaches ``10^{-7}``, on a +fixed `nx`-cell grid. `:ray` takes the series data where the residual reaches `smax_tol` = +``10^{-9}`` ([`pick_smax`](@ref InnerLayer.pick_smax)), marches them inward to ``s_m``, where every +exponential family has decayed 45 e-folds from the origin, and resolves the local WKB rate in +every cell on ``[0, s_m]``. Both radii are computed per solve; `S`, `smax_tol` and `sm_fac` +override them, and [`delta_convergence`](@ref InnerLayer.delta_convergence) perturbs them to give +an error bar. + +`:galerkin` takes a few milliseconds per solve, 3–60 times less than `:ray`, but its accuracy +depends on the surface and on ``\arg Q`` even at small ``|Q|``. Against a tightened `:ray` +reference, on the GW2020 Eq. (55) surface it errs by ``7\times10^{-5}`` at ``Q = 0.1`` and +``3\times10^{-3}`` at ``Q = 0.3``; near ``\arg Q = 157.5^\circ`` it errs by order unity at +``|Q| = 0.1`` on that surface and by 0.4 at ``|Q| = 0.3`` on the ``q = 4`` surface. `:ray` held +``2\times10^{-6}`` on all three surfaces tested, for ``|Q|`` from 0.01 to 1. + +The backends take disjoint numerical-knob keywords: pass `GGJModel(solver=:galerkin)` explicitly +when supplying the Galerkin `nx`/`xfac` knobs; passing them to `:ray` throws a `MethodError`. ```julia using GeneralizedPerturbedEquilibrium.InnerLayer -p = q4_surface_benchmark() # GGJParameters for the q=4 benchmark surface +p = GGJ.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 +Δ = solve_inner(GGJModel(), p, γ) # Δ.interchange, Δ.tearing with the default :ray backend -# Full result (raw Δ, contour, mesh, nodal solution, residuals) via solve_ray: +# Full result (raw Δ, contour, mesh, nodal solution, boundary conditioning) via solve_ray: res = solve_ray(p, GGJ.inner_Q(p, γ)) -res.Δ, res.resid, res.bc_cond +res.Δ, res.breaks, res.bc_cond # res.Δ = (Δ₊, Δ₋) = (interchange, tearing) ``` ## API Reference diff --git a/examples/LAR_resistive_match_test/gpec.toml b/examples/LAR_resistive_match_test/gpec.toml index 98d9d56a2..1e2bb7ed6 100644 --- a/examples/LAR_resistive_match_test/gpec.toml +++ b/examples/LAR_resistive_match_test/gpec.toml @@ -77,7 +77,7 @@ gal_rpec_flag = true # Append mpert coil-response columns (unit edge s # 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_inner_solver = "ray" # Inner-layer Delta backend: "ray" (rotated contour, 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 diff --git a/regression-harness/cases/ggj_ray_q500i.toml b/regression-harness/cases/ggj_ray_q500i.toml index f5722ee4d..0ccb2ac23 100644 --- a/regression-harness/cases/ggj_ray_q500i.toml +++ b/regression-harness/cases/ggj_ray_q500i.toml @@ -1,7 +1,7 @@ -# Regression case: GGJ inner-layer rotated-ray collocation backend at Q = 500i. +# Regression case: GGJ inner-layer rotated-ray 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 +# :galerkin backend — and tracks the interchange (Δ₊) / tearing (Δ₋) 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. @@ -10,40 +10,40 @@ 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₀). +# interchange (Δ₊) / tearing (Δ₋) 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" +[quantities.delta_interchange_real] +h5path = "ggj/delta_interchange_real" type = "real_scalar" extract = "value" -label = "Re(Δ_odd)" +label = "Re(Δ₊ interchange)" noise_threshold = 1e-4 order = 10 -[quantities.delta_odd_imag] -h5path = "ggj/delta_odd_imag" +[quantities.delta_interchange_imag] +h5path = "ggj/delta_interchange_imag" type = "real_scalar" extract = "value" -label = "Im(Δ_odd)" +label = "Im(Δ₊ interchange)" noise_threshold = 1e-4 order = 11 -[quantities.delta_even_real] -h5path = "ggj/delta_even_real" +[quantities.delta_tearing_real] +h5path = "ggj/delta_tearing_real" type = "real_scalar" extract = "value" -label = "Re(Δ_even)" +label = "Re(Δ₋ tearing)" noise_threshold = 1e-4 order = 12 -[quantities.delta_even_imag] -h5path = "ggj/delta_even_imag" +[quantities.delta_tearing_imag] +h5path = "ggj/delta_tearing_imag" type = "real_scalar" extract = "value" -label = "Im(Δ_even)" +label = "Im(Δ₋ tearing)" noise_threshold = 1e-4 order = 13 diff --git a/regression-harness/cases/ggj_reference.toml b/regression-harness/cases/ggj_reference.toml index b35ca20a7..bcf1b76e3 100644 --- a/regression-harness/cases/ggj_reference.toml +++ b/regression-harness/cases/ggj_reference.toml @@ -1,42 +1,42 @@ # Regression case: GGJ inner-layer Galerkin solver reference (Glasser & Wang 2020, Eq. 55). # A "computed" case (kind = "computed", no example_dir) that calls solve_inner directly and -# tracks the Δ_odd / Δ_even inner-layer coefficients. Each [quantities.*] block names an +# tracks the interchange (Δ₊) / tearing (Δ₋) inner-layer coefficients. Each [quantities.*] block names an # HDF5 path in the computed output, how to extract it, and the noise floor for diffs. [case] name = "ggj_reference" -description = "GGJ inner-layer Galerkin solver, Glasser & Wang 2020 Eq. 55, γ = 1+i" +description = "GGJ inner-layer Galerkin solver, Glasser & Wang 2020 Eq. 55 with physical τ_A, τ_R, 1 kHz rotation (Q ≈ 0.22i)" kind = "computed" -# Δ_odd / Δ_even from solve_inner(GGJModel(:galerkin), glasser_wang_2020_eq55(), 1+1im) -[quantities.delta_odd_real] -h5path = "ggj/delta_odd_real" +# interchange (Δ₊) / tearing (Δ₋) from solve_inner(GGJModel(:galerkin), glasser_wang_2020_eq55(; taua=2.11226e-7, taur=0.968219), 2πi·10³) +[quantities.delta_interchange_real] +h5path = "ggj/delta_interchange_real" type = "real_scalar" extract = "value" -label = "Re(Δ_odd)" +label = "Re(Δ₊ interchange)" noise_threshold = 1e-10 order = 10 -[quantities.delta_odd_imag] -h5path = "ggj/delta_odd_imag" +[quantities.delta_interchange_imag] +h5path = "ggj/delta_interchange_imag" type = "real_scalar" extract = "value" -label = "Im(Δ_odd)" +label = "Im(Δ₊ interchange)" noise_threshold = 1e-10 order = 11 -[quantities.delta_even_real] -h5path = "ggj/delta_even_real" +[quantities.delta_tearing_real] +h5path = "ggj/delta_tearing_real" type = "real_scalar" extract = "value" -label = "Re(Δ_even)" +label = "Re(Δ₋ tearing)" noise_threshold = 1e-10 order = 12 -[quantities.delta_even_imag] -h5path = "ggj/delta_even_imag" +[quantities.delta_tearing_imag] +h5path = "ggj/delta_tearing_imag" type = "real_scalar" extract = "value" -label = "Im(Δ_even)" +label = "Im(Δ₋ tearing)" noise_threshold = 1e-10 order = 13 diff --git a/regression-harness/src/runner.jl b/regression-harness/src/runner.jl index bd4a10617..2e7ec5415 100644 --- a/regression-harness/src/runner.jl +++ b/regression-harness/src/runner.jl @@ -98,28 +98,27 @@ elapsed = time() - t_start """ # Self-contained computation for the GGJ inner-layer reference benchmark. -# Runs the Galerkin solver on the Glasser & Wang 2020 Eq. 55 parameter set -# at γ = 1 + i and writes (Δ_odd, Δ_even) into a small HDF5 file at ARGS[1]. +# Runs the Galerkin solver on the Glasser & Wang 2020 Eq. 55 coefficients with physical +# time scales (the q = 4 benchmark surface's τ_A, τ_R) at a 1 kHz rotation, γ = 2πi·10³ s⁻¹ +# (Q ≈ 0.22i, where the stock Galerkin matching is accurate), and writes the interchange (Δ₊) +# and tearing (Δ₋) matching data into a small HDF5 file at ARGS[1]. # Runtime and environment are written to ARGS[2] by the shared run-info epilogue. -# `solve_inner` returns the named-field `InnerLayerResponse`; the historical -# delta_odd / delta_even slots map to interchange / tearing respectively, which -# preserves the numeric identity of each tracked quantity. const COMPUTED_GGJ_SCRIPT_TEMPLATE = """ using Pkg %INSTANTIATE% using GeneralizedPerturbedEquilibrium using GeneralizedPerturbedEquilibrium.InnerLayer using HDF5 -γ = ComplexF64(1.0, 1.0) -p = glasser_wang_2020_eq55() +γ = 2π * 1e3 * im +p = InnerLayer.GGJ.glasser_wang_2020_eq55(; taua=2.11226e-7, taur=0.968219) t_start = time() Δ = solve_inner(GGJModel(solver=:galerkin), p, γ) elapsed = time() - t_start h5open(ARGS[1], "w") do fid - fid["ggj/delta_odd_real"] = real(Δ.interchange) - fid["ggj/delta_odd_imag"] = imag(Δ.interchange) - fid["ggj/delta_even_real"] = real(Δ.tearing) - fid["ggj/delta_even_imag"] = imag(Δ.tearing) + fid["ggj/delta_interchange_real"] = real(Δ.interchange) + fid["ggj/delta_interchange_imag"] = imag(Δ.interchange) + fid["ggj/delta_tearing_real"] = real(Δ.tearing) + fid["ggj/delta_tearing_imag"] = imag(Δ.tearing) end %RUNINFO% """ @@ -128,23 +127,22 @@ end # :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 and environment are # written to ARGS[2] by the shared run-info epilogue. -# As in the galerkin template, the delta_odd / delta_even slots map to interchange / tearing. const COMPUTED_GGJ_RAY_SCRIPT_TEMPLATE = """ using Pkg %INSTANTIATE% using GeneralizedPerturbedEquilibrium using GeneralizedPerturbedEquilibrium.InnerLayer using HDF5 -p = q4_surface_benchmark() +p = InnerLayer.GGJ.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(Δ.interchange) - fid["ggj/delta_odd_imag"] = imag(Δ.interchange) - fid["ggj/delta_even_real"] = real(Δ.tearing) - fid["ggj/delta_even_imag"] = imag(Δ.tearing) + fid["ggj/delta_interchange_real"] = real(Δ.interchange) + fid["ggj/delta_interchange_imag"] = imag(Δ.interchange) + fid["ggj/delta_tearing_real"] = real(Δ.tearing) + fid["ggj/delta_tearing_imag"] = imag(Δ.tearing) end %RUNINFO% """ diff --git a/src/ForceFreeStates/CoreTypes.jl b/src/ForceFreeStates/CoreTypes.jl index c321fa98f..ef384af7b 100644 --- a/src/ForceFreeStates/CoreTypes.jl +++ b/src/ForceFreeStates/CoreTypes.jl @@ -205,7 +205,7 @@ gpec.toml. # --- 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) + gal_inner_solver::String = "ray" # inner-layer Δ backend for the match: "ray" (rotated contour + inward march; recommended for every Q) or "galerkin" (real-axis Hermite FEM; accuracy depends on the surface and arg Q) # --- Inner-layer "galerkin" backend overrides (gal_inner_solver = "galerkin"; unset keeps the backend default) --- gal_inner_xfac::Union{Nothing,Float64} = nothing # asymptotic-matching radius multiplier gal_inner_nx::Union{Nothing,Int} = nothing # inner-layer grid cells diff --git a/src/ForceFreeStates/Matching/MatchProblem.jl b/src/ForceFreeStates/Matching/MatchProblem.jl index b2aadd42a..5e9727e2d 100644 --- a/src/ForceFreeStates/Matching/MatchProblem.jl +++ b/src/ForceFreeStates/Matching/MatchProblem.jl @@ -12,13 +12,15 @@ field `bpen`. `ffs` needs the Δ′ coil block, from `Galerkin(; rpec_flag=true) with `vac_flag=true`. # Keywords -- `eta`, `rho`, `rotation`: resistivity in Ω·m, mass density in kg/m³ and rotation frequency - in Hz, a scalar for every surface or one value per `ffs.surfaces` entry (core to edge). - Required unless `ideal=true`. -- `gamma`: ratio of specific heats. -- `ideal`: skip the inner layer and return the perfectly shielding solution. + + - `eta`, `rho`, `rotation`: resistivity in Ω·m, mass density in kg/m³ and rotation frequency + in Hz, a scalar for every surface or one value per `ffs.surfaces` entry (core to edge). + Required unless `ideal=true`. + - `gamma`: ratio of specific heats. + - `ideal`: skip the inner layer and return the perfectly shielding solution. ## Fields + - `ffs` - The outer solve being matched; held by reference, not copied. - `eta`, `rho`, `rotation` - Per-surface η in Ω·m, ρ in kg/m³ and rotation frequency in Hz, following `ffs.surfaces`. - `gamma` - Ratio of specific heats. @@ -126,45 +128,57 @@ function _compute_match(prob::MatchProblem, model::InnerLayer.GGJModel) # Layer-center field weight per parity (match.f intotsol_b). pen = zeros(ComplexF64, msing, 2) inner_psi = Vector{Vector{Float64}}(undef, msing) # ψ_s ± X·dψdx, ascending - inner_odd = Vector{Vector{ComplexF64}}(undef, msing) # Ξ₁, odd about ψ_s - inner_even = Vector{Vector{ComplexF64}}(undef, msing) # Ξ₂, even - inner_bodd = Vector{Vector{ComplexF64}}(undef, msing) # b^ψ from Ψ₁, even - inner_beven = Vector{Vector{ComplexF64}}(undef, msing) # b^ψ from Ψ₂, odd + inner_odd = Vector{Vector{ComplexF64}}(undef, msing) # Ξ₋, odd about ψ_s + inner_even = Vector{Vector{ComplexF64}}(undef, msing) # Ξ₊, even + inner_bodd = Vector{Vector{ComplexF64}}(undef, msing) # b^ψ from Ψ₋, even + inner_beven = Vector{Vector{ComplexF64}}(undef, msing) # b^ψ from Ψ₊, odd inner_params = Vector{InnerLayer.GGJParameters}(undef, msing) + profiles = true for i in 1:msing params = ggj_parameters(sings[i], equil; eta=prob.eta[i], rho=prob.rho[i], gamma=prob.gamma, ising=i) inner_params[i] = params γ = 2π * im * nn * prob.rotation[i] # forced eigenvalue, f in Hz rpec_eig[i] = γ - inner = InnerLayer.solve_inner_profile(model, params, γ) + inner = InnerLayer.solve_inner_profile(model, params, γ) # parity columns (Δ₊, Δ₋) deltar[i, 1] = inner.Δ[1] deltar[i, 2] = inner.Δ[2] profΨ, profΞ, xg = inner.Ψ, inner.Ξ, inner.x resc = inner.rescale # inner big-branch amplitude → outer δψ normalization - # b_m = 2πi·χ₁·(m−nq)·ξ_m with m−nq = −n·q′·dψdx·X and Ψ = XΞ (GWP2016 Eq. 16) ⇒ b_m = scale·resc·Ψ. + # b_m = 2πi·χ₁·(m−nq)·ξ_m with m−nq = −n·q′·dψdx·X and Ψ = XΞ at leading order (GWP2016 Eqs. 16, 27; Ψ normalization Eq. A17) ⇒ b_m = scale·resc·Ψ. scale = -2π * chi1 * im * nn * sings[i].q1 * inner.dψdx - pen[i, 1] = scale * profΨ[1, 1] * resc # X = 0, parity 1 - pen[i, 2] = scale * profΨ[1, 2] * resc # parity 2 (Ψ(0) = 0) + pen[i, 1] = scale * profΨ[1, 1] * resc # X = 0, even parity (Ψ(0) = 0) + pen[i, 2] = scale * profΨ[1, 2] * resc # odd parity + # A backend returning only layer-center values (real-ψ profiles not requested) skips the inner profiles. + if length(xg) == 1 + profiles = false + continue + end xvar = xg .* inner.dψdx # X → ψ distance inner_psi[i] = vcat(reverse(sings[i].psifac .- xvar), sings[i].psifac .+ xvar) - inner_odd[i] = resc .* vcat(reverse(.-profΞ[:, 1]), profΞ[:, 1]) - inner_even[i] = resc .* vcat(reverse(profΞ[:, 2]), profΞ[:, 2]) - inner_bodd[i] = (scale * resc) .* vcat(reverse(profΨ[:, 1]), profΨ[:, 1]) - inner_beven[i] = (scale * resc) .* vcat(reverse(.-profΨ[:, 2]), profΨ[:, 2]) + inner_odd[i] = resc .* vcat(reverse(.-profΞ[:, 2]), profΞ[:, 2]) + inner_even[i] = resc .* vcat(reverse(profΞ[:, 1]), profΞ[:, 1]) + inner_bodd[i] = (scale * resc) .* vcat(reverse(profΨ[:, 2]), profΨ[:, 2]) + inner_beven[i] = (scale * resc) .* vcat(reverse(.-profΨ[:, 1]), profΨ[:, 1]) end cout, cin, residual = _match_system(dp.raw, dp.coil, deltar) - # Penetrated field at each layer center (match.f intotsol_b). + # Penetrated field at each layer center (match.f intotsol_b); cin rows per surface are (d₊, d₋). bpen = zeros(ComplexF64, msing, mcoil) for i in 1:msing, j in 1:mcoil - bpen[i, j] = pen[i, 1] * cin[2i, j] + pen[i, 2] * cin[2i-1, j] + bpen[i, j] = pen[i, 1] * cin[2i-1, j] + pen[i, 2] * cin[2i, j] end # Matched inner ξ_ψ and b^ψ per coil: odd parity × cin[2i], even × cin[2i-1] (match.f intotsol). - inner_xi = [inner_odd[i] * transpose(cin[2i, :]) .+ inner_even[i] * transpose(cin[2i-1, :]) for i in 1:msing] - inner_b = [inner_bodd[i] * transpose(cin[2i, :]) .+ inner_beven[i] * transpose(cin[2i-1, :]) for i in 1:msing] + if profiles + inner_xi = [inner_odd[i] * transpose(cin[2i, :]) .+ inner_even[i] * transpose(cin[2i-1, :]) for i in 1:msing] + inner_b = [inner_bodd[i] * transpose(cin[2i, :]) .+ inner_beven[i] * transpose(cin[2i-1, :]) for i in 1:msing] + else + inner_psi = Vector{Float64}[] + inner_xi = Matrix{ComplexF64}[] + inner_b = Matrix{ComplexF64}[] + end end reconnected_flux = Matrix{ComplexF64}(dp.coil .+ transpose(dp.raw) * cout) @@ -173,7 +187,9 @@ function _compute_match(prob::MatchProblem, model::InnerLayer.GGJModel) if gal_sol === nothing xi = Array{ComplexF64,3}(undef, 0, 0, 0) xi_deriv = Array{ComplexF64,3}(undef, 0, 0, 0) - empty!(inner_psi); empty!(inner_xi); empty!(inner_b) + empty!(inner_psi) + empty!(inner_xi) + empty!(inner_b) else sols = gal_sol.xi # (mpert, ngrid, nsol) sols_deriv = gal_sol.xi_deriv @@ -197,7 +213,9 @@ function _compute_match(prob::MatchProblem, model::InnerLayer.GGJModel) if isempty(sols_cut) @warn "MatchProblem: the solve retained no cut solution (set cut_solution=true on Galerkin), " * "so the composite inner-region profiles are skipped; bpen and the matched outer ξ are unaffected" maxlog=1 - empty!(inner_psi); empty!(inner_xi); empty!(inner_b) + empty!(inner_psi) + empty!(inner_xi) + empty!(inner_b) else cut_range = gal_sol.cut_range keep = .!gal_sol.issing diff --git a/src/ForceFreeStates/Matching/ResonantMatch.jl b/src/ForceFreeStates/Matching/ResonantMatch.jl index 2f14fb03a..8fc93cf65 100644 --- a/src/ForceFreeStates/Matching/ResonantMatch.jl +++ b/src/ForceFreeStates/Matching/ResonantMatch.jl @@ -1,6 +1,6 @@ # ResonantMatch.jl # -# Matching result and the 4·msing outer↔inner system (rmatch match_rpec; Wang et al. 2020 PoP 27 122509, Eq. 11). +# Matching result and the 4·msing outer↔inner system (rmatch match_rpec; Wang et al., Phys. Plasmas 27, 122503 (2020), Eq. 11). """ MatchResult @@ -9,14 +9,16 @@ Matching coefficients and matched profiles, found in `result.match`. The profile fields are empty when the solve kept no Galerkin basis; the inner fields are empty with `ideal`. ## Fields + - `cout`, `cin`: `(2msing, ncoil)` outer and inner coefficients. - - `deltar`: `(msing, 2)` inner-layer Δ per surface and parity. + - `deltar`: `(msing, 2)` inner-layer Δ per surface, columns (Δ₊, Δ₋) = (interchange, tearing). - `rpec_eig`: forced eigenvalue γ = 2πi·n·f per surface. - `residual`: relative residual of the linear solve. - `bpen`: `(msing, ncoil)` penetrated resonant field at each layer center. - `reconnected_flux`: `(2msing, ncoil)` `Δ_coil + Δ_outᵀ·cout`. - `xi`, `xi_deriv`: `(mpert, ngrid, ncoil)` matched outer ξ and ξ′, one column per coil. - - `inner_psi`, `inner_xi`, `inner_b`: per surface, the inner ψ grid and the composite ξ_ψ and b^ψ on it. + - `inner_psi`, `inner_xi`, `inner_b`: per surface, the inner ψ grid and the composite ξ_ψ and b^ψ on it; empty unless the + inner model returns real-ψ profiles (`GGJModel(; solver=:galerkin)` or `GGJModel(; real_profiles=true)`). - `inner_params`: per-surface `GGJParameters` used. - `model`: the inner-layer model solved, with its backend options. """ @@ -50,8 +52,8 @@ function _match_system(dp_raw::AbstractMatrix, dp_coil::AbstractMatrix, deltar:: idx2 = 2ising idx3 = idx1 + 2msing idx4 = idx2 + 2msing - delta1 = deltar[ising, 1] - delta2 = deltar[ising, 2] + delta1 = deltar[ising, 1] # Δ₊; unknowns idx3, idx4 (cin rows 2i−1, 2i) are d₊, d₋ (GWP2016 Eqs. 36–37) + delta2 = deltar[ising, 2] # Δ₋ mat[idx1, idx1] = 1 mat[idx2, idx2] = 1 mat[idx1, idx3] = -1 @@ -111,7 +113,7 @@ const MATCH_H5_ANNOTATIONS = [ "SingularSurfaces/Match/cin" => (; long_name="inner-region matching coefficients per surface parity, per unit coil-response column", dims=("surface_parity", "coil")), "SingularSurfaces/Match/Delta_r" => - (; long_name="inner-layer Δ per surface and parity (deltac convention)", dims=("surface", "parity")), + (; long_name="inner-layer Δ per surface and parity, (Δ₊ interchange, Δ₋ tearing)", dims=("surface", "parity")), "SingularSurfaces/Match/bpen" => (; long_name="penetrated resonant field b^ψ at each layer center, per unit coil-response column", dims=("surface", "coil")), "SingularSurfaces/Match/reconnected_flux" => diff --git a/src/InnerLayer/GGJ/GGJ.jl b/src/InnerLayer/GGJ/GGJ.jl index cdd6e37d6..fcb299164 100644 --- a/src/InnerLayer/GGJ/GGJ.jl +++ b/src/InnerLayer/GGJ/GGJ.jl @@ -1,40 +1,34 @@ # GGJ.jl # -# Glasser–Greene–Johnson resistive inner-layer model. Provides three -# interchangeable solvers selected via the `solver` type-parameter of -# `GGJModel`: +# Glasser–Greene–Johnson resistive inner-layer model, with three solvers selected by the +# `solver` type parameter of `GGJModel`: # -# - `: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 +# - `:shooting` – backward shoot from X_max → 0; valid only for |Q| ≪ 1 (tests and comparisons) +# - `:galerkin` – Hermite-cubic finite elements with resonant/extension matching cells +# (GWP2016 Sec. III); accuracy depends on the surface and arg Q +# - `:ray` – the same Hermite assembly on the rotated contour (Ray.jl), with the far-field +# condition marched in from where the series converges; robust to |Q| ~ 500 # -# They return the parity-projected matching data -# `(Δ_odd, Δ_even)` of GWP2016 Eqs. (34)–(35) in the same (deltac) convention. +# They return the matching data Δ± of GWP2016 Eqs. (34)–(35). Parity is that of the displacement Ξ: +# Δ₊ (even) is the interchange channel, Δ₋ (odd) the tearing channel. Positional pairs are always +# ordered (Δ₊, Δ₋), as d₊, d₋ in GWP2016 Eq. (38) and Pletzer, Bondeson & Dewar (1994) Eqs. 83, 85. # -# Equation references throughout this module use two source papers: +# Main source papers, cited by equation number throughout: # # GWP2016 — A. H. Glasser, Z. R. Wang & J.-K. Park, "Computation of resistive # instabilities by matched asymptotic expansions", Phys. Plasmas 23, -# 112506 (2016). Inner-region equations (Eq. 11), matrix form -# A Ψ'' + B Ψ' + C Ψ = 0 (Eqs. 12–15), singular-Galerkin weak form -# (Eq. 32), grid packing (Eq. 33), matching data (Eqs. 34–35), -# dimensionless parameters / scale factors (Appendix, Eqs. A8–A15). +# 112506 (2016). Inner-region equations (Eq. 11), singular-Galerkin method +# (Sec. III), matching data (Eqs. 34–35), scale factors (Appendix, Eqs. A8–A16). # # GW2020 — A. H. Glasser & Z. R. Wang, "Asymptotic solutions and convergence # studies of the resistive inner region equations", Phys. Plasmas 27, # 012506 (2020). Wasow construction of the large-x asymptotic basis # (Eqs. 1–55); implemented in InnerAsymptotics.jl. -# -# The inner-region equations are identical in both: GW2020 Eq. (1) ≡ GWP2016 Eq. (11). module GGJ using LinearAlgebra using StaticArrays -using SparseArrays using Random using Printf using DoubleFloats: Double64 @@ -47,10 +41,10 @@ import ..solve_inner, ..solve_inner_profile GGJModel(; solver=:ray, options...) 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). Backend keywords given at construction, e.g. -`GGJModel(; solver=:galerkin, nx=1280)`, apply to every solve. +backend: `:ray` (default; rotated contour with an inward-marched far-field condition, +robust at large |Q|), `:galerkin` (real-axis Hermite FEM with resonant/extension matching +cells; accuracy depends on the surface and arg Q, even at small |Q|), or `:shooting` (|Q| ≪ 1 only). Backend keywords given at +construction, e.g. `GGJModel(; solver=:galerkin, nx=1280)`, apply to every solve. ## Fields @@ -76,7 +70,6 @@ end include("GGJParameters.jl") include("InnerAsymptotics.jl") -include("RayAsymptotics.jl") include("Reference.jl") include("Shooting.jl") include("Galerkin.jl") @@ -84,10 +77,6 @@ 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 +export solve_ray, pick_smax, delta_convergence, solution_profile, asymptotic_profile end # module GGJ diff --git a/src/InnerLayer/GGJ/GGJParameters.jl b/src/InnerLayer/GGJ/GGJParameters.jl index 4acca7fef..3d5f59e8b 100644 --- a/src/InnerLayer/GGJ/GGJParameters.jl +++ b/src/InnerLayer/GGJ/GGJParameters.jl @@ -11,9 +11,9 @@ 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). +GWP2016 Eq. (A8); E, F, G, H, K enter the inner-region equations (Eq. 11). -Fields: +## Fields | field | meaning | |:------- |:-------------------------------------------------------------- | @@ -22,10 +22,10 @@ Fields: | `G` | Coupling coefficient (curvature × pressure gradient) | | `H` | Pfirsch–Schlüter coefficient | | `K` | Glasser parameter | -| `M` | Mercier-related auxiliary parameter (held but not used here) | +| `M` | Inertia factor (A8; enters K and Q₀, X₀); held, not used here | | `taua` | Local Alfvén time at the rational surface | | `taur` | Local resistive time at the rational surface | -| `v1` | Linear scale factor used in the V₁ rescaling | +| `v1` | Coordinate scale: inner X = v₁·δψ/X₀; enters `rescale_delta` | | `ising` | Index of the singular surface (traceability only) | The growth rate `γ` is not stored here; it is a separate argument to `solve_inner`. @@ -79,29 +79,30 @@ sfac(p::GGJParameters) = p.taur / p.taua """ x0(p::GGJParameters) -> Float64 -Inner-layer length scale `X₀ = S^(−1/3)` (GWP2016 Eq. A14); physical `x = X₀ X`. +Inner-layer length scale `X₀ = S^(−1/3)` (GWP2016 Eq. A16); physical `x = X₀ X` (Eq. A11). """ x0(p::GGJParameters) = sfac(p)^(-1.0 / 3.0) """ q0(p::GGJParameters) -> Float64 -Growth-rate scale `Q₀ = X₀ / τ_A` (GWP2016 Eq. A15); physical rate `= Q₀ Q`. +Growth-rate scale `Q₀ = X₀ / τ_A` (GWP2016 Eq. A16); physical rate `= Q₀ Q` (Eq. A11). """ q0(p::GGJParameters) = x0(p) / p.taua """ inner_Q(p::GGJParameters, γ::Number) -> ComplexF64 -Scaled inner-layer growth rate `Q = γ / Q₀` (GWP2016 Eq. A15). +Scaled inner-layer growth rate `Q = γ / Q₀` (GWP2016 Eq. A11). """ inner_Q(p::GGJParameters, γ::Number) = ComplexF64(γ) / q0(p) """ rescale_delta(Δ, p::GGJParameters) -> SVector{2,ComplexF64} -Rescale the matching data to physical Δ by `S^(2√(−D_I)/3) · v₁^(2√(−D_I))` -(GWP2016 Sec. IV), element-wise on `(Δ_odd, Δ_even)`. +Rescale the matching data to physical Δ by `S^(2√(−D_I)/3) · v₁^(2√(−D_I))`, element-wise +on the pair (Δ₊, Δ₋). The factor X₀^(−2√(−D_I)) = S^(2√(−D_I)/3) is Wang et al., Phys. Plasmas +27, 122503 (2020), text after Eq. 11. """ function rescale_delta(Δ::AbstractVector, p::GGJParameters) s = sfac(p) @@ -110,8 +111,6 @@ function rescale_delta(Δ::AbstractVector, p::GGJParameters) 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 for solve_inner_profile: δψ per unit inner coordinate (X = v₁·δψ/X₀) and the +# large-solution amplitude rescale (v₁/X₀)^(1/2+p₁) to the outer δψ normalization. _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 a6c287438..10fe8c052 100644 --- a/src/InnerLayer/GGJ/Galerkin.jl +++ b/src/InnerLayer/GGJ/Galerkin.jl @@ -1,85 +1,23 @@ # Galerkin.jl # -# Hermite-cubic finite element (Galerkin) solver for the GGJ inner-layer -# model, in the "resonant + noexp + inps" configuration. The half-domain -# problem (x ∈ [0, xmax]) is solved with parity boundary conditions at -# x = 0 and asymptotic matching at x = xmax. -# -# The solved field is the 3-component inner-layer solution u = (Ψ, Ξ, Υ) -# obeying the GGJ inner-region equations -# GWP2016 Eq. (11) ≡ GW2020 Eq. (1): -# -# Ψ_xx − H Υ_x − Q(Ψ − x Ξ) = 0, -# Q² Ξ_xx − Q x² Ξ + Q x Ψ + (E + F) Υ + H Ψ_x = 0, -# Q Υ_xx − x² Υ + x Ψ + Q²[G(Ξ − Υ) − K(E Ξ + F Υ + H Ψ_x)] = 0, -# -# recast in matrix form A Ψ'' + B Ψ' + C Ψ = 0 (GWP2016 Eqs. 12–15) and -# assembled here in the weak-form coefficient layout `I·u'' − V·u' − U·u = 0` -# with (I, V, U) = (A, −B, −C). See `_physical_uv` for the row-by-row map. -# -# This module ports the GWP2016 inner-region singular-Galerkin method (Sec. III): -# the half-domain weak form (Eq. 32) on a packed Hermite-cubic grid (Eq. 33), -# with resonant/extension cells (Fig. 1) supplying the matching data (Eqs. 34–35). -# -# Configuration assumed: -# gal_method = "resonant" method = true noexp = true -# basis_type = 0 (Hermite) fulldomain = 0 inps_type = "inps" -# mpert = 3 (Ψ, Ξ, Υ) np = 3 (Hermite cubic → 4 DOFs/node) -# -# References: GWP2016 (Phys. Plasmas 23, 112506) for the Galerkin method and -# matching; GW2020 (Phys. Plasmas 27, 012506) for the large-x asymptotic basis. +# Hermite-cubic finite-element (Galerkin) solver for the GGJ inner layer, porting the GWP2016 +# Sec. III singular-Galerkin method: the half-domain weak form (Eq. 32) on a packed grid +# (Eq. 33), parity conditions at x = 0, and resonant/extension cells at x_max supplying the +# matching data (Eqs. 34–35). Fortran deltac configuration: "resonant" + noexp + inps, +# Hermite basis, mpert = 3 (Ψ, Ξ, Υ). `_assemble_and_solve_point!` reuses the assembly on the +# rotated ray with a point far-field condition (Ray.jl). using LinearAlgebra using FastGaussQuadrature: gausslobatto using QuadGK: quadgk -# ----------------------------------------------------------------------- -# Physical-variable helpers. -# ----------------------------------------------------------------------- - """ -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. 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 (Ψ, Ξ, Υ). -""" -function _physical_ua_dua(cache::InnerAsymptoticsCache, x::Real) - U, dU = evaluate_asymptotics(cache, x; 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] / x - ua[2, j] = U[2, j] - ua[3, j] = U[3, j] - dua[1, j] = U[4, j] - U[1, j] / (x * x) - dua[2, j] = U[5, j] * x - dua[3, j] = U[6, j] * x - end - return ua, dua -end + _physical_uv(params, Q, x) -> (I, U, V) +3×3 coefficients of `I·u'' − V·u' − U·u = 0` for u = (Ψ, Ξ, Υ) at real or complex `x`, with +(I, V, U) = (A, −B, −C) of GWP2016 Eqs. 12–15; each row is one equation of GWP2016 Eq. 11. """ -Build the (I, U, V) coefficient matrices of the second-order system -`I·u'' − V·u' − U·u = 0` for u = (Ψ, Ξ, Υ) at coordinate `x`. -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)) -""" -function _physical_uv(params::GGJParameters, Q::ComplexF64, x::Real) +function _physical_uv(params::GGJParameters, Q::ComplexF64, x::Number) e = ComplexF64(params.E) f = ComplexF64(params.F) h = ComplexF64(params.H) @@ -89,26 +27,15 @@ function _physical_uv(params::GGJParameters, Q::ComplexF64, x::Real) q2 = q * q x2 = x * x - # I = A of GWP2016 Eq. (14): diag(1, Q², Q), coefficients of (Ψ_xx, Ξ_xx, Υ_xx). Imat = @SMatrix ComplexF64[1 0 0; 0 q2 0; 0 0 q] - # U = −C of GWP2016 Eq. (15): coefficients of −u in each equation - # (stored pre-scaled: row 2 ×Q², row 3 ×Q, matching - # the I-row normalization). - # row 1 (Ψ): (Q, −Q x, 0) - # row 2 (Ξ): (−Q x, Q x², −(E+F)) - # row 3 (Υ): (−x, −Q²(G−KE), x²+Q²(G+KF)) + # U = −C (GWP2016 Eq. 15) and V = −B (Eq. 14), rows scaled by diag(1, Q², Q) like `Imat`. #! 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 ] - - # V = −B of GWP2016 Eq. (14): coefficients of −u' in each equation (same row scaling as U): - # row 1 (Ψ): (0, 0, H) ⇒ −H Υ_x - # 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 @@ -119,12 +46,7 @@ function _physical_uv(params::GGJParameters, Q::ComplexF64, x::Real) return Imat, U, V end -# ----------------------------------------------------------------------- -# Hermite cubic basis functions. -# Returns (pb[1:4], qb[1:4]) where pb = values, qb = derivatives, -# indexed 1:4 → physical 0:3. -# ----------------------------------------------------------------------- - +# Hermite cubic basis on [x0, x1]: values pb and derivatives qb, slots 1:4 ↔ (u₀, u₀', u₁, u₁'). function _hermite(x::Real, x0::Real, x1::Real) dx = x1 - x0 t0 = (x - x0) / dx @@ -146,14 +68,8 @@ function _hermite(x::Real, x0::Real, x1::Real) return pb, qb end -# ----------------------------------------------------------------------- -# Grid packing. The pfac > 1 branch is the GWP2016 -# inner-region packing X(ξ,λ) = ln[(1+λξ)/(1−λξ)] / ln[(1+λ)/(1−λ)] -# (Eq. 33), which concentrates nodes near X = 0; pfac = 1 gives a uniform -# grid, and pfac in (0,1) the complementary edge-packed map. The grid-density -# ratio is P(λ) = 1 − λ² (GWP2016 between Eqs. 5 and 33). -# ----------------------------------------------------------------------- - +# Grid packing: pfac > 1 concentrates nodes near X = 0 (GWP2016 Eq. 33), pfac = 1 is uniform, +# 0 < pfac < 1 packs toward the edge (Eq. 5). function _pack(nx::Int, pfac::Float64, side::String) xi = zeros(Float64, 2nx + 1) if side == "right" @@ -192,24 +108,14 @@ function _pack(nx::Int, pfac::Float64, side::String) return x end -# ----------------------------------------------------------------------- -# Three-level xmax sweep. The cutoff -# x_max is the position where the asymptotic-basis residual Δ (GW2020 Eq. 54, -# computed by `asymptotic_residual`) drops below a target tolerance — the -# choice studied in GW2020 Sec. III and Fig. 3, where x_max is taken where the -# error falls to ~1e−7. Three tolerance levels set x_max and the two thin -# matching cells (dx1, dx2) just inside it. -# Returns (xmax, dx1, dx2, cache). -# ----------------------------------------------------------------------- - +# x_max and the two thin matching cells inside it, from where the series residual (GW2020 Eq. 54) +# crosses three tolerances. Returns (xmax, dx1, dx2, cache). function _xmax_3level(params::GGJParameters, Q::ComplexF64; kmax::Int=8, xfac::Float64=1.0, eps_vec::NTuple{3,Float64}=(1e-2, 5e-7, 1e-7)) cache = build_asymptotics(params, Q; kmax=kmax) - # Bisection-based root finder: find the exact x where max(residual) = eps, - # for each of the 3 tolerance levels. First do a coarse sweep to bracket, - # then bisect to machine precision. + # Bracket each crossing on a coarse log sweep, then bisect. dxfac = 10.0^0.01 xmax_vals = zeros(Float64, 3) brackets = Vector{Tuple{Float64,Float64}}(undef, 3) @@ -256,51 +162,70 @@ function _xmax_3level(params::GGJParameters, Q::ComplexF64; return (; xmax, dx1, dx2, cache) end -# ----------------------------------------------------------------------- -# Cell types and grid construction. These are the grid cells of GWP2016 -# Fig. 1, laid out moving toward the cutoff x_max: ordinary "normal" cells -# (CT_NONE) carrying 4 Hermite-cubic DOFs/node, "extension" cells (CT_EXT2, -# CT_EXT1, CT_EXT) where the small resonant power series is blended into the -# Hermite basis, and the "resonant" cell (CT_RES) at x_max holding the large -# resonant solution whose coefficient provides the matching data (Eqs. 34–35). -# Here the half-domain has a single resonant cell at x_max (the X = 0 end is -# the parity boundary, not a singular surface). -# ----------------------------------------------------------------------- - +# Cell types toward x_max (GWP2016 Sec. III text after Eq. 33; the outer-region analogue is Fig. 1): +# normal Hermite cells (CT_NONE), extension cells blending the asymptotic solutions into the Hermite +# basis (CT_EXT2, CT_EXT1, CT_EXT), and the resonant cell (CT_RES) whose large-solution coefficient is +# the matching data. x = 0 carries the parity conditions. @enum CellType CT_NONE CT_EXT2 CT_EXT1 CT_EXT CT_RES +""" + GalerkinCell + +One finite element of the Galerkin grid. + +## Fields + + - `etype` - Cell type (`CT_NONE`, `CT_EXT2`, `CT_EXT1`, `CT_EXT`, `CT_RES`). + - `x_left`, `x_right` - Cell end points. + - `x_lsode` - Upper limit of the resonant-cell integral and the x_max surface term (`CT_RES` only). + - `np` - Highest local Hermite slot per component: 3 for Hermite cells, 1 for `CT_EXT` (left node only), −1 for `CT_RES`. + - `map` - `(mpert, np+1)` local-to-global DOF map; `CT_EXT` appends a column of asymptotic DOFs, `CT_RES` holds only those. + - `emap` - Asymptotic DOF indices shared by `CT_EXT` and `CT_RES`; only `emap[1]` is active (noexp). +""" struct GalerkinCell etype::CellType x_left::Float64 x_right::Float64 - x_lsode::Float64 # right integration limit for res cell - np::Int # Hermite DOF count per component (-1 for res) - map::Matrix{Int} # (mpert, 0:np_eff) local-to-global DOF map - emap::Vector{Int} # extra (asymptotic) DOF indices (length 1 for noexp) + x_lsode::Float64 + np::Int + map::Matrix{Int} + emap::Vector{Int} end +""" + GalerkinWorkspace + +Grid, banded system and scratch buffers of one Galerkin solve. + +## Fields + + - `cells` - Grid cells from x = 0 outward. + - `ndim` - Number of global unknowns. + - `nx` - Number of cells. + - `kl` - Banded half-bandwidth (`ku = kl`). + - `mat` - `(2kl+ku+1, ndim, 2)` LAPACK banded matrices, one per parity (1 = even/Δ₊, 2 = odd/Δ₋). + - `rhs`, `sol` - `(ndim, 2)` right-hand sides and solutions per parity. + - `cell_mat_buf`, `cell_mat_ext_buf` - `(3, 3, 4, 4)` per-cell scratch for the Hermite and extension terms. + - `cell_rhs_ext_buf` - `(3, 4)` per-cell extension RHS scratch. + - `ab_buf`, `rhs_buf` - Scratch copies for the banded LU solve. +""" struct GalerkinWorkspace cells::Vector{GalerkinCell} ndim::Int nx::Int kl::Int - mat::Array{ComplexF64,3} # (ldab, ndim, 2) banded storage - rhs::Matrix{ComplexF64} # (ndim, 2) - sol::Matrix{ComplexF64} # (ndim, 2) - # Reusable scratch buffers, zeroed per-cell via `fill!`. Eliminates the - # per-cell `zeros(...)` that otherwise allocates thousands of MiB over a - # full dispersion scan. - cell_mat_buf::Array{ComplexF64,4} # (mpert=3, mpert, np+1=4, np+1=4) - cell_mat_ext_buf::Array{ComplexF64,4} # (3, 3, 4, 4) max over CT_EXT/EXT1/EXT2 - cell_rhs_ext_buf::Matrix{ComplexF64} # (3, 4) - ab_buf::Matrix{ComplexF64} # (ldab, ndim) scratch for banded LU - rhs_buf::Vector{ComplexF64} # (ndim,) scratch for banded solve + mat::Array{ComplexF64,3} + rhs::Matrix{ComplexF64} + sol::Matrix{ComplexF64} + cell_mat_buf::Array{ComplexF64,4} + cell_mat_ext_buf::Array{ComplexF64,4} + cell_rhs_ext_buf::Matrix{ComplexF64} + ab_buf::Matrix{ComplexF64} + rhs_buf::Vector{ComplexF64} end function _build_grid_and_workspace(nx::Int, xmax::Float64, dx1::Float64, dx2::Float64, - pfac::Float64, cutoff::Int, nq::Int) - mpert = 3 - np = 3 + pfac::Float64, cutoff::Int) side = pfac < 1 ? "left" : "right" # Cell types @@ -332,6 +257,18 @@ function _build_grid_and_workspace(nx::Int, xmax::Float64, dx1::Float64, dx2::Fl x_nodes[1] = 0.0 x_nodes[ixmax+1] = x1_packed + return _galerkin_workspace(x_nodes, etypes) +end + +# Cells, local-to-global Hermite DOF maps, bandwidth and buffers for given nodes and cell types. +# The stock layout ends in an EXT/RES pair; an all-CT_NONE layout is the plain Hermite grid +# used with an external far-field condition (`_assemble_and_solve_point!`). +function _galerkin_workspace(x_nodes::Vector{Float64}, etypes::Vector{CellType}) + mpert = 3 + np = 3 + nx = length(etypes) + has_res = etypes[end] == CT_RES + # Build cells and local-to-global mapping cells = Vector{GalerkinCell}(undef, nx) imap = 1 @@ -343,7 +280,7 @@ function _build_grid_and_workspace(nx::Int, xmax::Float64, dx1::Float64, dx2::Fl cell_np = if et == CT_NONE || et == CT_EXT1 || et == CT_EXT2 np elseif et == CT_EXT - 1 # method=true: left Hermite pair only + 1 extra + 1 # left-node Hermite pair only, plus one asymptotic DOF else # CT_RES -1 # no Hermite DOFs end @@ -385,25 +322,25 @@ function _build_grid_and_workspace(nx::Int, xmax::Float64, dx1::Float64, dx2::Fl end # Set resonant cell mapping: shares emap with ext cell - ext_cell = cells[nx-1] - emap_vals = ext_cell.emap - res_map = reshape(copy(emap_vals), mpert, 1) - cells[nx] = GalerkinCell(CT_RES, x_nodes[nx], x_nodes[nx+1], x_nodes[nx+1], - -1, res_map, copy(emap_vals)) + if has_res + ext_cell = cells[nx-1] + emap_vals = ext_cell.emap + res_map = reshape(copy(emap_vals), mpert, 1) + cells[nx] = GalerkinCell(CT_RES, x_nodes[nx], x_nodes[nx+1], x_nodes[nx+1], + -1, res_map, copy(emap_vals)) + end - ndim = imap # noexp → ndim = imap (only emap[1] ≤ ndim) + ndim = has_res ? imap : imap - 1 # noexp → ndim = imap (only emap[1] ≤ ndim); plain grid: Hermite DOFs only # Bandwidth - kl = mpert * (np + 1) + 1 - 1 # resonant method with noexp; +1-1 kept for indexing agreement + kl = mpert * (np + 1) # one cell's mpert·(np+1) DOFs ku = kl ldab = 2kl + ku + 1 mat = zeros(ComplexF64, ldab, ndim, 2) rhs = zeros(ComplexF64, ndim, 2) sol = zeros(ComplexF64, ndim, 2) - # Preallocate per-cell scratch buffers sized to the max case (np+1=4). - # Smaller cells (e.g. CT_EXT with cell.np=1) use a (2×2) sub-slice and - # rely on fill!(buf, 0) to keep the remainder zero. + # Scratch buffers sized for np + 1 = 4; smaller cells use a zeroed sub-slice. cell_mat_buf = zeros(ComplexF64, mpert, mpert, np + 1, np + 1) cell_mat_ext_buf = zeros(ComplexF64, mpert, mpert, np + 1, np + 1) cell_rhs_ext_buf = zeros(ComplexF64, mpert, np + 1) @@ -415,19 +352,11 @@ function _build_grid_and_workspace(nx::Int, xmax::Float64, dx1::Float64, dx2::Fl ab_buf, rhs_buf) end -# ----------------------------------------------------------------------- -# Local assembly: Gauss quadrature for "none" cells. -# Evaluates the GWP2016 Eq. (32) bilinear form on one cell. After integrating -# the A-term by parts, the symmetric weak form is -# Σⱼ [−(αᵢ', A αⱼ') + (αᵢ, B αⱼ') + (αᵢ, C αⱼ)] uⱼ , -# which in this code's (I, V, U) = (A, −B, −C) convention is assembled below as -# I·(αᵢ' αⱼ') + V·(αᵢ αⱼ') + U·(αᵢ αⱼ) with αᵢ the Hermite-cubic basis -# (pb = values, qb = derivatives). Scalar products use Gauss–Lobatto quadrature. -# ----------------------------------------------------------------------- - +# One cell's GWP2016 Eq. 32 bilinear form after integrating the A-term by parts, +# I·(αᵢ' αⱼ') + V·(αᵢ αⱼ') + U·(αᵢ αⱼ), by Gauss–Lobatto quadrature. function _gauss_quad!(cell_mat::Array{ComplexF64,4}, cell::GalerkinCell, quad_nodes::Vector{Float64}, quad_weights::Vector{Float64}, - params::GGJParameters, Q::ComplexF64) + params::GGJParameters, Q::ComplexF64; θ::Float64=0.0) np_cell = cell.np nq = length(quad_nodes) x0 = (cell.x_right + cell.x_left) / 2 @@ -437,11 +366,16 @@ function _gauss_quad!(cell_mat::Array{ComplexF64,4}, cell::GalerkinCell, @inbounds for iq in 1:nq x = x0 + dx * quad_nodes[iq] w = dx * quad_weights[iq] - Imat, Umat, Vmat = _physical_uv(params, Q, x) + if iszero(θ) + Imat, Umat, Vmat = _physical_uv(params, Q, x) + else + # Rotated contour x = e^{iθ}s: d/dx = e^{-iθ} d/ds; times e^{2iθ}, I is unchanged, V gains e^{iθ}, U gains e^{2iθ}. + Imat, U0, V0 = _physical_uv(params, Q, cis(θ) * x) + Umat, Vmat = cis(2θ) * U0, cis(θ) * V0 + end pb, qb = _hermite(x, cell.x_left, cell.x_right) for ip in 0:np_cell, jp in 0:np_cell - # cell_mat[:,:,ip+1,jp+1] += w*(I*qb[ip]*qb[jp] + V*pb[ip]*qb[jp] + U*pb[ip]*pb[jp]) for jpert in 1:3, ipert in 1:3 cell_mat[ipert, jpert, ip+1, jp+1] += w * ( @@ -454,15 +388,8 @@ function _gauss_quad!(cell_mat::Array{ComplexF64,4}, cell::GalerkinCell, end end -# ----------------------------------------------------------------------- -# Extension assembly. Handles ext, ext1, ext2 cell -# types: the GWP2016 extension cells (Fig. 1) where the small resonant power -# series (column 2 of the asymptotic basis) drives the response and the large -# resonant solution (column 1) is blended into the Hermite cubics with -# C¹ continuity. This couples the bulk Hermite DOFs to the asymptotic emap DOF -# via the same Eq. (32) bilinear form; the small-solution terms form the RHS. -# ----------------------------------------------------------------------- - +# Extension cells (GWP2016 Sec. III): the large resonant solution (column 1) is blended into the +# Hermite cubics with C¹ continuity as an extra DOF; the small one (column 2) drives the RHS. function _extension!(cell_mat::Array{ComplexF64,4}, cell_rhs::Matrix{ComplexF64}, cell::GalerkinCell, quad_nodes::Vector{Float64}, quad_weights::Vector{Float64}, @@ -475,23 +402,19 @@ function _extension!(cell_mat::Array{ComplexF64,4}, cell_rhs::Matrix{ComplexF64} if cell.etype == CT_EXT # Asymptotic basis at right boundary of ext cell - ua_bdy, dua_bdy = _physical_ua_dua(cache, cell.x_right) + ua_bdy, dua_bdy = physical_ua_dua(cache, cell.x_right) ep = np_cell + 1 # index of the extra DOF (1-based → ep+1 in array) @inbounds for iq in 1:nq x = x0c + dxc * quad_nodes[iq] w = dxc * quad_weights[iq] - uax, duax = _physical_ua_dua(cache, x) + uax, duax = physical_ua_dua(cache, x) Imat, Umat, Vmat = _physical_uv(params, Q, x) pb, qb = _hermite(x, cell.x_left, cell.x_right) - # The asymptotic basis function in the ext cell uses Hermite blending - # at the right node: phi_j(x) = ua_j * pb[3] + dua_j * pb[4] - # (pb[3] = H_10, pb[4] = H_11 at right node) - - # (h_i, L*phi_j) for i=0:np_cell, j=1 (only large algebraic, noexp) - # and (phi_j^T, L*h_i) for the transpose - jp = 1 # only tid(1) = column 1 = large algebraic (with noexp) + # Extra DOF φ = ua₁(x_r)·H₁₀ + ua₁'(x_r)·H₁₁, the large solution Hermite-blended at the right + # node; assembles (hᵢ, Lφ), (φ, Lhᵢ) and (φ, Lφ). + jp = 1 # column 1 = large solution (noexp) ua1 = ua_bdy[:, jp] * pb[3] + dua_bdy[:, jp] * pb[4] dua1 = ua_bdy[:, jp] * qb[3] + dua_bdy[:, jp] * qb[4] @@ -511,7 +434,7 @@ function _extension!(cell_mat::Array{ComplexF64,4}, cell_rhs::Matrix{ComplexF64} term = transpose(dua2_s) * Imat * dua1 + transpose(ua2_s) * Vmat * dua1 + transpose(ua2_s) * Umat * ua1 cell_mat[1, 1, ep+1, ep+1] += w * term - # RHS: driving from small algebraic (column 2 = tid(4)) + # RHS: driven by the small solution (column 2) ua_drv = uax[:, 2] dua_drv = duax[:, 2] for ip in 0:np_cell @@ -532,11 +455,11 @@ function _extension!(cell_mat::Array{ComplexF64,4}, cell_rhs::Matrix{ComplexF64} pb, qb = _hermite(x, cell.x_left, cell.x_right) if cell.etype == CT_EXT1 - uax, duax = _physical_ua_dua(cache, x) + uax, duax = physical_ua_dua(cache, x) ua_drv = uax[:, 2] dua_drv = duax[:, 2] else # CT_EXT2 - ua_bdy, dua_bdy = _physical_ua_dua(cache, cell.x_right) + ua_bdy, dua_bdy = physical_ua_dua(cache, cell.x_right) ua_drv = ua_bdy[:, 2] * pb[3] + dua_bdy[:, 2] * pb[4] dua_drv = ua_bdy[:, 2] * qb[3] + dua_bdy[:, 2] * qb[4] end @@ -549,28 +472,16 @@ function _extension!(cell_mat::Array{ComplexF64,4}, cell_rhs::Matrix{ComplexF64} end end -# ----------------------------------------------------------------------- -# Resonant integral, evaluated with QuadGK. Evaluates the -# GWP2016 Eq. (32) bilinear form over the resonant cell using the asymptotic -# (power-series) solutions as both test and trial functions — the scalar -# products that GWP2016 notes "are non-analytic and therefore require an -# adaptive integration" rather than Gauss quadrature. -# Computes ∫_{x_left}^{x_lsode} ua_1^T L ua_j dx for j=1 and j=2. -# With noexp, only (1,1) and (1,2) entries of the 1×2 result matter -# (ip=1 only for test function, jp=1 for stiffness, jp=2 for RHS). -# ----------------------------------------------------------------------- - +# Resonant-cell scalar products ∫ ua₁ᵀ L ua_j dx (j = 1 stiffness, j = 2 RHS) with the asymptotic +# solutions as test and trial functions; non-analytic, so adaptive QuadGK (GWP2016 Sec. III). function _resonant_integral(cell::GalerkinCell, params::GGJParameters, Q::ComplexF64, cache::InnerAsymptoticsCache; tol::Float64=1e-5) x_left = cell.x_left x_right = cell.x_lsode - # u_res[1,1] = ∫ ua_1^T L ua_1 dx (stiffness) - # u_res[1,2] = ∫ ua_1^T L ua_2 dx (RHS/driving) - # Weak form after IBP: (test') I (trial') + (test) V (trial') + (test) U (trial) function integrand_11(x) - ua, dua = _physical_ua_dua(cache, x) + ua, dua = physical_ua_dua(cache, x) Imat, Umat, Vmat = _physical_uv(params, Q, x) ua1 = ua[:, 1] dua1 = dua[:, 1] @@ -578,7 +489,7 @@ function _resonant_integral(cell::GalerkinCell, params::GGJParameters, end function integrand_12(x) - ua, dua = _physical_ua_dua(cache, x) + ua, dua = physical_ua_dua(cache, x) Imat, Umat, Vmat = _physical_uv(params, Q, x) dua1 = dua[:, 1] ua1 = ua[:, 1] @@ -609,43 +520,23 @@ function _assemble_and_solve!(ws::GalerkinWorkspace, fill!(ws.mat, 0) fill!(ws.rhs, 0) - # Per-cell assembly — reuse the preallocated scratch buffers, zeroing - # only the sub-slice actually used by this cell's np_eff. cell_mat = ws.cell_mat_buf cell_mat_ext = ws.cell_mat_ext_buf cell_rhs_ext = ws.cell_rhs_ext_buf for ix in 1:ws.nx cell = ws.cells[ix] - # Gauss quadrature for Hermite contribution (all cell types) + # Gauss quadrature for the Hermite contribution (all cells but CT_RES) if cell.np >= 0 np_eff = cell.np fill!(cell_mat, 0) _gauss_quad!(cell_mat, cell, quad_nodes, quad_weights, params, Q) - - # Assemble into global banded matrix (both parities use same base matrix) - for ip in 0:np_eff, ipert in 1:mpert - i = cell.map[ipert, ip+1] - if i > ws.ndim - continue - end - for jp in 0:np_eff, jpert in 1:mpert - j = cell.map[jpert, jp+1] - if j > ws.ndim - continue - end - ws.mat[offset+i-j, j, 1] += cell_mat[ipert, jpert, ip+1, jp+1] - end - end + _scatter_cell!(ws, cell, cell_mat, np_eff) end # Extension terms if cell.etype in (CT_EXT, CT_EXT1, CT_EXT2) - # np_eff matches the semantic size: CT_EXT has cell.np=1 → ext slot - # at index cell.np+1=2 (using 0-based; +1 in Julia), so the array - # used by the current code is (3,3,cell.np+2,cell.np+2)=(3,3,3,3). - # For CT_EXT1/EXT2 it's (3,3,cell.np+1,cell.np+1)=(3,3,4,4). - # Either way npp = cell.etype == CT_EXT ? cell.np + 1 : cell.np. + # CT_EXT carries one extra (asymptotic) slot after its Hermite pair. np_eff = cell.etype == CT_EXT ? cell.np + 1 : cell.np fill!(cell_mat_ext, 0) fill!(cell_rhs_ext, 0) @@ -688,12 +579,10 @@ function _assemble_and_solve!(ws::GalerkinWorkspace, end end - # Surface term at xmax for "resonant" method: the boundary contribution - # [αᵢ A uⱼ']_{x_max} from the IBP of the A-term in GWP2016 Eq. (32), kept - # here (not dropped) because the resonant cell carries the matching data. + # IBP surface term [αᵢ A uⱼ'] at x_max, dropped in GWP2016 Eq. 32 and kept here because the resonant cell carries Δ. res_cell = ws.cells[ws.nx] x_lsode = res_cell.x_lsode - ua_xmax, dua_xmax = _physical_ua_dua(cache, x_lsode) + ua_xmax, dua_xmax = physical_ua_dua(cache, x_lsode) Imat_xmax, _, _ = _physical_uv(params, Q, x_lsode) emap1 = res_cell.emap[1] if emap1 <= ws.ndim @@ -705,6 +594,34 @@ function _assemble_and_solve!(ws::GalerkinWorkspace, ws.rhs[emap1, 1] += surf_12 end + _parity_bcs_and_solve!(ws, params, Q, ws.ndim) +end + +# Add one cell's Hermite block to the base (isol = 1) banded matrix; both parities share it. +function _scatter_cell!(ws::GalerkinWorkspace, cell::GalerkinCell, cell_mat::Array{ComplexF64,4}, np_eff::Int) + mpert = 3 + offset = ws.kl + ws.kl + 1 + for ip in 0:np_eff, ipert in 1:mpert + i = cell.map[ipert, ip+1] + if i > ws.ndim + continue + end + for jp in 0:np_eff, jpert in 1:mpert + j = cell.map[jpert, jp+1] + if j > ws.ndim + continue + end + ws.mat[offset+i-j, j, 1] += cell_mat[ipert, jpert, ip+1, jp+1] + end + end +end + +# Copy the base system to both parities, impose the x = 0 parity conditions, and solve the +# leading n×n banded system of each parity by LAPACK banded LU. +function _parity_bcs_and_solve!(ws::GalerkinWorkspace, params::GGJParameters, Q::ComplexF64, n::Int) + mpert = 3 + offset = ws.kl + ws.kl + 1 + # Copy base matrix/rhs for both parities ws.mat[:, :, 2] .= ws.mat[:, :, 1] ws.rhs[:, 2] .= ws.rhs[:, 1] @@ -713,7 +630,7 @@ function _assemble_and_solve!(ws::GalerkinWorkspace, cell1 = ws.cells[1] Imat0, _, _ = _physical_uv(params, Q, cell1.x_left) - # Restore surface term at x=0 (IBP boundary) + # IBP surface term at x = 0 for ipert in 1:mpert, jpert in 1:mpert i = cell1.map[ipert, 1] # ip=0 j = cell1.map[jpert, 2] # ip=1 @@ -723,49 +640,36 @@ function _assemble_and_solve!(ws::GalerkinWorkspace, end end - # Apply parity BCs for each solution (isol=1: odd, isol=2: even). The - # inner equations have reflection symmetry about X=0 (GWP2016 Sec. III): - # one solution has even Ξ, Υ and odd Ψ, the other odd Ξ, Υ and even Ψ. - # Imposing both at X=0 gives the two parities that each contribute a Δ for - # outer-region matching (the Δ_odd, Δ_even of GWP2016 Eqs. 34–35). - # For each isol, build a modified local - # matrix for ip=0..1 of cell 1, then write it into the global matrix. + # Parity conditions at X = 0 (GWP2016 Sec. III), parity of the displacement Ξ (Eq. 34): + # isol = 1 even (Δ₊, interchange), isol = 2 odd (Δ₋, tearing). The ip = 0 rows are replaced. for isol in 1:2 # Zero out ip=0 rows in the global matrix for ipert in 1:mpert i = cell1.map[ipert, 1] # ip=0 DOFs if i > ws.ndim - continue end for jj in max(1, i-ws.kl):min(ws.ndim, i+ws.kl) ws.mat[offset+i-jj, jj, isol] = 0 end end - # Odd parity (isol=1): Ψ'(0)=0, Ξ(0)=0, Υ(0)=0 - # → row=Ψ(ip=0), col=Ψ(ip=1): A[map[1,1], map[1,2]] = 1 - # → row=Ξ(ip=0), col=Ξ(ip=0): A[map[2,1], map[2,1]] = 1 - # → row=Υ(ip=0), col=Υ(ip=0): A[map[3,1], map[3,1]] = 1 - # Even parity (isol=2): Ψ(0)=0, Ξ'(0)=0, Υ'(0)=0 - # → row=Ψ(ip=0), col=Ψ(ip=0): A[map[1,1], map[1,1]] = 1 - # → 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 + # Even: Ψ(0) = Ξ'(0) = Υ'(0) = 0; odd: Ψ'(0) = Ξ(0) = Υ(0) = 0 (ip = 1 slots are derivatives). if isol == 1 i = cell1.map[1, 1] - j = cell1.map[1, 2] + j = cell1.map[1, 1] ws.mat[offset+i-j, j, isol] = 1 for ipert in 2:3 i = cell1.map[ipert, 1] - j = cell1.map[ipert, 1] + j = cell1.map[ipert, 2] ws.mat[offset+i-j, j, isol] = 1 end else i = cell1.map[1, 1] - j = 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] - j = cell1.map[ipert, 2] + j = cell1.map[ipert, 1] ws.mat[offset+i-j, j, isol] = 1 end end @@ -778,35 +682,99 @@ function _assemble_and_solve!(ws::GalerkinWorkspace, end # Solve for each parity using LAPACK banded LU (gbtrf! + gbtrs!) - n = ws.ndim kl = ws.kl ku = kl + ab = @view ws.ab_buf[:, 1:n] + rb = @view ws.rhs_buf[1:n] for isol in 1:2 - copyto!(ws.ab_buf, @view(ws.mat[:, :, isol])) - copyto!(ws.rhs_buf, @view(ws.rhs[:, isol])) - _, ipiv = LinearAlgebra.LAPACK.gbtrf!(kl, ku, n, ws.ab_buf) - LinearAlgebra.LAPACK.gbtrs!('N', kl, ku, n, ws.ab_buf, ipiv, ws.rhs_buf) - ws.sol[:, isol] .= ws.rhs_buf + copyto!(ab, @view(ws.mat[:, 1:n, isol])) + copyto!(rb, @view(ws.rhs[1:n, isol])) + _, ipiv = LinearAlgebra.LAPACK.gbtrf!(kl, ku, n, ab) + LinearAlgebra.LAPACK.gbtrs!('N', kl, ku, n, ab, ipiv, rb) + ws.sol[1:n, isol] .= rb end end -# ----------------------------------------------------------------------- -# Top-level Galerkin solver. -# ----------------------------------------------------------------------- +""" + _assemble_and_solve_point!(ws, params, Q, θ, Us, Ub, E; nq=6) -> (Δraw, cexp) -# ----------------------------------------------------------------------- -# Diagnostic: reconstruct the physical inner-layer field profiles (Ψ, Ξ, Υ) -# from the solved Hermite DOFs over the bulk (Hermite) cells. The thin -# CT_EXT / CT_RES matching cells near xmax are skipped (their solution is -# carried by the asymptotic emap DOF, not the Hermite basis). Returns the -# inner-coordinate grid `x` and the three fields for both parities -# (columns: 1 = odd, 2 = even). Used to diagnose the resonant-layer -# structure vs resistivity. Not part of the production matching path. -# ----------------------------------------------------------------------- +Hermite Galerkin solve on the all-`CT_NONE` grid of `ws` along the ray x = e^{iθ}s, with the +far-field condition imposed strongly at the last node: its state is v = Us + Δ·Ub + c₁E₁ + c₂E₂, +with Us, Ub, E given as (Ψ, Ξ, Υ, Ψ', Ξ', Υ') in x-derivatives. The last node's six DOFs are +replaced by (Δ, c₁, c₂), tested in the same three directions. Returns the raw Δ and decaying-pair coefficients per parity +(isol = 1 even/Δ₊, 2 odd/Δ₋); `ws.sol` holds the full nodal solution afterwards. +""" +function _assemble_and_solve_point!(ws::GalerkinWorkspace, params::GGJParameters, Q::ComplexF64, + θ::Float64, Us::Vector{ComplexF64}, Ub::Vector{ComplexF64}, E::Matrix{ComplexF64}; nq::Int=6) + quad_nodes, quad_weights = gausslobatto(nq + 1) + offset = ws.kl + ws.kl + 1 + N = ws.nx + fill!(ws.mat, 0) + fill!(ws.rhs, 0) + cell_mat = ws.cell_mat_buf + for ix in 1:(N-1) + _gauss_quad!(cell_mat, ws.cells[ix], quad_nodes, quad_weights, params, Q; θ=θ) + _scatter_cell!(ws, ws.cells[ix], cell_mat, ws.cells[ix].np) + end + + # Last cell as a dense 12×12 block, local DOF index 3·ip + ipert (ip = 0..3 Hermite slots); the + # far-field state is imposed at its right node. + cell = ws.cells[N] + _gauss_quad!(cell_mat, cell, quad_nodes, quad_weights, params, Q; θ=θ) + Kc = zeros(ComplexF64, 12, 12) + for ip in 0:3, jp in 0:3, ipert in 1:3, jpert in 1:3 + Kc[3ip+ipert, 3jp+jpert] = cell_mat[ipert, jpert, ip+1, jp+1] + end + # IBP boundary term at the far end, mirroring the x = 0 term of `_parity_bcs_and_solve!` (I is x-independent). + Imat, _, _ = _physical_uv(params, Q, 0.0) + for ipert in 1:3, jpert in 1:3 + Kc[6+ipert, 9+jpert] -= Imat[ipert, jpert] + end + + # Far-field trace in DOF units (values, s-derivatives = e^{iθ}·x-derivatives); Ub normalized. + dof(v) = ComplexF64[v[1], v[2], v[3], cis(θ)*v[4], cis(θ)*v[5], cis(θ)*v[6]] + βb = norm(Ub) + Bd = hcat(dof(Ub) ./ βb, dof(E[:, 1]), dof(E[:, 2])) + usd = dof(Us) + # Tl maps the 9 reduced unknowns (left node's 6 DOFs, Δ̃, c₁, c₂) to the cell's 12 Hermite DOFs; test space = Tlᵀ. + Tl = zeros(ComplexF64, 12, 9) + for k in 1:6 + Tl[k, k] = 1 + end + Tl[7:12, 7:9] .= Bd + K9 = transpose(Tl) * Kc * Tl + r9 = -(transpose(Tl) * (Kc[:, 7:12] * usd)) + # (Δ̃, c₁, c₂) occupy the last node's value slots; its derivative slots drop out of the reduced system. + yslots = cell.map[:, 3] + glob = vcat(cell.map[:, 1], cell.map[:, 2], yslots) + for a in 1:9 + i = glob[a] + ws.rhs[i, 1] += r9[a] + for b in 1:9 + j = glob[b] + ws.mat[offset+i-j, j, 1] += K9[a, b] + end + end + + _parity_bcs_and_solve!(ws, params, Q, yslots[end]) + + Δraw = MVector{2,ComplexF64}(0, 0) + cexp = zeros(MMatrix{2,2,ComplexF64}) + for isol in 1:2 + y = ws.sol[yslots, isol] + Δraw[isol] = y[1] / βb + cexp[1, isol] = y[2] + cexp[2, isol] = y[3] + ws.sol[vcat(cell.map[:, 3], cell.map[:, 4]), isol] .= usd .+ Bd * y # reconstruct the last node for profiles + end + return SVector(Δraw), SMatrix(cexp) +end + +# Profiles (Ψ, Ξ, Υ) from the Hermite DOFs, columns 1 = even/Δ₊, 2 = odd/Δ₋; CT_EXT/CT_RES cells are +# skipped because their solution lives in the asymptotic emap DOF. function _solution_profile(ws::GalerkinWorkspace; npc::Int=10) mpert = 3 - # Preallocate the output arrays (every qualifying cell contributes npc points) instead of growing - # them with vcat inside the loop, which reallocates and copies the whole array each iteration (O(n²)). + # Every Hermite cell contributes npc points. npts = 0 for cell in ws.cells (cell.etype == CT_RES || cell.etype == CT_EXT) && continue @@ -844,25 +812,18 @@ function _solution_profile(ws::GalerkinWorkspace; npc::Int=10) return (; x=xs[p], Ψ=Ψ[p, :], Ξ=Ξ[p, :], Υ=Υ[p, :]) end -""" - solve_inner_profile(params::GGJParameters, γ::Number; kwargs...) - -Diagnostic variant of `solve_inner(GGJModel(:galerkin), ...)` that also returns the reconstructed -inner-layer field profiles. Returns `(Δ, Q, profile)` where `profile = (x, Ψ, Ξ, Υ)` (Ψ = inner ξ, -columns odd/even parity). Same numerics/kwargs as `solve_inner`; for investigating resonant-layer -structure vs resistivity, not the production matching path. -""" -function solve_inner_profile(params::GGJParameters, γ::Number; +# Stock Galerkin solve shared by `solve_inner` and `solve_inner_profile`. The matching data are the +# resonant cell's emap DOF, the coefficient of the large resonant solution (GWP2016 Eqs. 34–35), per +# parity isol = 1 (even, Δ₊), 2 (odd, Δ₋); returns them rescaled, with the solved workspace. +function _galerkin_solve(params::GGJParameters, γ::Number; kmax::Int=8, nx::Int=512, nq::Int=4, pfac::Float64=1.0, cutoff::Int=5, xfac::Float64=1.0, tol_res::Float64=1e-5) Q = inner_Q(params, γ) xmax_info = _xmax_3level(params, Q; kmax=kmax, xfac=xfac) - ws = _build_grid_and_workspace(nx, xmax_info.xmax, xmax_info.dx1, xmax_info.dx2, pfac, cutoff, nq) + ws = _build_grid_and_workspace(nx, xmax_info.xmax, xmax_info.dx1, xmax_info.dx2, pfac, cutoff) _assemble_and_solve!(ws, params, Q, xmax_info.cache; nq=nq, tol_res=tol_res) - res_cell = ws.cells[ws.nx] - emap1 = res_cell.emap[1] - Δ = rescale_delta(SVector{2,ComplexF64}(ws.sol[emap1, 2], ws.sol[emap1, 1]), params) - return Δ, Q, _solution_profile(ws), (; xmax=xmax_info.xmax, x0=x0(params), sfac=sfac(params)) + emap1 = ws.cells[ws.nx].emap[1] + return rescale_delta(SVector{2,ComplexF64}(ws.sol[emap1, 1], ws.sol[emap1, 2]), params), ws end """ @@ -874,7 +835,8 @@ 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(::_BareGGJ{:galerkin}, params::GGJParameters, γ::Number; kwargs...) - Δ, _, prof, _ = solve_inner_profile(params, γ; kwargs...) + Δ, ws = _galerkin_solve(params, γ; kwargs...) + prof = _solution_profile(ws) return (; Δ=Δ, x=prof.x, Ψ=prof.Ψ, Ξ=prof.Ξ, _profile_conversions(params)...) end @@ -888,91 +850,15 @@ Solve the GGJ inner-layer matching problem using the Hermite-cubic finite element (Galerkin) method (GWP2016 Sec. III), in the "resonant + noexp + inps" configuration. -Returns the parity-projected matching data (GWP2016 Eqs. 34–35) with the -`X₀^{2√(−D_I)}` physical rescaling applied, as an `InnerLayerResponse` whose -`tearing`/`interchange` fields are the isol=1/isol=2 element solutions -respectively — no parity swap (see the boundary-condition block above for the -parity derivation). -""" -function solve_inner(::_BareGGJ{:galerkin}, params::GGJParameters, γ::Number; - kmax::Int=8, nx::Int=512, nq::Int=4, pfac::Float64=1.0, - cutoff::Int=5, xfac::Float64=1.0, tol_res::Float64=1e-5) - Q = inner_Q(params, γ) +Returns the parity-projected matching data (GWP2016 Eqs. 34–35), rescaled by +[`rescale_delta`](@ref), as an `InnerLayerResponse`: `interchange` is the even-parity Δ₊ +(isol = 1), `tearing` the odd-parity Δ₋ (isol = 2). - # Build asymptotic cache and determine grid parameters - xmax_info = _xmax_3level(params, Q; kmax=kmax, xfac=xfac) - cache = xmax_info.cache - - # Build grid and workspace - ws = _build_grid_and_workspace(nx, xmax_info.xmax, xmax_info.dx1, xmax_info.dx2, - pfac, cutoff, nq) - - # Assemble and solve - _assemble_and_solve!(ws, params, Q, cache; nq=nq, tol_res=tol_res) - - # Extract the matching data from the resonant cell's emap DOF: the solved - # coefficient of the large resonant solution is Δ for each parity - # (GWP2016 Eqs. 34–35), here columns isol=1 (odd) and isol=2 (even). - res_cell = ws.cells[ws.nx] - emap1 = res_cell.emap[1] - Δ_raw = SVector{2,ComplexF64}(ws.sol[emap1, 1], ws.sol[emap1, 2]) - - # Rescaling is linear & diagonal; apply to the (tearing, interchange) - # pair directly, no parity swap. - Δ_rescaled = rescale_delta(Δ_raw, params) - - return InnerLayerResponse(Δ_rescaled[1], Δ_rescaled[2]) -end - - -# THIS IS AN ATTEMPT AT A CONVERGENCE-GUARDED INNER-LAYER SOLVE. IT PRESENTLY HAS SOME -# 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, max_levels=6, kwargs...) - -> (; delta, converged, err, kmax, xfac, nx, nlevels) - -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. +Knobs: `kmax` series order; `nx` cell count; `nq` Gauss–Lobatto order; `pfac` grid packing; +`cutoff` extension-cell span (cells nx−cutoff−1 … nx−1); `xfac` scale on the two inner matching radii; `tol_res` +resonant-cell quadrature tolerance. """ -function solve_inner_converged(model::GGJModel{:galerkin}, params::GGJParameters, γ::Number; - rtol::Float64=1e-2, kmax0::Int=12, xfac0::Float64=1.5, - cells_per_unit::Float64=3.0, nx_min::Int=1024, nx_max::Int=8192, - kmax_step::Int=2, kmax_max::Int=28, xfac_growth::Float64=1.5, - max_levels::Int=6, nq::Int=6, pfac::Float64=1.0, cutoff::Int=8, tol_res::Float64=1e-4) - Q = inner_Q(params, γ) - # nx sized so cells-per-unit-x ≈ cells_per_unit over [0, xmax]; gal_nx must be even. - nx_for(km, xf) = clamp(2 * cld(ceil(Int, cells_per_unit * _xmax_3level(params, Q; kmax=km, xfac=xf).xmax), 2), nx_min, nx_max) - function level(km, xf) - nx = nx_for(km, xf) - return solve_inner(model, params, γ; kmax=km, nx=nx, nq=nq, pfac=pfac, cutoff=cutoff, xfac=xf, tol_res=tol_res), nx - end - # Per real/imag component relative change with a significance floor (5% of |Δᵢ|). - function comperr(a, b) - e = 0.0 - for i in 1:2 - fl = 0.05 * abs(b[i]) - e = max(e, abs(real(a[i]) - real(b[i])) / (abs(real(b[i])) + fl)) - e = max(e, abs(imag(a[i]) - imag(b[i])) / (abs(imag(b[i])) + fl)) - end - return e - end - - kmax = kmax0 - xfac = xfac0 - Δ, nx = level(kmax, xfac) - err = Inf - for ℓ in 1:max_levels - kmax = min(kmax + kmax_step, kmax_max) - xfac *= xfac_growth - Δnew, nx = level(kmax, xfac) - err = comperr(Δnew, Δ) - Δ = Δnew - err < rtol && return (; delta=Δ, converged=true, err=err, kmax=kmax, xfac=xfac, nx=nx, nlevels=ℓ) - end - return (; delta=Δ, converged=false, err=err, kmax=kmax, xfac=xfac, nx=nx, nlevels=max_levels) +function solve_inner(::_BareGGJ{:galerkin}, params::GGJParameters, γ::Number; kwargs...) + Δ, _ = _galerkin_solve(params, γ; kwargs...) + return InnerLayerResponse(; interchange=Δ[1], tearing=Δ[2]) end diff --git a/src/InnerLayer/GGJ/InnerAsymptotics.jl b/src/InnerLayer/GGJ/InnerAsymptotics.jl index 1c752bf8d..3abcaced9 100644 --- a/src/InnerLayer/GGJ/InnerAsymptotics.jl +++ b/src/InnerLayer/GGJ/InnerAsymptotics.jl @@ -3,10 +3,10 @@ # Wasow asymptotic basis ("inps" basis) for the GGJ inner-layer system. # # Reference: A. H. Glasser & Z. R. Wang, Phys. Plasmas 27, 012506 (2020) -# ("GW2020"), Sections II.A–II.G (Eqs. 3–55). Notation map (all equation +# ("GW2020"), Secs. II.A–II.G (Eqs. 1–53) and III (Eqs. 54–55). Notation map (all equation # numbers are GW2020): # -# A0/A1/A2 -> Eqs. 3–5 (physical-system coefficient matrices A = A₀+x⁻²A₁+x⁻⁴A₂) +# A0/A1/A2 -> Eqs. 3–5 (A = A₀+x⁻²A₁+x⁻⁴A₂ (Eq. 3); A₀, A₂ Eq. 4; A₁ Eq. 5) # T, Tinv -> Eqs. 6–8 (Jordan basis of A₀; u = Tv, v' = xJv, λ = Q^{-1/2}) # J0/J1/J2 -> Eqs. 9–10 (Jₖ = T^{-1} Aₖ T) # P_k, B_k -> Eqs. 11–24 (splitting transformation V = PW; Lyapunov solve Eq. 22) @@ -15,17 +15,11 @@ # Y0, R -> Eqs. 48–49 (lowest-order solution; Mercier Frobenius exponents r±) # Y_k, Z_k -> Eqs. 50–52 (Y-series Y = Y₀Z in shifted-exponent form) # U(x) -> Eq. 53 (final asymptotic solution U = TPQSY at large x) -# residual -> Eq. 54 (convergence measure Δ± used to pick x_max) +# residual -> Eq. 54 (convergence measure Δ±, used to pick the matching radius) using LinearAlgebra using StaticArrays -# Index ranges that partition the 6-component first-order system into the -# 2-dimensional algebraic ("small") subspace and the 4-dimensional -# exponential ("large") subspace. -const _R1 = SVector(1, 2) -const _R2 = SVector(3, 4, 5, 6) - """ InnerAsymptoticsCache @@ -33,49 +27,29 @@ 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: +## Fields - `params, Q, kmax` — input parameters and series truncation order. - - `λ = 1/√Q` — complex scale factor used by the Wasow split. - `R = (r₊, r₋)` — Mercier-shifted Frobenius exponents at infinity (Eq. 49). - - `T, Tinv` — 6×6 eigenvector basis of A₀ (Eq. 7–8). - - `J` — `(J₀, J₁, J₂)`, the J-rotated coefficient matrices (Eq. 9–10). - - `P, B` — splitting matrices, k = 0..kmax+2 (Eqs. 16, 22). - - `K2` — 2×2 inner working matrices, k = 0..kmax+2 (Eq. 32; entry k=0 unused). - - `Qmat, Cmat` — 2×2 inner-block transformation matrices, k = 0..kmax+2 (Eqs. 32–39). - - `Dmat` — 2×2 shearing matrices, k = 0..kmax (Eq. 43). - - `Y0, Y0inv` — lowest-order Y matrix and its inverse (Eq. 48; exponents from Eq. 49). - - `Y, Z` — 2×2 series matrices, k = 0..kmax (Eq. 52). + - `T` — 6×6 eigenvector basis of A₀ (Eq. 7). + - `J` — `(J₀, J₁, J₂)`, the J-rotated coefficient matrices (Eqs. 9–10). + - `P` — splitting matrices, k = 0..kmax+2 (Eqs. 16, 22). + - `Qmat` — 2×2 inner-block transformation matrices, k = 0..kmax+2 (Eqs. 32–39). + - `Y` — 2×2 series matrices Y_k = Y₀Z_k, k = 0..kmax (Eqs. 50–52). """ struct InnerAsymptoticsCache params::GGJParameters Q::ComplexF64 kmax::Int - λ::ComplexF64 R::SVector{2,Float64} T::SMatrix{6,6,ComplexF64,36} - Tinv::SMatrix{6,6,ComplexF64,36} J::NTuple{3,SMatrix{6,6,ComplexF64,36}} - P::Vector{SMatrix{6,6,ComplexF64,36}} # k = 0..kmax+2 - B::Vector{SMatrix{6,6,ComplexF64,36}} # k = 0..kmax+2 - K2::Vector{SMatrix{2,2,ComplexF64,4}} # k = 0..kmax+2 - Qmat::Vector{SMatrix{2,2,ComplexF64,4}} # k = 0..kmax+2 - Cmat::Vector{SMatrix{2,2,ComplexF64,4}} # k = 0..kmax+2 - Dmat::Vector{SMatrix{2,2,ComplexF64,4}} # k = 0..kmax - Y0::SMatrix{2,2,ComplexF64,4} - Y0inv::SMatrix{2,2,ComplexF64,4} - Y::Vector{SMatrix{2,2,ComplexF64,4}} # k = 0..kmax - Z::Vector{SMatrix{2,2,ComplexF64,4}} # k = 0..kmax + P::Vector{SMatrix{6,6,ComplexF64,36}} + Qmat::Vector{SMatrix{2,2,ComplexF64,4}} + Y::Vector{SMatrix{2,2,ComplexF64,4}} end -# Helpers for cache indexing using physical (zero-based) order k. -@inline _at(v::Vector, k::Int) = v[k+1] - -# ----------------------------------------------------------------------- -# Build T, Tinv, A_0..A_2, J_0..J_2. The coefficient matrices A₀,A₁,A₂ are -# GW2020 Eqs. (3)–(5); the Jordan basis T, T⁻¹ are Eqs. (7)–(8); λ = Q^{-1/2} -# (Eq. 6); and J_i = T⁻¹A_iT are Eqs. (9)–(10). -# ----------------------------------------------------------------------- +# T, T⁻¹, λ = Q^{-1/2} and J_i = T⁻¹A_iT from the coefficient matrices A₀–A₂ (GW2020 Eqs. 3–10). function _build_tjmat(p::GGJParameters, Q::ComplexF64) e = ComplexF64(p.E) @@ -87,8 +61,7 @@ function _build_tjmat(p::GGJParameters, Q::ComplexF64) q2 = q * q λ = 1 / sqrt(q) - # T (Eq. 7) — built by column-major reshape; the listing below is - # row-major Julia order matching that layout. + # T (Eq. 7) and T⁻¹ (Eq. 8), listed row by row. #! format: off T = @SMatrix ComplexF64[ 1 0 h*q q2/λ h*q -q2/λ @@ -109,7 +82,7 @@ function _build_tjmat(p::GGJParameters, Q::ComplexF64) ] #! format: on - # A_0, A_1, A_2 — GW2020 Eqs. (4), (5), and the A₂ matrix of Eq. (3). Build mutable then freeze. + # A₀ and A₂: GW2020 Eq. (4); A₁: Eq. (5). A0 = zeros(ComplexF64, 6, 6) A1 = zeros(ComplexF64, 6, 6) A2 = zeros(ComplexF64, 6, 6) @@ -141,7 +114,6 @@ function _build_tjmat(p::GGJParameters, Q::ComplexF64) A2[5, 1] = h / q2 A2[6, 1] = -h * k * q - # Freeze A_i and build J_i = T^{-1} A_i T. A0s = SMatrix{6,6,ComplexF64}(A0) A1s = SMatrix{6,6,ComplexF64}(A1) A2s = SMatrix{6,6,ComplexF64}(A2) @@ -150,21 +122,11 @@ function _build_tjmat(p::GGJParameters, Q::ComplexF64) J1 = Tinv * A1s * T J2 = Tinv * A2s * T - return T, Tinv, (J0, J1, J2), λ + return T, (J0, J1, J2), λ end -# ----------------------------------------------------------------------- -# Closed-form Lyapunov solve for the splitting transformation — GW2020 Eq. (22), -# J₀²² P_k²¹ − P_k²¹ J₀¹¹ = −K_k²¹ and J₀¹¹ P_k¹² − P_k¹² J₀²² = −K_k¹², with the -# block-diagonal part B_k = K_k from Eq. (21). -# -# Given a 6×6 K, returns (B, P) such that: -# - B is block-diagonal with B[r1,r1] = K[r1,r1] and B[r2,r2] = K[r2,r2] -# - P has zero diagonal blocks, P[r1,r2] and P[r2,r1] given by the -# closed-form expressions below (which exploit the special structure -# of the lower-right block of J_0). -# ----------------------------------------------------------------------- - +# Closed-form solve of the GW2020 Eq. (22) Lyapunov equations for the off-diagonal blocks of P_k, +# using the structure of J₀'s lower-right block; B_k is the block-diagonal part of K_k (Eq. 21). function _lyap_solve(K::SMatrix{6,6,ComplexF64}, λ::ComplexF64) Bm = zeros(ComplexF64, 6, 6) Pm = zeros(ComplexF64, 6, 6) @@ -205,27 +167,19 @@ function _lyap_solve(K::SMatrix{6,6,ComplexF64}, λ::ComplexF64) return SMatrix{6,6,ComplexF64}(Bm), SMatrix{6,6,ComplexF64}(Pm) end -# ----------------------------------------------------------------------- -# Split recurrence — GW2020 Eqs. (14)–(19). Seeks descending power series -# J,P,B = Σ (·)_k x^{-2k} (Eq. 14) with B₀ = J₀, P₀ = I (Eq. 17), and solves the -# order-by-order recurrence Eq. (16)/(18) with K_k of Eq. (19) via the Lyapunov -# solve above. Builds B_k, P_k, K6_k for k = 0..kmax+2. -# ----------------------------------------------------------------------- - +# Splitting recurrence (GW2020 Eqs. 14–19): P_k, B_k for k = 0..kmax+2, each order by `_lyap_solve` of K_k. function _split_recurrence(J::NTuple{3,SMatrix{6,6,ComplexF64}}, λ::ComplexF64, kmax::Int) N = kmax + 3 # store k = 0..kmax+2 P = Vector{SMatrix{6,6,ComplexF64,36}}(undef, N) B = Vector{SMatrix{6,6,ComplexF64,36}}(undef, N) - K6 = Vector{SMatrix{6,6,ComplexF64,36}}(undef, N) Z6 = zero(SMatrix{6,6,ComplexF64}) P[1] = SMatrix{6,6,ComplexF64}(I) # P_0 = I B[1] = J[1] # B_0 = J_0 - K6[1] = J[1] for k in 1:(kmax+2) - kk = k + 1 # 1-based slot for K6/B/P at index k - # K6_k starts as J_k (only k = 1, 2 contribute), else 0 + kk = k + 1 # 1-based slot for B/P at index k + # K_k starts as J_k (only k = 1, 2 contribute), else 0 Kacc = (k <= 2) ? J[k+1] : Z6 # +2*(k-1)*P_{k-1} (k > 1 only) if k > 1 @@ -239,22 +193,15 @@ function _split_recurrence(J::NTuple{3,SMatrix{6,6,ComplexF64}}, λ::ComplexF64, end Kacc = Kacc - P[l+1] * B[kml+1] end - K6[kk] = Kacc Bk, Pk = _lyap_solve(Kacc, λ) B[kk] = Bk P[kk] = Pk end - return P, B, K6 + return P, B end -# ----------------------------------------------------------------------- -# Coefs recurrence — the 2×2 reduced-equation construction, GW2020 Eqs. (25)–(52). -# Builds the inner-block diagonalization Q_k, C_k via K2_k (Eqs. 31–32, 36), -# the shearing matrices D_k (Eqs. 42–43), the lowest-order Y₀ and Mercier -# exponents r± (Eqs. 48–49), and the shifted-exponent series Y_k = Y₀Z_k with -# Z_k from Eqs. (50)–(52). -# ----------------------------------------------------------------------- - +# Reduced 2×2 construction (GW2020 Eqs. 25–52): Q_k, C_k via K2_k, shearing D_k, lowest-order Y₀ +# with exponents r±, and the series Y_k = Y₀Z_k. function _coefs(B::Vector{SMatrix{6,6,ComplexF64,36}}, R::SVector{2,Float64}, kmax::Int) @@ -281,21 +228,18 @@ function _coefs(B::Vector{SMatrix{6,6,ComplexF64,36}}, end K2[k+1] = K2acc - # Q_k = [[0, 0], [-K2[1,1], -K2[1,2]]] (column-major reshape gives this layout) Qm[k+1] = @SMatrix ComplexF64[ 0 0 -K2acc[1, 1] -K2acc[1, 2] ] - # C_k = [[0, 0], [K2[2,1], K2[1,1] + K2[2,2]]] Cm[k+1] = @SMatrix ComplexF64[ 0 0 - K2acc[2, 1] K2acc[1, 1] + K2acc[2, 2] + K2acc[2, 1] K2acc[1, 1]+K2acc[2, 2] ] end # Build D_k for k = 0..kmax — GW2020 Eq. (43); c_2^21 = −(E+F+H+2) is Eq. (39). D = Vector{SMatrix{2,2,ComplexF64,4}}(undef, kmax + 1) - # D_0 = [[0, 1], [C_2[2,1], 3]] D[1] = @SMatrix ComplexF64[ 0 1 Cm[3][2, 1] 3 @@ -319,20 +263,17 @@ function _coefs(B::Vector{SMatrix{6,6,ComplexF64,36}}, r1 -1 ] - # E_k, Z_k, Y_k recurrence — GW2020 Eqs. (50)–(52): E_k = Y₀⁻¹D_kY₀ (Eq. 51), - # Z_k(R−2kI) − RZ_k = Σ E_l Z_{k−l} (Eq. 52), Y_k = Y₀Z_k (Eq. 50). + # Z_k and Y_k = Y₀Z_k from the GW2020 Eqs. (50)–(52) recurrence, with E_l = Y₀⁻¹D_lY₀ (Eq. 51). Y = Vector{SMatrix{2,2,ComplexF64,4}}(undef, kmax + 1) Z = Vector{SMatrix{2,2,ComplexF64,4}}(undef, kmax + 1) Z[1] = I2 Y[1] = Y0 for k in 1:kmax - # Build Z_k = sum_{l=1..k} E_l * Z_{k-l}, with E_l = Y0inv * D_l * Y0 Zacc = Z2 for l in 1:k El = Y0inv * D[l+1] * Y0 Zacc = Zacc + El * Z[k-l+1] end - # Divide each entry by (R[j] - R[i] - 2k) Zk = MMatrix{2,2,ComplexF64}(Zacc) for i in 1:2, j in 1:2 Zk[i, j] = Zk[i, j] / (R[j] - R[i] - 2 * k) @@ -341,51 +282,33 @@ function _coefs(B::Vector{SMatrix{6,6,ComplexF64,36}}, Y[k+1] = Y0 * Z[k+1] end - return K2, Qm, Cm, D, Y0, Y0inv, Y, Z + return Qm, Y end -# ----------------------------------------------------------------------- -# Public builder. -# ----------------------------------------------------------------------- - """ build_asymptotics(params::GGJParameters, Q::ComplexF64; kmax::Int=8) -> InnerAsymptoticsCache 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). +pick a matching radius with [`pick_smax`](@ref). """ function build_asymptotics(params::GGJParameters, Q::ComplexF64; kmax::Int=8) p1v = p1(params) # Mercier Frobenius exponents r± = 3/2 ± √(−D_I) — GW2020 Eq. (49). R = SVector{2,Float64}(p1v + 1.5, -p1v + 1.5) - T, Tinv, J, λ = _build_tjmat(params, Q) - P, B, _ = _split_recurrence(J, λ, kmax) - K2, Qm, Cm, D, Y0, Y0inv, Y, Z = _coefs(B, R, kmax) - - return InnerAsymptoticsCache( - params, Q, kmax, λ, R, T, Tinv, J, - P, B, K2, Qm, Cm, D, Y0, Y0inv, Y, Z - ) + T, J, λ = _build_tjmat(params, Q) + P, B = _split_recurrence(J, λ, kmax) + Qm, Y = _coefs(B, R, kmax) + return InnerAsymptoticsCache(params, Q, kmax, R, T, J, P, Qm, Y) end build_asymptotics(params::GGJParameters, Q::Number; kmax::Int=8) = build_asymptotics(params, ComplexF64(Q); kmax=kmax) -# ----------------------------------------------------------------------- -# Horner evaluator with optional fractional-power prefactor. Sums the -# descending power series in x^{-2k} (GW2020 Eq. 14 form) for the Y, Q, and -# P[r2,r1] coefficient blocks, with the x^{R/2} Mercier prefactor. -# -# Computes y[i] = (Σ_{k=0..n} c[i,k] * x^k) * x^rvec[i] -# and dy[i] = d/dx of the above. -# -# `c` is an n_components × (n+1) matrix; the second axis is k = 0..n, -# 1-indexed in Julia. -# ----------------------------------------------------------------------- - -function _horner(x::Real, c::AbstractMatrix{ComplexF64}; +# Horner sum y[i] = (Σ_{k=0..n} c[i,k+1] zᵏ)·z^rvec[i] and its z-derivative, z = 1/x² (GW2020 Eq. 14 form); +# rvec = −R/2 gives the x^R Mercier prefactor. c is n_components × (n+1); complex z uses the principal branch. +function _horner(x::Number, c::AbstractMatrix{ComplexF64}; rvec::Union{Nothing,AbstractVector{<:Real}}=nothing, derivative::Bool=false) nrows, ncols = size(c) @@ -439,9 +362,7 @@ function _horner(x::Real, c::AbstractMatrix{ComplexF64}; return y, dy end -# Pack the cache's Y, Qmat, and the [r2, r1] block of P into the row-major -# reshape order used by `_horner`. Returned matrices are -# n_components × (kmax+1). +# Pack Y, and Qmat with P[3:6, 1:2], column-major into `_horner`'s n_components × (kmax+1) layout. function _pack_y_coefs(cache::InnerAsymptoticsCache) kmax = cache.kmax cc = Matrix{ComplexF64}(undef, 4, kmax + 1) @@ -488,19 +409,19 @@ end # ----------------------------------------------------------------------- """ - evaluate_asymptotics(cache::InnerAsymptoticsCache, x::Real; derivative=true, apply_T=true) -> (U, dU) + evaluate_asymptotics(cache::InnerAsymptoticsCache, x::Number; derivative=true, apply_T=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 +Evaluate the inps asymptotic basis at `x`, real (`x > 0`) or complex on a +rotated contour (principal branch of `x^r`; the series is in 1/x², so valid for |arg x| < π/2). +Returns the 6×2 complex matrix `U` whose two columns are the power-like (Mercier) solutions of +the GGJ system (column 1 large, column 2 small), 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 `asymptotic_residual` for residual checks). The default `apply_T=true` -returns the solutions in the original 6-component first-order-system -basis used by `_physical_uv` and the shooting / Galerkin solvers. +If `apply_T=false`, the result is left in the J-rotated basis (used by `asymptotic_residual`). +The default `apply_T=true` returns it in the GW2020 Eq. (2) state basis, which +[`physical_ua_dua`](@ref) converts to (Ψ, Ξ, Υ). """ -function evaluate_asymptotics(cache::InnerAsymptoticsCache, x::Real; +function evaluate_asymptotics(cache::InnerAsymptoticsCache, x::Number; derivative::Bool=true, apply_T::Bool=true) xfac = 1.0 / (x * x) R = cache.R @@ -569,22 +490,42 @@ function evaluate_asymptotics(cache::InnerAsymptoticsCache, x::Real; end end -# ----------------------------------------------------------------------- -# Residual `delta` and adaptive xmax — the convergence measure Δ± of GW2020 -# Eq. (54), ‖u' − xJu‖/max(‖u'‖,‖xJu‖) in the J-rotated basis, used to choose -# x_max (GW2020 Sec. III, Fig. 3). -# ----------------------------------------------------------------------- +""" + physical_ua_dua(cache, x::Number) -> (ua, dua) +Convert the inps 6×2 basis at real or complex `x` to physical `(Ψ, Ξ, Υ)` +values and x-derivatives, each 3×2 (column 1 = large power solution, 2 = small). +The inps columns are in the GW2020 Eq. (2) state `u = (xΨ, Ξ, Υ, (xΨ)'/x, Ξ'/x, Υ'/x)`, +hence the /x and ·x factors; same map as the deltac `inpso_get_ua/dua` convention. """ - asymptotic_residual(cache::InnerAsymptoticsCache, x::Real) -> SVector{2,Float64} +function physical_ua_dua(cache::InnerAsymptoticsCache, x::Number) + U, dU = evaluate_asymptotics(cache, x; 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] / x + ua[2, j] = U[2, j] + ua[3, j] = U[3, j] + dua[1, j] = U[4, j] - U[1, j] / (x * x) + dua[2, j] = U[5, j] * x + dua[3, j] = U[6, j] * x + end + return ua, dua +end + +# Series residual, the analogue of GW2020 Eq. (54) for v' = xJv (Eq. 6), used to choose the matching +# radius (`pick_smax`, `_xmax_3level`; GW2020 Sec. III, Fig. 3). -Compute the convergence measure `Δ±` of the asymptotic basis at `x` for each -of the two algebraic columns (GW2020 Eq. 54). Returns +""" + asymptotic_residual(cache::InnerAsymptoticsCache, x::Number) -> SVector{2,Float64} + +Convergence measure of the asymptotic basis at real or complex `x` for each of the two algebraic +columns, the analogue of GW2020 Eq. 54 in the J-rotated basis. 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). """ -function asymptotic_residual(cache::InnerAsymptoticsCache, x::Real) +function asymptotic_residual(cache::InnerAsymptoticsCache, x::Number) U, dU = evaluate_asymptotics(cache, x; derivative=true, apply_T=false) xfac = 1.0 / (x * x) @@ -596,8 +537,6 @@ function asymptotic_residual(cache::InnerAsymptoticsCache, x::Real) M = M + xfac * xfac * cache.J[3] end - # matvec(:,:,1) = dU; matvec(:,:,2) = -x*M*U; matvec(:,:,0) = sum. - # delta(j) = ||matvec(:,j,0)||∞ / max(||matvec(:,j,1)||∞, ||matvec(:,j,2)||∞). matvec1 = dU matvec2 = -x * (M * U) matvec0 = matvec1 + matvec2 @@ -619,35 +558,36 @@ function asymptotic_residual(cache::InnerAsymptoticsCache, x::Real) end """ - pick_xmax(params, Q; eps=1e-7, kmax=8, xlogmin=-1.0, xlogmax=4.0, dxlog=0.01) -> (x_max, cache) + 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) -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. - -Throws an `ErrorException` if no `x` in the sweep range achieves the -target tolerance. +Sweep `s` log-uniformly upward and return the smallest `s` at which the series residual +(GW2020 Eq. 54) along `x = e^{iθ}s` drops below `eps`. If never reached, returns the +residual-minimizing `s` with `achieved = false`. """ -function 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) - cache = build_asymptotics(params, Q; kmax=kmax) - dxfac = 10.0^dxlog - xlog = xlogmin - x = 10.0^xlog - while xlog <= xlogmax - delta = asymptotic_residual(cache, x) - if maximum(delta) < eps - return x, cache +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 - xlog += dxlog - x *= dxfac + slog += dslog end - error("pick_xmax: failed to reach residual < $eps for x in [10^$xlogmin, 10^$xlogmax]") + return best_s, c, false end -pick_xmax(params::GGJParameters, Q::Number; kwargs...) = - pick_xmax(params, ComplexF64(Q); kwargs...) +pick_smax(params::GGJParameters, Q::Number; kwargs...) = + pick_smax(params, ComplexF64(Q); kwargs...) diff --git a/src/InnerLayer/GGJ/Ray.jl b/src/InnerLayer/GGJ/Ray.jl index 8440432f6..f7e1804d0 100644 --- a/src/InnerLayer/GGJ/Ray.jl +++ b/src/InnerLayer/GGJ/Ray.jl @@ -1,58 +1,16 @@ # 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). +# Rotated-contour backend (`GGJModel{:ray}`): the GGJ equations continued onto x = e^{iθ}s, +# θ = arg(Q)/4. Far-field series data at the radius S are marched inward to s_m modulo the +# decaying exponential pair, then the Galerkin.jl Hermite assembly solves [0, s_m] with +# v(s_m) = Us + Δ·Ub + c₁E₁ + c₂E₂. Method and rationale: docs/src/inner_layer.md. """ 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`; +Coefficient matrix `M(x)` of GWP2016 Eq. 11 as `dv/dx = M v`, `v = (Ψ, Ξ, Υ, Ψ', Ξ', Υ')`, for +real or complex `x`. `M` is polynomial in `x`, so x = 0 is an ordinary point and the system +continues analytically to complex `x`. `CT` selects the element type (default `ComplexF64`; `Complex{Double64}` for the extended-precision march). """ ode_matrix(p::GGJParameters, Q::ComplexF64, x::Number) = @@ -93,186 +51,44 @@ function ode_matrix(::Type{CT}, p::GGJParameters, Q::Number, x::Number) where {C 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} + ray_breaks(params, Q, θ, S; cfl=0.2, min_cells=1024) -> 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. +Cell boundaries `0 = s₀ < s₁ < … = S` for the Hermite BVP: each cell resolves the full +local WKB rate μ + ρs (fast Ohm and Υ-family pair μ = |√Q| + √(|Q|·|G+KF|), parabolic pair +ρ = Re(e^{2iθ}/√Q)) whether or not the mode is present, since unresolved dormant modes +alias into spurious discrete ones. `min_cells` floors the resolution where the local rates +are small (|Q| ≪ 1). """ -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) +function ray_breaks(params::GGJParameters, Q::ComplexF64, θ::Float64, S::Float64; + cfl::Float64=0.2, min_cells::Int=1024) sq = sqrt(Q) - ρ = real(cis(2θ) / sq) # parabolic pair: exponent ρ s²/2 + ρ = real(cis(2θ) / sq) ρ > 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 + hmax = S / min_cells 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) + h = min(cfl / (μ + ρ * breaks[end]), hmax) + snew = breaks[end] + h + push!(breaks, snew >= S - 0.25 * h ? S : snew) 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. + +# Far-field boundary data. The admissible solutions at large s are the power pair (u_small, u_big) +# from the inps series and the forward-decaying exponential pair E. The series is trusted only at +# S, which can lie far beyond s_m, so the power pair is marched S → s_m by 2-stage Radau IIA, +# R(z) = 2(3+z)/(z² − 4z + 6). Unresolved steps (z > 8, clear of |R| > 1 on 0 < z < 6) damp the +# exponential directions, which is harmless because Δ is defined modulo E. + +# Small and large power-like solutions at x as 6-component states (Ψ, Ξ, Υ, Ψ', Ξ', Υ'). +function _power_states(cache::InnerAsymptoticsCache, x::Number) + ua, dua = physical_ua_dua(cache, x) + return vcat(ua[:, 2], dua[:, 2]), vcat(ua[:, 1], dua[:, 1]) +end + +# Condition number of the column-normalized boundary basis [Us, Ub, E]. +_basis_cond(Us, Ub, E) = cond(hcat(Us ./ norm(Us), Ub ./ norm(Ub), E)) """ decaying_pair(params, Q, θ, S; growth=30.0, seed=20260707, dp_res=20.0) -> E (6×2) @@ -287,13 +103,7 @@ function decaying_pair(params::GGJParameters, Q::ComplexF64, θ::Float64, S::Flo ρ = 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. + # Window ΔS: the slowest backward grower separates from the power pair by `growth` e-folds. g2 = Inf gmax_loc = 0.0 for λ in eigvals(Matrix(ph * ode_matrix(params, Q, ph * S))) @@ -304,11 +114,7 @@ function decaying_pair(params::GGJParameters, Q::ComplexF64, θ::Float64, S::Flo 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). + # Resolve the true local spectrum; ρs + |√Q| under-counts the Υ-family rate √|Q(G+KF)|. for λ in eigvals(Matrix(ph * ode_matrix(params, Q, ph * (S + ΔS)))) gmax_loc = max(gmax_loc, abs(real(λ))) end @@ -338,14 +144,15 @@ function decaying_pair(params::GGJParameters, Q::ComplexF64, θ::Float64, S::Flo end """ - march_boundary(params, Q, θ, S, s_m, Us, Ub; rtol=1e-9, growth=30.0, seed=20260707) + march_boundary(params, Q, θ, S, s_m, Us, Ub; rtol=1e-9, ztrk=0.3, purge_budget=2.5, + growth=30.0, seed=20260707, dp_res=20.0) -> (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). +which Δ is defined): 2-stage Radau IIA with purges against a carried decaying frame. +`rtol` sets the step factor c = (rtol/100)^{1/3}; `ztrk` caps resolved-band steps at ztrk/g_max; +`purge_budget` is the e-folds between purges; `growth`, `seed`, `dp_res` go to [`decaying_pair`](@ref). """ function march_boundary(params::GGJParameters, Q::ComplexF64, θ::Float64, S::Float64, s_m::Float64, Us::Vector{ComplexF64}, Ub::Vector{ComplexF64}; @@ -361,67 +168,35 @@ function march_boundary(params::GGJParameters, Q::ComplexF64, θ::Float64, 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, :])) + # One Radau IIA step in the precision of W (`Mfun` returns the matrix in that precision): + # A = [5/12 -1/12; 3/4 1/4], c = (1/3, 1), b = A[2,:]. + function radau_step(Mfun, s, W, h) + CT = eltype(W) + hc = real(CT)(h) + Id = Matrix{CT}(I, 6, 6) + M1 = Mfun(s + h / 3) + M2 = Mfun(s + h) + Ablk = [Id-(5hc/12)*M1 (hc/12)*M1; -(3hc / 4)*M2 Id-(hc/4)*M2] + K = Ablk \ [M1 * W; M2 * W] + return W + hc * ((CT(3) / 4) * @view(K[1:6, :]) + (CT(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). + # Damped-zone steps have cond ~ z ≫ 1; their structured LU backward error mixes u_small into u_big + # by ~eps·ρS²/2·‖u_big‖/‖u_small‖, so that zone runs in Double64. 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. + # Step law, from the local growth rates g of 𝔐(s): + # damped zone (every grower has g·c·s > 8): h = −c·s with c = (rtol/100)^{1/3}, so the step + # count grows only like ln(S/s_m); + # resolved band: h = −min(c·s, ztrk/g_max), with the decaying frame F carried by the same + # steps and purged from W every `purge_budget` e-folds. + # On the real axis with imaginary Q the band entry lands on the turning point x² = −Q²(G+KF), + # where F cannot be identified; the default θ keeps it off the contour. 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. + # Damped zone, in Double64 from the first purge at S to the resolved-band entry. Wd = CTd.(W) let Floc = decaying_pair(params, Q, θ, S; growth=growth, seed=seed, dp_res=dp_res) Fd = CTd.(Floc) @@ -438,21 +213,19 @@ function march_boundary(params::GGJParameters, Q::ComplexF64, θ::Float64, 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) + 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). --- + # Resolved band, in Float64 with the carried frame. Φ = 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. + # Backward-growing rates (Re λ < 0); rates ≲ 10/s belong to the algebraic power pair. λs = eigvals(Matrix(𝔐(s))) gmin, gmax = Inf, 0.0 for λ in λs @@ -470,9 +243,9 @@ function march_boundary(params::GGJParameters, Q::ComplexF64, θ::Float64, have_F = true Φ = 0.0 end - W = radau_step(s, W, h) + W = radau_step(𝔐, s, W, h) if have_F - F = radau_step(s, F, h) + F = radau_step(𝔐, s, F, h) F = Matrix(qr(F).Q) Φ += gmax * abs(h) if Φ > purge_budget @@ -484,16 +257,9 @@ function march_boundary(params::GGJParameters, Q::ComplexF64, θ::Float64, 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). + # E must span the solutions that decay at infinity, which differ from the locally decaying pair at + # s_m once the mode families morph (s ~ |Q|); the carried F is that subspace. Local extraction is + # the fallback when the march never entered the resolved band. if have_F W .-= F * (F' * W) else @@ -503,11 +269,11 @@ function march_boundary(params::GGJParameters, Q::ComplexF64, θ::Float64, end """ - boundary_basis(params, Q, cache, θ, S; growth=30.0) + boundary_basis(params, Q, cache, θ, S; growth=30.0, seed=20260707, dp_res=20.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 +Boundary data at s = S taken directly from the inps series, used when the BVP edge is S itself +(s_m ≳ S): `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. @@ -515,272 +281,70 @@ 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]] + Us, Ub = _power_states(cache, cis(θ) * S) 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 + return Us, Ub, E, _basis_cond(Us, Ub, E) 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. +# Driver: series radius S, BVP edge s_m, march S → s_m, Hermite BVP on [0, s_m]. """ 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`. +Output of [`solve_ray`](@ref). + +## Fields + + - `Δ`: rescaled matching data (Δ₊, Δ₋) = (interchange, tearing), the even (isol = 1) and + odd (isol = 2) parities of the displacement Ξ. + - `Δraw`: the same pair in raw inps normalization, before `rescale_delta`. + - `cexp`: decaying-pair coefficients (c₁, c₂) at s_m, one column per parity. + - `Q`, `θ`: scaled eigenvalue and contour angle. + - `S`: series radius where the far-field data were taken. + - `series_ok`: whether `pick_smax` reached its residual tolerance. + - `bc_cond`: condition number of the column-normalized boundary basis [Us, Ub, E]. + - `breaks`: Hermite cell boundaries on [0, s_m]. + - `cfl`, `min_cells`: mesh resolution parameters used by [`ray_breaks`](@ref). + - `ws`: the full Galerkin workspace of the BVP (system buffers included); the nodal solution is `ws.sol`, both parities. """ struct RaySolveResult Δ::SVector{2,ComplexF64} Δraw::SVector{2,ComplexF64} - cexp::SMatrix{2,2,ComplexF64,4} # decaying-pair coefficients (col = isol) + cexp::SMatrix{2,2,ComplexF64,4} Q::ComplexF64 θ::Float64 S::Float64 - series_ok::Bool # pick_smax reached its tolerance - bc_cond::Float64 # boundary-basis condition number + series_ok::Bool + bc_cond::Float64 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 + cfl::Float64 + min_cells::Int + ws::GalerkinWorkspace 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) + solve_ray(params::GGJParameters, Q; θ=angle(Q)/4, kmax=12, S=nothing, smax_tol=1e-9, + cfl=0.2, min_cells=1024, nq=6, verbose=false, …) -> RaySolveResult - # 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) +Solve the GGJ inner-layer problem on the rotated ray x = e^{iθ}s. `Δ` is the rescaled pair +(Δ₊, Δ₋) = (interchange, tearing), directly comparable to the `:galerkin` backend; θ = 0 is a +real-axis solve. - # 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`). +Both boundary points are computed per solve. `S` defaults to the smallest radius where the series +residual is below `smax_tol` with order `kmax` ([`pick_smax`](@ref); it warns and uses the best +point if the tolerance is not reached). The BVP edge s_m is where every exponential family has +decayed 45 e-folds from the origin, from the local spectrum, scaled by `sm_fac`. `cfl`, +`min_cells` set the Hermite grid ([`ray_breaks`](@ref)) and `nq` its quadrature. `march_rtol` +is [`march_boundary`](@ref)'s `rtol`; `growth`, `seed`, `ztrk`, `purge_budget`, `dp_res` go to +`march_boundary`/[`decaying_pair`](@ref); `handoff` (diagnostic) supplies (Us, Ub) at S. """ function solve_ray(params::GGJParameters, Q::ComplexF64; - θ::Float64=angle(Q) / 4, kmax::Int=12, p::Int=12, + θ::Float64=angle(Q) / 4, kmax::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, + growth::Float64=30.0, march_rtol::Float64=1e-9, sm_fac::Float64=1.0, + cfl::Float64=0.2, min_cells::Int=1024, nq::Int=6, 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) @@ -792,17 +356,8 @@ function solve_ray(params::GGJParameters, Q::ComplexF64; 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: where every exponential family has decayed 45 e-folds from the origin, from the local spectrum. + # When s_m ≳ S (small |Q|) the series data are applied at S directly. s_m = let acc = 0.0, sprev = 0.0, out = Sv for sk in exp10.(range(-2, log10(Sv); length=400)) gdec = Inf @@ -821,28 +376,14 @@ function solve_ray(params::GGJParameters, Q::ComplexF64; 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 + if s_m < 0.98 * Sv + # `handoff` is a diagnostic: externally supplied series data (Us, Ub) at x = e^{iθ}S. + UsS, UbS = handoff === nothing ? _power_states(cache, cis(θ) * Sv) : handoff 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)) + bc_cond = _basis_cond(Us, 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; @@ -851,72 +392,15 @@ function solve_ray(params::GGJParameters, Q::ComplexF64; 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 + breaks = ray_breaks(params, Q, θ, Sbvp; cfl=cfl, min_cells=min_cells) + ws = _galerkin_workspace(breaks, fill(CT_NONE, length(breaks) - 1)) + Δraw, cexp = _assemble_and_solve_point!(ws, params, Q, θ, Us, Ub, E; nq=nq) + verbose && @printf(" BVP: %d Hermite cells on [0, %.4g], Δraw[1] = %.6e %+.6eim\n", + length(breaks) - 1, Sbvp, real(Δraw[1]), imag(Δraw[1])) - Δswapped = SVector{2,ComplexF64}(Δraw[2], Δraw[1]) # deltac_solve convention - Δ = rescale_delta(Δswapped, params) + Δ = rescale_delta(Δraw, params) - return RaySolveResult(Δ, SVector(Δraw), SMatrix(cexp), Q, θ, Sv, series_ok, - bc_cond, breaks, p, sols, resid, δΔ, round_used) + return RaySolveResult(Δ, Δraw, cexp, Q, θ, Sv, series_ok, bc_cond, breaks, cfl, min_cells, ws) end solve_ray(params::GGJParameters, Q::Number; kwargs...) = @@ -927,117 +411,57 @@ solve_ray(params::GGJParameters, Q::Number; kwargs...) = -> InnerLayerResponse Rotated-ray backend: converts γ via `inner_Q`, forwards keywords to -[`solve_ray`](@ref), and returns the parity-projected matching data. The -`solve_ray` pair is ordered `(Δ₁, Δ₂) = (interchange, tearing)`, so the two -channels are swapped into the named fields of [`InnerLayerResponse`](@ref). +[`solve_ray`](@ref), and returns the parity-projected matching data (Δ₊, Δ₋) as the +`interchange` and `tearing` fields of [`InnerLayerResponse`](@ref). Preferred for |Q| ≳ 1 and near the imaginary axis; use `solve_ray` directly -when the full [`RaySolveResult`](@ref) is wanted. +when the full [`RaySolveResult`](@ref) is wanted. `real_profiles` is a profile-only model +option (see `solve_inner_profile`) and is ignored here. """ -function solve_inner(::_BareGGJ{:ray}, params::GGJParameters, γ::Number; kwargs...) +function solve_inner(::_BareGGJ{:ray}, params::GGJParameters, γ::Number; real_profiles::Bool=false, kwargs...) res = solve_ray(params, inner_Q(params, γ); kwargs...) - return InnerLayerResponse(res.Δ[2], res.Δ[1]) + return InnerLayerResponse(; interchange=res.Δ[1], tearing=res.Δ[2]) end """ solve_inner_profile(::GGJModel{:ray}, params::GGJParameters, γ::Number; - npc=8, certify_rtol=1e-3, kwargs...) + real_profiles=false, 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). +Rotated-ray implementation of the [`solve_inner_profile`](@ref) interface. `Δ` and the +layer-center fields come from the solve at θ = arg(Q)/4; x = 0 lies on every contour, so +Ψ(0), Ξ(0) are the real-axis values, and by default only they are returned (`x = [0.0]`, +`certΔ = NaN`). With `real_profiles = true` the fields across the layer (x > 0) come from an +extra θ = 0 solve (`npc` points per cell), a diagnostic; the x = 0 row stays the rotated +solve's, and `certΔ` is the largest relative difference of the two solves' Δ₊, Δ₋, warned +above `certify_rtol`. Other keywords forward to [`solve_ray`](@ref), except θ. """ function solve_inner_profile(::_BareGGJ{:ray}, params::GGJParameters, γ::Number; - npc::Int=8, certify_rtol::Float64=1e-3, kwargs...) + real_profiles::Bool=false, npc::Int=8, certify_rtol::Float64=1e-3, kwargs...) Q = inner_Q(params, γ) - rr = solve_ray(params, Q; kwargs...) # certified Δ (optimal θ) + rr = solve_ray(params, Q; kwargs...) # Δ and layer-center fields (optimal θ) + c = solution_profile(rr; npc=1) # first point is s = 0 + if !real_profiles + return (; Δ=rr.Δ, x=[0.0], Ψ=c.Ψ[1:1, :], Ξ=c.Ξ[1:1, :], _profile_conversions(params)..., certΔ=NaN) + end 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Δ = maximum(abs.(r0.Δ .- rr.Δ) ./ max.(abs.(rr.Δ), 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) + prof.Ψ[1, :] .= c.Ψ[1, :] # x = 0 from the rotated solve, as without profiles + prof.Ξ[1, :] .= c.Ξ[1, :] 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. +Reconstruct (Ψ, Ξ, Υ) on the solution contour from the Hermite nodal solution, `npc` +points per cell; each field is `npts × 2` with columns (even, odd) as in `Δ`, and `x = e^{iθ}s`. +For rotated solves the fields live on the ray, not the real axis. """ 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, Ψ=Ψ, Ξ=Ξ, Υ=Υ) + prof = _solution_profile(res.ws; npc=npc) + return (; s=prof.x, x=cis(res.θ) .* prof.x, Ψ=prof.Ψ, Ξ=prof.Ξ, Υ=prof.Υ) end """ @@ -1067,28 +491,19 @@ function asymptotic_profile(params::GGJParameters, res::RaySolveResult, end """ - delta_convergence(params, Q; verbose=true, refine_tol=1e-9, - max_rounds=16, max_cells=12000, kwargs...) - -> (; Δ, spread, base, table) + delta_convergence(params, Q; verbose=true, 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. +Error bar for the matching data: re-solve with each numerical knob perturbed in turn (contour +angle, mesh, series order and radius, march tolerance, BVP edge s_m, decaying-pair separation) and +report the relative change of (Δ₊, Δ₋); `spread` is the worst case per parity. θ moves outward +only, away from the pseudo-resonance x² = −Q²(G+KF). """ -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...) +function delta_convergence(params::GGJParameters, Q::ComplexF64; verbose::Bool=true, kwargs...) + base = solve_ray(params, Q; kwargs...) variations = [ ("θ → 1.2θ", (; θ=(base.θ == 0 ? 0.1 : 1.2 * base.θ))), - ("p → p+4", (; p=base.p + 4)), + ("mesh ×2", (; cfl=(base.cfl / 2), min_cells=2 * base.min_cells)), ("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)) @@ -1096,19 +511,18 @@ function delta_convergence(params::GGJParameters, Q::ComplexF64; 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", + 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...) + r = solve_ray(params, Q; 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) + verbose && @printf(" %-16s δΔ₊ = %.2e δΔ₋ = %.2e\n", name, d1, d2) end - verbose && @printf(" %-16s δΔ₁ = %.2e δΔ₂ = %.2e\n", "WORST CASE", spread[1], spread[2]) + verbose && @printf(" %-16s δΔ₊ = %.2e δΔ₋ = %.2e\n", "WORST CASE", spread[1], spread[2]) return (; Δ=base.Δ, spread=SVector(spread), base=base, table=table) end diff --git a/src/InnerLayer/GGJ/RayAsymptotics.jl b/src/InnerLayer/GGJ/RayAsymptotics.jl deleted file mode 100644 index 630d5c397..000000000 --- a/src/InnerLayer/GGJ/RayAsymptotics.jl +++ /dev/null @@ -1,237 +0,0 @@ -# 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 db5bfefde..13168d3d9 100644 --- a/src/InnerLayer/GGJ/Reference.jl +++ b/src/InnerLayer/GGJ/Reference.jl @@ -7,18 +7,12 @@ 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. +**27**, 012506 (2020), Eq. 55 (coefficients transcribed verbatim), the benchmark of their +inps Wasow-basis convergence study (Figs. 1–3). At its operating point `Q = 1.234e-1` +(Fig. 1) the GGJ backends are cross-checked against the Fortran `rmatch deltac`; the +paper's `Δ_±` (Eq. 54) is a series-convergence norm, not matching data. -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 an -independent inner-layer reference, not this paper alone. - -Timescale parameters (taua, taur, v1) are set to canonical normalization; -callers should override them for physical cases. +`taua`, `taur`, `v1` default to a unit normalization; override them for physical cases. """ function glasser_wang_2020_eq55(; taua::Float64=1.0, taur::Float64=1e6, v1::Float64=1.0) return GGJParameters(; @@ -31,9 +25,8 @@ 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). +D_I ≈ −0.31166, p₁ = √(−D_I) ≈ 0.5583). Primary validation point for the +rotated-ray backend on the imaginary-Q axis (pinned at Q = 500i in the test suite). """ function q4_surface_benchmark() return GGJParameters(; diff --git a/src/InnerLayer/GGJ/Shooting.jl b/src/InnerLayer/GGJ/Shooting.jl index b83b2c369..aebe029da 100644 --- a/src/InnerLayer/GGJ/Shooting.jl +++ b/src/InnerLayer/GGJ/Shooting.jl @@ -1,27 +1,12 @@ # Shooting.jl # -# NOTE: This solver is NOT used in production and is not to be used. It was a -# porting exercise retained for reference only. The Galerkin solver -# (Galerkin.jl) is the sole supported GGJ inner-layer solver; shooting underflows -# at large |γ| and shares none of the production matching path. Do not call it. -# -# Stable backward shooting solver for the GGJ inner-layer model. -# Integrates the 4×4 origin-diagonalized -# resistive-layer ODE from `tmax` (large-x asymptotic regime) backward to -# a small `tmin` (origin Frobenius regime), then projects onto the local -# Frobenius basis at the origin and reads off the parity-projected -# matching data. -# -# This solves the same inner-region equations as the production Galerkin path -# (GWP2016 Eq. 11 ≡ GW2020 Eq. 1) but by a different numerical method: a -# backward shoot in the origin-diagonalized basis rather than the matched -# singular-Galerkin weak form. The 4 origin Frobenius exponents are -# `al0 = (p1, 1/2, −p1, −1/2)` where `p1 = √(−D_I)` (Mercier-stable required; -# D_I = E+F+H−1/4, GWP2016 Eq. A9). The two "small" (singular) modes are -# columns 3 and 4 (`x^{−p1}` and `x^{−1/2}`); the two "large" (smooth) modes -# are columns 1 and 2 (`x^{p1}` and `x^{1/2}`). -# -# Reference: A. H. Glasser, Phys. Plasmas **23**, 072505 (2016). +# Backward shooting solver for the GGJ inner layer, valid only for |Q| ≪ 1 (the large-t initial +# condition underflows at large |Q|); used for tests and backend comparisons. It follows the Fourier- +# transform method of Glasser, Jardin & Tesauro, Phys. Fluids 27, 1225 (1984), Sec. V: the 4×4 system +# in the transform variable t is diagonalized at the origin (exponents `al0 = (p1, 1/2, −p1, −1/2)`, +# p1 = √(−D_I), GJT Eqs. 49, 56), integrated from `tmax` on the decaying large-t solutions back to +# `tmin`, projected onto the origin Frobenius basis, and Δ± read off from the coefficient ratios +# (GJT Eqs. 90–91). using LinearAlgebra using SpecialFunctions: gamma @@ -32,23 +17,35 @@ using OrdinaryDiffEq Precomputed origin- and infinity-side arrays for the backward shoot. Construct via [`_build_shooting_system`](@ref). + +## Fields + + - `params`, `Q` - Layer parameters and scaled growth rate. + - `p1`, `pplus` - p₁ = √(−D_I) and p₁ − 1/2. + - `al0` - Origin Frobenius exponents (p₁, 1/2, −p₁, −1/2), the diagonal of the origin exponent matrix. + - `al1` - First-order coupling in the origin-diagonalized basis. + - `d0`, `d0inv` - Origin diagonalizing matrix and its inverse. + - `ups` - `(4, 4, nps)` origin power-series coefficients in t², `ups[:, :, 1] = I`. + - `nps` - Origin series length. + - `d1` - Diagonalizing matrix at large t. + - `pexp` - Logarithmic exponents of the large-t asymptotic basis. + - `bl1` - Coefficients ±√Q of t²/2 in the large-t exponentials. + - `tmax`, `tmin` - Shoot start (large t) and end (small t). """ struct GGJShootingSystem params::GGJParameters Q::ComplexF64 p1::Float64 pplus::Float64 - # origin side - al0::SVector{4,Float64} # diagonal of origin exponent matrix - al1::SMatrix{4,4,ComplexF64,16} # first-order coupling (already in al0-diag basis) - d0::SMatrix{4,4,ComplexF64,16} # diagonalizing matrix at origin + al0::SVector{4,Float64} + al1::SMatrix{4,4,ComplexF64,16} + d0::SMatrix{4,4,ComplexF64,16} d0inv::SMatrix{4,4,ComplexF64,16} - ups::Array{ComplexF64,3} # (4, 4, nps) origin power-series coeffs + ups::Array{ComplexF64,3} nps::Int - # infinity side - d1::SMatrix{4,4,ComplexF64,16} # diagonalizing matrix at infinity - pexp::SVector{4,ComplexF64} # log-exponents of asymptotic basis - bl1::SVector{4,ComplexF64} # ±λ exponents of t² in asymptotic basis + d1::SMatrix{4,4,ComplexF64,16} + pexp::SVector{4,ComplexF64} + bl1::SVector{4,ComplexF64} tmax::Float64 tmin::Float64 end @@ -149,10 +146,7 @@ function _build_origin_arrays(p::GGJParameters, Q::ComplexF64; nps::Int=8, rtol: return (; al0, al1=al1_s, d0=d0_s, d0inv=d0inv_s, ups, nps, p1v, pplus, tmin) end -# ----------------------------------------------------------------------- -# Infinity arrays — large-t asymptotic basis. -# With nxps = 1 the higher-order vps recurrence is dead code; vps[:,:,1] = I. -# ----------------------------------------------------------------------- +# Large-t asymptotic basis, leading order only. function _build_infinity_arrays(p::GGJParameters, Q::ComplexF64, d0inv::SMatrix{4,4,ComplexF64}; @@ -217,7 +211,7 @@ function _build_infinity_arrays(p::GGJParameters, Q::ComplexF64, d1 = SMatrix{4,4,ComplexF64}(d1_m) bl1 = SVector{4,ComplexF64}(-lamda, -lamda, lamda, lamda) - # tmax (ntmax = 1 default). + # tmax: the smaller of the asymptotic-validity bound and 6/|Q|. p1v = sqrt(-mercier_di(p)) tmax_inner = sqrt((max(0.5, p1v)^2 - log(rtol)) / abs(lamda)) tmax_outer = 6.0 / abs(q) @@ -272,15 +266,7 @@ function _origin_basis(sys::GGJShootingSystem, t::Real) return SMatrix{4,4,ComplexF64}(u) end -# ----------------------------------------------------------------------- -# Initial condition at tmax. -# With nxps = 1, vps[:,:,1] = I, so the asymptotic-basis -# evaluator returns a diagonal v. The IC is then -# y(tmax) = d1 * v(:, 1:2) -# which extracts the first two columns of d1 scaled by the algebraic-mode -# exponential factors. This selects the two "small" decaying modes of the -# large-t asymptotic basis. -# ----------------------------------------------------------------------- +# Initial condition at tmax: the two decaying modes of the large-t basis, d1[:, 1:2]·exp(bl1·t²/2 + pexp·ln t). function _initial_condition(sys::GGJShootingSystem) t = sys.tmax @@ -345,13 +331,8 @@ end Solve the GGJ inner-layer matching problem by stable backward shooting in the origin-diagonalized 4×4 basis. -Returns an `InnerLayerResponse(tearing, interchange)` with rescaling -applied. `_delta_from_c0` returns a `(Δ₁, Δ₂)` pair where `Δ₁` is the -**interchange** (anti-symmetric / W-odd) channel and `Δ₂` is the -**tearing** (symmetric / W-even) channel. We therefore map -`Δ₂ → tearing` and `Δ₁ → interchange` into the named fields, matching the -physics channel labels used by the Galerkin solver and by the -`InnerLayerResponse` docstring. +Returns an `InnerLayerResponse` rescaled by [`rescale_delta`](@ref), with `interchange` = Δ₊ +(Ξ even) and `tearing` = Δ₋ (Ξ odd), the same convention as the other GGJ backends. Tolerances `reltol`/`abstol` are the integrator tolerances; `rtol_origin` controls the truncation error of the origin Frobenius series and the @@ -377,8 +358,7 @@ function solve_inner(::_BareGGJ{:shooting}, params::GGJParameters, γ::Number; Δ_raw = _delta_from_c0(c0, sys) Δ_rescaled = rescale_delta(Δ_raw, params) - # Δ_rescaled ≡ (Δ₁, Δ₂) = (interchange, tearing). - return InnerLayerResponse(Δ_rescaled[2], Δ_rescaled[1]) + return InnerLayerResponse(; interchange=Δ_rescaled[1], tearing=Δ_rescaled[2]) end solve_inner(::_BareGGJ{:shooting}, params::GGJParameters, γ::Real; kwargs...) = diff --git a/src/InnerLayer/InnerLayer.jl b/src/InnerLayer/InnerLayer.jl index 7f47632e6..f8c3a30c3 100644 --- a/src/InnerLayer/InnerLayer.jl +++ b/src/InnerLayer/InnerLayer.jl @@ -1,9 +1,8 @@ # InnerLayer.jl # # Resistive inner-layer models for matched-asymptotic resistive/extended MHD stability. -# Provides an abstract `InnerLayerModel` interface and the GGJ (Glasser–Greene– -# Johnson) submodule. Future submodules (SLAYER, other inner layers) will plug in via the -# same interface. +# Provides an abstract `InnerLayerModel` interface and its implementations: the GGJ +# (Glasser–Greene–Johnson) and SLAYER submodules. module InnerLayer @@ -16,12 +15,8 @@ include("InnerLayerInterface.jl") include("GGJ/GGJ.jl") include("SLAYER/SLAYER.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 -import .GGJ: 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 +import .GGJ: GGJModel, GGJParameters, mercier_di, mercier_dr, inner_Q, rescale_delta +import .GGJ: solve_ray, pick_smax, delta_convergence, solution_profile, asymptotic_profile import .SLAYER: SLAYERModel, SLAYERParameters, slayer_parameters, r_based_shear import .SLAYER: riccati_del_s, slayer_layer_thickness, LayerWidths @@ -37,11 +32,8 @@ model_name(::SLAYERModel{S}) where {S} = "slayer_$S" export InnerLayerModel, InnerLayerParameters, InnerLayerResponse, 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 +export solve_ray, pick_smax, delta_convergence, solution_profile, asymptotic_profile export SLAYER, SLAYERModel, SLAYERParameters, slayer_parameters, r_based_shear export riccati_del_s, slayer_layer_thickness, LayerWidths diff --git a/src/InnerLayer/InnerLayerInterface.jl b/src/InnerLayer/InnerLayerInterface.jl index 52e1b36ea..86885d993 100644 --- a/src/InnerLayer/InnerLayerInterface.jl +++ b/src/InnerLayer/InnerLayerInterface.jl @@ -1,7 +1,7 @@ # InnerLayerInterface.jl # -# Abstract interface for resistive inner-layer models. Concrete models -# (e.g. GGJ, SLAYER, kinetic) live in submodules and specialize `solve_inner`. +# Abstract interface for resistive inner-layer models. Concrete models (GGJ, SLAYER) live in +# submodules and specialize `solve_inner`. """ InnerLayerModel @@ -35,25 +35,18 @@ the `Δ_{j,±}(γ)` of Glasser, Wang & Park, Phys. Plasmas **23**, 112506 # Fields - - `tearing` — the **odd-parity** matching coefficient (GWP Δ_+, the - "odd mode"). Corresponds to a flux perturbation W that is EVEN in x and - a velocity/temperature perturbation that is ODD — i.e., the - reconnecting mode with a current sheet at the rational surface. This is - the tearing drive that appears as Δ' in the classical constant-ψ - tearing equation. Must be populated by every resistive inner-layer model. - - - `interchange` — the **even-parity** matching coefficient (GWP Δ_−, the - "even mode"). Corresponds to W odd, N and Θ even — i.e., the - non-reconnecting interchange/ballooning channel. Its dissipative piece - in toroidal geometry is the Glasser, Greene & Johnson stabilization - term that opposes tearing growth (Glasser, Greene & Johnson 1975; - Lütjens-Bondeson-Roy 1993). Pressureless inner-layer models (e.g. - SLAYER's Fitzpatrick Riccati) set this identically zero. - -The naming follows the physics channel rather than a mathematical parity -label because `odd/even` carries different meanings across the literature -depending on whether you label by the parity of W (GWP paper convention) -or the parity of (N, Θ). Using `tearing` and `interchange` avoids ambiguity. + - `tearing` — odd-parity matching coefficient Δ₋: Ξ odd, Ψ even, the reconnecting channel + whose Δ is the constant-ψ tearing Δ′ (Pletzer, Bondeson & Dewar, J. Comput. Phys. 115, 530 + (1994), Eqs. 77, 86–87, 89). Every resistive inner-layer model populates it. + - `interchange` — even-parity matching coefficient Δ₊: Ξ even, Ψ odd, the non-reconnecting + interchange channel carrying the Glasser, Greene & Johnson (1975) stabilization. + Pressureless models (SLAYER) set it to zero. + +Parity is that of the displacement Ξ, matched to the outer ξ (GWP2016 Eq. 34; Pletzer, +Bondeson & Dewar Eqs. 32–33), with + even and − odd. The fields are named by the physics +channel because the parity of the flux Ψ is the opposite. Positional +pairs elsewhere in `InnerLayer` (e.g. the `Δ` of [`solve_inner_profile`](@ref)) are +ordered (Δ₊, Δ₋) = (interchange, tearing). """ struct InnerLayerResponse tearing::ComplexF64 @@ -66,13 +59,14 @@ InnerLayerResponse(; tearing::Number=0, interchange::Number=0) = """ solve_inner(model::InnerLayerModel, params, γ::Number; kwargs...) -> InnerLayerResponse -Compute the parity-projected matching data `(Δ_tearing, Δ_interchange)` for +Compute the parity-projected matching data (the `tearing` and `interchange` fields of +`InnerLayerResponse`) for the given inner-layer `model`, physical parameters `params`, and complex growth rate `γ`. Concrete models specialize this function. See `InnerLayerResponse` for the physics-oriented field definitions. Pressureless models (SLAYER) populate only `tearing` and leave -`interchange` at zero; two-fluid / finite-β models (GGJ) populate both. +`interchange` at zero; finite-pressure models (GGJ) populate both. """ function solve_inner end @@ -85,15 +79,16 @@ 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) + - `Δ` — the same matching data as [`solve_inner`](@ref), as the pair + (Δ₊, Δ₋) = (interchange, tearing) - `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 + `x ≥ 0` with the rational surface at `x = 0` + - `Ψ`, `Ξ` — `length(x) × 2` profiles, columns in the order of `Δ`, 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) + 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. diff --git a/src/InnerLayer/SLAYER/LayerInputs.jl b/src/InnerLayer/SLAYER/LayerInputs.jl index 949f02963..7e7fb4175 100644 --- a/src/InnerLayer/SLAYER/LayerInputs.jl +++ b/src/InnerLayer/SLAYER/LayerInputs.jl @@ -23,8 +23,8 @@ using FastInterpolations: DerivOp, integrate, cubic_interp, cumulative_integrate Minor radius at normalized flux `psi` and poloidal angle `theta`, computed from `equil.rzphi_rsquared` as `√((R − R₀)² + (Z − Z₀)²)`. -`theta = 0.0` (outboard midplane) is the default; pass `θ = π` to measure -the inboard side if you want an average. +`theta` is normalized to [0, 1): `0.0` (outboard midplane) is the default and +`0.5` is the inboard midplane. """ function surface_minor_radius(equil, psi::Real; theta::Real=0.0) r_sq = equil.rzphi_rsquared((Float64(psi), Float64(theta))) @@ -183,12 +183,16 @@ profiles, without an intermediate file round-trip. physical field, so passing it as a scalar is almost always wrong. - `mu_i` -- ion mass in proton-mass units (default `2.0` for D). + - `zeff` -- effective charge. Scalar or a callable of `psi`, evaluated per surface (default `1.0`). + - `chi_perp` -- perpendicular heat diffusivity [m²/s]. Scalar or a callable of `psi` (default `1.0`). + - `chi_tor` -- toroidal heat diffusivity [m²/s]. Scalar or a callable of `psi` (default `1.0`). + - `dr_val` -- resistive interchange index `D_R = E + F + H²` (Glasser-Greene-Johnson 1975) feeding the critical-Δ formulas (`:lar`, `:fitzpatrick`, `:toroidal`). When `nothing` (default), Julia @@ -203,6 +207,7 @@ profiles, without an intermediate file round-trip. NOT the Mercier index `D_I = E + F + H − 1/4`. The two differ by `(H − 1/2)²`, which is non-trivial on shaped equilibria (~factor 3 on DIII-D); this code uses the physically correct `D_R`. + - `dgeo_val` -- Connor 2015 (PPCF 57 065001) Eq. 59 geometric factor used by `dc_type=:toroidal`. When `nothing` (default), an error is raised if `dc_type=:toroidal` is also requested — the auto-derived @@ -210,21 +215,26 @@ profiles, without an intermediate file round-trip. doesn't currently expose. Pass a scalar / vector / callable to use a prescribed value. (For `dc_type=:fitzpatrick` and `:lar`, dgeo_val is not consulted.) + - `dc_type` -- `:none` (default), `:lar`, `:fitzpatrick`, or `:toroidal`. + - `rs_method` -- radial label defining `r_s` for the whole layer stack: `:midplane` (default), `:halfwidth`, `:fsa`, `:volume`, or `:flux`. See [`radial_label`](@ref) for the definitions. S, the r-based shear, W_d, and the Δ' reference-length factor `k_ref` all follow the choice together, so every option is self-consistent. Not exposed via TOML — programmatic use only. + - `theta` -- poloidal angle at which to measure minor radius (default `0.0`, outboard midplane). + - `resistivity_model` -- `SauterNeoModel()` (default), `RedlNeoModel()`, `SpitzerModel()`, or `SpitzerHarmModel()` (legacy Fitzpatrick σ_∥). Sets the η entering τ_R = μ₀r_s²/η. With a neoclassical model, `f_trap` and ν*_e are taken from the surface's `ResistGeometry` if populated (via `ForceFreeStates.resist_eval_all!`), otherwise fall back to the ε-only Lin-Liu-Miller form and `rs/R_0` aspect ratio. + - `lnLambda_form` -- Coulomb-log form passed through to `slayer_parameters` (default `:nrl`; `:wesson` + `SpitzerHarmModel()` reproduces legacy SLAYER exactly). @@ -370,9 +380,10 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; 1.0 end alpha_k = if rg === nothing - @warn("build_slayer_inputs: sing.restype not populated; using the " * - "slab Mercier exponent α = 1/2 for the Δ' reference-length " * - "conversion at all such surfaces.", maxlog=1) + @warn( + "build_slayer_inputs: sing.restype not populated; using the " * + "slab Mercier exponent α = 1/2 for the Δ' reference-length " * + "conversion at all such surfaces.", maxlog=1) 0.5 else sqrt(max(-(rg.E + rg.F + rg.H - 0.25), 0.0)) diff --git a/src/InnerLayer/SLAYER/Riccati.jl b/src/InnerLayer/SLAYER/Riccati.jl index 87781259b..8f4e518c6 100644 --- a/src/InnerLayer/SLAYER/Riccati.jl +++ b/src/InnerLayer/SLAYER/Riccati.jl @@ -261,7 +261,7 @@ function solve_inner(::SLAYERModel{:fitzpatrick}, Δ::ComplexF64 = π / dW_end # Fitzpatrick / pressureless SLAYER has no interchange channel - # (the Δ_− / even-parity matching quantity is identically zero in + # (the Δ₊ / even-parity matching quantity is identically zero in # the pressureless limit), so populate only the tearing field. return InnerLayerResponse(Δ, zero(ComplexF64)) end diff --git a/src/Tearing/Dispersion/CoupledFullMatch.jl b/src/Tearing/Dispersion/CoupledFullMatch.jl index 9b0f136c9..3ae55d515 100644 --- a/src/Tearing/Dispersion/CoupledFullMatch.jl +++ b/src/Tearing/Dispersion/CoupledFullMatch.jl @@ -4,7 +4,7 @@ # the m inner-layer resonances decoupled via the matching-identity rows # # C^j_L = d^j_+ − d^j_- -# C^j_R = -(d^j_+ + d^j_-) +# C^j_R = d^j_+ + d^j_- # # (see Wang-Glasser-Brennan-Liu-Park 2020, Phys. Plasmas **27**, 122503, # Eq. (11a)-(11d) and Glasser-Wang-Park 2016, Phys. Plasmas **23**, 112506, @@ -155,7 +155,7 @@ function (mc::MultiSurfaceCouplingFull)(Q::Number) # Allocate the matching matrix and fill the lower-left 2m × 2m block # with transpose(dp_raw[1:s2, 1:s2]). mat = zeros(ComplexF64, s4, s4) - @views mat[s2+1:s4, 1:s2] .= transpose(mc.dp_raw[1:s2, 1:s2]) + @views mat[(s2+1):s4, 1:s2] .= transpose(mc.dp_raw[1:s2, 1:s2]) # Per-surface inner-layer assembly @inbounds for k in 1:m @@ -171,7 +171,7 @@ function (mc::MultiSurfaceCouplingFull)(Q::Number) Q_k = Qc * (ref_tauk / sc.tauk) + 1im * mc.ntor * mc.rotation[k] resp = solve_inner(sc.model, sc.params, Q_k; mc.inner_kwargs...) - # delta1 = interchange (parity −), delta2 = tearing (parity +); named + # delta1 = interchange (parity +), delta2 = tearing (parity −); named # fields expose the two channels directly. sc.scale converts inner-basis # Δ to outer units (1.0 for GGJ since rescale_delta is applied inside # solve_inner; S^(1/3) for SLAYER). @@ -195,15 +195,15 @@ function (mc::MultiSurfaceCouplingFull)(Q::Number) # --- Upper-right 2×2 block: matching identity --- # C^k_L = d^k_+ − d^k_- ⇒ mat[idx1,idx3]=-1, mat[idx1,idx4]=+1 - # C^k_R = -(d^k_+ + d^k_-) ⇒ mat[idx2,idx3]=-1, mat[idx2,idx4]=-1 + # C^k_R = d^k_+ + d^k_- ⇒ mat[idx2,idx3]=-1, mat[idx2,idx4]=-1 mat[idx1, idx3] = -1 mat[idx1, idx4] = 1 mat[idx2, idx3] = -1 mat[idx2, idx4] = -1 # --- Lower-right 2×2 block: inner Δ matching --- - # d^k_+ eqn: -Δ_int·d^k_+ + Δ_tear·d^k_- + (outer D' terms) = 0 - # d^k_- eqn: -Δ_int·d^k_+ - Δ_tear·d^k_- + (outer D' terms) = 0 + # L-side row: -Δ_int·d^k_+ + Δ_tear·d^k_- + (outer D' terms) = 0 (GWP2016 Eq. 37) + # R-side row: -Δ_int·d^k_+ - Δ_tear·d^k_- + (outer D' terms) = 0 mat[idx3, idx3] = -delta1 mat[idx3, idx4] = delta2 mat[idx4, idx3] = -delta1 diff --git a/src/Tearing/Runner/run_slayer.jl b/src/Tearing/Runner/run_slayer.jl index f50f4101e..93f67163c 100644 --- a/src/Tearing/Runner/run_slayer.jl +++ b/src/Tearing/Runner/run_slayer.jl @@ -178,7 +178,7 @@ the `S^(1/3)` Δ(Q) scale, both built on r_s). Near surface `k` the outer solution combines the large and small Frobenius solutions with the Mercier exponent `α = √(−D_I)`, `D_I = E + F + H − 1/4` -(Glasser, Greene & Johnson 1975, Eq. 48; STRIDE's `alpha`). The displacement +(Glasser, Greene & Johnson 1975, Eqs. 16, 24; STRIDE's `alpha`). The displacement scales as `|x|^(−1/2 ± α)` (Glasser, Wang & Park 2016, Eq. 26), so the flux-like variable `ψ ∝ x·ξ` of the slab Δ' definition scales as `A_L·|x|^(1/2−α) + A_S·|x|^(1/2+α)`. Rescaling the radial variable @@ -387,7 +387,7 @@ end Evaluate both parity channels of the GGJ inner-layer matching data at the complex normalized growth rate `Q` for each rational surface. Returns a vector of `(ising, tearing, interchange)` named tuples, where `tearing` -(GWP Δ_+) is the reconnecting channel and `interchange` (GWP Δ_−) the +(GWP Δ₋, odd) is the reconnecting channel and `interchange` (GWP Δ₊, even) the non-reconnecting Glasser-stabilization channel (see `InnerLayerResponse`). This is the supported GGJ diagnostic: both parity Δ's are physical and diff --git a/test/runtests_dispersion_coupled.jl b/test/runtests_dispersion_coupled.jl index 1b950016e..19ab7c99d 100644 --- a/test/runtests_dispersion_coupled.jl +++ b/test/runtests_dispersion_coupled.jl @@ -1,6 +1,7 @@ @testset "Dispersion coupled determinant" begin using GeneralizedPerturbedEquilibrium.InnerLayer using GeneralizedPerturbedEquilibrium.InnerLayer: InnerLayerModel, solve_inner + using GeneralizedPerturbedEquilibrium.InnerLayer.GGJ: glasser_wang_2020_eq55 using GeneralizedPerturbedEquilibrium.Dispersion using LinearAlgebra using StaticArrays diff --git a/test/runtests_dispersion_residual.jl b/test/runtests_dispersion_residual.jl index 30ebf24ff..69d9cec08 100644 --- a/test/runtests_dispersion_residual.jl +++ b/test/runtests_dispersion_residual.jl @@ -1,6 +1,7 @@ @testset "Dispersion residual (SurfaceCoupling)" begin using GeneralizedPerturbedEquilibrium.InnerLayer using GeneralizedPerturbedEquilibrium.InnerLayer: InnerLayerModel, solve_inner + using GeneralizedPerturbedEquilibrium.InnerLayer.GGJ: glasser_wang_2020_eq55 using GeneralizedPerturbedEquilibrium.Dispersion using StaticArrays diff --git a/test/runtests_innerlayer.jl b/test/runtests_innerlayer.jl index 9aec5395b..6cda60780 100644 --- a/test/runtests_innerlayer.jl +++ b/test/runtests_innerlayer.jl @@ -19,7 +19,7 @@ const IL = GeneralizedPerturbedEquilibrium.InnerLayer const GGJ = IL.GGJ @testset "InnerLayer GGJ (Glasser & Wang 2020, Eq. 55)" begin - p = IL.glasser_wang_2020_eq55() + p = IL.GGJ.glasser_wang_2020_eq55() @testset "Mercier index and matching powers (Eq. 49)" begin # D_I = E + F + H − 1/4 from the verbatim Eq. 55 coefficients. @@ -52,11 +52,10 @@ const GGJ = IL.GGJ 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 "InnerLayer GGJ :ray backend (rotated contour, Hermite BVP)" begin + # These tests pin the implementation; the method is validated against an independent + # strict reference over the Q plane outside the test suite. + p = IL.GGJ.glasser_wang_2020_eq55() @testset "agrees with :galerkin at the paper point Q = 0.1234" begin γ = 0.1234 * GGJ.q0(p) @@ -71,7 +70,7 @@ end end @testset "q4 physical benchmark at Q = 500i (regime beyond :galerkin)" begin - q4 = IL.q4_surface_benchmark() + q4 = IL.GGJ.q4_surface_benchmark() γ = 500.0im * GGJ.q0(q4) Δ = IL.solve_inner(IL.GGJModel(), q4, γ) # Pins from the pre-port validation suite (post extended-precision @@ -90,15 +89,14 @@ 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 + p = IL.GGJ.q4_surface_benchmark() + + @testset "Hermite BVP resolution floor at small |Q|" begin + # At |Q| ≪ 1 the local WKB rates are small; refining the grid (cfl/8, 2× min_cells) must not move Δ. + Q = 0.01 + 0.0im + r = IL.solve_ray(p, Q) + r2 = IL.solve_ray(p, Q; cfl=0.025, min_cells=2048) + @test maximum(abs.(r2.Δraw .- r.Δraw) ./ abs.(r.Δraw)) < 1e-6 end @testset "ode_matrix: ordinary point and type-generic build" begin @@ -112,11 +110,6 @@ end @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 @@ -140,7 +133,7 @@ end Q = 5.0im conv = IL.delta_convergence(p, Q; verbose=false) Δ = IL.solve_inner(IL.GGJModel(; solver=:ray), p, Q * GGJ.q0(p)) - # conv.Δ is the raw solve_ray pair (Δ₁, Δ₂) = (interchange, tearing). + # conv.Δ is the solve_ray pair (Δ₊, Δ₋) = (interchange, tearing). @test conv.Δ[1] ≈ Δ.interchange rtol = 1e-6 # baseline == the plain solve @test conv.Δ[2] ≈ Δ.tearing rtol = 1e-6 @test maximum(conv.spread) < 1e-4 # honest error bar is small here @@ -148,12 +141,12 @@ end end @testset "solve_inner_profile interface (matching-driver contract)" begin - p = IL.glasser_wang_2020_eq55() + p = IL.GGJ.glasser_wang_2020_eq55() γ = 0.1234 * GGJ.q0(p) - for model in (IL.GGJModel(; solver=:ray), IL.GGJModel(; solver=:galerkin)) + for model in (IL.GGJModel(; solver=:ray, real_profiles=true), IL.GGJModel(; solver=:galerkin)) prof = IL.solve_inner_profile(model, p, γ) # Δ agrees with the plain matching solve of the same backend (identical solve path). - # prof.Δ stays the positional pair (Δ₁, Δ₂) = (interchange, tearing) that feeds deltar. + # prof.Δ is the positional pair (Δ₊, Δ₋) = (interchange, tearing) that feeds deltar. Δ_plain = IL.solve_inner(model, p, γ) @test prof.Δ[1] ≈ Δ_plain.interchange rtol = 1e-12 @test prof.Δ[2] ≈ Δ_plain.tearing rtol = 1e-12 @@ -162,13 +155,25 @@ end @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]) + # Parity at the layer center: Ψ(0) = 0 in the even (Δ₊) column, Ψ(0) ≠ 0 in the odd (Δ₋) one. + @test abs(prof.Ψ[1, 1]) < 1e-6 * abs(prof.Ψ[1, 2]) # 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, γ) + ray = IL.solve_inner_profile(IL.GGJModel(; solver=:ray, real_profiles=true), p, γ) @test ray.certΔ < 1e-12 + + # Default :ray returns only the layer-center fields of the rotated solve. With real profiles the + # x = 0 row is still the rotated solve's, so the penetrated field does not depend on the option, + # and the θ = 0 profile next to it is continuous with it. + γi = 5.0im * GGJ.q0(p) + c = IL.solve_inner_profile(IL.GGJModel(; solver=:ray), p, γi) + full = IL.solve_inner_profile(IL.GGJModel(; solver=:ray, real_profiles=true), p, γi) + @test c.x == [0.0] && size(c.Ψ) == (1, 2) && isnan(c.certΔ) + @test c.Δ == full.Δ + @test full.x[1] == 0 && full.Ψ[1, :] == c.Ψ[1, :] && full.Ξ[1, :] == c.Ξ[1, :] + @test full.Ψ[2, 2] ≈ c.Ψ[1, 2] rtol = 1e-2 # odd Ψ is even, so smooth through x = 0 + @test full.certΔ < 1e-5 end diff --git a/test/runtests_matching_models.jl b/test/runtests_matching_models.jl index 920e03d4c..d7410d618 100644 --- a/test/runtests_matching_models.jl +++ b/test/runtests_matching_models.jl @@ -11,7 +11,7 @@ @test !FFS.closure_capable(IL.GGJModel(; solver=:shooting)) # no layer profiles @test IL.GGJModel(; solver=:galerkin, nx=640).options == (nx=640,) # Carried options reach the backend exactly as call-site keywords do. - p = IL.glasser_wang_2020_eq55() + p = IL.GGJ.glasser_wang_2020_eq55() γ = 1e-3im @test IL.solve_inner(IL.GGJModel(; solver=:galerkin, nx=256), p, γ) == IL.solve_inner(IL.GGJModel(; solver=:galerkin), p, γ; nx=256) @test IL.solve_inner(IL.GGJModel(; solver=:galerkin, nx=256), p, γ; nx=512) == IL.solve_inner(IL.GGJModel(; solver=:galerkin), p, γ) diff --git a/test/runtests_solve_api.jl b/test/runtests_solve_api.jl index 18d4e6082..3de6f581f 100644 --- a/test/runtests_solve_api.jl +++ b/test/runtests_solve_api.jl @@ -244,6 +244,29 @@ using TOML end end + @testset "matched inner profiles: opt-in, and parity columns agree across inner models" begin + mktempdir() do dir + gal = solve(equil, Galerkin(; nx=32, rpec_flag=true, cut_solution=true); nn=1, dir_path=dir, ffs_kwargs...) + prob = MatchProblem(gal; eta=1e-7, rho=1e-7, rotation=1.0) + center = solve(prob, GGJModel()).match + ray = solve(prob, GGJModel(; real_profiles=true)).match + galm = solve(prob, GGJModel(; solver=:galerkin)).match + # Δ and the penetrated field come from the rotated solve whether or not profiles are requested. + @test ray.deltar == center.deltar && ray.bpen == center.bpen + @test isempty(center.inner_psi) && isempty(center.inner_xi) && isempty(center.inner_b) + @test !isempty(ray.inner_xi) && !isempty(galm.inner_xi) + # At ψ_s the composite ξ is the even (Δ₊) column and b the odd (Δ₋) one, so a parity swap in the + # profile assembly makes the two backends disagree here. + for (i, sing) in enumerate(gal.surfaces) + ir = findfirst(==(sing.psifac), ray.inner_psi[i]) + ig = findfirst(==(sing.psifac), galm.inner_psi[i]) + @test ir !== nothing && ig !== nothing + @test isapprox(ray.inner_xi[i][ir, :], galm.inner_xi[i][ig, :]; rtol=1e-4) + @test isapprox(ray.inner_b[i][ir, :], galm.inner_b[i][ig, :]; rtol=1e-4) + end + end + end + @testset "a resistive match keeps its coefficients and drops the ideal δW" begin mktempdir() do dir ric = solve(equil, Riccati(); nn=1, dir_path=dir, ffs_kwargs...)