diff --git a/benchmarks/benchmark_diiid_ideal_ntv_torque.jl b/benchmarks/benchmark_diiid_ideal_ntv_torque.jl index 3b2ae96fe4..c3724c3035 100644 --- a/benchmarks/benchmark_diiid_ideal_ntv_torque.jl +++ b/benchmarks/benchmark_diiid_ideal_ntv_torque.jl @@ -145,7 +145,6 @@ singfac_min = 1e-4 [Riccati] eulerlagrange_tolerance = 1e-7 save_interval = 3 -renorm_threshold = 1e4 """ ) end diff --git a/benchmarks/benchmark_riccati_der.jl b/benchmarks/benchmark_riccati_der.jl deleted file mode 100644 index 9a1d6411b5..0000000000 --- a/benchmarks/benchmark_riccati_der.jl +++ /dev/null @@ -1,117 +0,0 @@ -# Sanity check: riccati_der! correctly evaluates the explicit Riccati ODE. -# -# riccati_der! implements [Glasser 2018 Phys. Plasmas 25, 032507, Eq. 19]: -# dS/dψ = w†·F̄⁻¹·w - S·Ḡ·S, w = Q - K̄·S -# -# where Q = diag(1/(m - n·q)), F̄ = L·L† (Cholesky), K̄ and Ḡ are the MHD -# metric matrices evaluated at ψ. -# -# NOTE: The identity between this Riccati ODE and the EL chain rule -# dS/dψ = dU₁·U₂⁻¹ - S·dU₂·U₂⁻¹ -# holds ONLY for Hermitian S (physical states evolved from the axis, where -# S†=S is preserved by the EL symmetry). For arbitrary non-Hermitian (U₁, U₂), -# the two expressions differ — so this script compares riccati_der! against the -# explicit formula rather than against sing_der!. -# -# Usage (from the repository root): -# julia --project=. benchmarks/benchmark_riccati_der.jl - -using LinearAlgebra, Random, Printf, TOML -using GeneralizedPerturbedEquilibrium - -const FFS = GeneralizedPerturbedEquilibrium.ForceFreeStates - -function setup_solovev() - GPEC = GeneralizedPerturbedEquilibrium - ex = joinpath(@__DIR__, "..", "test", "test_data", "regression_solovev_ideal_example") - inputs = TOML.parsefile(joinpath(ex, "gpec.toml")) - equil = GPEC.Equilibrium.PlasmaEquilibrium(GPEC.equilibrium_config(inputs, ex), GPEC.analytic_config(inputs)) - prob = GPEC.ForceBalanceProblem(equil; GPEC.force_balance_keywords(inputs, ex)...) - ctrl, intr, equil, _, mats = GPEC._prepare_solve(prob, GPEC.force_balance_algorithm(inputs)) - return ctrl, equil, mats, intr -end - -# Evaluate the Riccati RHS explicitly from splines: dS = w†·F̄⁻¹·w - S·Ḡ·S -function riccati_rhs_manual(S, psi, equil, mats, intr) - N = intr.numpert_total - L = zeros(ComplexF64, N, N) - Kmat = zeros(ComplexF64, N, N) - Gmat = zeros(ComplexF64, N, N) - mats.ideal.F_spline_lower(vec(L), psi; hint=mats._hint) - mats.ideal.K_spline(vec(Kmat), psi; hint=mats._hint) - mats.ideal.G_spline(vec(Gmat), psi; hint=mats._hint) - - q = equil.profiles.q_spline(psi) - singfac = vec(1.0 ./ ((intr.mlow:intr.mhigh) .- q .* (intr.nlow:intr.nhigh)')) - - # w = Q - K̄·S (Q is diagonal; add only the diagonal entries) - w = -Kmat * S - for i in 1:N - w[i, i] += singfac[i] - end - - # v = F̄⁻¹·w via stored Cholesky factor L (L·L† = F̄) - v = copy(w) - ldiv!(LowerTriangular(L), v) - ldiv!(UpperTriangular(L'), v) - - return adjoint(w) * v - S * Gmat * S -end - -println("\n=== riccati_der! formula verification ===") -println("Verifies riccati_der! output matches manual evaluation of Glasser 2018 Eq. 19.") -println("Test state: Hermitian S (physical constraint). Expected error: ~machine epsilon.\n") - -ctrl, equil, mats, intr = setup_solovev() -N = intr.numpert_total - -odet = FFS.OdeState(N, ctrl.algorithm.numsteps_init, ctrl.algorithm.numunorms_init, intr.msing) -FFS.initialize_el_at_axis!(odet, ctrl, mats, equil.profiles, intr) -chunks = FFS.chunk_el_integration_bounds(odet, ctrl, intr) - -# 30% into each chunk: well inside the interval, away from singularities at psi_end -test_psis = [c.psi_start + 0.3 * (c.psi_end - c.psi_start) for c in chunks] - -println(" N=$N modes, $(length(test_psis)) test ψ points (30% into each chunk)\n") -@printf(" %8s %14s %14s %12s\n", "ψ", "‖dS_manual‖", "‖dS_ric‖", "rel error") -println(" " * "-"^54) - -rng = Random.MersenneTwister(42) -threshold = 1e-10 - -max_err = let max_err = 0.0 - for psi in test_psis - # Hermitian S: physical Riccati matrix is Hermitian (preserved by EL symmetry) - A = randn(rng, ComplexF64, N, N) - S = (A + A') / 2 # Hermitian by construction - - # Manual RHS - dS_manual = riccati_rhs_manual(S, psi, equil, mats, intr) - - # riccati_der! RHS - u_ric = zeros(ComplexF64, N, N, 2) - du_ric = zeros(ComplexF64, N, N, 2) - u_ric[:, :, 1] .= S - u_ric[:, :, 2] .= Matrix{ComplexF64}(I, N, N) - dummy_chunk = FFS.IntegrationChunk(psi, psi, false, 0, 1) - params = (ctrl, equil, mats, intr, odet, dummy_chunk) - FFS.riccati_der!(du_ric, u_ric, params, psi) - dS_ric = du_ric[:, :, 1] - - ref = max(norm(dS_manual), 1e-10) - err = norm(dS_ric - dS_manual) / ref - max_err = max(max_err, err) - status = err < threshold ? "" : " ← FAIL" - @printf(" %8.4f %14.4e %14.4e %12.4e%s\n", psi, norm(dS_manual), norm(dS_ric), err, status) - end - max_err -end - -println() -if max_err < threshold - @printf("PASSED — max rel error = %.2e (threshold %.0e)\n", max_err, threshold) -else - @printf("FAILED — max rel error = %.2e exceeds threshold %.0e\n", max_err, threshold) - exit(1) -end -println() diff --git a/docs/development/references.md b/docs/development/references.md index 9dbb3c4112..e1ca1cda8e 100644 --- a/docs/development/references.md +++ b/docs/development/references.md @@ -112,7 +112,7 @@ GPEC will eventually implement resistive MHD stability analysis based on: - **Glasser (2018)**: "A robust solution for the resistive MHD toroidal Δ′ matrix in near real-time" - Location: `docs/resources/2018-Glasser-A robust solution for the resistive MHD toroidal Delta-prime matrix in near real-time.pdf` - - Published: Physics of Plasmas **25**, 032501 (2018) + - Published: Physics of Plasmas **25**, 082502 (2018) - Describes: Fast computation of Δ' matrix for resistive stability - **Wang et al. (2020)**: "Modeling of resistive plasma response in toroidal geometry using an asymptotic matching approach" diff --git a/docs/development/regression-harness.md b/docs/development/regression-harness.md index c6da88b19a..22a715e67f 100644 --- a/docs/development/regression-harness.md +++ b/docs/development/regression-harness.md @@ -115,7 +115,7 @@ regress --cases solovev_n1 --ref-range develop~10..develop - `--jobs N` — run N cases at a time (see *Parallel lanes* above) - `--isolate` — run each case in its own fresh subprocess instead of a warm worker -GPEC workers run with `-t auto` (all cores) so GPEC's threaded kernels are active; set `GPEC_REGRESS_THREADS=1` to force single-threaded runs. Tracked quantities are not bit-identical across thread counts. Most move only in the last few bits, inside their noise floors; the parallel Riccati integrator splits its domain by thread count, so its energies move by ~5e-9 and its step counts by ~2 between 1 and 16 threads, and the Riccati cases' floors are set above that spread. The count each run actually used is recorded in its environment fingerprint (shown in the report's `env:` lines), and the report flags two compared runs with different counts. Thread count is not part of the cache key; a `CHANGED` row in a comparison across thread counts that sits just above its floor is worth re-checking at a matched `GPEC_REGRESS_THREADS`. +GPEC workers run with `-t auto` (all cores) so GPEC's threaded kernels are active; set `GPEC_REGRESS_THREADS=1` to force single-threaded runs. Tracked quantities are not bit-identical across thread counts. Most move only in the last few bits, inside their noise floors; the Riccati integrator's chunking does not depend on the thread count. The count each run actually used is recorded in its environment fingerprint (shown in the report's `env:` lines), and the report flags two compared runs with different counts. Thread count is not part of the cache key; a `CHANGED` row in a comparison across thread counts that sits just above its floor is worth re-checking at a matched `GPEC_REGRESS_THREADS`. On a machine you share with other people or with your own parallel sessions, `-t auto` is antisocial: set `GPEC_REGRESS_THREADS` to a bounded count, and `JULIA_NUM_PRECOMPILE_TASKS` alongside it, since `Pkg.instantiate()` otherwise fans out to `Sys.CPU_THREADS + 1` precompile workers at the start of every ref. If the machine has a batch scheduler, submit the run through it and pin both variables to the allocation, with `--jobs` × `GPEC_REGRESS_THREADS` equal to the CPUs requested. diff --git a/docs/src/stability.md b/docs/src/stability.md index aa86df2ecb..be8eaebff7 100644 --- a/docs/src/stability.md +++ b/docs/src/stability.md @@ -34,14 +34,14 @@ where ``\mathbf{L}`` is built from the Fourier-space metric matrices ``\bar{F}`` The resonance factor vanishes at each rational surface, so the solutions there behave as ``|\psi - \psi_s|^{-1/2 \pm \sqrt{-D_I}}``, a large and a small solution governed by the Mercier index ``D_I``; GPEC assembles the right-hand side from reduced matrices that carry explicit ``(m - nq)`` factors rather than their inverses, the convention used on every page. The Newcomb criterion states that the plasma is fixed-boundary stable if and only if the solution that is regular at the axis remains nonsingular across every rational surface and up to the edge. -The Riccati integrator works with the matrix ``S = U_1 U_2^{-1}`` ([Glasser, Kolemen and Glasser (2018)](citations.md#ForceFreeStates), Eq. 19), which obeys a quadratic matrix ODE and stays bounded where ``U_1`` and ``U_2`` grow exponentially, so that the rational surfaces can be crossed without the Gaussian reduction of the forward integration. +The Riccati integrator works with the matrix ``S = U_1 U_2^{-1}`` ([Glasser, Kolemen and Glasser (2018)](citations.md#ForceFreeStates), Eq. 36), which obeys a quadratic matrix ODE and stays bounded where ``U_1`` and ``U_2`` grow exponentially, so that the rational surfaces can be crossed without the Gaussian reduction of the forward integration. Its global boundary-value problem over all surfaces gives the Δ′ matrix of [Glasser and Kolemen (2018)](citations.md#Tearing), stored in the PEST-3 convention: the diagonal entry of each surface is its outer-region Δ′, positive when the outer region drives the tearing mode, and the off-diagonal entries couple the surfaces. ## Integrators `[ForceFreeStates] integrator` selects the radial method; its numerical settings sit in a section of the same name. -- `"Riccati"` (the default; `[Riccati]`) splits the domain into independent chunks integrated in parallel, crosses each rational surface in the Riccati basis and assembles the Δ′ matrix. It stores the chunk endpoints only, so it yields no dense ``\xi`` profiles. +- `"Riccati"` (the default; `[Riccati]`) splits the domain into independent chunks integrated in parallel, crosses each rational surface in its asymptotic basis and assembles the Δ′ matrix. It stores the solution at every solver step, like Forward, so its result carries ``\xi`` profiles; `plasma_response` does not accept a Riccati result yet. - `"Forward"` (`[Forward]`) integrates the system directly in ``(U_1, U_2)`` with Gaussian reduction at each rational surface. It returns dense displacement profiles, and is the only integrator whose result `plasma_response` accepts (the response operators are built from the dense ``\xi`` solution; a wall is refused too) and the only one that supports the kinetic solve. - `"Galerkin"` (`[Galerkin]`) solves the outer region variationally on a finite-element basis, giving the resistive Δ′ matrix and the PEST-3 blocks but no free-boundary energies; see [Galerkin Solver](galerkin.md). @@ -69,7 +69,7 @@ save_interval = 1 # save every Nth ODE step; steps near rational surfaces renorm_threshold = 1e3 # renormalize the solution basis above this column norm ``` -`examples/DIIID-like_riccati_deltaprime_example` sets the same three keys under `[Riccati]`. +`examples/DIIID-like_riccati_deltaprime_example` sets `eulerlagrange_tolerance` and `save_interval` under `[Riccati]`; `renorm_threshold` belongs to `[Forward]` only. The same solve from a script takes the deck keys as keywords of the problem and the integrator: @@ -136,10 +136,10 @@ The perpendicular diffusivities that set the layer widths are described with the ## Numerical notes -- The Riccati integrator splits the domain into chunks whose count is derived from the number of singular surfaces (`[Riccati] nchunks` overrides the target, `0` is automatic); every chunk integrates from identity initial conditions, so the result is independent of the thread count. -- The chunk ending at each rational surface is integrated backward, from the surface inward, because the forward propagator is ill-conditioned there; the Δ′ boundary-value problem is assembled from the split forward and backward propagators for the same reason. -- Both integrators use the adaptive ninth-order `Vern9` solver at `eulerlagrange_tolerance`; the Riccati basis is recovered by renormalizing whenever the solution exceeds `renorm_threshold`. -- The two crossing conventions differ slightly, so Riccati and Forward free-boundary energies agree to within a few percent, not to machine precision. +- The Riccati integrator splits the domain into chunks whose count is derived from the number of singular surfaces (`[Riccati] nchunks` overrides the target, `0` is automatic); every chunk integrates from its own start (the identity, or the asymptotic solutions of [Glasser and Kolemen (2018)](citations.md#Tearing) for the two chunks next to a rational surface), so the result is independent of the thread count. +- The chunk ending at each rational surface is integrated backward, from the surface inward, because the forward propagator is ill-conditioned there; the Δ′ boundary-value problem is assembled from products of the same propagators, matched midway between the surfaces. +- Both integrators use the adaptive ninth-order `Vern9` solver at `eulerlagrange_tolerance`; the Riccati integrator reduces its state to ``(S, I)`` after every chunk, and Forward applies Gaussian reduction whenever a solution norm grows past `renorm_threshold`. +- The two integrators cross a rational surface differently. Riccati carries the nonresonant asymptotic coefficients across it (Glasser and Kolemen (2018), Eq. 36), so its energies do not depend on `singfac_min`. Forward takes a trapezoidal step between the two sides, with an error of first order in `singfac_min`, so its energies converge to Riccati's as `singfac_min` shrinks. Local stability (Mercier, resistive interchange and ballooning) is a separate stage on the same equilibrium; see [Local Stability](ballooning.md). diff --git a/examples/DIIID-like_SLAYER_example/gpec.toml b/examples/DIIID-like_SLAYER_example/gpec.toml index fa41a0dd50..0f5f28d085 100644 --- a/examples/DIIID-like_SLAYER_example/gpec.toml +++ b/examples/DIIID-like_SLAYER_example/gpec.toml @@ -56,8 +56,7 @@ rational_offset_fraction = 0.2 # Truncate integration at (last_rational_q + rati # Numerical settings of the Riccati integrator [Riccati] eulerlagrange_tolerance = 1e-10# Relative tolerance for ODE integration of Euler-Lagrange equations -save_interval = 1 # Save every Nth ODE step (1=all). Always saves near rational surfaces. -renorm_threshold = 1e3 # Column-norm threshold that triggers solution renormalization +save_interval = 1 # Store every Nth solver step of each chunk (1=all); chunk ends and the edge-dW scan band are always stored # SLAYER tearing-mode growth-rate analysis on the DIII-D-like H-mode equilibrium. # Runs off the ForceFreeStates result (no PE). The plasma (profiles incl. the χ⊥/χ_φ diff --git a/examples/DIIID-like_riccati_deltaprime_example/gpec.toml b/examples/DIIID-like_riccati_deltaprime_example/gpec.toml index 31457af0d4..ddc2ee098a 100644 --- a/examples/DIIID-like_riccati_deltaprime_example/gpec.toml +++ b/examples/DIIID-like_riccati_deltaprime_example/gpec.toml @@ -55,5 +55,4 @@ rational_offset_fraction = 0.2 # Truncate integration at (last_rational_q + rati # Numerical settings of the Riccati integrator [Riccati] eulerlagrange_tolerance = 1e-10 # Relative tolerance for ODE integration of Euler-Lagrange equations -save_interval = 1 # Save every Nth ODE step (1=all). Always saves near rational surfaces. -renorm_threshold = 1e3 # Column-norm threshold that triggers solution renormalization +save_interval = 1 # Store every Nth solver step of each chunk (1=all); chunk ends and the edge-dW scan band are always stored diff --git a/examples/LAR_beta_scan/gpec.toml b/examples/LAR_beta_scan/gpec.toml index d40e82a476..5906bf4c4c 100644 --- a/examples/LAR_beta_scan/gpec.toml +++ b/examples/LAR_beta_scan/gpec.toml @@ -69,5 +69,4 @@ rational_offset_fraction = 0.2 # Truncate integration at (last_rational_q + rati # Numerical settings of the Riccati integrator [Riccati] eulerlagrange_tolerance = 1e-12 # Relative tolerance for ODE integration of Euler-Lagrange equations -renorm_threshold = 1e3 # Column-norm threshold that triggers solution renormalization -save_interval = 1 # Save every Nth ODE step (1=all). Always saves near rational surfaces. +save_interval = 1 # Store every Nth solver step of each chunk (1=all); chunk ends and the edge-dW scan band are always stored diff --git a/examples/LAR_epsilon_scan/gpec.toml b/examples/LAR_epsilon_scan/gpec.toml index 119bb894a7..a5a19e9ac9 100644 --- a/examples/LAR_epsilon_scan/gpec.toml +++ b/examples/LAR_epsilon_scan/gpec.toml @@ -70,5 +70,4 @@ rational_offset_fraction = 0.2 # Truncate integration at (last_rational_q + rati # Numerical settings of the Riccati integrator [Riccati] eulerlagrange_tolerance = 1e-12 # Relative tolerance for ODE integration of Euler-Lagrange equations -renorm_threshold = 1e3 # Column-norm threshold that triggers solution renormalization -save_interval = 1 # Save every Nth ODE step (1=all). Always saves near rational surfaces. +save_interval = 1 # Store every Nth solver step of each chunk (1=all); chunk ends and the edge-dW scan band are always stored diff --git a/examples/Solovev_ideal_example_3D/gpec.toml b/examples/Solovev_ideal_example_3D/gpec.toml index a7b68e8d62..e5a6681c68 100644 --- a/examples/Solovev_ideal_example_3D/gpec.toml +++ b/examples/Solovev_ideal_example_3D/gpec.toml @@ -51,8 +51,7 @@ rational_offset_fraction = 0.2 # Truncate integration at (last_rational_q + rati # Numerical settings of the Riccati integrator [Riccati] eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations -renorm_threshold = 1e3 # Column-norm threshold that triggers solution renormalization -save_interval = 1 # Save every Nth ODE step (1=all). Always saves near rational surfaces. +save_interval = 1 # Store every Nth solver step of each chunk (1=all); chunk ends and the edge-dW scan band are always stored [Wall] # See examples/Solovev_ideal_example/gpec.toml for the rationale behind the close conformal wall. diff --git a/examples/Solovev_ideal_example_multi_n/gpec.toml b/examples/Solovev_ideal_example_multi_n/gpec.toml index c8e51b8c78..569280fc47 100644 --- a/examples/Solovev_ideal_example_multi_n/gpec.toml +++ b/examples/Solovev_ideal_example_multi_n/gpec.toml @@ -58,7 +58,6 @@ rational_offset_fraction = 0.2 # Truncate integration at (last_rational_q + rati # Numerical settings of the Riccati integrator [Riccati] eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations -renorm_threshold = 1e3 # Column-norm threshold that triggers solution renormalization # Solovev analytic equilibrium parameters; this section selects the equilibrium. See SolovevConfig in src/Equilibrium. [Solovev] diff --git a/examples/Solovev_ideal_example_multi_n/single_n_1/gpec.toml b/examples/Solovev_ideal_example_multi_n/single_n_1/gpec.toml index 0c272a9a0a..8e611f4956 100644 --- a/examples/Solovev_ideal_example_multi_n/single_n_1/gpec.toml +++ b/examples/Solovev_ideal_example_multi_n/single_n_1/gpec.toml @@ -52,7 +52,6 @@ singfac_min = 1e-4 # Fractional distance from rational q at which ide # Numerical settings of the Riccati integrator [Riccati] eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations -renorm_threshold = 1e3 # Column-norm threshold that triggers solution renormalization # Solovev analytic equilibrium parameters; this section selects the equilibrium. See SolovevConfig in src/Equilibrium. [Solovev] diff --git a/examples/Solovev_ideal_example_multi_n/single_n_2/gpec.toml b/examples/Solovev_ideal_example_multi_n/single_n_2/gpec.toml index f200df8c2a..0aa031281d 100644 --- a/examples/Solovev_ideal_example_multi_n/single_n_2/gpec.toml +++ b/examples/Solovev_ideal_example_multi_n/single_n_2/gpec.toml @@ -52,7 +52,6 @@ singfac_min = 1e-4 # Fractional distance from rational q at which ide # Numerical settings of the Riccati integrator [Riccati] eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations -renorm_threshold = 1e3 # Column-norm threshold that triggers solution renormalization # Solovev analytic equilibrium parameters; this section selects the equilibrium. See SolovevConfig in src/Equilibrium. [Solovev] diff --git a/regression-harness/cases/diiid_n1_riccati.toml b/regression-harness/cases/diiid_n1_riccati.toml index 4c2642145e..3dbbb2d9b4 100644 --- a/regression-harness/cases/diiid_n1_riccati.toml +++ b/regression-harness/cases/diiid_n1_riccati.toml @@ -36,7 +36,6 @@ label = "edge coil response delta_coil" noise_threshold = 1e-8 order = 12 -# The parallel Riccati chunking depends on the thread count, which moves energies by ~5e-9 (spectrum ~1e-6 absolute) and step counts by ~2 between 1 and 16 threads; the floors below sit well above that. # Energies — leading generalized (W,N) pencil eigenvalues at the final truncation (psilim): power-normalized (⟨|ξ|²⟩ = 1 metric) and invariant to the working-coordinate Jacobian. [quantities.et_real] h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_energies" @@ -95,7 +94,7 @@ label = "singular q values" noise_threshold = 1e-8 order = 32 -# Integration — chunk-endpoint step counts move when the chunk sizing changes. +# Integration — stored and total solver-step counts move when the chunk layout or the tolerance changes. [quantities.nstep] h5path = "ForceFreeStates/Solutions/ForwardIntegration/nstep" type = "int_scalar" diff --git a/regression-harness/cases/diiid_n1_riccati_q1.toml b/regression-harness/cases/diiid_n1_riccati_q1.toml new file mode 100644 index 0000000000..bf43935d32 --- /dev/null +++ b/regression-harness/cases/diiid_n1_riccati_q1.toml @@ -0,0 +1,12 @@ +# Regression case: the diiid_n1_riccati deck with the q = 1 surface inside the integration domain. +# The on-axis q is lowered to 0.9 and the integration starts inside q = 1, so the m = 1 surface is +# crossed like the others. It tracks the energies and Δ' where the singular-surface crossing +# matters most: the base deck starts outside q = 1.02 and never crosses that surface. +[case] +name = "diiid_n1_riccati_q1" +base = "diiid_n1_riccati" # Inherit example_dir, [quantities.*] and [overrides] from this case +description = "DIII-D-like equilibrium, n=1, Riccati integrator with the q = 1 surface inside the domain" + +[overrides] +"Equilibrium.newq0" = 0.9 # On-axis q below 1, so a q = 1 surface exists +"ForceFreeStates.qlow" = 0.1 # Start the integration inside q = 1 diff --git a/regression-harness/cases/solovev_multi_n.toml b/regression-harness/cases/solovev_multi_n.toml index 7f366331d1..caadd39091 100644 --- a/regression-harness/cases/solovev_multi_n.toml +++ b/regression-harness/cases/solovev_multi_n.toml @@ -8,7 +8,6 @@ description = "Solovev analytical equilibrium, multi-n, ideal stability" example_dir = "examples/Solovev_ideal_example_multi_n" # Energies — leading generalized (W,N) pencil eigenvalues at the final truncation (psilim): power-normalized (⟨|ξ|²⟩ = 1 metric) and invariant to the working-coordinate Jacobian. -# The parallel Riccati chunking depends on the thread count, which moves energies by ~2e-9 (spectra ~1e-6 absolute) and step counts by ~2 between 1 and 16 threads; the floors below sit well above that. [quantities.et_real] h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_energies" type = "complex_vector" diff --git a/src/Analysis/Registry.jl b/src/Analysis/Registry.jl index 03940daa32..159531f0d1 100644 --- a/src/Analysis/Registry.jl +++ b/src/Analysis/Registry.jl @@ -44,7 +44,7 @@ What each capability needs, worded as what to change in the run to get it. const REASONS = Dict{Symbol,String}( :equilibrium => "needs an equilibrium", :singular_surfaces => "needs rational surfaces in the integration domain", - :solution_profiles => "needs stored ξ profiles: use the Forward integrator or a matched Galerkin solve (Riccati and unmatched Galerkin store none)", + :solution_profiles => "needs stored ξ profiles: use the Forward or Riccati integrator, or a matched Galerkin solve (unmatched Galerkin stores none)", :fixed_boundary_criterion => "needs the critical determinant scan: use the Forward or Riccati integrator (Galerkin computes none)", :free_boundary => "needs the free-boundary energies: run with vac_flag = true and the Forward or Riccati integrator", :edge_scan => "needs the edge stability scan: set dW_edge_scan_start below psilim", diff --git a/src/Deck.jl b/src/Deck.jl index ea0f4b17ce..c4403b60b1 100644 --- a/src/Deck.jl +++ b/src/Deck.jl @@ -20,7 +20,8 @@ const _RUN_SECTION_KEYS = ("verbose", "HDF5_filename", "output_filename", "write const _RUN_KEY_SECTIONS = ("Equilibrium", "LocalStability", "ForceFreeStates", "PerturbedEquilibrium", "KineticForces", "ErrorFields") # Keys that moved to another section, each with where it lives now. -const _ODE_DECK_KEYS = ("eulerlagrange_tolerance", "renorm_threshold", "numsteps_init", "numunorms_init", "save_interval") +const _ODE_DECK_KEYS = ("eulerlagrange_tolerance", "save_interval") +const _FORWARD_DECK_KEYS = ("renorm_threshold", "numsteps_init", "numunorms_init") const _SEARCH_DECK_KEYS = ("scan_mode", "Q_re_range", "Q_im_range", "nre", "nim", "amr_passes", "amr_max_cells", "boxes", "multi_box_prescreen_n", "pole_threshold", "pole_threshold_adaptive", "filter_above_poles", "filter_outside_re", "gamma_gap_threshold", "polish_roots", "validity_rtol", "root_classification", "store_scan") @@ -35,6 +36,7 @@ const _MOVED_KEYS = Dict{String,Dict{String,String}}( "gal_ideal_flag" => "moved to [Match] as ideal", "gal_inner_solver" => "moved to [GGJ] as solver", (k => "moved to the integrator's section, [Forward] or [Riccati]" for k in _ODE_DECK_KEYS)..., + (k => "moved to [Forward]" for k in _FORWARD_DECK_KEYS)..., ("gal_$k" => "moved to [Match] as $k" for k in ("eta", "rho", "rotation", "gamma"))..., ("gal_inner_$k" => "moved to [GGJ] as $k" for k in ("xfac", "nx", "nq", "cutoff", "kmax"))..., ("gal_$k" => "moved to [Galerkin] as $k" for k in ("nx", "nq", "pfac", "dx0", "dx1", "dx2", "cutoff", "tol", "gnstep", diff --git a/src/Display.jl b/src/Display.jl index 662cbbd476..f2bb509569 100644 --- a/src/Display.jl +++ b/src/Display.jl @@ -180,13 +180,12 @@ end function _solution_note(r::ForceFreeStates.ForceFreeStatesResult) r.solution === nothing || return "ξ profiles stored ($(r.solution.step) ψ nodes, $(r.solution.basis) basis)" - return r.integrator === :riccati ? "not computed (Riccati stores no ξ profiles)" : "not computed ($(r.integrator) without a match)" + return "not computed ($(r.integrator) without a match)" end # The crit sign-change count is not yet a validated stability verdict, so it is reported raw. function _fixed_boundary_note(r::ForceFreeStates.ForceFreeStatesResult) r.diagnostics === nothing && return "not computed by the $(r.integrator) integrator" - r.integrator === :riccati && return "not reported (the Riccati (S, I) history cannot show crossings)" r.control.kinetic === nothing || return "not reported (crit assumes a Hermitian W)" return "$(_plural(r.diagnostics.nzero, "crit sign change")) (unvalidated, not a stability verdict)" end @@ -254,6 +253,7 @@ _resonances(r::PlasmaResponse) = ["$(m)/$(n)" for (m, n) in zip(r.rational_m_res function _missing_operators_note(r::PlasmaResponse) r.ffs.free_boundary === nothing && return "not computed (needs a free-boundary result)" r.ffs.solution === nothing && return "not computed (needs stored ξ profiles)" + r.ffs.integrator === :riccati && return "not computed (a Riccati result is not accepted yet)" return "not computed" end diff --git a/src/ForceFreeStates/CoreTypes.jl b/src/ForceFreeStates/CoreTypes.jl index 1cb794f65c..c1d0e8f17e 100644 --- a/src/ForceFreeStates/CoreTypes.jl +++ b/src/ForceFreeStates/CoreTypes.jl @@ -11,6 +11,22 @@ Supertype for objects that carry the resolved (m, n) mode space — `mlow`, `mhi """ abstract type ModeSpace end +""" + ModeIndices(modes::ModeSpace) + +Immutable copy of a mode space's indexing, for code that must not see the rest of its source. +""" +struct ModeIndices <: ModeSpace + mlow::Int + mhigh::Int + mpert::Int + nlow::Int + nhigh::Int + npert::Int + numpert_total::Int +end +ModeIndices(m::ModeSpace) = ModeIndices(m.mlow, m.mhigh, m.mpert, m.nlow, m.nhigh, m.npert, m.numpert_total) + """ ForceFreeStatesInternal @@ -59,26 +75,6 @@ Solve-time scratch state for a force-balance solve: the resolved mode space, sin q1lim::Float64 = 0.0 wall_settings::Union{Nothing,Vacuum.WallSettings} = nothing verbose::Bool = false - """ - Inter-surface Δ' matrix, `msing × msing`, PEST3 convention: the odd-parity projection of - `delta_prime_raw`, set by `compute_delta_prime_matrix!` (Riccati path only). - """ - delta_prime_matrix::Matrix{ComplexF64} = Matrix{ComplexF64}(undef, 0, 0) - - """ - Edge coil-response matrix, `2msing × numpert_total`: column k is the resonant small-solution - response at each surface side to a unit source on edge mode k (Eq. 37 rpec edge condition, - `_solve_bvp_edge_coil`). Empty unless the S-axis vacuum-edge BVP was assembled. - """ - delta_coil::Matrix{ComplexF64} = Matrix{ComplexF64}(undef, 0, 0) - - """ - Raw `2msing × 2msing` outer-region matching matrix D' (Pletzer–Dewar 1991, before parity - rotation), side-major ordering `[L_s1, R_s1, …, L_sm, R_sm]`. `pest3_decompose` recovers the - A', B', Γ', Δ' blocks. No ½ prefactor (matches Fortran rdcon; the paper applies ½). - Empty unless the Riccati integrator was used. - """ - delta_prime_raw::Matrix{ComplexF64} = Matrix{ComplexF64}(undef, 0, 0) end """ diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl index d407000592..ef90625414 100644 --- a/src/ForceFreeStates/EulerLagrange.jl +++ b/src/ForceFreeStates/EulerLagrange.jl @@ -53,8 +53,6 @@ records needed to undo it, and singular-surface asymptotic coefficients. - `du_store::Array{ComplexF64,3}` - du₁/dψ at each stored step, shape `(numpert_total, numpert_total, step)`; empty until `materialize_derivative_stores!` fills it (the galerkin-matched path supplies it at construction). - `xi_s_store::Array{ComplexF64,3}` - Clebsch displacement Ξ_s (Glasser 2016 Eq. 18) at each stored step; empty until materialized, like `du_store`. - - `u_store_el_basis::Bool` - True when `u_store` holds the Euler-Lagrange state; false on the sparse parallel path (Riccati-basis chunk endpoints), where derivative - stores cannot be materialized. - `du_store_populated::Bool` - True once `du_store`/`xi_s_store` hold valid data in the final basis. - `crit_store::Vector{Float64}` - Critical eigenvalue (smallest eigenvalue of W⁻¹) at each stored step. - `ca_r::Array{ComplexF64,4}` - Asymptotic coefficients right of each singular surface, shape `(numpert_total, numpert_total, 2, msing)`. @@ -94,7 +92,6 @@ records needed to undo it, and singular-surface asymptotic coefficients. u_store::Array{ComplexF64,4} = Array{ComplexF64}(undef, numpert_total, numpert_total, 2, numsteps_init) du_store::Array{ComplexF64,3} = Array{ComplexF64}(undef, numpert_total, numpert_total, 0) xi_s_store::Array{ComplexF64,3} = Array{ComplexF64}(undef, numpert_total, numpert_total, 0) - u_store_el_basis::Bool = true du_store_populated::Bool = false crit_store::Vector{Float64} = Vector{Float64}(undef, numsteps_init) ca_r::Array{ComplexF64,4} = zeros(ComplexF64, numpert_total, numpert_total, 2, msing) @@ -133,9 +130,9 @@ OdeState(numpert_total::Int, numsteps_init::Int, numunorms_init::Int, msing::Int """ compute_delta_prime_from_ca!(odet, intr, equil) -**Stub, not physically valid.** Per-surface Δ' estimate `(ca_r[i,i,2,s] - ca_l[i,i,2,s]) / (4π²·psio)` -from the asymptotic coefficients. It ignores inter-surface coupling and the vacuum, so it does not -agree with the physical `ForceFreeStatesInternal.delta_prime_matrix`. Used only by tests. +**Stub, not physically valid.** Fill each surface's `delta_prime` with the per-surface estimate +`(ca_r[i,i,2,s] - ca_l[i,i,2,s]) / (4π²·psio)` from the asymptotic coefficients. It ignores inter-surface +coupling and the vacuum, so it does not agree with the physical `ForceFreeStatesResult.delta_prime`. """ function compute_delta_prime_from_ca!(odet::OdeState, intr::ForceFreeStatesInternal, equil::Equilibrium.PlasmaEquilibrium) denom = (2π)^2 * equil.psio # = twopi * chi1 in SingularCoupling.jl @@ -182,80 +179,60 @@ end """ balance_integration_chunks(chunks, ctrl, intr) -> Vector{IntegrationChunk} -Split the highest-`ode_itime_cost` chunk at its equal-cost midpoint until the count reaches -`ctrl.algorithm.nchunks` (`0` = auto). The target never depends on the thread count, so Riccati -output is thread-count independent. Only the last sub-chunk of a split keeps its crossing. +Split every crossing chunk once, then the costliest chunk (`ode_itime_cost`) until there are +`ctrl.algorithm.nchunks` (`0` = auto). Splits are at the equal-cost point; the upper part keeps the crossing. """ function balance_integration_chunks(chunks::Vector{IntegrationChunk}, ctrl::ForceFreeStatesControl, intr::ForceFreeStatesInternal) - min_chunks = 2 * intr.msing + 3 - # Enough sub-chunks per segment that no single long-span propagator is too ill-conditioned for the BVP - min_bvp_intervals = 8 * (intr.msing + 1) + intr.msing - if ctrl.algorithm.nchunks > 0 - if ctrl.algorithm.nchunks < min_chunks - @warn "nchunks = $(ctrl.algorithm.nchunks) is below the $min_chunks chunks required by $(intr.msing) singular surfaces; clamping up." - end - target_n = max(ctrl.algorithm.nchunks, min_chunks) - else - target_n = max(min_chunks, min_bvp_intervals) - end - - result = collect(chunks) - - while length(result) < target_n - # Find the highest-cost splittable chunk - best_idx = 0 - best_cost = -Inf - for (i, chunk) in enumerate(result) - width = chunk.psi_end - chunk.psi_start - if width > 1e-8 - c = ode_itime_cost(chunk.psi_start, chunk.psi_end, intr) - if c > best_cost - best_cost = c - best_idx = i - end - end - end - - best_idx == 0 && break # No more splittable chunks - - chunk = result[best_idx] - total_cost = best_cost - target_cost = total_cost / 2.0 - - # Bisect for ψ_mid with cost(psi_start, ψ_mid) ≈ target_cost + function split(chunk::IntegrationChunk) + half_cost = ode_itime_cost(chunk.psi_start, chunk.psi_end, intr) / 2 lo, hi = chunk.psi_start, chunk.psi_end for _ in 1:50 - mid = (lo + hi) / 2.0 - if ode_itime_cost(chunk.psi_start, mid, intr) < target_cost + mid = (lo + hi) / 2 + if ode_itime_cost(chunk.psi_start, mid, intr) < half_cost lo = mid else hi = mid end end - psi_mid = (lo + hi) / 2.0 - - left = IntegrationChunk(; psi_start=chunk.psi_start, psi_end=psi_mid, - needs_crossing=false, ising=0, direction=1) - right = IntegrationChunk(; psi_start=psi_mid, psi_end=chunk.psi_end, - needs_crossing=chunk.needs_crossing, ising=chunk.ising, - direction=chunk.direction) - splice!(result, best_idx, [left, right]) + psi_mid = (lo + hi) / 2 + return [IntegrationChunk(; psi_start=chunk.psi_start, psi_end=psi_mid, needs_crossing=false), + IntegrationChunk(; psi_start=psi_mid, psi_end=chunk.psi_end, chunk.needs_crossing, chunk.ising, chunk.direction)] + end + + result = IntegrationChunk[] + for chunk in chunks + append!(result, chunk.needs_crossing ? split(chunk) : [chunk]) + end + + target_n = length(result) + if ctrl.algorithm.nchunks > 0 + if ctrl.algorithm.nchunks < target_n + @warn "nchunks = $(ctrl.algorithm.nchunks) is below the $target_n chunks required by the singular surfaces crossed; clamping up." + end + target_n = max(ctrl.algorithm.nchunks, target_n) + else + # Enough sub-chunks per segment that no single long-span propagator is too ill-conditioned for the BVP + target_n = max(target_n, 8 * (intr.msing + 1) + intr.msing) + end + + while length(result) < target_n + costs = [c.psi_end - c.psi_start > 1e-8 ? ode_itime_cost(c.psi_start, c.psi_end, intr) : -Inf for c in result] + cost, costliest = findmax(costs) + isfinite(cost) || break # nothing left wide enough to split + splice!(result, costliest, split(result[costliest])) end return result end """ - eulerlagrange_integration(ctrl, equil, mats, intr) -> (odet, propagators, chunks, S_left) + eulerlagrange_integration(ctrl, equil, mats, intr) -> (odet, riccati) Integrate the Euler-Lagrange equations from the axis to `intr.psilim`, crossing each singular surface on the way. Dispatches on `ctrl.algorithm` to [`riccati_eulerlagrange_integration`](@ref) (the chunked propagator BVP) or -[`forward_eulerlagrange_integration`](@ref). - -Only the Riccati branch populates `propagators` / `chunks` / `S_left`, which -`compute_delta_prime_matrix!` consumes for the Δ' BVP; the forward branch returns `nothing` -for all three. +[`forward_eulerlagrange_integration`](@ref). `riccati` holds the inputs of +`compute_delta_prime_matrix`; it is `nothing` on the forward branch. """ function eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, intr::ForceFreeStatesInternal) @@ -272,13 +249,12 @@ function eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil::Equilibr end """ - forward_eulerlagrange_integration(ctrl, equil, mats, intr; verbose=intr.verbose) -> (odet, nothing, nothing, nothing) + forward_eulerlagrange_integration(ctrl, equil, mats, intr; verbose=intr.verbose) -> (odet, nothing) Forward branch of [`eulerlagrange_integration`](@ref): integrates chunk by chunk from the axis, applying Gaussian reduction whenever a solution norm ratio exceeds `ctrl.algorithm.renorm_threshold` and undoing it -via `transform_u!` at the end, so `odet.u_store` comes back dense in the axis basis. Call -directly to force this branch regardless of `ctrl.algorithm`; `verbose` overrides -`intr.verbose` for progress logging. +via `transform_u!` at the end, so `odet.u_store` comes back dense in the axis basis. `ctrl.algorithm` +must be a `Forward`, whose settings it reads; `verbose` overrides `intr.verbose` for progress logging. """ function forward_eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, intr::ForceFreeStatesInternal; verbose::Bool=intr.verbose) @@ -322,7 +298,7 @@ function forward_eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil:: # Undo Gaussian reduction to get true solution vectors transform_u!(odet, intr) - return (odet, nothing, nothing, nothing) + return odet, nothing end """ @@ -630,7 +606,7 @@ function cross_ideal_singular_surf!( odet.ca_r[:, :, :, ising] .= sing_get_ca(odet.u, ua, intr) odet.ca_populated = true - # No Δ' here: it needs the Riccati gauge (U₂ = I); see compute_delta_prime_matrix! + # No Δ' here: it needs the Riccati propagators; see compute_delta_prime_matrix odet.q = equil.profiles.q_spline(odet.psifac; hint=odet.spline_hint) store_ode_data!(odet, odet.psifac, odet.u) end @@ -697,14 +673,16 @@ end const COLUMN_ABSTOL_CALLBACK = DiscreteCallback((u, t, integrator) -> true, refresh_column_abstol!; save_positions=(false, false)) """ - solve_el(prob, rtol; callback=COLUMN_ABSTOL_CALLBACK) + solve_el(prob, rtol; callback=COLUMN_ABSTOL_CALLBACK, save_everystep=false) -Vern9 solve of an Euler-Lagrange problem returning the end state, with the per-column abstol of `column_abstol!`. +Vern9 solve of an Euler-Lagrange problem with the per-column abstol of `column_abstol!`, returning the end state +or, with `save_everystep`, the state at every solver step. A custom `callback` must refresh that abstol itself, after any reduction or renormalization of the state. Errors when the solver stops short of the interval's end (e.g. the step size underflows), naming the ψ reached. """ -function solve_el(prob, rtol; callback=COLUMN_ABSTOL_CALLBACK) - sol = solve(prob, Vern9(); reltol=rtol, abstol=column_abstol(prob.u0, rtol), callback, save_everystep=false, save_end=true) +function solve_el(prob, rtol; callback=COLUMN_ABSTOL_CALLBACK, save_everystep::Bool=false) + # dense=false: saved steps keep the state only, not the interpolation stages. + sol = solve(prob, Vern9(); reltol=rtol, abstol=column_abstol(prob.u0, rtol), callback, save_everystep, dense=false, save_end=true) SciMLBase.successful_retcode(sol) || error( "Euler-Lagrange integration aborted at ψ = $(@sprintf("%.6f", sol.t[end])) before ψ = $(@sprintf("%.6f", prob.tspan[2])) " * @@ -1011,7 +989,7 @@ end In-place Euler-Lagrange derivative du/dψ at `psieval` [Glasser Phys. Plasmas 2016 112506 Eq. 24], ideal or kinetic, in the `(du, u, p, t)` form the ODE solver calls; `du` and `u` have shape -`(numpert_total, numpert_total, 2)`. Ξ_s is not computed here; see [`compute_node_xi_s!`](@ref). +`(numpert_total, ncolumns, 2)`. Ξ_s is not computed here; see [`compute_node_xi_s!`](@ref). """ function sing_der!(du::Array{ComplexF64,3}, u::Array{ComplexF64,3}, params::Tuple{ForceFreeStatesControl,Equilibrium.PlasmaEquilibrium, @@ -1039,7 +1017,7 @@ end el_derivatives!(du, u, kinetic, equil, mats, intr, psieval, spline_hint, mats_hint) -> q Euler-Lagrange (or, when `kinetic` is true, FKG) derivative kernel: writes du₁/dψ and du₂/dψ at -`psieval` into `du` and returns q there. Holds no state of its own — the two hints are the +`psieval` into `du` and returns q there, for any number of solution columns. Holds no state of its own — the two hints are the caller's interval-search accelerators, so concurrent callers just pass their own. """ @with_pool pool function el_derivatives!(du::Array{ComplexF64,3}, u::Array{ComplexF64,3}, @@ -1054,7 +1032,7 @@ caller's interval-search accelerators, so concurrent callers just pass their own fmat_lower = acquire!(pool, ComplexF64, Npert, Npert) kmat = similar!(pool, fmat_lower) gmat = similar!(pool, fmat_lower) - tmp_mat = similar!(pool, fmat_lower) + tmp_mat = acquire!(pool, ComplexF64, Npert, size(u, 2)) fill!(tmp_mat, zero(ComplexF64)) u1 = @view(u[:, :, 1]) diff --git a/src/ForceFreeStates/Integrators.jl b/src/ForceFreeStates/Integrators.jl index 8210f4f8cc..3705a5e76c 100644 --- a/src/ForceFreeStates/Integrators.jl +++ b/src/ForceFreeStates/Integrators.jl @@ -14,7 +14,7 @@ Serial forward integration with Gaussian reduction and dense solution profiles. ## Fields - `eulerlagrange_tolerance::Float64` - Relative tolerance for ODE integration. - - `renorm_threshold::Float64` - Column-norm growth threshold for Gaussian reduction or Riccati renormalization. + - `renorm_threshold::Float64` - Column-norm growth threshold for Gaussian reduction. - `numsteps_init::Int` - Initial capacity for stored ODE solution steps. - `numunorms_init::Int` - Initial capacity for solution-normalization records. - `save_interval::Int` - Save every Nth ODE step; always retain steps near rational surfaces. @@ -35,21 +35,15 @@ Chunked Riccati integration and the singular-surface boundary-value solve. ## Fields - `eulerlagrange_tolerance::Float64` - Relative tolerance for ODE integration. - - `renorm_threshold::Float64` - Column-norm growth threshold for Gaussian reduction or Riccati renormalization. - - `numsteps_init::Int` - Initial capacity for stored ODE solution steps. - - `numunorms_init::Int` - Initial capacity for solution-normalization records. - - `save_interval::Int` - Save every Nth ODE step; always retain steps near rational surfaces. + - `save_interval::Int` - Store every Nth solver step of each chunk; chunk ends and the edge-dW scan band are always stored. - `nchunks::Int` - Target chunk count; zero derives the count from the singular surfaces. - - `extended_precision_bvp::Bool` - Use Complex{Double64} for the Delta-prime BVP and PEST3 combination to reduce cancellation error. + - `extended_precision_bvp::Bool` - Solve the Delta-prime BVP in Complex{Double64} instead of ComplexF64; far slower, and needed only beyond about 150 poloidal harmonics. """ @kwdef struct Riccati <: AbstractIntegrator eulerlagrange_tolerance::Float64 = 1e-8 - renorm_threshold::Float64 = 1e3 - numsteps_init::Int = 4000 - numunorms_init::Int = 100 save_interval::Int = 1 nchunks::Int = 0 - extended_precision_bvp::Bool = true + extended_precision_bvp::Bool = false end """ diff --git a/src/ForceFreeStates/Matching/MatchProblem.jl b/src/ForceFreeStates/Matching/MatchProblem.jl index f80f52eb46..522409c46d 100644 --- a/src/ForceFreeStates/Matching/MatchProblem.jl +++ b/src/ForceFreeStates/Matching/MatchProblem.jl @@ -276,8 +276,9 @@ end # Copy of ffs with the match applied; the ideal reference keeps closure :ideal, its zero bpen and δW. function _matched_result(ffs::ForceFreeStatesResult, match::MatchResult, ideal::Bool) gal = ffs.galerkin + # Without a matched ξ, the ideal solution stands only for the ideal reference. solution = (gal !== nothing && gal.solution !== nothing && !isempty(match.xi)) ? - _matched_gal_profiles(gal, match, ffs.mats, ffs) : ffs.solution + _matched_gal_profiles(gal, match, ffs.mats, ffs) : (ideal ? ffs.solution : nothing) ideal && return _with(ffs; match=match, solution=solution) # The ideal δW does not describe the matched plasma; PE gates on its absence. return _with(ffs; closure=:matched, bpen=match.bpen, match=match, solution=solution, diff --git a/src/ForceFreeStates/Result.jl b/src/ForceFreeStates/Result.jl index 2d2affa457..2251f1ee7a 100644 --- a/src/ForceFreeStates/Result.jl +++ b/src/ForceFreeStates/Result.jl @@ -63,7 +63,7 @@ Every final result records how its basis crosses the rational surfaces, so bpen inner layer, per surface (rows follow `surfaces`) and driving mode; all zeros under `:ideal` and `:regular` closure. - `match::Union{Nothing,MatchResult}` - The [`MatchResult`](@ref) of a `MatchProblem` solve; `nothing` otherwise. - `solution::Union{Nothing,SolutionProfiles}` - THE solve's ξ solution; `nothing` when the - formalism produced none (Riccati, and Galerkin without a match). + formalism produced none (Galerkin without a match). - `diagnostics::Union{Nothing,OdeState}` - The integrator's raw ODE state (ψ trace, `crit`, edge scan, asymptotic coefficients). Written to HDF5; not a consumable solution. - `wp::Union{Nothing,Matrix{ComplexF64}}` - Fixed-boundary plasma energy matrix `W_p = (U₂·U₁⁻¹)/ψ₀²` at @@ -138,7 +138,7 @@ with a message naming what a solution takes to produce. """ function require_solution(result::ForceFreeStatesResult, calc::AbstractString) result.solution === nothing || return true - @warn "Skipping $calc: no ξ solution — dense profiles require a Forward (or matched Galerkin) run; " * + @warn "Skipping $calc: no ξ solution — dense profiles require a Forward, Riccati or matched Galerkin run; " * "this result came from the $(result.integrator) integrator" return false end @@ -175,15 +175,15 @@ function _matched_gal_profiles(gal_result::GalerkinResult, m::MatchResult, mats: end """ - build_result(integrator, ctrl, equil, intr, metric, mats, odet, free_energies, gal_data, gal_dp) + build_result(integrator, ctrl, equil, intr, metric, mats, odet, free_energies, gal_data, delta_prime) -> ForceFreeStatesResult Assemble the published result once the solve is finished, with the ideal or regular closure (a `MatchProblem` solve transforms it afterwards). Materializes the forward path's derivative stores and forms `W_p` when the free-boundary stage did not run; everything else is copied or aliased. -`odet` is the ODE state (`nothing` for Galerkin); `gal_data`/`gal_dp` are the Galerkin internals and -its Δ′ payload (`nothing` otherwise). +`odet` is the ODE state (`nothing` for Galerkin), `gal_data` the Galerkin internals (`nothing` otherwise), +and `delta_prime` the Δ′ payload of a Riccati or Galerkin solve (`nothing` otherwise). """ function build_result( integrator::Symbol, @@ -195,11 +195,11 @@ function build_result( odet::Union{Nothing,OdeState}, free_energies::Union{Nothing,FreeBoundaryResult}, gal_data::Union{Nothing,GalerkinResult}, - gal_dp::Union{Nothing,DeltaPrimeData} + delta_prime::Union{Nothing,DeltaPrimeData} ) - # The forward sweep is the only formalism whose stores need materializing; doing it here - # keeps `SolutionProfiles.du_store`/`xi_s_store` populated by construction. - solution = if integrator === :forward && odet !== nothing + # The ODE integrators' stores need materializing; doing it here keeps + # `SolutionProfiles.du_store`/`xi_s_store` populated by construction. + solution = if odet !== nothing materialize_derivative_stores!(odet, equil, mats, intr) SolutionProfiles(:el_axis, odet.step, odet.psi_store, odet.q_store, odet.u_store, odet.du_store, odet.xi_s_store) @@ -207,15 +207,6 @@ function build_result( nothing end - # One Δ′ payload whichever formalism produced it; Riccati leaves the PEST-3 blocks unfilled - # because it persists only the raw D′ and recovers them on demand via `pest3_decompose`. - delta_prime = if gal_dp !== nothing - gal_dp - elseif !isempty(intr.delta_prime_matrix) - DeltaPrimeData(intr.delta_prime_matrix, intr.delta_prime_raw, intr.delta_coil, nothing, nothing, nothing) - else - nothing - end kinetic = (kmsing=intr.kmsing, kinsing=intr.kinsing, scan_psi=intr.kinsing_scan_psi, scan_cond=intr.kinsing_scan_cond, scan_threshold=intr.kinsing_scan_threshold) diff --git a/src/ForceFreeStates/Riccati/Crossings.jl b/src/ForceFreeStates/Riccati/Crossings.jl index d0b157d4aa..50ec5e4d76 100644 --- a/src/ForceFreeStates/Riccati/Crossings.jl +++ b/src/ForceFreeStates/Riccati/Crossings.jl @@ -1,131 +1,87 @@ -# Singular-surface crossing algorithms for the Riccati/fundamental-matrix integration. +# Singular-surface crossing of the Riccati integrator, in the asymptotic basis. """ - riccati_cross_ideal_singular_surf!(odet, ctrl, equil, mats, intr, ising) + asymptotic_bases(chunks, ctrl, equil, mats, intr) -> Dict{Int,NamedTuple} -Cross a singular surface on the Riccati path: the Riccati counterpart of `cross_ideal_singular_surf!`, -without Gaussian reduction (S is bounded) and zeroing the resonant columns `ipert_res` directly. +Asymptotic basis `(; left, right)`, each N × 2N × 2, on both sides of every crossed surface, where its crossing +chunk ends. Keyed by the surface's index in `intr.sing`. +""" +function asymptotic_bases(chunks::Vector{IntegrationChunk}, ctrl::ForceFreeStatesControl, + equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, intr::ForceFreeStatesInternal) + bases = Dict{Int,@NamedTuple{left::Array{ComplexF64,3}, right::Array{ComplexF64,3}}}() + for chunk in chunks + chunk.needs_crossing || continue + singp = intr.sing[chunk.ising] + dpsi = singp.psifac - chunk.psi_end + bases[chunk.ising] = (; left=sing_get_ua(compute_sing_asymptotics(singp, ctrl, equil, mats, intr; sig=-1.0), dpsi), + right=sing_get_ua(compute_sing_asymptotics(singp, ctrl, equil, mats, intr; sig=1.0), dpsi)) + end + return bases +end -Expects `odet.u` in the bounded (U₁, U₂) form from `riccati_integrate_chunk!` with `needs_crossing=true`. -`ca_l`/`ca_r` are captured before `renormalize_riccati!` restores the (S, I) form, so the introduced -small solution has `ca_r[ipert_res, ipert_res, 2] = 1`. """ -function riccati_cross_ideal_singular_surf!( - odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, - mats::MatrixSplines, intr::ForceFreeStatesInternal, ising::Int -) - # Skip Gaussian reduction — S is bounded so no large-norm columns exist. + cross_singular_surface!(odet, intr, ising) -> Matrix{ComplexF64} + +Cross singular surface `ising` [Glasser & Kolemen 2018, Phys. Plasmas 25, 082502, Eq. 36]. `odet.u` holds +asymptotic coefficients, in the left basis on entry and the right basis on return, with no big solution on +either side; both sets are recorded in `odet.ca_l` and `odet.ca_r`. Returns the N × N matrix that takes the +columns on return to their continuation left of the surface, in the columns on entry. +""" +function cross_singular_surface!(odet::OdeState, intr::ForceFreeStatesInternal, ising::Int) singp = intr.sing[ising] - dpsi = singp.psifac - odet.psifac # ψ_res - ψ_current (positive) ipert_res = 1 .+ singp.m .- intr.mlow .+ (singp.n .- intr.nlow) .* intr.mpert - - sing_asymp_left, sing_asymp_right = _two_sided_singular_asymptotics(singp, ctrl, equil, mats, intr) - _log_riccati_crossing_diagnostics(odet, intr, ising, singp, dpsi, sing_asymp_left, sing_asymp_right) - - _capture_left_crossing_data!(odet, singp, sing_asymp_left, dpsi, intr, ising) - _predict_across_singular_surface!(odet, ctrl, equil, mats, intr, ising, ipert_res, dpsi, sing_asymp_right) - _capture_right_crossing_data!(odet, singp, sing_asymp_right, dpsi, intr, ising, ipert_res, ctrl) - - _stash_per_surface_delta_prime_stub!(odet, intr, ising, ipert_res, sing_asymp_right, equil, ctrl) - _store_crossing_step!(odet) - - # Restore canonical (S_new, I) form before continuing integration. - renormalize_riccati!(odet, intr) + continuation = eliminate_big_solutions!(odet.u, ipert_res) + odet.ca_l[:, :, :, ising] .= odet.u + start_small_solutions!(odet.u, ipert_res) + odet.ca_r[:, :, :, ising] .= odet.u + odet.ca_populated = true + # The small solutions introduced on the right do not exist on the left. + continuation[:, ipert_res] .= 0 + return continuation end """ - _two_sided_singular_asymptotics(singp, ctrl, equil, mats, intr) -> (left, right) + eliminate_big_solutions!(ca, ipert_res) -> Matrix{ComplexF64} -Left (`sig=-1`) and right (`sig=+1`) singular asymptotics, each from its own resonant block. +Column-reduce the asymptotic coefficients `ca` (N × N × 2: basis solution, solution, big or small half of the +basis) so that only column `ipert_res[i]` carries the big solution of resonance `i`. Returns the N × N column +transformation applied. """ -function _two_sided_singular_asymptotics(singp::SingType, ctrl::ForceFreeStatesControl, - equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, - intr::ForceFreeStatesInternal) - sing_asymp_right = compute_sing_asymptotics(singp, ctrl, equil, mats, intr; sig=1.0) - sing_asymp_left = compute_sing_asymptotics(singp, ctrl, equil, mats, intr; sig=-1.0) - return sing_asymp_left, sing_asymp_right -end - -# @debug-only per-crossing diagnostics. Enable via JULIA_DEBUG=GeneralizedPerturbedEquilibrium. -function _log_riccati_crossing_diagnostics(odet, intr, ising, singp, dpsi, sing_asymp_left, sing_asymp_right) - @debug begin - ipert_res_diag = 1 .+ singp.m .- intr.mlow .+ (singp.n .- intr.nlow) .* intr.mpert - msg = " ising=$ising: psi_sing=$(@sprintf("%.10f", singp.psifac)), psi_eval=$(@sprintf("%.10f", odet.psifac)), dpsi=$(@sprintf("%.10e", dpsi))\n" - msg *= " alpha_L = $(sing_asymp_left.alpha), alpha_R = $(sing_asymp_right.alpha)\n" - for ip in ipert_res_diag - msg *= " vmatL[0] big: vmat[$ip,$ip,1,1]=$(@sprintf("%.8e", real(sing_asymp_left.vmat[ip,ip,1,1]))), vmat[$ip,$ip,2,1]=$(@sprintf("%.8e", real(sing_asymp_left.vmat[ip,ip,2,1])))\n" - msg *= " vmatR[0] big: vmat[$ip,$ip,1,1]=$(@sprintf("%.8e", real(sing_asymp_right.vmat[ip,ip,1,1]))), vmat[$ip,$ip,2,1]=$(@sprintf("%.8e", real(sing_asymp_right.vmat[ip,ip,2,1])))\n" +function eliminate_big_solutions!(ca::Array{ComplexF64,3}, ipert_res) + transform = Matrix{ComplexF64}(I, size(ca, 2), size(ca, 2)) + for (i, ipert) in enumerate(ipert_res) + candidates = setdiff(axes(ca, 2), @view ipert_res[1:(i-1)]) + pivot = candidates[argmax(abs.(@view ca[ipert, candidates, 1]))] + for col in axes(ca, 2) + col == pivot && continue + factor = ca[ipert, col, 1] / ca[ipert, pivot, 1] + @views ca[:, col, :] .-= factor .* ca[:, pivot, :] + @views transform[:, col] .-= factor .* transform[:, pivot] + ca[ipert, col, 1] = 0 end - msg - end -end - -# Capture left-side asymptotic data into odet.ca_l and singp.ua_left/psi_ua_left. -function _capture_left_crossing_data!(odet::OdeState, singp::SingType, sing_asymp_left, - dpsi::Float64, intr::ForceFreeStatesInternal, ising::Int) - ua = sing_get_ua(sing_asymp_left, dpsi) - singp.ua_left = copy(ua) - singp.psi_ua_left = odet.psifac - odet.ca_l[:, :, :, ising] .= sing_get_ca(odet.u, ua, intr) -end - -# Trapezoidal predictor across the singular surface with the resonant columns zeroed (they stay zero). -function _predict_across_singular_surface!(odet::OdeState, ctrl::ForceFreeStatesControl, - equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, - intr::ForceFreeStatesInternal, ising::Int, - ipert_res, dpsi::Float64, sing_asymp_right) - if ctrl.kinetic === nothing - for i in eachindex(sing_asymp_right.r1) - odet.u[:, ipert_res[i], :] .= 0 + if pivot != ipert + swapped = ca[:, pivot, :] + @views ca[:, pivot, :] .= ca[:, ipert, :] + ca[:, ipert, :] .= swapped + transform[:, [pivot, ipert]] .= transform[:, [ipert, pivot]] end end - params = (ctrl, equil, mats, intr, odet, IntegrationChunk(0.0, 0.0, false, ising, 1)) - du1 = zeros(ComplexF64, intr.numpert_total, intr.numpert_total, 2) - du2 = zeros(ComplexF64, intr.numpert_total, intr.numpert_total, 2) - sing_der!(du1, odet.u, params, odet.psifac) - odet.psifac += 2 * dpsi # jump to other side of singular surface - sing_der!(du2, odet.u, params, odet.psifac) - odet.u .+= (du1 .+ du2) .* dpsi + return transform end -# Inject the right-side small asymptotic into the resonant columns, then capture odet.ca_r and singp.ua_right. -function _capture_right_crossing_data!(odet::OdeState, singp::SingType, sing_asymp_right, - dpsi::Float64, intr::ForceFreeStatesInternal, ising::Int, - ipert_res, ctrl::ForceFreeStatesControl) - ua = sing_get_ua(sing_asymp_right, dpsi) - singp.ua_right = copy(ua) - singp.psi_ua_right = odet.psifac - if ctrl.kinetic === nothing - for i in eachindex(sing_asymp_right.r1) - odet.u[ipert_res[i], :, :] .= 0 - odet.u[:, ipert_res[i], :] .= ua[:, ipert_res[i]+intr.numpert_total, :] - end - end - odet.ca_r[:, :, :, ising] .= sing_get_ca(odet.u, ua, intr) - odet.ca_populated = true -end +""" + start_small_solutions!(ca, ipert_res) -# STUB: per-surface ca-based Δ' (not physically valid; see SingType.delta_prime docstring). -# The canonical Δ' is intr.delta_prime_matrix from compute_delta_prime_matrix!. -function _stash_per_surface_delta_prime_stub!(odet::OdeState, intr::ForceFreeStatesInternal, - ising::Int, ipert_res, sing_asymp_right, - equil::Equilibrium.PlasmaEquilibrium, - ctrl::ForceFreeStatesControl) - ctrl.kinetic === nothing || return - denom = (2π)^2 * equil.psio - n_res = length(sing_asymp_right.r1) - N = intr.numpert_total - resize!(intr.sing[ising].delta_prime, n_res) - intr.sing[ising].delta_prime_col = zeros(ComplexF64, N, n_res) - for i in eachindex(sing_asymp_right.r1) - Δca_col = (odet.ca_r[:, ipert_res[i], 2, ising] - odet.ca_l[:, ipert_res[i], 2, ising]) / denom - intr.sing[ising].delta_prime_col[:, i] .= Δca_col - intr.sing[ising].delta_prime[i] = Δca_col[ipert_res[i]] +Turn reduced left-side coefficients `ca` into right-side ones: the resonant coefficients vanish and column +`ipert_res[i]` becomes the small solution of resonance `i`. +""" +function start_small_solutions!(ca::Array{ComplexF64,3}, ipert_res) + for ipert in ipert_res + ca[ipert, :, :] .= 0 + ca[:, ipert, :] .= 0 end -end - -# Store (U₁_new, U₂_new) into u_store before renormalization so that -# evaluate_stability_criterion! can recover S_new = U₁_new / U₂_new via compute_smallest_eigenvalue. -function _store_crossing_step!(odet::OdeState) - store_ode_data!(odet, odet.psifac, odet.u) + for ipert in ipert_res + ca[ipert, ipert, 2] = 1 + end + return ca end diff --git a/src/ForceFreeStates/Riccati/DeltaPrimeBVP.jl b/src/ForceFreeStates/Riccati/DeltaPrimeBVP.jl index fcac7b1151..d2269e0c07 100644 --- a/src/ForceFreeStates/Riccati/DeltaPrimeBVP.jl +++ b/src/ForceFreeStates/Riccati/DeltaPrimeBVP.jl @@ -1,123 +1,69 @@ # STRIDE global boundary-value problem: Delta-prime matrix assembly, solve, PEST-3 decomposition. -""" - compute_delta_prime_matrix!(intr, propagators, chunks; wv, psio, debug, S_at_surface_left, ctrl, equil, mats) - -Compute the inter-surface tearing stability matrix (msing × msing) using the -STRIDE global BVP formulation [Glasser 2018 Phys. Plasmas 25, 032501, Sec. III.B]. - -The BVP encodes the full plasma response with unknowns at each surface boundary: - -``` - x_axis (N): free IC parameters at the axis (U₁ = 0 regular solutions) - x_left[j] (2N): state at left inner-layer boundary of surface j - x_right[j] (2N): state at right inner-layer boundary of surface j - x_edge (N): free IC parameters at the edge - Total unknowns: nMat = (2 + 4·msing)·N -``` - -## Boundary conditions - -With `wv` (the singfac-scaled vacuum response) the edge BC is `U₁ = c, U₂ = -wv·ψ₀²·c`, i.e. -`wp + wv = 0`; with `wv = nothing` it is a conducting wall (`U₁ = 0`). When `S_at_surface_left` -is given (one per crossed surface), the axis is matched through a QR-conditioned axis propagator and -ua-initialized shooting propagators; otherwise through a Gaussian-reduced axis FM (STRIDE `ode_fixup`), -which requires forward crossing chunks. +# Relative difference between the two plasma response matrices above which `compute_delta_prime_matrix` warns. +const DELTA_PRIME_RESPONSE_MISMATCH = 1e-6 -## Output: PEST3-convention Δ' (deltap) - -The raw BVP solution is a 2·msing × 2·msing matrix `dp` with left/right -sub-indices at each surface, stored in `intr.delta_prime_raw`. The PEST3-convention Δ' is -[Chance, PPPL-2527]: +""" + compute_delta_prime_matrix(intr, propagators, chunks, axis_state; wv=nothing, wp=nothing, psio=0.0, extended_precision=false) -> Union{Nothing,DeltaPrimeData} -``` - deltap(i,j) = dp(2i,2j) - dp(2i,2j-1) - dp(2i-1,2j) + dp(2i-1,2j-1) -``` +Solve the STRIDE boundary-value problem for the inter-surface Δ' matrix +[Glasser & Kolemen 2018, Phys. Plasmas 25, 082502, Eq. 37]. -stored in `intr.delta_prime_matrix` (msing × msing). +### Arguments -## Limitations + - `intr`: Read for its mode space and singular surfaces. + - `propagators`, `chunks`, `axis_state`: From `riccati_eulerlagrange_integration`. + - `wv`: Singfac-scaled vacuum response for the edge condition `wp + wv = 0`; `nothing` puts a conducting wall there. + - `wp`: Plasma response from `free_run`. A warning is issued when the BVP's own differs from it by more than + `DELTA_PRIME_RESPONSE_MISMATCH`. + - `psio`: Poloidal flux normalization. + - `extended_precision`: Solve in `Complex{Double64}`. -Assumes one resonant mode per surface: if any surface carries several (a multi-`n` run where one -q satisfies two `(m, n)`, e.g. q = 2 with (2, 1) and (4, 2)), it warns and skips the BVP. +### Returns -`intr.delta_prime_matrix` is the only physically valid Δ'; the per-surface ca-based -`intr.sing[*].delta_prime` is a diagnostic placeholder and need not agree with it. +A [`DeltaPrimeData`](@ref) whose `coil` is Δ_EF (empty without `wv`), or `nothing` when no surface is crossed +or, with a warning, when a surface carries more than one resonant mode. """ -function compute_delta_prime_matrix!( +function compute_delta_prime_matrix( intr::ForceFreeStatesInternal, propagators::Vector{ChunkPropagator}, - chunks::Vector{IntegrationChunk}; + chunks::Vector{IntegrationChunk}, + axis_state::Matrix{ComplexF64}; wv::Union{Nothing,Matrix{ComplexF64}}=nothing, + wp::Union{Nothing,Matrix{ComplexF64}}=nothing, psio::Float64=0.0, - debug::Bool=false, - S_at_surface_left::Union{Nothing,Vector{Matrix{ComplexF64}}}=nothing, - ctrl::Union{Nothing,ForceFreeStatesControl}=nothing, - equil::Union{Nothing,Equilibrium.PlasmaEquilibrium}=nothing, - mats::Union{Nothing,MatrixSplines}=nothing + extended_precision::Bool=false ) - intr.msing == 0 && return - _has_unsupported_multi_resonance(intr) && return + intr.msing == 0 && return nothing + _has_unsupported_multi_resonance(intr) && return nothing sing, i_crossings, msing = _select_active_surfaces(intr, chunks) - msing == 0 && return + msing == 0 && return nothing N = intr.numpert_total - - use_S_axis = S_at_surface_left !== nothing && length(S_at_surface_left) == msing - - # The FM-axis fallback uses crossing-chunk propagators as forward maps; backward chunks would give a wrong Δ'. - if !use_S_axis - for ic in i_crossings - chunks[ic].direction == 1 || - error( - "compute_delta_prime_matrix!: FM-axis fallback (use_S_axis=false) requires forward crossing chunks; " * - "chunk $ic has direction=$(chunks[ic].direction). Either provide S_at_surface_left or use bidirectional=false." - ) - end - end - - Phi_L_mats, Phi_R_mats, Phi_R_halves = _assemble_segment_propagators( - propagators, chunks, i_crossings, msing, N, use_S_axis) - ipert_all = [1 + sing[j].m[1] - intr.mlow + (sing[j].n[1] - intr.nlow) * intr.mpert for j in 1:msing] - has_ua = all(j -> !isempty(sing[j].ua_left), 1:msing) - T_left_mats, T_right_mats, T_left_inv, T_right_inv = - _build_asymptotic_basis_matrices(sing, has_ua, N, msing) - - debug && _log_bvp_setup(chunks, sing, S_at_surface_left, use_S_axis, has_ua, - Phi_L_mats, Phi_R_mats, Phi_R_halves, ipert_all, wv, psio, N, msing) - - if use_S_axis - uShootR, uShootL, uAxis = _build_S_axis_shooting_propagators( - propagators, chunks, i_crossings, sing, msing, N, - T_left_mats, T_right_mats, has_ua, ctrl, equil, mats, intr, debug) - debug && _log_S_axis_shooting_propagators(uShootR, uShootL, uAxis, - S_at_surface_left, T_left_mats, - ipert_all, has_ua, msing, N) - M, nMat, col_edge = _assemble_bvp_S_axis( - uShootR, uShootL, uAxis, ipert_all, msing, N, wv, psio) - else - M, nMat, col_edge = _assemble_bvp_FM_axis( - Phi_L_mats, Phi_R_mats, ipert_all, msing, N, - T_left_inv, T_right_inv, has_ua, wv, psio) - end - if debug - @info "Δ' BVP: nMat=$nMat, rank(M)=$(rank(M)), cond(M)=$(@sprintf("%.2e", cond(M)))" - end - - # rpec coil-response block: needs the S-axis row layout that `_solve_bvp_edge_coil` assumes - # for its edge rows, and a vacuum edge in the assembled matrix (col_edge in junc_rows). - if use_S_axis && wv !== nothing - intr.delta_coil = _solve_bvp_edge_coil(M, col_edge, msing, N, ipert_all) + uShootR, uShootL = _shooting_matrices(propagators, chunks, i_crossings) + uAxis = Matrix(qr(axis_state).Q)[:, 1:N] + M, col_edge = _assemble_bvp(uShootR, uShootL, uAxis, ipert_all, msing, N, wv, psio) + + # Δ_EF: needs a vacuum edge in the assembled matrix. + delta_ef = Matrix{ComplexF64}(undef, 0, 0) + if wv !== nothing + delta_ef, edge_coefficients = _solve_bvp_delta_ef(M, col_edge, msing, N, ipert_all) + if wp !== nothing + # Unit edge displacement, so the edge U₂ is the plasma response times ψ₀². + wp_bvp = uShootR[msing][(N+1):(2N), :] * edge_coefficients ./ psio^2 + mismatch = norm(wp_bvp - wp) / norm(wp) + if mismatch > DELTA_PRIME_RESPONSE_MISMATCH + @warn "Riccati: the plasma response from the Δ' system differs from the axis-to-edge one by $(@sprintf("%.1e", mismatch)) " * + "(relative), so the chunk propagators have lost accuracy. Raise [Riccati] nchunks (now $(length(chunks)) chunks); " * + "beyond about 150 poloidal harmonics also set extended_precision_bvp." + end + end end - - deltap, dp_raw_persisted = _solve_bvp_and_combine_pest3( - M, msing, N, nMat, use_S_axis, ipert_all, col_edge, ctrl, debug) - - # The raw side-major D' (Fortran gal_write_delta layout) feeds `pest3_decompose`. - intr.delta_prime_matrix = deltap - intr.delta_prime_raw = dp_raw_persisted + # The raw side-major D' (Fortran gal_write_delta layout) feeds `pest3_decompose`; Riccati leaves the blocks unfilled. + delta_prime, delta_prime_raw = _solve_bvp_and_combine_pest3(M, msing, N, ipert_all, extended_precision) + return DeltaPrimeData(delta_prime, delta_prime_raw, delta_ef, nothing, nothing, nothing) end # Column index helpers for the BVP matrix. j is the 1-based singular-surface index, @@ -131,7 +77,7 @@ function _has_unsupported_multi_resonance(intr::ForceFreeStatesInternal) n_res_per_surface = [length(intr.sing[j].m) for j in 1:msing] any(>(1), n_res_per_surface) || return false offenders = [(j, intr.sing[j].m, intr.sing[j].n) for j in 1:msing if n_res_per_surface[j] > 1] - @warn "compute_delta_prime_matrix!: skipping inter-surface Δ' BVP because some surfaces carry more than one resonant mode " * + @warn "compute_delta_prime_matrix: skipping inter-surface Δ' BVP because some surfaces carry more than one resonant mode " * "(multi-n collision; generalization tracked as follow-up). " * "Per-surface Δ' is unaffected. Multi-resonance surfaces: $offenders" return true @@ -148,120 +94,29 @@ function _select_active_surfaces(intr::ForceFreeStatesInternal, chunks::Vector{I if msing_active < msing excluded = setdiff(1:msing, sing_indices) excluded_ms = [intr.sing[j].m for j in excluded] - @debug "compute_delta_prime_matrix!: $msing singular surfaces, $msing_active crossed (excluded: m=$excluded_ms)" + @debug "compute_delta_prime_matrix: $msing singular surfaces, $msing_active crossed (excluded: m=$excluded_ms)" end sing = [intr.sing[si] for si in sing_indices] return sing, i_crossings, msing_active end -# Assemble per-surface crossing FMs (Phi_L), inter-surface and edge FMs (Phi_R; Phi_R[1] only for the -# FM-axis fallback), and midpoint-split halves (Phi_R_halves) used only by the debug diagnostics. -function _assemble_segment_propagators(propagators::Vector{ChunkPropagator}, - chunks::Vector{IntegrationChunk}, - i_crossings::Vector{Int}, msing::Int, N::Int, - use_S_axis::Bool) - Phi_L_mats = [assemble_fm_matrix(propagators, i_crossings[j]:i_crossings[j]) for j in 1:msing] - Phi_R_mats = Vector{Matrix{ComplexF64}}(undef, msing + 1) - if !use_S_axis - Phi_R_mats[1] = assemble_fm_matrix(propagators, 1:(i_crossings[1]-1); condition=true) - end - for j in 2:msing - Phi_R_mats[j] = assemble_fm_matrix(propagators, (i_crossings[j-1]+1):(i_crossings[j]-1)) - end - Phi_R_mats[msing+1] = assemble_fm_matrix(propagators, (i_crossings[msing]+1):length(chunks)) - - Phi_R_halves = Vector{Tuple{Matrix{ComplexF64},Matrix{ComplexF64}}}(undef, msing - 1) - for j in 1:(msing-1) - chunk_start = i_crossings[j] + 1 - chunk_end = i_crossings[j+1] - 1 - n_chunks = chunk_end - chunk_start + 1 - if n_chunks >= 2 - i_mid = chunk_start + div(n_chunks, 2) - 1 - Phi_left_half = assemble_fm_matrix(propagators, chunk_start:i_mid) - Phi_right_half = assemble_fm_matrix(propagators, (i_mid+1):chunk_end) - Phi_R_halves[j] = (Phi_left_half, Phi_right_half) - else - Phi_R_halves[j] = (Matrix{ComplexF64}(I, 2N, 2N), Phi_R_mats[j+1]) - end - end - return Phi_L_mats, Phi_R_mats, Phi_R_halves -end - -# Asymptotic-basis transformation T = [ua[:,:,1]; ua[:,:,2]] maps (big, small) coefficients to the raw -# state; ua columns 1:N are big (z^{-α}) and N+1:2N small (z^{+α}) solutions. -function _build_asymptotic_basis_matrices(sing::Vector{SingType}, has_ua::Bool, N::Int, msing::Int) - T_left_mats = Vector{Matrix{ComplexF64}}(undef, msing) - T_right_mats = Vector{Matrix{ComplexF64}}(undef, msing) - T_left_inv = Vector{Matrix{ComplexF64}}(undef, msing) - T_right_inv = Vector{Matrix{ComplexF64}}(undef, msing) - if has_ua - for j in 1:msing - sp = sing[j] - T_left_mats[j] = [sp.ua_left[:, :, 1]; sp.ua_left[:, :, 2]] - T_right_mats[j] = [sp.ua_right[:, :, 1]; sp.ua_right[:, :, 2]] - T_left_inv[j] = inv(T_left_mats[j]) - T_right_inv[j] = inv(T_right_mats[j]) - end - end - return T_left_mats, T_right_mats, T_left_inv, T_right_inv -end - -# Build the S-axis shooting propagators uShootR (surface j right → midpoint), uShootL (surface j -# left → midpoint, backward), and the QR-conditioned axis propagator uAxis that uShootL[1] meets. -function _build_S_axis_shooting_propagators( - propagators::Vector{ChunkPropagator}, chunks::Vector{IntegrationChunk}, - i_crossings::Vector{Int}, sing::Vector{SingType}, msing::Int, N::Int, - T_left_mats::Vector{Matrix{ComplexF64}}, T_right_mats::Vector{Matrix{ComplexF64}}, - has_ua::Bool, ctrl, equil, mats, intr::ForceFreeStatesInternal, debug::Bool) - - can_reintegrate = has_ua && ctrl !== nothing && equil !== nothing && mats !== nothing +# Shooting matrices of the BVP [Glasser & Kolemen 2018, Phys. Plasmas 25, 082502, Eqs. 33-35]: `uShootR[j]` and +# `uShootL[j]` map the asymptotic coefficients right and left of surface j to the state at the matching point +# toward the next and the previous surface. +function _shooting_matrices(propagators::Vector{ChunkPropagator}, chunks::Vector{IntegrationChunk}, i_crossings::Vector{Int}) + msing = length(i_crossings) + matrix(i) = [propagators[i].u[:, :, 1]; propagators[i].u[:, :, 2]] uShootR = Vector{Matrix{ComplexF64}}(undef, msing) - uShootL = Vector{Matrix{ComplexF64}}(undef, msing) # uShootL[1] handled separately below - + uShootL = Vector{Matrix{ComplexF64}}(undef, msing) for j in 1:msing - shoot_range_R = _midpoint_shoot_range(chunks, i_crossings, j, msing; side=:right) - if debug && !isempty(shoot_range_R) - psi_surf_R = chunks[first(shoot_range_R)].psi_start - psi_mid_R = chunks[last(shoot_range_R)].psi_end - psi_ua_R = sing[j].psi_ua_right - @info " uShootR[$j]: shoot_range=$(shoot_range_R), psi_chunk=$(@sprintf("%.6f", psi_surf_R)), psi_ua=$(@sprintf("%.6f", psi_ua_R)), psi_mid=$(@sprintf("%.6f", psi_mid_R)), Δψ_fix=$(@sprintf("%.6e", psi_ua_R - psi_surf_R))" - end - if can_reintegrate && !isempty(shoot_range_R) - uShootR[j] = integrate_fm_with_ua_ic(chunks, shoot_range_R, sing[j].ua_right, - ctrl, equil, mats, intr; backward=false, psi_ua=sing[j].psi_ua_right) - else - T_init = has_ua ? T_right_mats[j] : nothing - uShootR[j] = assemble_fm_matrix(propagators, shoot_range_R; T_init=T_init) - end - - # uShootL[j>=2]: backward from surface j left to midpoint. uShootL[1] handled below. - j == 1 && continue - shoot_range_L = _midpoint_shoot_range(chunks, i_crossings, j, msing; side=:left) - if debug - psi_mid = chunks[first(shoot_range_L)].psi_start - psi_surf = chunks[last(shoot_range_L)].psi_end - psi_ua_L = sing[j].psi_ua_left - @info " uShootL[$j]: shoot_range=$(shoot_range_L), psi_mid=$(@sprintf("%.6f", psi_mid)), psi_chunk=$(@sprintf("%.6f", psi_surf)), psi_ua=$(@sprintf("%.6f", psi_ua_L)), Δψ_fix=$(@sprintf("%.6e", psi_ua_L - psi_surf))" - end - if can_reintegrate && !isempty(shoot_range_L) - uShootL[j] = integrate_fm_with_ua_ic(chunks, shoot_range_L, sing[j].ua_left, - ctrl, equil, mats, intr; backward=true, psi_ua=sing[j].psi_ua_left) - else - T_init = has_ua ? T_left_mats[j] : nothing - uShootL[j] = assemble_fm_matrix(propagators, shoot_range_L; T_init=T_init) - end - end - - uAxis, i_axis_mid = _build_conditioned_axis_propagator(propagators, i_crossings, N) - uShootL[1] = _build_uShootL_first(propagators, chunks, i_crossings, sing, - T_left_mats, has_ua, can_reintegrate, i_axis_mid, - ctrl, equil, mats, intr, N) - if debug - shoot_range_L1 = (i_axis_mid+1):(i_crossings[1]-1) - @info " Axis propagator: $(i_axis_mid) chunks, cond=$(@sprintf("%.2e", cond(uAxis)))" - @info " uShootL[1]: range=$(shoot_range_L1), cond=$(@sprintf("%.2e", cond(uShootL[1])))" - end - return uShootR, uShootL, uAxis + # The chunk after the surface starts from the right basis; the rest multiply on. + range_R = _midpoint_shoot_range(chunks, i_crossings, j, msing; side=:right) + uShootR[j] = foldl((X, i) -> matrix(i) * X, range_R[2:end]; init=matrix(first(range_R))) + # The crossing chunk carries the left basis back to its start; forward chunks before it are inverted. + range_L = j == 1 ? (1:0) : _midpoint_shoot_range(chunks, i_crossings, j, msing; side=:left) + uShootL[j] = foldl((X, i) -> matrix(i) \ X, reverse(range_L); init=matrix(i_crossings[j])) + end + return uShootR, uShootL end # Chunk range from surface j to the ψ midpoint toward the next surface (`:right`; the edge for @@ -289,59 +144,14 @@ function _midpoint_shoot_range(chunks::Vector{IntegrationChunk}, i_crossings::Ve return side === :right ? (chunk_start:i_mid_inter) : ((i_mid_inter+1):chunk_end) end -# Propagate [0; I] from the axis with a QR fixup after each chunk (Fortran ode_fixup), stopping one -# chunk before the first surface so uShootL[1] spans a single well-conditioned chunk. -function _build_conditioned_axis_propagator(propagators::Vector{ChunkPropagator}, - i_crossings::Vector{Int}, N::Int) - n_pre_cross = i_crossings[1] - 1 - i_axis_mid = max(1, n_pre_cross - 1) - uAxis = zeros(ComplexF64, 2N, N) - for i in 1:N - uAxis[N+i, i] = 1 - end - for ic in 1:i_axis_mid - prop = propagators[ic] - upper_old = uAxis[1:N, :] - lower_old = uAxis[(N+1):2N, :] - uAxis[1:N, :] .= prop.block_upper_ic[:, :, 1] * upper_old .+ prop.block_lower_ic[:, :, 1] * lower_old - uAxis[(N+1):2N, :] .= prop.block_upper_ic[:, :, 2] * upper_old .+ prop.block_lower_ic[:, :, 2] * lower_old - Q, _ = qr(uAxis) - uAxis .= Matrix(Q)[:, 1:N] - end - for j in 1:N - uAxis[:, j] ./= norm(@view uAxis[:, j]) - end - return uAxis, i_axis_mid -end - -# uShootL[1]: backward propagator from surface 1's left boundary to the axis midpoint; T_left_mats[1] -# (or the identity without ua) when that range is empty. -function _build_uShootL_first(propagators::Vector{ChunkPropagator}, - chunks::Vector{IntegrationChunk}, i_crossings::Vector{Int}, - sing::Vector{SingType}, T_left_mats::Vector{Matrix{ComplexF64}}, - has_ua::Bool, can_reintegrate::Bool, i_axis_mid::Int, - ctrl, equil, mats, intr::ForceFreeStatesInternal, N::Int) - shoot_range_L1 = (i_axis_mid+1):(i_crossings[1]-1) - if can_reintegrate && !isempty(shoot_range_L1) - return integrate_fm_with_ua_ic(chunks, shoot_range_L1, sing[1].ua_left, - ctrl, equil, mats, intr; - backward=true, psi_ua=sing[1].psi_ua_left) - elseif !isempty(shoot_range_L1) - return assemble_fm_matrix(propagators, shoot_range_L1; - T_init=has_ua ? T_left_mats[1] : nothing) - else - return has_ua ? T_left_mats[1] : Matrix{ComplexF64}(I, 2N, 2N) - end -end - -# Assemble the BVP matrix M with the axis matched to the conditioned axis propagator, avoiding the -# ill-conditioned full axis FM. Fortran-matched layout, nMat = (2 + 4·msing)·N. Returns (M, nMat, col_edge). -function _assemble_bvp_S_axis(uShootR::Vector{Matrix{ComplexF64}}, +# Assemble the BVP matrix M with the axis matched to the axis-regular solutions `uAxis` at the start of the first +# crossing chunk. Fortran-matched layout, nMat = (2 + 4·msing)·N. Returns (M, col_edge). +function _assemble_bvp(uShootR::Vector{Matrix{ComplexF64}}, uShootL::Vector{Matrix{ComplexF64}}, uAxis::Matrix{ComplexF64}, ipert_all::Vector{Int}, msing::Int, N::Int, wv::Union{Nothing,Matrix{ComplexF64}}, psio::Float64) - # STRIDE global BVP block structure [Glasser-Kolemen 2018 PoP 25, 032501 Eq. 37]. + # STRIDE global BVP block structure [Glasser & Kolemen 2018, Phys. Plasmas 25, 082502, Eq. 37]. nMat = (2 + 4 * msing) * N col_axis = 1:N col_edge = (nMat-N+1):nMat @@ -390,90 +200,35 @@ function _assemble_bvp_S_axis(uShootR::Vector{Matrix{ComplexF64}}, M[row_offset, _col_right(j, N)[ipert_j]] = 1 end @assert row_offset == nMat "Row count mismatch: expected $nMat, got $row_offset" - return M, nMat, col_edge + return M, col_edge end -# Coil-response block for the Eq. (37) edge [Glasser-Kolemen 2018 PoP 25, 032501]: impose the rpec edge BC -# (RDCON `gal_set_boundary`, unit source per poloidal mode) and read the small-solution coefficient at each -# surface side, solving all N modes with one factorization. Returns delta_coil (2·msing × N). -function _solve_bvp_edge_coil(M::Matrix{ComplexF64}, col_edge, msing::Int, N::Int, ipert_all::Vector{Int}) +# Δ_EF (2·msing × N): the small-solution coefficient at each surface side per unit displacement of each edge mode, +# from the BVP with its edge fixed. Also returns the coefficients right of the last surface (2N × N). +function _solve_bvp_delta_ef(M::Matrix{ComplexF64}, col_edge, msing::Int, N::Int, ipert_all::Vector{Int}) nMat = size(M, 1) bot = (nMat-2msing-N+1):(nMat-2msing) # Eq. (38) bottom rows, carrying the W_V block Mc = copy(M) Mc[bot, :] .= 0 Mc[bot, col_edge] .= I(N) # Dirichlet edge: the edge coefficients equal the source B = zeros(ComplexF64, nMat, N) - B[bot, :] .= I(N) # unit drive, one column per edge poloidal mode + B[bot, :] .= I(N) # unit displacement, one column per edge poloidal mode X = lu(Mc) \ B - delta_coil = zeros(ComplexF64, 2msing, N) + delta_ef = zeros(ComplexF64, 2msing, N) for j in 1:msing row_left = _col_left(j, N)[ipert_all[j]+N] row_right = _col_right(j, N)[ipert_all[j]+N] - @views delta_coil[2j-1, :] .= X[row_left, :] - @views delta_coil[2j, :] .= X[row_right, :] + @views delta_ef[2j-1, :] .= X[row_left, :] + @views delta_ef[2j, :] .= X[row_right, :] end - return delta_coil + return delta_ef, X[_col_right(msing, N), :] end -# Fallback BVP assembly with FM-based axis BC (used when no Riccati S matrices are available). -# Uses the conditioned axis propagator Phi_R[1][:,N+1:2N] in place of S-axis matching. -function _assemble_bvp_FM_axis(Phi_L_mats::Vector{Matrix{ComplexF64}}, - Phi_R_mats::Vector{Matrix{ComplexF64}}, ipert_all::Vector{Int}, - msing::Int, N::Int, - T_left_inv::Vector{Matrix{ComplexF64}}, - T_right_inv::Vector{Matrix{ComplexF64}}, has_ua::Bool, - wv::Union{Nothing,Matrix{ComplexF64}}, psio::Float64) - nMat = (2 + 4 * msing) * N - col_axis = 1:N - col_edge = (N+4N*msing+1):nMat - M = zeros(ComplexF64, nMat, nMat) - - M[1:2N, (N+1):(N+2N)] .= Phi_L_mats[1] - M[1:2N, col_axis] .= -view(Phi_R_mats[1], :, (N+1):2N) - - row_drive_base = 2N + (4N - 2) * msing - for j in 1:msing - ipert_j = ipert_all[j] - cl = _col_left(j, N) - cr = _col_right(j, N) - row_cont = 2N + (4N - 2) * (j - 1) - for i in 1:2N - if i != ipert_j && i != ipert_j + N - row_cont += 1 - M[row_cont, cl[i]] = 1 - M[row_cont, cr[i]] = -1 - end - end - junc_rows = (row_cont+1):(2N+(4N-2)*j) - if j < msing - M[junc_rows, cr] .= Phi_R_mats[j+1] - M[junc_rows, _col_left(j + 1, N)] .= -Phi_L_mats[j+1] - else - M[junc_rows, cr] .= Phi_R_mats[msing+1] - if wv !== nothing - M[junc_rows[1:N], col_edge] .= -I(N) - M[junc_rows[(N+1):end], col_edge] .= wv .* psio^2 - else - M[junc_rows[(N+1):end], col_edge] .= -I(N) - end - end - if has_ua - M[row_drive_base+2j-1, cl] .= T_left_inv[j][ipert_j, :] - M[row_drive_base+2j, cr] .= T_right_inv[j][ipert_j, :] - else - M[row_drive_base+2j-1, cl[ipert_j]] = 1 - M[row_drive_base+2j, cr[ipert_j]] = 1 - end - end - return M, nMat, col_edge -end - -# Solve the BVP and apply the PEST3 four-term combination. -function _solve_bvp_and_combine_pest3(M::Matrix{ComplexF64}, msing::Int, N::Int, nMat::Int, - use_S_axis::Bool, ipert_all::Vector{Int}, col_edge, - ctrl, debug::Bool) +# Solve the BVP and apply the PEST3 four-term combination, in Complex{Double64} when `extended_precision`. +function _solve_bvp_and_combine_pest3(M::Matrix{ComplexF64}, msing::Int, N::Int, ipert_all::Vector{Int}, extended_precision::Bool) + nMat = size(M, 1) s2 = 2 * msing - Tc = (ctrl === nothing || ctrl.algorithm.extended_precision_bvp) ? Complex{Double64} : ComplexF64 + Tc = extended_precision ? Complex{Double64} : ComplexF64 M_solve = Tc.(M) M_lu = lu(M_solve; check=false) use_lu = issuccess(M_lu) @@ -482,18 +237,13 @@ function _solve_bvp_and_combine_pest3(M::Matrix{ComplexF64}, msing::Int, N::Int, @warn "Δ' BVP: LU factorization singular (rank $(rank(M))/$nMat), using pseudo-inverse fallback" end + # One solve per surface side, driving that side's big-solution coefficient. dp_raw = zeros(Tc, s2, s2) b = zeros(Tc, nMat) - for jsing in 1:msing, side in 1:2 - dRow = 2jsing - (2 - side) + for dRow in 1:s2 fill!(b, 0) - drive_row = use_S_axis ? (nMat - s2 + dRow) : (2N + (4N - 2) * msing + dRow) - b[drive_row] = 1 + b[nMat-s2+dRow] = 1 x = use_lu ? (M_lu \ b) : (M_pinv * b) - - debug && _log_bvp_solve(x, b, M_solve, jsing, side, dRow, msing, N, - ipert_all, col_edge, use_S_axis) - for ksing in 1:msing ipert_k = ipert_all[ksing] dp_raw[dRow, 2ksing-1] = x[_col_left(ksing, N)[ipert_k+N]] @@ -501,123 +251,14 @@ function _solve_bvp_and_combine_pest3(M::Matrix{ComplexF64}, msing::Int, N::Int, end end - # PEST3 four-term combination [Chance PPPL-2527; Glasser-Kolemen 2018 PoP 25, 032501 Eq. 31]. - # Δ'[i,j] = (NW − NE − SW + SE) on each 2×2 block of dp_raw. + # PEST3 four-term combination [Chance PPPL-2527]: Δ'[i,j] = (NW − NE − SW + SE) on each 2×2 block of dp_raw. deltap = zeros(Tc, msing, msing) for i in 1:msing, j in 1:msing deltap[i, j] = dp_raw[2i, 2j] - dp_raw[2i, 2j-1] - dp_raw[2i-1, 2j] + dp_raw[2i-1, 2j-1] end - - debug && _log_bvp_pest3(dp_raw, deltap, s2, msing) return ComplexF64.(deltap), ComplexF64.(dp_raw) end -# Logging helpers for `compute_delta_prime_matrix!`. Called only when debug=true. -function _log_bvp_setup(chunks, sing, S_at_surface_left, use_S_axis, has_ua, - Phi_L_mats, Phi_R_mats, Phi_R_halves, ipert_all, wv, psio, N, msing) - @info "Δ' BVP: $(length(chunks)) chunks, $msing surfaces, N=$N" - @info "Δ' BVP: Axis BC: $(use_S_axis ? "S-based (Riccati)" : "FM-based (conditioned)")" - @info "Δ' BVP: Asymptotic basis: $(has_ua ? "available" : "NOT available (raw basis driving)")" - if use_S_axis - for j in 1:msing - @info " S_left[$j]: max=$(@sprintf("%.2e", maximum(abs, S_at_surface_left[j]))), cond=$(@sprintf("%.2e", cond(S_at_surface_left[j])))" - end - end - if has_ua - for j in 1:msing - sp = sing[j] - T_l = [sp.ua_left[:, :, 1]; sp.ua_left[:, :, 2]] - T_r = [sp.ua_right[:, :, 1]; sp.ua_right[:, :, 2]] - @info " Surface $j: cond(T_left)=$(@sprintf("%.2e", cond(T_l))), cond(T_right)=$(@sprintf("%.2e", cond(T_r)))" - ipert_j = ipert_all[j] - @info " Surface $j ua_left (ipert=$ipert_j, psi_ua_left=$(@sprintf("%.8f", sp.psi_ua_left))):" - for i in 1:min(5, N) - @info " ua($i,$ipert_j,1)=$(@sprintf("%16.8e %16.8e", real(sp.ua_left[i,ipert_j,1]), imag(sp.ua_left[i,ipert_j,1]))) ua($i,$ipert_j,2)=$(@sprintf("%16.8e %16.8e", real(sp.ua_left[i,ipert_j,2]), imag(sp.ua_left[i,ipert_j,2])))" - end - @info " small: ua(1,$(ipert_j+N),1)=$(@sprintf("%16.8e %16.8e", real(sp.ua_left[1,ipert_j+N,1]), imag(sp.ua_left[1,ipert_j+N,1])))" - end - end - for j in 1:(msing-1) - Phi_L_h, Phi_R_h = Phi_R_halves[j] - @info " Inter-surface $j→$(j+1): half_L cond=$(@sprintf("%.2e",cond(Phi_L_h))), half_R cond=$(@sprintf("%.2e",cond(Phi_R_h))), full cond=$(@sprintf("%.2e",cond(Phi_R_mats[j+1])))" - end - @info " Phi_R[$(msing+1)] (edge): cond=$(@sprintf("%.2e",cond(Phi_R_mats[msing+1])))" - for j in 1:msing - @info " Surface $j (m=$(sing[j].m[1])): ipert=$(ipert_all[j]), cond(Phi_L)=$(@sprintf("%.2e", cond(Phi_L_mats[j])))" - end - @info "Δ' BVP: Vacuum BC $(wv === nothing ? "off (conducting wall)" : "on (psio=$psio)")" - for j in 1:msing - if !isempty(sing[j].delta_prime) - @info " Surface $j ca-based Δ' = $(@sprintf("%.6f%+.6fi", real(sing[j].delta_prime[1]), imag(sing[j].delta_prime[1])))" - end - end -end - -function _log_S_axis_shooting_propagators(uShootR, uShootL, uAxis, S_at_surface_left, - T_left_mats, ipert_all, has_ua, msing, N) - @info " Shooting propagators (S-based axis BC, no axis unknowns):" - for j in 1:msing - shoot_R_str = @sprintf("%.2e", cond(uShootR[j])) - shoot_L_str = j >= 2 ? @sprintf("%.2e", cond(uShootL[j])) : "N/A (S axis BC)" - @info " uShootL[$j]: cond=$shoot_L_str, uShootR[$j]: cond=$shoot_R_str" - end - S1 = S_at_surface_left[1] - if has_ua - T1 = T_left_mats[1] - axis_BC = T1[1:N, :] - S1 * T1[(N+1):2N, :] - @info " S-axis BC matrix: cond=$(@sprintf("%.2e", cond(axis_BC)))" - end - for j in 1:msing - ipert_j = ipert_all[j] - col_norms_R = [norm(view(uShootR[j], :, k)) for k in 1:2N] - @info " uShootR[$j] column norms: min=$(@sprintf("%.2e", minimum(col_norms_R))), max=$(@sprintf("%.2e", maximum(col_norms_R)))" - @info " uShootR[$j] col ipert=$ipert_j norm=$(@sprintf("%.2e", col_norms_R[ipert_j])), col ipert+N=$(ipert_j+N) norm=$(@sprintf("%.2e", col_norms_R[ipert_j+N]))" - if j >= 2 - col_norms_L = [norm(view(uShootL[j], :, k)) for k in 1:2N] - @info " uShootL[$j] column norms: min=$(@sprintf("%.2e", minimum(col_norms_L))), max=$(@sprintf("%.2e", maximum(col_norms_L)))" - @info " uShootL[$j] col ipert=$ipert_j norm=$(@sprintf("%.2e", col_norms_L[ipert_j])), col ipert+N=$(ipert_j+N) norm=$(@sprintf("%.2e", col_norms_L[ipert_j+N]))" - end - end - for j in 1:(msing-1) - mid_block = hcat(uShootR[j], -uShootL[j+1]) - @info " Midpoint $j→$(j+1): cond([uShootR[$j] | -uShootL[$(j+1)]]) = $(@sprintf("%.2e", cond(mid_block)))" - col_norms_Ljp1 = [norm(view(uShootL[j+1], :, k)) for k in 1:2N] - @info " uShootL[$(j+1)] all col norms: $([(@sprintf("%.2e", c)) for c in col_norms_Ljp1])" - end -end - -function _log_bvp_solve(x, b, M_solve, jsing, side, dRow, msing, N, - ipert_all, col_edge, use_S_axis) - residual = norm(ComplexF64.(M_solve * x - b)) - side_str = side == 1 ? "left" : "right" - @info " BVP solve: jsing=$jsing side=$side_str (dRow=$dRow): ||Mx-b||=$(@sprintf("%.2e", residual)), ||x||=$(@sprintf("%.2e", Float64(norm(x))))" - for ks in 1:msing - ipert_ks = ipert_all[ks] - cl = _col_left(ks, N) - cr = _col_right(ks, N) - xl_big = ComplexF64(x[cl[ipert_ks]]) - xl_small = ComplexF64(x[cl[ipert_ks+N]]) - xr_big = ComplexF64(x[cr[ipert_ks]]) - xr_small = ComplexF64(x[cr[ipert_ks+N]]) - @info " surf $ks: x_left[big]=$(@sprintf("%+.4e%+.4ei", real(xl_big), imag(xl_big))), x_left[small]=$(@sprintf("%+.4e%+.4ei", real(xl_small), imag(xl_small)))" - @info " surf $ks: x_right[big]=$(@sprintf("%+.4e%+.4ei", real(xr_big), imag(xr_big))), x_right[small]=$(@sprintf("%+.4e%+.4ei", real(xr_small), imag(xr_small)))" - @info " surf $ks: ||x_left||=$(@sprintf("%.2e", Float64(norm(x[cl])))), ||x_right||=$(@sprintf("%.2e", Float64(norm(x[cr]))))" - end - if use_S_axis - @info " ||x_edge||=$(@sprintf("%.2e", Float64(norm(x[col_edge]))))" - end -end - -function _log_bvp_pest3(dp_raw, deltap, s2, msing) - @info "Δ' BVP: Full dp_raw matrix ($(s2)×$(s2)):" - for i in 1:s2 - row_str = join([@sprintf("%+.6e", Float64(real(dp_raw[i, j]))) for j in 1:s2], " ") - @info " dp_raw[$i,:] = $row_str" - end - @info "Δ' BVP: Raw dp diagonal = $([@sprintf("%.4f%+.4fi", Float64(real(dp_raw[i,i])), Float64(imag(dp_raw[i,i]))) for i in 1:s2])" - @info "Δ' BVP: deltap diagonal = $([@sprintf("%.4f%+.4fi", real(deltap[i,i]), imag(deltap[i,i])) for i in 1:msing])" -end - """ pest3_decompose(dp_raw::AbstractMatrix) -> NamedTuple{(:A, :B, :Γ, :Δ)} @@ -638,12 +279,11 @@ where `RR = dp_raw[2i, 2j]`, `RL = dp_raw[2i, 2j−1]`, No Pletzer–Dewar ½ prefactor, matching Fortran RDCON; it cancels in `det(D' − D(γ)) = 0`. -The Δ' block returned here equals `intr.delta_prime_matrix` (the m×m PEST3 -tearing projection computed inside `compute_delta_prime_matrix!`). +The Δ' block returned here equals `DeltaPrimeData.matrix`. # Arguments - - `dp_raw` — 2m×2m complex matrix (typically `intr.delta_prime_raw`). + - `dp_raw` — 2m×2m complex matrix (typically `DeltaPrimeData.raw`). # Returns diff --git a/src/ForceFreeStates/Riccati/Driver.jl b/src/ForceFreeStates/Riccati/Driver.jl index ef85b7bfa4..fb55eef23b 100644 --- a/src/ForceFreeStates/Riccati/Driver.jl +++ b/src/ForceFreeStates/Riccati/Driver.jl @@ -1,228 +1,196 @@ """ - Riccati/ - Dual Riccati reformulation of the Euler-Lagrange ODE + riccati_eulerlagrange_integration(ctrl, equil, mats, intr) -> (odet, riccati) -Tracks the dual Riccati matrix S = U₁·U₂⁻¹ = P⁻¹ (P = U₂·U₁⁻¹ the plasma response matrix), which stays -bounded near singular surfaces where U₁ and U₂ grow exponentially. Reference: Glasser (2018) Phys. Plasmas -25, 032507 — Eq. 19 (dual form), dS/dψ = w†·F̄⁻¹·w − S·Ḡ·S with w = Q − K̄·S. +Integrate the Euler-Lagrange equations from the axis to `intr.psilim` with parallel chunk propagators +[Glasser, Kolemen & Glasser 2018, Phys. Plasmas 25, 032507; Glasser & Kolemen 2018, Phys. Plasmas 25, 082502]. -The explicit Riccati ODE is not integrated directly: its quadratic term defeats relative error control when -|S| is large. Instead the linear EL ODE (`sing_der!`) is integrated in (U₁, U₂), and the exact change of -variables `renormalize_riccati_inplace!` resets the state to (S, I) whenever max(|U₁|) or max(|U₂|) exceeds -`renorm_threshold`, the Riccati analogue of Gaussian reduction. - -## Storage convention - - - During a chunk: u[:,:,1] = U₁ (starts as S), u[:,:,2] = U₂ (starts as I). - - After renormalization: u[:,:,1] = S, u[:,:,2] = I, so downstream U₁/U₂ ratios (wp in Free.jl, the - fixed-boundary criterion) read S directly. - - Axis: S(ψ₀) = 0 (U₁ = 0, U₂ = I), unless a positive `ctrl.frobenius_psi_max` at or above the start - opts into the diagonal Frobenius init (U₂ = I, per-mode U₁). - -Unlike the forward path, `riccati_integrator_callback!` renormalizes instead of Gaussian reduction, -`riccati_cross_ideal_singular_surf!` crosses without Gaussian reduction, and `transform_u!` is a no-op. -""" - -""" - riccati_eulerlagrange_integration(ctrl, equil, mats, intr) -> (odet, propagators, chunks, S_left) - -The Riccati/STRIDE integrator: solves the same system as [`forward_eulerlagrange_integration`](@ref) -by integrating chunk fundamental matrices (propagators) concurrently with `Threads.@threads`, -assembling them serially with Riccati crossings (no Gaussian reduction), and re-integrating the -outer plasma with `riccati_integrate_chunk!` for dense, bounded checkpoints. Crossing chunks are -integrated backward and inverted by LU, which keeps them well-conditioned at large N (STRIDE, -Glasser 2018 Phys. Plasmas 25, 032501). The trailing three return values feed the Δ' BVP in -`compute_delta_prime_matrix!`; this is the only branch that produces them. - -Select with `Riccati()`, or `integrator = "Riccati"` in a deck's `[ForceFreeStates]`. Requires -`singfac_min != 0`. Results do not depend on the thread count; `ctrl.algorithm.nchunks` is the only -tunable. `odet.u_store` holds chunk-endpoint Riccati states, not the EL basis -(`odet.u_store_el_basis == false`), so PerturbedEquilibrium and the HDF5 ξ datasets require `Forward()`. +Returns `odet`, whose stores hold the solution at every stored solver step in columns that are continuous in ψ +and equal (S, I), S = U₁U₂⁻¹, at the edge, and `riccati`, the inputs of `compute_delta_prime_matrix`: +`(; propagators, chunks, columns)`, with `columns` from `_assemble_propagators_serially!`. """ function riccati_eulerlagrange_integration( ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, intr::ForceFreeStatesInternal ) odet = _initialize_parallel_odet(ctrl, equil, mats, intr) - chunks, propagators, odet_proxies = _setup_parallel_chunks_and_proxies(odet, ctrl, intr) - intr.verbose && _log_parallel_start(odet, equil, chunks) - - _run_parallel_bvp_phase!(propagators, chunks, ctrl, equil, mats, intr, odet_proxies) - - # Harvest thread-local solver-step counts; the outer re-integration counts via its own callback. - odet.total_steps += sum(p.total_steps for p in odet_proxies) - - S_at_surface_left, last_crossing_step = - _assemble_propagators_serially!(odet, propagators, chunks, ctrl, equil, mats, intr) - - _reintegrate_outer_plasma!(odet, last_crossing_step, ctrl, equil, mats, intr) + chunks = balance_integration_chunks(chunk_el_integration_bounds(odet, ctrl, intr; bidirectional=true), ctrl, intr) + bases = asymptotic_bases(chunks, ctrl, equil, mats, intr) + if intr.verbose + @info " ψ = $((@sprintf "%.3f" odet.psifac)), q = $((@sprintf "%.3f" equil.profiles.q_spline(odet.psifac)))" + @info " Riccati FM: $(length(chunks)) chunks over $(Threads.nthreads()) thread$(Threads.nthreads() == 1 ? "" : "s")" + end - chunks, propagators = _handle_edge_dW_scan!(odet, chunks, propagators, ctrl, equil, mats, intr) - resize!(S_at_surface_left, intr.msing) # one per surviving surface, as the Δ' BVP requires + rtol = ctrl.algorithm.eulerlagrange_tolerance + propagators = chunk_propagators(chunks, bases, rtol, equil, mats, ModeIndices(intr)) + odet.total_steps += sum(p -> length(p.psi) - 1, propagators) + columns = _assemble_propagators_serially!(odet, propagators, chunks, ctrl, equil, intr) + + if scan_edge_dW!(odet, ctrl, equil, mats, intr) + # The edge moved in: shorten the chunk list and multiply the propagators again over the new domain. + rebuild = truncate_chunks!(chunks, intr.psilim, intr.msing) + resize!(propagators, length(chunks)) + if rebuild !== nothing + start = propagator_start(chunks, rebuild, bases, intr.numpert_total) + propagators[rebuild] = chunk_propagator(chunks[rebuild], start, rtol, equil, mats, ModeIndices(intr)) + end + fill!(odet.u, 0) + initialize_el_at_axis!(odet, ctrl, mats, equil.profiles, intr) + columns = _assemble_propagators_serially!(odet, propagators, chunks, ctrl, equil, intr) + end + compute_delta_prime_from_ca!(odet, intr, equil) - # compute_delta_prime_matrix! runs later, after free_run supplies the vacuum edge response wv. if intr.verbose @info "Evaluating fixed-boundary stability criterion" end odet.nzero = evaluate_stability_criterion!(odet, equil.profiles) - transform_u!(odet, intr) # no-op when ifix=0 (no Gaussian reduction) + restore_solution_columns!(odet, propagators, chunks, columns) - return odet, propagators, chunks, S_at_surface_left + return odet, (; propagators, chunks, columns) end -# Build odet and initialize at the magnetic axis. Same path as serial eulerlagrange_integration. +# OdeState at the axis start, with empty stores. function _initialize_parallel_odet(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, intr::ForceFreeStatesInternal) - odet = OdeState(intr.numpert_total, ctrl.algorithm.numsteps_init, ctrl.algorithm.numunorms_init, intr.msing) + odet = OdeState(intr.numpert_total, 1, 1, intr.msing) initialize_el_at_axis!(odet, ctrl, mats, equil.profiles, intr) - # No Gaussian reduction on the Riccati path. odet.new = false fill!(odet.unorm0, 1.0) return odet end -# Build the bidirectional chunk list, per-chunk propagators, and per-thread proxy OdeStates -# (sized by maxthreadid(), since threadid() can exceed nthreads() with an interactive pool). -function _setup_parallel_chunks_and_proxies(odet::OdeState, ctrl::ForceFreeStatesControl, - intr::ForceFreeStatesInternal) - # Crossing chunks get direction=-1: growing EL solutions decay backward, so Φ_bwd is well-conditioned. - base_chunks = chunk_el_integration_bounds(odet, ctrl, intr; bidirectional=true) - chunks = balance_integration_chunks(base_chunks, ctrl, intr) - N = intr.numpert_total - propagators = [ChunkPropagator(N) for _ in chunks] - odet_proxies = [OdeState(N, 1, 1, 0) for _ in 1:Threads.maxthreadid()] - return chunks, propagators, odet_proxies -end - -function _log_parallel_start(odet::OdeState, equil::Equilibrium.PlasmaEquilibrium, chunks::Vector{IntegrationChunk}) - @info " ψ = $((@sprintf "%.3f" odet.psifac)), q = $((@sprintf "%.3f" equil.profiles.q_spline(odet.psifac)))" - @info " Riccati FM: $(length(chunks)) chunks over $(Threads.nthreads()) thread$(Threads.nthreads() == 1 ? "" : "s")" -end - -# Integrate each chunk's propagator from identity IC in parallel; chunks are independent, so the -# result does not depend on the thread layout. :static makes threadid() a stable proxy index. -function _run_parallel_bvp_phase!(propagators::Vector{ChunkPropagator}, - chunks::Vector{IntegrationChunk}, - ctrl::ForceFreeStatesControl, - equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, - intr::ForceFreeStatesInternal, - odet_proxies::Vector{OdeState}) +# Every chunk's propagator, integrated across threads from inputs none of them can modify. +function chunk_propagators(chunks::Vector{IntegrationChunk}, bases::AbstractDict, rtol::Float64, + equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, modes::ModeIndices) + propagators = Vector{ChunkPropagator}(undef, length(chunks)) nblas = BLAS.get_num_threads() BLAS.set_num_threads(1) try - Threads.@threads :static for i in eachindex(chunks) - integrate_propagator_chunk!(propagators[i], chunks[i], ctrl, equil, mats, intr, - odet_proxies[Threads.threadid()]) + Threads.@threads for i in eachindex(chunks) + propagators[i] = chunk_propagator(chunks[i], propagator_start(chunks, i, bases, modes.numpert_total), rtol, equil, mats, modes) end finally BLAS.set_num_threads(nblas) end + return propagators +end + +# Solver steps of a chunk whose states are stored, in increasing ψ: every `save_interval`-th, all from +# `dW_edge_scan_start` on, and the upper end. `after_crossing` adds the first, the state right of the surface. +function stored_step_indices(prop::ChunkPropagator, chunk::IntegrationChunk, after_crossing::Bool, ctrl::ForceFreeStatesControl) + nstep = length(prop.psi) + order = chunk.direction == 1 ? ((after_crossing ? 1 : 2):nstep) : ((nstep-1):-1:1) + return [k for (i, k) in enumerate(order) if + i == length(order) || i % ctrl.algorithm.save_interval == 0 || prop.psi[k] >= ctrl.dW_edge_scan_start] end -# Apply chunk propagators serially, renormalizing to (S, I) after each (STRIDE's ode_fixup role: -# keeps conditioning at O(cond_per_chunk) instead of its K-th power). Backward chunks are applied -# via apply_propagator_inverse!. S_at_surface_left is the Riccati S at each surface's left boundary, -# the Δ' BVP axis BC. Returns (S_at_surface_left, last_crossing_step). +# Multiply the chunk propagators from the axis outward, crossing each singular surface, and fill `odet`'s stores +# with the (S, I) state at every stored step. Returns what later stages need to rebuild solutions from the +# propagators: `axis_state` (2N × N, the state entering the first crossing chunk), `starts` (per chunk, the +# 2N × N coefficients its propagator multiplies), `stored` and `slots` (per chunk, the stored solver steps and the +# store slot before its first), and `continuations` (per crossing chunk, from `cross_singular_surface!`). function _assemble_propagators_serially!(odet::OdeState, propagators::Vector{ChunkPropagator}, - chunks::Vector{IntegrationChunk}, - ctrl::ForceFreeStatesControl, - equil::Equilibrium.PlasmaEquilibrium, - mats::MatrixSplines, intr::ForceFreeStatesInternal) + chunks::Vector{IntegrationChunk}, ctrl::ForceFreeStatesControl, + equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal) N = intr.numpert_total - S_at_surface_left = Matrix{ComplexF64}[] - last_crossing_step = 1 + stored = [stored_step_indices(propagators[i], chunks[i], i > 1 && chunks[i-1].needs_crossing, ctrl) for i in eachindex(chunks)] + slots = cumsum([0; length.(stored)]) + odet.step = slots[end] + odet.psi_store = Vector{Float64}(undef, odet.step) + odet.q_store = Vector{Float64}(undef, odet.step) + odet.u_store = Array{ComplexF64}(undef, N, N, 2, odet.step) + + starts = Vector{Matrix{ComplexF64}}(undef, length(chunks)) + continuations = Dict{Int,Matrix{ComplexF64}}() + axis_state = Matrix{ComplexF64}(undef, 0, 0) for (i, chunk) in enumerate(chunks) - if chunk.direction == -1 + if chunk.needs_crossing + isempty(axis_state) && (axis_state = [odet.u[:, :, 1]; odet.u[:, :, 2]]) apply_propagator_inverse!(odet, propagators[i]) + starts[i] = [odet.u[:, :, 1]; odet.u[:, :, 2]] + continuations[i] = cross_singular_surface!(odet, intr, chunk.ising) + odet.psifac = chunks[i+1].psi_start else + starts[i] = [odet.u[:, :, 1]; odet.u[:, :, 2]] apply_propagator!(odet, propagators[i]) + renormalize_riccati_inplace!(odet.u, N) + odet.psifac = chunk.psi_end end - renormalize_riccati_inplace!(odet.u, N) - odet.psifac = chunk.psi_end - odet.q = equil.profiles.q_spline(odet.psifac) - if intr.verbose - @info " ψ = $((@sprintf "%.3f" odet.psifac)), q= $((@sprintf "%.3f" odet.q)), max(S) = $((@sprintf "%.2e" maximum(abs, odet.u[:,:,1]))), steps = $(odet.step-1)" + @info " ψ = $((@sprintf "%.3f" odet.psifac)), q= $((@sprintf "%.3f" equil.profiles.q_spline(odet.psifac))), steps = $(slots[i+1])" end + end + odet.q = equil.profiles.q_spline(odet.psifac) + for (i, prop) in enumerate(propagators), (j, k) in enumerate(stored[i]) + odet.psi_store[slots[i]+j] = prop.psi[k] + odet.q_store[slots[i]+j] = equil.profiles.q_spline(prop.psi[k]) + end + store_propagated!(odet.u_store, propagators, starts, stored, slots; reduce=true) + return (; axis_state, starts, stored, slots, continuations) +end - if chunk.needs_crossing - ctrl.kinetic !== nothing && error("Riccati() does not support kinetic runs; use Forward().") - # State is (S, I) from the renorm above — well-conditioned at the surface's left boundary. - push!(S_at_surface_left, copy(odet.u[:, :, 1])) - riccati_cross_ideal_singular_surf!(odet, ctrl, equil, mats, intr, chunk.ising) - last_crossing_step = odet.step - 1 - else - # Store the chunk-endpoint state; these are Riccati states, not the EL basis. - odet.u_store_el_basis = false - if odet.step >= size(odet.u_store, 4) - resize_storage!(odet) +# Fill `u_store` with each chunk's propagator at its stored steps times that chunk's `coefficients` (2N × N), +# reduced to (S, I) when `reduce`. +function store_propagated!(u_store::Array{ComplexF64,4}, propagators::Vector{ChunkPropagator}, + coefficients::Vector{Matrix{ComplexF64}}, stored::Vector{Vector{Int}}, slots::Vector{Int}; reduce::Bool) + N = size(u_store, 1) + nblas = BLAS.get_num_threads() + BLAS.set_num_threads(1) + try + Threads.@threads for i in eachindex(propagators) + u = zeros(ComplexF64, N, N, 2) + for (j, k) in enumerate(stored[i]) + @views for half in 1:2 + mul!(u[:, :, half], propagators[i].steps[k][:, :, half], coefficients[i]) + end + reduce && renormalize_riccati_inplace!(u, N) + @views u_store[:, :, :, slots[i]+j] .= u end - odet.psi_store[odet.step] = odet.psifac - odet.q_store[odet.step] = odet.q - @views odet.u_store[:, :, :, odet.step] .= odet.u - odet.step += 1 end + finally + BLAS.set_num_threads(nblas) end - return S_at_surface_left, last_crossing_step -end - -# Re-integrate last rational surface → psilim with Riccati: FM propagation loses precision there at -# high N, and findmax_dW_edge! needs dense checkpoints. Starts from the pre-renorm crossing state. -function _reintegrate_outer_plasma!(odet::OdeState, last_crossing_step::Int, - ctrl::ForceFreeStatesControl, - equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, - intr::ForceFreeStatesInternal) - N = intr.numpert_total - odet.u .= odet.u_store[:, :, :, last_crossing_step] - odet.psifac = odet.psi_store[last_crossing_step] - odet.q = odet.q_store[last_crossing_step] - odet.step = last_crossing_step + 1 - renormalize_riccati_inplace!(odet.u, N) - outer_chunk = IntegrationChunk(; psi_start=odet.psifac, psi_end=intr.psilim * (1 - eps), - needs_crossing=false, ising=0) - riccati_integrate_chunk!(odet, ctrl, equil, mats, intr, outer_chunk) - # Post: odet.u is in (S, I) form; odet.step points to next empty slot. end -# Edge-dW scan; on truncation, also cut the Δ' chunks and propagators back to the new edge. -function _handle_edge_dW_scan!(odet::OdeState, chunks::Vector{IntegrationChunk}, - propagators::Vector{ChunkPropagator}, - ctrl::ForceFreeStatesControl, - equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, - intr::ForceFreeStatesInternal) - N = intr.numpert_total - odet.step -= 1 - trim_storage!(odet) - scan_edge_dW!(odet, ctrl, equil, mats, intr) || return chunks, propagators - - renormalize_riccati_inplace!(odet.u, N) # stored snapshot may be pre-renorm - rebuild = truncate_chunks!(chunks, intr.psilim, intr.msing) - resize!(propagators, length(chunks)) - if rebuild !== nothing - integrate_propagator_chunk!(propagators[rebuild], chunks[rebuild], ctrl, equil, mats, intr, OdeState(N, 1, 1, 0)) +# Rewrite `odet.u_store` from (S, I) into solution columns that are continuous in ψ and equal (S, I) at the edge: +# what `transform_u!` leaves for the forward integrator. +function restore_solution_columns!(odet::OdeState, propagators::Vector{ChunkPropagator}, + chunks::Vector{IntegrationChunk}, columns) + (; starts, stored, slots, continuations) = columns + N = size(odet.u, 1) + edge_columns = Matrix{ComplexF64}(I, N, N) # each chunk's start columns in terms of the edge columns + coefficients = similar(starts) + for i in reverse(eachindex(chunks)) + if chunks[i].needs_crossing + edge_columns = continuations[i] * edge_columns + else + chunk_end = [propagators[i].u[:, :, 1]; propagators[i].u[:, :, 2]] * starts[i] + edge_columns = chunk_end[(N+1):(2N), :] \ edge_columns + end + coefficients[i] = starts[i] * edge_columns end - return chunks, propagators + store_propagated!(odet.u_store, propagators, coefficients, stored, slots; reduce=false) end """ truncate_chunks!(chunks, psi_edge, msing) -> Union{Int,Nothing} -Cut `chunks` back to an edge at `psi_edge`: drop those past it, end the straddling one at it, and -clear any crossing onto a surface no longer in the plasma. Returns the index of a chunk that was -shortened and so must be re-integrated. +Cut `chunks` back to an edge at `psi_edge`: drop those past it, end the straddling one at it, and turn a +crossing chunk whose surface is no longer in the plasma into a forward chunk. An edge at the right boundary of +a surface keeps the chunk after its crossing, at zero width. Returns the index of a chunk that changed and so +must be re-integrated. """ function truncate_chunks!(chunks::Vector{IntegrationChunk}, psi_edge::Real, msing::Int) - last = findlast(c -> c.psi_start < psi_edge, chunks) + last = findlast(i -> chunks[i].psi_start < psi_edge || (i > 1 && chunks[i-1].needs_crossing && chunks[i].psi_start == psi_edge), + eachindex(chunks)) last === nothing && error("truncate_chunks!: edge ψ=$psi_edge lies before all chunk starts") resize!(chunks, last) c = chunks[last] - shortened = c.psi_end > psi_edge - # ising <= msing means the surface survived: edge truncation drops surfaces from the outside only. + # Truncation drops surfaces from the outside, so ising <= msing means the surface survived. crosses = c.needs_crossing && c.ising <= msing - if shortened || crosses != c.needs_crossing - chunks[last] = IntegrationChunk(; psi_start=c.psi_start, psi_end=min(c.psi_end, psi_edge), - needs_crossing=crosses, ising=crosses ? c.ising : 0, direction=c.direction) - end - return shortened ? last : nothing + (c.psi_end > psi_edge || crosses != c.needs_crossing) || return nothing + chunks[last] = IntegrationChunk(; psi_start=c.psi_start, psi_end=min(c.psi_end, psi_edge), + needs_crossing=crosses, ising=crosses ? c.ising : 0, direction=crosses ? c.direction : 1) + return last end diff --git a/src/ForceFreeStates/Riccati/Propagators.jl b/src/ForceFreeStates/Riccati/Propagators.jl index 81bd3d831c..65662fd779 100644 --- a/src/ForceFreeStates/Riccati/Propagators.jl +++ b/src/ForceFreeStates/Riccati/Propagators.jl @@ -1,249 +1,9 @@ -# Chunk-propagator integration: EL-ODE fundamental matrices, Riccati renormalization, assembly. - -# `riccati_integrator_callback!` saves every step within this fraction of the span (or absolute ψ -# floor, for short chunks) of a segment's right end; elsewhere every `save_interval`-th step. -const SAVE_NEAR_END_FRAC = 0.05 -const SAVE_NEAR_END_PSI = 1e-4 - -""" - assemble_fm_matrix(propagators, idx_range; condition=false, T_init=nothing) -> Matrix{ComplexF64} - -Assemble the 2N×2N fundamental matrix Φ_end * ... * Φ_start * T_init over chunks `idx_range`, -mapping the IC at the start of `idx_range[1]` to the state at the end of `idx_range[end]` -(`T_init` defaults to the identity). - -Each `ChunkPropagator` stores the 2N columns of Φ split into two N×N×2 blocks: - -``` - block_upper_ic[:,:,1:2] ↔ Φ[:,1:N] (result from IC=(I,0)) - block_lower_ic[:,:,1:2] ↔ Φ[:,N+1:2N] (result from IC=(0,I)) -``` - -`condition=true` applies `condition_propagator!` after each product (STRIDE's `ode_fixup`), so only -the N U₂ columns survive. Use it only for the axis→first-surface segment, where the axis BC (U₁ = 0) -needs only U₂ ICs; never for inter-surface segments, where both halves carry information. -""" -function assemble_fm_matrix(propagators::Vector{ChunkPropagator}, idx_range; - condition::Bool=false, - T_init::Union{Nothing,Matrix{ComplexF64}}=nothing) - # Take N from T_init when given, so an empty idx_range (or propagator list) still works. - N = if T_init !== nothing - size(T_init, 1) ÷ 2 - elseif !isempty(idx_range) - size(propagators[first(idx_range)].block_upper_ic, 1) - else - @assert !isempty(propagators) "assemble_fm_matrix: cannot infer N from empty propagators with no T_init" - size(propagators[1].block_upper_ic, 1) - end - Phi = T_init !== nothing ? copy(T_init) : Matrix{ComplexF64}(I, 2N, 2N) - isempty(idx_range) && return Phi - for i in idx_range - p = propagators[i] - #! format: off - Phi_i = [p.block_upper_ic[:,:,1] p.block_lower_ic[:,:,1]; - p.block_upper_ic[:,:,2] p.block_lower_ic[:,:,2]] - #! format: on - Phi = Phi_i * Phi - if condition - condition_propagator!(Phi, N) - end - end - return Phi -end - -""" - condition_propagator!(Phi, N) - -Apply Gaussian reduction to the U₂-columns (columns N+1:2N) of a 2N×2N propagator -matrix in-place, following STRIDE's `ode_fixup` convention. Triangularizes the U₁ -(upper N rows) subblock by pivoted elimination, improving the condition number so -the propagator can be used in a BVP without losing numerical rank. - -After conditioning, only the U₂ columns carry meaningful information; the U₁ columns -(1:N) are zeroed. The BVP axis block uses `Phi[:, N+1:2N]` (the conditioned half). -""" -function condition_propagator!(Phi::Matrix{ComplexF64}, N::Int) - # Work on the right half: columns N+1:2N (U₂ initial conditions) - cols = view(Phi, :, N+1:2N) - - # Sort columns by norm of the U₁ (upper N) block — largest first - norms = [norm(view(cols, 1:N, k)) for k in 1:N] - order = sortperm(norms; rev=true) - - mask_col = trues(N) # which columns remain to process - mask_row = trues(N) # which pivot rows remain available - - for isol in 1:N - kcol = order[isol] - mask_col[kcol] = false - - # Find best pivot row (largest |element| among unmasked rows) - best_row = 0 - best_val = 0.0 - for r in 1:N - if mask_row[r] && abs(cols[r, kcol]) > best_val - best_val = abs(cols[r, kcol]) - best_row = r - end - end - if best_row == 0 || best_val == 0 - continue - end - mask_row[best_row] = false - - # Eliminate this pivot from all other unmasked columns - pivot = cols[best_row, kcol] - for jcol in 1:N - if mask_col[jcol] - factor = -cols[best_row, jcol] / pivot - @views cols[:, jcol] .+= factor .* cols[:, kcol] - cols[best_row, jcol] = 0 # exact zero - end - end - end - - # Zero the U₁ columns (left half) — they are no longer meaningful - Phi[:, 1:N] .= 0 - return Phi -end - -""" - riccati_der!(du, u, params, psieval) - -Evaluate the explicit dual Riccati ODE right-hand side: -dS/dψ = w†·F̄⁻¹·w - S·Ḡ·S, w = Q - K̄·S - -where Q = diag(1/(m - n·q)) is the diagonal singular factor matrix. -The identity slice u[:,:,2] = I does not evolve (du[:,:,2] = 0). - -Reference implementation, not called in production: the quadratic S·Ḡ·S term makes the explicit -Riccati ODE unstable for explicit solvers, so the driver integrates `sing_der!` with -`renormalize_riccati_inplace!` instead. Exercised only by unit tests. - -See: Glasser (2018) Phys. Plasmas 25, 032507 — Eq. 19 (dual Riccati form) -""" -@with_pool pool function riccati_der!( - du::Array{ComplexF64,3}, - u::Array{ComplexF64,3}, - params::Tuple{ForceFreeStatesControl,Equilibrium.PlasmaEquilibrium, - MatrixSplines,ForceFreeStatesInternal,OdeState,IntegrationChunk}, - psieval::Float64 -) - - _, equil, mats, intr, odet, _ = params - - Npert = intr.numpert_total - S = @view u[:, :, 1] - dS = @view du[:, :, 1] - @view(du[:, :, 2]) .= 0 # identity does not evolve - - # Compute singfac = 1/(m - n·q) as column vector Q = diag(singfac_vec) - # [Glasser 2016 eq. 24] - singfac_vec = acquire!(pool, Float64, Npert) - singfac_mat = reshape(singfac_vec, intr.mpert, intr.npert) - odet.q = equil.profiles.q_spline(psieval; hint=odet.spline_hint) - singfac_mat .= 1.0 ./ ((intr.mlow:intr.mhigh) .- odet.q .* (intr.nlow:intr.nhigh)') - - fmat_lower = acquire!(pool, ComplexF64, Npert, Npert) - kmat = similar!(pool, fmat_lower) - gmat = similar!(pool, fmat_lower) - w = similar!(pool, fmat_lower) # w = Q - K̄·S - v = similar!(pool, fmat_lower) # v = F̄⁻¹·w (then reused for S·Ḡ·S) - tmp = similar!(pool, fmat_lower) # scratch - - # Evaluate F̄ (Cholesky factor), K̄, Ḡ splines at current ψ - mats.ideal.F_spline_lower(vec(fmat_lower), psieval; hint=mats._hint) - mats.ideal.K_spline(vec(kmat), psieval; hint=mats._hint) - mats.ideal.G_spline(vec(gmat), psieval; hint=mats._hint) - - # w = Q - K̄·S; Q is diagonal, so singfac_vec is added on the diagonal only. - mul!(w, kmat, S) # w = K̄·S - @. w = -w # w = -K̄·S - for i in 1:Npert - @inbounds w[i, i] += singfac_vec[i] # add diagonal Q: w = Q - K̄·S - end - - # v = F̄⁻¹·w (in-place Cholesky solve with stored lower-triangular factor) - v .= w - ldiv!(LowerTriangular(fmat_lower), v) - ldiv!(UpperTriangular(fmat_lower'), v) - - # dS = w†·v - S·Ḡ·S [Glasser 2018 eq. 19, dual Riccati] - mul!(dS, adjoint(w), v) # dS = w†·v - - # Subtract S·Ḡ·S (reuse v and tmp to avoid extra allocation) - mul!(tmp, gmat, S) # tmp = Ḡ·S - mul!(v, S, tmp) # v = S·Ḡ·S - dS .-= v -end - -""" - riccati_integrator_callback!(integrator) - -Per-step callback for `riccati_integrate_chunk!`: counts the step, applies -`renormalize_riccati_inplace!` when max(|U₁|) or max(|U₂|) exceeds `ctrl.algorithm.renorm_threshold` -(the Riccati analogue of Gaussian reduction), refreshes the per-column abstol, and stores selected steps. -""" -function riccati_integrator_callback!(integrator) - - ctrl, _, _, intr, odet, chunk = integrator.p - - odet.total_steps += 1 # count every accepted solver step (saved or not), as segment_callback! does - - # Renormalize to (S, I) when either block exceeds renorm_threshold. - if maximum(abs, @view(integrator.u[:, :, 1])) > ctrl.algorithm.renorm_threshold || - maximum(abs, @view(integrator.u[:, :, 2])) > ctrl.algorithm.renorm_threshold - renormalize_riccati_inplace!(integrator.u, intr.numpert_total) - end - - # Refresh after any renorm so the tolerance matches the current state. - column_abstol!(integrator.opts.abstol, integrator.u, ctrl.algorithm.eulerlagrange_tolerance) - - # Save the first two steps, every step near the right end, and every save_interval-th step otherwise. - psi_range = abs(integrator.sol.prob.tspan[2] - integrator.sol.prob.tspan[1]) - psi_remaining = abs(integrator.sol.prob.tspan[2] - integrator.t) - near_end = psi_remaining < SAVE_NEAR_END_FRAC * psi_range || psi_remaining < SAVE_NEAR_END_PSI - steps_in_segment = length(integrator.sol.t) - near_start = steps_in_segment <= 2 - should_save = near_start || near_end || (odet.step % ctrl.algorithm.save_interval == 0) - - if should_save - store_ode_data!(odet, integrator.t, integrator.u) - end -end - -""" - riccati_integrate_chunk!(odet, ctrl, equil, mats, intr, chunk) - -Integrate the Riccati state from `chunk.psi_start` to `chunk.psi_end` with `sing_der!` and -`riccati_integrator_callback!`, starting from u = (S, I). The end state is renormalized to (S, I) -unless `chunk.needs_crossing`, in which case the bounded (U₁, U₂) is left for the crossing. -""" -function riccati_integrate_chunk!( - odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, - mats::MatrixSplines, intr::ForceFreeStatesInternal, chunk::IntegrationChunk -) - cb = DiscreteCallback((u, t, integrator) -> true, riccati_integrator_callback!) - rtol = ctrl.algorithm.eulerlagrange_tolerance - prob = ODEProblem(sing_der!, odet.u, (chunk.psi_start, chunk.psi_end), - (ctrl, equil, mats, intr, odet, chunk)) - sol = solve_el(prob, rtol; callback=cb) - odet.u .= sol.u[end] - odet.psifac = sol.t[end] - # Before a crossing, skip the renorm so ca_l and ca_r share the pre-renorm normalization (needed for Δ'). - if !chunk.needs_crossing - renormalize_riccati_inplace!(odet.u, intr.numpert_total) - end -end +# Chunk propagators: their integration and application, and the reduction of a state to (S, I). """ renormalize_riccati_inplace!(u, N) -In-place Riccati renormalization on an arbitrary N×N×2 array: -u[:,:,1] = U₁ · U₂⁻¹ (new S) -u[:,:,2] = I - -The Riccati analogue of Gaussian reduction: keeps the inputs to `sing_der!` bounded. +Reduce the N × N × 2 state `u` = (U₁, U₂) in place to (S, I), S = U₁U₂⁻¹. """ function renormalize_riccati_inplace!(u::Array{ComplexF64,3}, N::Int) U2_copy = copy(@view u[:, :, 2]) @@ -255,176 +15,63 @@ function renormalize_riccati_inplace!(u::Array{ComplexF64,3}, N::Int) end """ - renormalize_riccati!(odet, intr) - -[`renormalize_riccati_inplace!`](@ref) on `odet.u`. After a singular surface crossing, which leaves -u[:,:,2] = U₂_new ≠ I, it restores the canonical storage u[:,:,1] = S_new = U₁_new·U₂_new⁻¹, u[:,:,2] = I -before the integration continues. The crossing's u_store entry keeps U₁_new and U₂_new, so -`compute_smallest_eigenvalue` is unaffected. -""" -renormalize_riccati!(odet::OdeState, intr::ForceFreeStatesInternal) = renormalize_riccati_inplace!(odet.u, intr.numpert_total) - -""" - integrate_propagator_chunk!(prop, chunk, ctrl, equil, mats, intr, odet_proxy) + propagator_start(chunks, i, bases, N) -> Array{ComplexF64,3} -Fill `prop` with the chunk's fundamental matrix by integrating the EL ODE from IC (I, 0) into -`prop.block_upper_ic` and from (0, I) into `prop.block_lower_ic`, without renormalization or storage. -`odet_proxy` holds the `sing_der!` side effects; concurrent calls need distinct proxies. +Start basis (N × 2N × 2) of chunk `i`'s propagator: the left basis of its surface for a crossing chunk, the +right basis for the chunk after one [Glasser & Kolemen 2018, Phys. Plasmas 25, 082502, Eq. 33], else the identity. """ -function integrate_propagator_chunk!( - prop::ChunkPropagator, - chunk::IntegrationChunk, - ctrl::ForceFreeStatesControl, - equil::Equilibrium.PlasmaEquilibrium, - mats::MatrixSplines, - intr::ForceFreeStatesInternal, - odet_proxy::OdeState -) - N = intr.numpert_total - # Backward chunks integrate psi_end → psi_start, where forward-growing solutions decay. - tspan = chunk.direction == 1 ? - (chunk.psi_start, chunk.psi_end) : - (chunk.psi_end, chunk.psi_start) - rtol = ctrl.algorithm.eulerlagrange_tolerance - params = (ctrl, equil, mats, intr, odet_proxy, chunk) - - # Upper block IC: U₁ = I, U₂ = 0 - u_upper = zeros(ComplexF64, N, N, 2) - for i in 1:N - u_upper[i, i, 1] = 1 - end - odet_proxy.spline_hint[] = 1 - odet_proxy.mats_hint[] = 1 - prob = ODEProblem(sing_der!, u_upper, tspan, params) - sol = solve_el(prob, rtol) - prop.block_upper_ic .= sol.u[end] - odet_proxy.total_steps += sol.stats.naccept # thread-local; summed into odet after the BVP barrier - - # Lower block IC: U₁ = 0, U₂ = I - u_lower = zeros(ComplexF64, N, N, 2) - for i in 1:N - u_lower[i, i, 2] = 1 +function propagator_start(chunks::Vector{IntegrationChunk}, i::Int, bases::AbstractDict, N::Int) + chunks[i].needs_crossing && return bases[chunks[i].ising].left + i > 1 && chunks[i-1].needs_crossing && return bases[chunks[i-1].ising].right + u0 = zeros(ComplexF64, N, 2N, 2) + for j in 1:N + u0[j, j, 1] = 1 + u0[j, N+j, 2] = 1 end - odet_proxy.spline_hint[] = 1 - odet_proxy.mats_hint[] = 1 - prob = ODEProblem(sing_der!, u_lower, tspan, params) - sol = solve_el(prob, rtol) - prop.block_lower_ic .= sol.u[end] - odet_proxy.total_steps += sol.stats.naccept + return u0 end """ - integrate_fm_with_ua_ic(chunks, chunk_range, ua, ctrl, equil, mats, intr; - backward=false, psi_ua=NaN) -> Matrix{ComplexF64} + chunk_propagator(chunk, u0, rtol, equil, mats, modes) -> ChunkPropagator -Re-integrate a span of chunks from the asymptotic solutions `ua` as initial conditions (Fortran -STRIDE `uFM_sing_init`). Column j of the 2N×2N result is the solution at the far end with IC -column j of T = [ua[:,:,1]; ua[:,:,2]]. With `backward=true` the IC sits at the span's right end -and the integration runs toward its start. `psi_ua`, when given, replaces the chunk boundary as the -IC location. Integrating from `ua` keeps the small solutions that multiplying an identity-IC -propagator by T loses to roundoff. +Integrate the ideal EL ODE over `chunk` from the start basis `u0` at relative tolerance `rtol`, keeping every +solver step. """ -function integrate_fm_with_ua_ic( - chunks::Vector{IntegrationChunk}, - chunk_range::UnitRange{Int}, - ua::Array{ComplexF64,3}, - ctrl::ForceFreeStatesControl, - equil::Equilibrium.PlasmaEquilibrium, - mats::MatrixSplines, - intr::ForceFreeStatesInternal; - backward::Bool=false, - psi_ua::Float64=NaN -) - N = intr.numpert_total - psi_start = chunks[first(chunk_range)].psi_start - psi_end = chunks[last(chunk_range)].psi_end - # ua lives at the inner-layer boundary psi_ua, which may differ slightly from the chunk boundary. - if backward && !isnan(psi_ua) - psi_end = psi_ua - elseif !backward && !isnan(psi_ua) - psi_start = psi_ua +function chunk_propagator(chunk::IntegrationChunk, u0::Array{ComplexF64,3}, rtol::Float64, + equil::Equilibrium.PlasmaEquilibrium, mats::MatrixSplines, modes::ModeIndices) + # A zero-width chunk only carries its start basis. + chunk.psi_start == chunk.psi_end && return ChunkPropagator(copy(u0), [chunk.psi_start], [copy(u0)]) + tspan = chunk.direction == 1 ? (chunk.psi_start, chunk.psi_end) : (chunk.psi_end, chunk.psi_start) + spline_hint, mats_hint = Ref(1), Ref(1) + function derivative!(du, u, _, psi) + el_derivatives!(du, u, false, equil, mats, modes, psi, spline_hint, mats_hint) + return nothing end - tspan = backward ? (psi_end, psi_start) : (psi_start, psi_end) - rtol = ctrl.algorithm.eulerlagrange_tolerance - - result = zeros(ComplexF64, 2N, 2N) - odet_proxy = OdeState(N, 1, 1, 0) - dummy_chunk = IntegrationChunk(psi_start, psi_end, false, 0, backward ? -1 : 1) - params = (ctrl, equil, mats, intr, odet_proxy, dummy_chunk) - - # Batch 1: columns 1:N of T (big solutions) - u0 = zeros(ComplexF64, N, N, 2) - u0[:, :, 1] .= ua[:, 1:N, 1] - u0[:, :, 2] .= ua[:, 1:N, 2] - odet_proxy.spline_hint[] = 1 - odet_proxy.mats_hint[] = 1 - prob = ODEProblem(sing_der!, u0, tspan, params) - sol = solve_el(prob, rtol) - result[1:N, 1:N] .= sol.u[end][:, :, 1] - result[N+1:2N, 1:N] .= sol.u[end][:, :, 2] - - # Batch 2: columns N+1:2N of T (small solutions) - u0[:, :, 1] .= ua[:, N+1:2N, 1] - u0[:, :, 2] .= ua[:, N+1:2N, 2] - odet_proxy.spline_hint[] = 1 - odet_proxy.mats_hint[] = 1 - prob = ODEProblem(sing_der!, u0, tspan, params) - sol = solve_el(prob, rtol) - result[1:N, N+1:2N] .= sol.u[end][:, :, 1] - result[N+1:2N, N+1:2N] .= sol.u[end][:, :, 2] - - return result + sol = solve_el(ODEProblem(derivative!, copy(u0), tspan), rtol; save_everystep=true) + return ChunkPropagator(sol.u[end], sol.t, sol.u) end """ apply_propagator!(odet, prop) -Apply the chunk propagator `prop` to any state `odet.u` (including the Riccati (S, I) form) in-place. -Implements the subpropagator composition Φ(ψ₂, ψ₀) = Φ(ψ₂, ψ₁) · Φ(ψ₁, ψ₀) of -Glasser-Kolemen (2018) Phys. Plasmas 25, 032501 Eq. 29. +Multiply `odet.u` by the chunk propagator `prop` [Glasser, Kolemen & Glasser 2018, Phys. Plasmas 25, 032507, Eq. 29]. """ function apply_propagator!(odet::OdeState, prop::ChunkPropagator) - U1_upper = @view prop.block_upper_ic[:, :, 1] - U2_upper = @view prop.block_upper_ic[:, :, 2] - U1_lower = @view prop.block_lower_ic[:, :, 1] - U2_lower = @view prop.block_lower_ic[:, :, 2] - - u1_prev = copy(@view odet.u[:, :, 1]) - u2_prev = copy(@view odet.u[:, :, 2]) - tmp = similar(u1_prev) - - # U₁_new = U1_upper · u1_prev + U1_lower · u2_prev - mul!(view(odet.u, :, :, 1), U1_upper, u1_prev) - mul!(tmp, U1_lower, u2_prev) - odet.u[:, :, 1] .+= tmp - - # U₂_new = U2_upper · u1_prev + U2_lower · u2_prev - mul!(view(odet.u, :, :, 2), U2_upper, u1_prev) - mul!(tmp, U2_lower, u2_prev) - odet.u[:, :, 2] .+= tmp + x = [odet.u[:, :, 1]; odet.u[:, :, 2]] + @views for k in 1:2 + mul!(odet.u[:, :, k], prop.u[:, :, k], x) + end end """ apply_propagator_inverse!(odet, prop) -Apply the *inverse* of the chunk propagator `prop` to the current state `odet.u` in-place. - -Used for backward chunks: the well-conditioned Φ_bwd maps `psi_end` → `psi_start`, so solving -Φ_bwd · x = u_old by LU advances the state `psi_start` → `psi_end` far more accurately than Φ_fwd. - -Implements the inverse subpropagator identity Φ(ψ₂, ψ₁) = Φ(ψ₁, ψ₂)⁻¹ of -Glasser-Kolemen (2018) Phys. Plasmas 25, 032501 Eq. 33. +Solve the chunk propagator `prop` against `odet.u`, which carries a state across a backward chunk +[Glasser, Kolemen & Glasser 2018, Phys. Plasmas 25, 032507, Eq. 33]. """ function apply_propagator_inverse!(odet::OdeState, prop::ChunkPropagator) N = size(odet.u, 1) - # Assemble 2N×2N backward FM Φ_bwd - #! format: off - Φ = [prop.block_upper_ic[:,:,1] prop.block_lower_ic[:,:,1]; - prop.block_upper_ic[:,:,2] prop.block_lower_ic[:,:,2]] - #! format: on - # x = Φ_bwd⁻¹ · [U₁_old; U₂_old] advances psi_start → psi_end. - u_old = [odet.u[:, :, 1]; odet.u[:, :, 2]] # 2N × N - u_new = Φ \ u_old # LU solve, 2N × N - odet.u[:, :, 1] .= u_new[1:N, :] - odet.u[:, :, 2] .= u_new[N+1:2N, :] + x = [prop.u[:, :, 1]; prop.u[:, :, 2]] \ [odet.u[:, :, 1]; odet.u[:, :, 2]] + odet.u[:, :, 1] .= x[1:N, :] + odet.u[:, :, 2] .= x[(N+1):(2N), :] end diff --git a/src/ForceFreeStates/Riccati/Types.jl b/src/ForceFreeStates/Riccati/Types.jl index 04fd25355a..8ea1b70495 100644 --- a/src/ForceFreeStates/Riccati/Types.jl +++ b/src/ForceFreeStates/Riccati/Types.jl @@ -25,22 +25,17 @@ end """ ChunkPropagator -Fundamental matrix for one integration chunk, stored as two N×N×2 solution blocks. -Represents the propagator Φ(ψ₂,ψ₁) computed by integrating the EL ODE from two -identity-block initial conditions: +Fundamental matrix of one integration chunk: column `j` is the solution started from column `j` of the chunk's +start basis (`propagator_start`). - - `block_upper_ic`: result of integrating with IC = (I_N, 0_N) (U₁ = I, U₂ = 0) - - `block_lower_ic`: result of integrating with IC = (0_N, I_N) (U₁ = 0, U₂ = I) - -Applying the propagator to the current state `u_prev`: - -u₁_new = block_upper_ic[:,:,1] · u₁_prev + block_lower_ic[:,:,1] · u₂_prev -u₂_new = block_upper_ic[:,:,2] · u₁_prev + block_lower_ic[:,:,2] · u₂_prev +## Fields -Each chunk starts from an identity IC rather than the accumulated state, so chunks integrate in parallel. + - `u::Array{ComplexF64,3}` - The matrix at the end of the integration, `(N, 2N, 2)`: U₁ in slice 1, U₂ in slice 2. + - `psi::Vector{Float64}` - ψ of every solver step, in integration order. + - `steps::Vector{Array{ComplexF64,3}}` - The matrix at each `psi`. """ struct ChunkPropagator - block_upper_ic::Array{ComplexF64,3} # shape (N, N, 2) — result from IC = (I, 0) - block_lower_ic::Array{ComplexF64,3} # shape (N, N, 2) — result from IC = (0, I) + u::Array{ComplexF64,3} + psi::Vector{Float64} + steps::Vector{Array{ComplexF64,3}} end -ChunkPropagator(N::Int) = ChunkPropagator(zeros(ComplexF64, N, N, 2), zeros(ComplexF64, N, N, 2)) diff --git a/src/ForceFreeStates/Run.jl b/src/ForceFreeStates/Run.jl index 7e440b1f59..7fee4f550b 100644 --- a/src/ForceFreeStates/Run.jl +++ b/src/ForceFreeStates/Run.jl @@ -61,20 +61,20 @@ function run_force_free_states( odet = nothing free_energies = nothing gal_data = nothing - gal_dp = nothing + delta_prime = nothing if ctrl.algorithm isa Galerkin ctrl.kinetic === nothing || error("Galerkin() does not support kinetic runs; use Forward().") gal_start = time() wv = ctrl.vac_flag ? first(compute_scaled_wv(ctrl, equil, intr)) : nothing - gal_data, gal_dp = galerkin_solve(ctrl, equil, mats, intr; wv=wv, cut_solution=cut_solution) + gal_data, delta_prime = galerkin_solve(ctrl, equil, mats, intr; wv=wv, cut_solution=cut_solution) gal_dt = time() - gal_start @info "Galerkin solve completed in $(@sprintf("%.3f", gal_dt)) s" else if intr.verbose @info "Integrating Euler-Lagrange equation" end - odet, fm_propagators, fm_chunks, fm_S_left = eulerlagrange_integration(ctrl, equil, mats, intr) + odet, riccati = eulerlagrange_integration(ctrl, equil, mats, intr) if odet.nzero > 0 && intr.verbose @warn "Fixed-boundary mode unstable for n = $nstring" end @@ -97,20 +97,18 @@ function run_force_free_states( end # Inter-surface Δ' (STRIDE BVP): needs Riccati() propagators and wv from free_run. - if ctrl.kinetic === nothing && intr.msing > 0 && fm_propagators !== nothing + if ctrl.kinetic === nothing && intr.msing > 0 && riccati !== nothing if intr.verbose @info "Computing Δ' matrix (STRIDE BVP with vacuum coupling)" end - compute_delta_prime_matrix!(intr, fm_propagators, fm_chunks; - wv=free_energies.wv, psio=equil.psio, debug=intr.verbose, - S_at_surface_left=fm_S_left, - ctrl=ctrl, equil=equil, mats=mats) + delta_prime = compute_delta_prime_matrix(intr, riccati.propagators, riccati.chunks, riccati.columns.axis_state; + wv=free_energies.wv, wp=free_energies.wp, psio=equil.psio, extended_precision=ctrl.algorithm.extended_precision_bvp) end end end # Publish the solve: from here on the downstream stages read the result, never `intr`. - result = build_result(_integrator_symbol(ctrl.algorithm), ctrl, equil, intr, metric, mats, odet, free_energies, gal_data, gal_dp) + result = build_result(_integrator_symbol(ctrl.algorithm), ctrl, equil, intr, metric, mats, odet, free_energies, gal_data, delta_prime) return result end diff --git a/src/ForceFreeStates/Surfaces/Types.jl b/src/ForceFreeStates/Surfaces/Types.jl index 0099590ab1..333a93ca6c 100644 --- a/src/ForceFreeStates/Surfaces/Types.jl +++ b/src/ForceFreeStates/Surfaces/Types.jl @@ -14,9 +14,7 @@ A mutable struct holding data related to the singular surfaces in the equilibriu - `q::Float64` - Safety factor (= m/n) - `q1::Float64` - Derivative of safety factor with respect to ψ - `delta_prime::Vector{ComplexF64}` - **Debug stub, not physically valid.** Per-surface ca-based Δ' estimate; - use `ForceFreeStatesInternal.delta_prime_matrix` (STRIDE global BVP, Glasser 2018 PoP 25, 032501) for tearing stability. - - `delta_prime_col::Matrix{ComplexF64}` - **Debug stub, not physically valid.** Per-surface ca-based Δ' columns, `(numpert_total, n_res_modes)`, - `delta_prime_col[j, i] = (ca_r[j,ipert_res_i,2] - ca_l[j,ipert_res_i,2]) / (4π²·psio)`; diagonal matches `delta_prime`. Riccati/parallel paths only. + use `ForceFreeStatesResult.delta_prime` (STRIDE global BVP, Glasser & Kolemen 2018, Phys. Plasmas 25, 082502) for tearing stability. """ @kwdef mutable struct SingType psifac::Float64 = 0.0 @@ -26,11 +24,6 @@ A mutable struct holding data related to the singular surfaces in the equilibriu q::Float64 = 0.0 q1::Float64 = 0.0 delta_prime::Vector{ComplexF64} = ComplexF64[] - delta_prime_col::Matrix{ComplexF64} = Matrix{ComplexF64}(undef, 0, 0) - ua_left::Array{ComplexF64,3} = Array{ComplexF64}(undef, 0, 0, 0) # asymptotic basis at left inner-layer boundary - ua_right::Array{ComplexF64,3} = Array{ComplexF64}(undef, 0, 0, 0) # asymptotic basis at right inner-layer boundary - psi_ua_left::Float64 = 0.0 # ψ where ua_left was evaluated (left inner-layer boundary) - psi_ua_right::Float64 = 0.0 # ψ where ua_right was evaluated (right inner-layer boundary) restype::Any = nothing # ResistGeometry from ResistEval.jl (populated by resist_eval_all!); typed `Any` to avoid a cross-file type reference end diff --git a/src/ForceFreeStates/Utils.jl b/src/ForceFreeStates/Utils.jl index b7948380af..f3fb7af636 100644 --- a/src/ForceFreeStates/Utils.jl +++ b/src/ForceFreeStates/Utils.jl @@ -83,9 +83,7 @@ fixup transforms and the free-boundary normalization is exact rather than merely Euler-Lagrange system is linear in `u`, and both operations right-multiply the solution by a mixing matrix `T`, so `du(ψ, u·T) = du(ψ, u)·T`. -Returns `false` without allocating when there is nothing to work from — no `mats`, no stored -steps, or a solution held in a basis the Euler-Lagrange kernel does not apply to (the sparse -parallel path, which stores chunk-endpoint Riccati matrices). +Returns `false` without allocating when there is nothing to work from: no `mats` or no stored steps. """ function materialize_derivative_stores!( odet::OdeState, @@ -94,7 +92,7 @@ function materialize_derivative_stores!( intr::ModeSpace ) odet.du_store_populated && return true - (isnothing(mats) || odet.step == 0 || isempty(odet.u_store) || !odet.u_store_el_basis) && return false + (isnothing(mats) || odet.step == 0 || isempty(odet.u_store)) && return false kinetic = mats.kinetic !== nothing nstep = min(odet.step, size(odet.u_store, 4)) diff --git a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl index 2beaa59163..b8d7fabe64 100644 --- a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl +++ b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl @@ -75,7 +75,10 @@ function plasma_response(ffs::ForceFreeStatesResult; compute_singular_coupling:: resp = PlasmaResponse(; ffs, m_modes, n_modes, compute_singular_coupling, compute_torque_response, odet_from_gal=solution !== nothing && solution.basis === :gal_native) - if ForceFreeStates.require(ffs, :free_boundary, "plasma response calculation") && + # A Riccati result carries ξ profiles, but the response operators are not validated on them yet. + riccati = ffs.integrator === :riccati + riccati && @warn "Skipping plasma response calculation: a Riccati result is not accepted yet; use the Forward integrator" + if !riccati && ForceFreeStates.require(ffs, :free_boundary, "plasma response calculation") && ForceFreeStates.require_solution(ffs, "plasma response calculation") build_response_operators!(resp; verbose) if compute_singular_coupling diff --git a/src/Precompile.jl b/src/Precompile.jl index 79091d88d5..6e00648c2e 100644 --- a/src/Precompile.jl +++ b/src/Precompile.jl @@ -22,7 +22,7 @@ import Logging for formalism in ("Riccati", "Forward") inputs["ForceFreeStates"]["integrator"] = formalism foreach(name -> delete!(inputs, name), ("Riccati", "Forward")) - inputs[formalism] = Dict{String,Any}("renorm_threshold" => 1e3, "eulerlagrange_tolerance" => 1e-7) + inputs[formalism] = Dict{String,Any}("eulerlagrange_tolerance" => 1e-7) mktempdir() do workdir open(io -> TOML.print(io, inputs), joinpath(workdir, "gpec.toml"), "w") Logging.with_logger(Logging.NullLogger()) do diff --git a/src/Tearing/Dispersion/CoupledFullMatch.jl b/src/Tearing/Dispersion/CoupledFullMatch.jl index 19e7e52ba0..578f1bcd93 100644 --- a/src/Tearing/Dispersion/CoupledFullMatch.jl +++ b/src/Tearing/Dispersion/CoupledFullMatch.jl @@ -11,7 +11,7 @@ # Eq. (36)-(40); original matching construction Pletzer & Dewar 1991, # J. Plasma Phys. **45**, 427). # -# The outer D' (`intr.delta_prime_raw`) is in the side-major L/R basis while the inner layer returns +# The outer D' (`DeltaPrimeData.raw`) is in the side-major L/R basis while the inner layer returns # Δ in the parity ± basis, so both amplitude sets are kept as unknowns and joined by the matching # identity rows, giving a 4m × 4m system instead of a basis-mixing 2m × 2m determinant. # @@ -64,9 +64,8 @@ end inner_kwargs=NamedTuple()) -> MultiSurfaceCouplingFull Construct the 4m × 4m dispersion matrix driver. `dp_raw` must be the -2m × 2m matrix in side-major ordering (the `intr.delta_prime_raw` -field populated by `ForceFreeStates.compute_delta_prime_matrix!` on the -Riccati path). Surface k's inner layer is evaluated at `Q·ratio_k` Doppler-shifted by its +2m × 2m matrix in side-major ordering (the `DeltaPrimeData.raw` +field of a Riccati solve). Surface k's inner layer is evaluated at `Q·ratio_k` Doppler-shifted by its E×B offset `q_shift_k`, exactly as in the reduced `multi_surface_coupling`: added to Q for a SLAYER surface, and as `i·q_shift_k` to the growth rate s for a GGJ surface (see `SurfaceCoupling`). diff --git a/test/runtests_analysis_plots.jl b/test/runtests_analysis_plots.jl index 6951228387..8b3bd4bb1c 100644 --- a/test/runtests_analysis_plots.jl +++ b/test/runtests_analysis_plots.jl @@ -140,17 +140,13 @@ using TOML end end - @testset "Riccati run offers Δ′ but no solution profiles" begin + @testset "Riccati run offers Δ′ and solution profiles" begin mktempdir() do dir ret, h5 = _run_fixture(dir, "regression_solovev_ideal_example"; - edit=d -> (d["ForceFreeStates"]["integrator"] = "Riccati"; d["Riccati"] = pop!(d, "Forward"))) + edit=d -> (d["ForceFreeStates"]["integrator"] = "Riccati"; d["Riccati"] = delete!(pop!(d, "Forward"), "renorm_threshold"))) _check_sources(ret.ffs, h5) - @test :delta_prime in A.capabilities(ret.ffs) && :solution_profiles ∉ A.capabilities(ret.ffs) - @test_throws A.REASONS[:solution_profiles] A.ForceFreeStates.plot_mode_displacement(ret.ffs) - @test occursin("✗ plot_mode_displacement", repr(MIME"text/plain"(), available_plots(ret.ffs))) - # The summary leaves out the ξ panel rather than failing, keeping the other four. - caps = A.capabilities(ret.ffs) - @test count(((f, _),) -> A._can(f, caps), A.ForceFreeStates._ffs_panels(ret.ffs)) == 4 + @test :delta_prime in A.capabilities(ret.ffs) && :solution_profiles in A.capabilities(ret.ffs) + @test A.ForceFreeStates.plot_mode_displacement(ret.ffs) isa Figure @test A.ForceFreeStates.plot_ffs_summary(ret.ffs) isa Figure # The Δ′ family draws any part of its complex entries and refuses others. @test A.ForceFreeStates.plot_delta_prime(ret.ffs; part=:abs) isa Figure diff --git a/test/runtests_eulerlagrange.jl b/test/runtests_eulerlagrange.jl index 2fe69f4417..75f0edc93d 100644 --- a/test/runtests_eulerlagrange.jl +++ b/test/runtests_eulerlagrange.jl @@ -135,18 +135,25 @@ end @testset "truncate_chunks! cuts the Riccati chunks back to the edge" begin F = GeneralizedPerturbedEquilibrium.ForceFreeStates - # Crossing chunks end just short of surfaces 1-3; the edge drops surface 3 (msing = 2). + # Backward crossing chunks end just short of surfaces 1-3; the edge drops surface 3 (msing = 2). chunks() = [ - F.IntegrationChunk(; psi_start=a, psi_end=b, needs_crossing=i <= 3, ising=i <= 3 ? i : 0) + F.IntegrationChunk(; psi_start=a, psi_end=b, needs_crossing=i <= 3, ising=i <= 3 ? i : 0, direction=i <= 3 ? -1 : 1) for (i, (a, b)) in enumerate(((0.10, 0.29), (0.31, 0.59), (0.61, 0.89), (0.91, 0.99))) ] cs = chunks() @test F.truncate_chunks!(cs, 0.75, 2) == 3 # straddling chunk shortened, so re-integrate it @test length(cs) == 3 && cs[3].psi_end == 0.75 && count(c -> c.needs_crossing, cs) == 2 + @test cs[3].direction == 1 # no surface left to cross: a forward chunk + # Edge exactly at a crossing chunk's end: its surface is dropped, so it turns forward and is re-integrated. cs = chunks() - @test F.truncate_chunks!(cs, 0.89, 2) === nothing # edge exactly at a chunk end: nothing to re-integrate - @test count(c -> c.needs_crossing, cs) == 2 # but its crossing onto the dropped surface is still cleared + @test F.truncate_chunks!(cs, 0.89, 2) == 3 + @test cs[3].psi_end == 0.89 && !cs[3].needs_crossing && cs[3].direction == 1 + + # Edge at the end of the last chunk with every surface kept: nothing changes. + cs = chunks() + @test F.truncate_chunks!(cs, 0.99, 3) === nothing + @test cs[3].needs_crossing && cs[3].direction == -1 && cs[4].psi_end == 0.99 end @testset "transform_u!" begin @@ -321,9 +328,8 @@ end @test FFS.column_abstol!(tol, u, 1e-6) === tol @test all(FFS.column_abstol(zeros(ComplexF64, 2, 2, 2), 1e-6) .> 0) # Pin the defaults: every shipped deck sets both, so the harness would not catch a revert - c = FFS.ForceFreeStatesControl() - @test c.algorithm.renorm_threshold == 1e3 - @test c.algorithm.save_interval == 1 + @test Forward().renorm_threshold == 1e3 + @test Forward().save_interval == 1 && Riccati().save_interval == 1 end @testset "compute_solution_norms!" begin @@ -449,7 +455,6 @@ end @test size(odet.du_store) == (numpert_total, numpert_total, 0) @test size(odet.xi_s_store) == (numpert_total, numpert_total, 0) @test odet.du_store_populated == false - @test odet.u_store_el_basis == true @test length(odet.psi_store) == numsteps_init @test length(odet.q_store) == numsteps_init @test size(odet.ca_r) == (numpert_total, numpert_total, 2, msing) @@ -666,7 +671,7 @@ end intr.numpert_total = intr.mpert * intr.npert metric = FFS.make_metric(equil, intr.mpert) mats = FFS.build_matrix_splines(equil, intr, metric) - odet, _, _, _ = FFS.eulerlagrange_integration(ctrl, equil, mats, intr) + odet = first(FFS.eulerlagrange_integration(ctrl, equil, mats, intr)) return odet, ctrl, equil, mats, intr end @@ -722,15 +727,6 @@ end @test isapprox(odet_t.xi_s_store[:, :, istep], odet.xi_s_store[:, :, istep] * T; rtol=1e-10) end end - - @testset "refuses a solution outside the Euler-Lagrange basis" begin - odet.du_store_populated = false - odet.du_store = Array{ComplexF64}(undef, intr.numpert_total, intr.numpert_total, 0) - odet.u_store_el_basis = false - @test !FFS.materialize_derivative_stores!(odet, equil, mats, intr) - @test isempty(odet.du_store) - @test !odet.du_store_populated - end end @testset "resonant block Frobenius exponents" begin diff --git a/test/runtests_parallel_integration.jl b/test/runtests_parallel_integration.jl index 7ddba18d9f..9a11dbcd1b 100644 --- a/test/runtests_parallel_integration.jl +++ b/test/runtests_parallel_integration.jl @@ -9,6 +9,7 @@ function _prepared_solovev(integrator::String="Forward"; nchunks=nothing) if integrator != "Forward" inputs["ForceFreeStates"]["integrator"] = integrator integrator == "Galerkin" ? delete!(inputs, "Forward") : (inputs[integrator] = pop!(inputs, "Forward")) + integrator == "Riccati" && delete!(inputs["Riccati"], "renorm_threshold") end nchunks === nothing || (inputs[integrator]["nchunks"] = nchunks) equil = GPEC.Equilibrium.PlasmaEquilibrium(GPEC.equilibrium_config(inputs, ex), GPEC.analytic_config(inputs)) @@ -17,7 +18,7 @@ function _prepared_solovev(integrator::String="Forward"; nchunks=nothing) return ctrl, intr, equil, metric, mats end -@testset "BVP without controls preserves extended precision" begin +@testset "BVP solve precision" begin FFS = GeneralizedPerturbedEquilibrium.ForceFreeStates M = ComplexF64[1 / (i + j - 1) + 0.02im / (i + j) for i in 1:6, j in 1:6] reference = setprecision(BigFloat, 256) do @@ -32,123 +33,71 @@ end delta = reshape([raw[2, 2] - raw[2, 1] - raw[1, 2] + raw[1, 1]], 1, 1) return ComplexF64.(delta), ComplexF64.(raw) end - default = FFS._solve_bvp_and_combine_pest3(M, 1, 1, 6, false, [1], nothing, nothing, false) - extended = FFS._solve_bvp_and_combine_pest3(M, 1, 1, 6, false, [1], nothing, FFS.ForceFreeStatesControl(), false) - ordinary = FFS._solve_bvp_and_combine_pest3(M, 1, 1, 6, false, [1], nothing, FFS.ForceFreeStatesControl(; algorithm=Riccati(; extended_precision_bvp=false)), false) - @test default == reference + extended = FFS._solve_bvp_and_combine_pest3(M, 1, 1, [1], true) + ordinary = FFS._solve_bvp_and_combine_pest3(M, 1, 1, [1], false) @test extended == reference @test ordinary != reference + @test all(isapprox.(ordinary, reference; rtol=1e-8)) + @test !Riccati().extended_precision_bvp end @testset "Riccati FM Integration Tests" begin - @testset "ChunkPropagator identity on trivial interval" begin - # Integrating over a zero-width interval should give the identity propagator. - # We test that apply_propagator! on an identity state preserves the state. + @testset "apply_propagator! and its inverse" begin + FFS = GeneralizedPerturbedEquilibrium.ForceFreeStates N = 3 - prop = GeneralizedPerturbedEquilibrium.ForceFreeStates.ChunkPropagator(N) - - # Set propagator to identity (block_upper_ic = (I, 0), block_lower_ic = (0, I)) - for i in 1:N - prop.block_upper_ic[i, i, 1] = 1 # U1 block from IC=(I,0) - prop.block_lower_ic[i, i, 2] = 1 # U2 block from IC=(0,I) - end - - # Apply identity propagator to an arbitrary state - odet = GeneralizedPerturbedEquilibrium.ForceFreeStates.OdeState(N, 10, 5, 0) - u1_in = [1.0+0.5im 0.2im 0.0; - 0.1+0.1im 1.2+0.1im 0.0; - 0.0im 0.0 0.9+0.3im] - u2_in = [0.8+0.1im 0.1im 0.0; - 0.0im 1.0+0.2im 0.1; - 0.1im 0.0 1.1+0.0im] - odet.u[:, :, 1] .= u1_in - odet.u[:, :, 2] .= u2_in - - GeneralizedPerturbedEquilibrium.ForceFreeStates.apply_propagator!(odet, prop) - - @test odet.u[:, :, 1] ≈ u1_in rtol = 1e-12 - @test odet.u[:, :, 2] ≈ u2_in rtol = 1e-12 - end - - @testset "apply_propagator! linearity" begin - # Verify that apply_propagator! applies the correct linear map. - N = 3 - prop = GeneralizedPerturbedEquilibrium.ForceFreeStates.ChunkPropagator(N) - - # Fill block_upper_ic and block_lower_ic with random data - rng_upper = [1.1+0.2im 0.1im 0.05; - 0.0im 0.9+0.3im 0.1; - 0.2+0.1im 0.0 1.0+0.1im] - rng_lower = [0.8+0.1im 0.1im 0.0; - 0.0im 1.2+0.2im 0.1; - 0.0im 0.1 0.9+0.1im] - prop.block_upper_ic[:, :, 1] .= rng_upper - prop.block_upper_ic[:, :, 2] .= 0.5 * rng_upper - prop.block_lower_ic[:, :, 1] .= 0.3 * rng_lower - prop.block_lower_ic[:, :, 2] .= rng_lower - - odet = GeneralizedPerturbedEquilibrium.ForceFreeStates.OdeState(N, 10, 5, 0) - u1_in = 0.5 * I(N) .+ 0.1im * ones(N, N) - u2_in = I(N) .+ 0.2im * ones(N, N) - odet.u[:, :, 1] .= u1_in - odet.u[:, :, 2] .= u2_in - - GeneralizedPerturbedEquilibrium.ForceFreeStates.apply_propagator!(odet, prop) - - # Manual computation of expected result - U1_upper = prop.block_upper_ic[:, :, 1] - U2_upper = prop.block_upper_ic[:, :, 2] - U1_lower = prop.block_lower_ic[:, :, 1] - U2_lower = prop.block_lower_ic[:, :, 2] - u1_expected = U1_upper * u1_in + U1_lower * u2_in - u2_expected = U2_upper * u1_in + U2_lower * u2_in - - @test odet.u[:, :, 1] ≈ u1_expected rtol = 1e-12 - @test odet.u[:, :, 2] ≈ u2_expected rtol = 1e-12 - end - - @testset "apply_propagator_inverse! is inverse of apply_propagator!" begin - # Verify that apply_propagator_inverse! is the algebraic inverse of apply_propagator!: - # applying inverse then forward should recover the original state exactly. - # This checks the LU-solve path: Φ \ (Φ * u) = u for an arbitrary invertible Φ. - N = 3 - prop = GeneralizedPerturbedEquilibrium.ForceFreeStates.ChunkPropagator(N) - - # Near-identity blocks guarantee the 2N×2N matrix [A B; C D] is invertible A = I(N) .+ 0.15 * [1.0+0.2im 0.1im 0.05; 0.0im 0.9+0.3im 0.1; 0.2+0.1im 0.0 1.0+0.1im] B = 0.1 * [0.8+0.1im 0.1im 0.0; 0.0im 1.2+0.2im 0.1; 0.0im 0.1 0.9+0.1im] C = 0.1 * [0.5+0.1im 0.0im 0.1; 0.1im 0.8+0.2im 0.0; 0.0im 0.0 0.7+0.1im] D = I(N) .+ 0.15 * [0.9+0.1im 0.0im 0.05; 0.0im 1.0+0.2im 0.0; 0.1+0.1im 0.0 0.95+0.1im] + # Fundamental matrix [A B; C D] in the propagator layout: U₁ rows in slice 1, U₂ rows in slice 2. + prop = FFS.ChunkPropagator(cat([A B], [C D]; dims=3), Float64[], Array{ComplexF64,3}[]) - prop.block_upper_ic[:, :, 1] .= A - prop.block_lower_ic[:, :, 1] .= B - prop.block_upper_ic[:, :, 2] .= C - prop.block_lower_ic[:, :, 2] .= D - - u1_in = [1.0+0.5im 0.2im 0.0; - 0.1+0.1im 1.2+0.1im 0.0; - 0.0im 0.0 0.9+0.3im] + u1_in = [1.0+0.5im 0.2im 0.0; 0.1+0.1im 1.2+0.1im 0.0; 0.0im 0.0 0.9+0.3im] u2_in = I(N) .+ 0.1im * ones(N, N) - - odet = GeneralizedPerturbedEquilibrium.ForceFreeStates.OdeState(N, 10, 5, 0) + odet = FFS.OdeState(N, 10, 5, 0) odet.u[:, :, 1] .= u1_in odet.u[:, :, 2] .= u2_in - # Round-trip: inverse then forward = identity - GeneralizedPerturbedEquilibrium.ForceFreeStates.apply_propagator_inverse!(odet, prop) - GeneralizedPerturbedEquilibrium.ForceFreeStates.apply_propagator!(odet, prop) + FFS.apply_propagator!(odet, prop) + @test odet.u[:, :, 1] ≈ A * u1_in + B * u2_in rtol = 1e-12 + @test odet.u[:, :, 2] ≈ C * u1_in + D * u2_in rtol = 1e-12 + FFS.apply_propagator_inverse!(odet, prop) @test odet.u[:, :, 1] ≈ u1_in rtol = 1e-12 @test odet.u[:, :, 2] ≈ u2_in rtol = 1e-12 end + @testset "propagator steps reproduce a direct integration" begin + # The fundamental matrix at a solver step, applied to a start state, is the solution from that state. + FFS = GeneralizedPerturbedEquilibrium.ForceFreeStates + ctrl, intr, equil, metric, mats = _prepared_solovev("Riccati") + N = intr.numpert_total + odet = FFS.OdeState(N, 1, 1, intr.msing) + FFS.initialize_el_at_axis!(odet, ctrl, mats, equil.profiles, intr) + chunks = FFS.balance_integration_chunks(FFS.chunk_el_integration_bounds(odet, ctrl, intr), ctrl, intr) + chunk = chunks[1] + start = FFS.propagator_start(chunks, 1, Dict{Int,Any}(), N) + prop = FFS.chunk_propagator(chunk, start, ctrl.algorithm.eulerlagrange_tolerance, equil, mats, FFS.ModeIndices(intr)) + @test length(prop.psi) == length(prop.steps) + @test prop.psi[1] == chunk.psi_start && prop.psi[end] == chunk.psi_end + @test prop.steps[1][:, :, 1] == [I(N) zeros(N, N)] && prop.steps[1][:, :, 2] == [zeros(N, N) I(N)] + @test prop.steps[end] == prop.u + + k = length(prop.psi) ÷ 2 + direct = FFS.solve_el(FFS.ODEProblem(FFS.sing_der!, copy(odet.u), (chunk.psi_start, prop.psi[k]), + (ctrl, equil, mats, intr, FFS.OdeState(N, 1, 1, 0), chunk)), ctrl.algorithm.eulerlagrange_tolerance).u[end] + x = [odet.u[:, :, 1]; odet.u[:, :, 2]] + @test prop.steps[k][:, :, 1] * x ≈ direct[:, :, 1] rtol = 1e-5 + @test prop.steps[k][:, :, 2] * x ≈ direct[:, :, 2] rtol = 1e-5 + end + @testset "balance_integration_chunks produces target count" begin - # Verify that balance_integration_chunks creates at least - # max(2*msing+3, 8*(msing+1)+msing) chunks from a small set of base chunks. + # Verify that balance_integration_chunks creates 8*(msing+1)+msing chunks from a small set of + # base chunks, and never fewer than the base chunks with every crossing chunk split once. ctrl, intr, equil, metric, mats = _prepared_solovev("Riccati") - odet = GeneralizedPerturbedEquilibrium.ForceFreeStates.OdeState(intr.numpert_total, ctrl.algorithm.numsteps_init, ctrl.algorithm.numunorms_init, intr.msing) + odet = GeneralizedPerturbedEquilibrium.ForceFreeStates.OdeState(intr.numpert_total, 1, 1, intr.msing) GeneralizedPerturbedEquilibrium.ForceFreeStates.initialize_el_at_axis!(odet, ctrl, mats, equil.profiles, intr) base_chunks = GeneralizedPerturbedEquilibrium.ForceFreeStates.chunk_el_integration_bounds(odet, ctrl, intr) @@ -158,7 +107,8 @@ end # (src/ForceFreeStates/EulerLagrange.jl). Keep this in sync. The formula reads only # intr.msing — no thread count enters it, which is what makes Riccati outputs # independent of how many threads `julia -t` provides. - target_n = max(2 * intr.msing + 3, 8 * (intr.msing + 1) + intr.msing) + min_chunks = length(base_chunks) + count(c -> c.needs_crossing, base_chunks) + target_n = max(min_chunks, 8 * (intr.msing + 1) + intr.msing) # After balancing, chunk count equals target_n: the while-loop adds exactly one # chunk per iteration (a bisection split) and exits when length(result) >= target_n, @@ -201,11 +151,12 @@ end balanced_more = GeneralizedPerturbedEquilibrium.ForceFreeStates.balance_integration_chunks(base_chunks, ctrl_more, intr) @test length(balanced_more) == target_n + 7 - # An nchunks below the singular-surface floor is clamped up, with a warning. - min_chunks = 2 * intr.msing + 3 + # An nchunks below the singular-surface floor is clamped up, with a warning. At the floor a + # forward chunk still separates every crossing chunk from the previous surface or the axis. ctrl_few = first(_prepared_solovev("Riccati"; nchunks=1)) balanced_few = @test_logs (:warn,) match_mode = :any GeneralizedPerturbedEquilibrium.ForceFreeStates.balance_integration_chunks(base_chunks, ctrl_few, intr) - @test length(balanced_few) == max(min_chunks, length(base_chunks)) + @test length(balanced_few) == min_chunks + @test all(i > 1 && !balanced_few[i-1].needs_crossing for i in eachindex(balanced_few) if balanced_few[i].needs_crossing) end @testset "chunk_el_integration_bounds direction field — bidirectional mode" begin @@ -214,7 +165,7 @@ end # the right sub-chunk inherits direction from the parent, the left sub-chunk is always +1. ctrl, intr, equil, metric, mats = _prepared_solovev("Riccati") - odet = GeneralizedPerturbedEquilibrium.ForceFreeStates.OdeState(intr.numpert_total, ctrl.algorithm.numsteps_init, ctrl.algorithm.numunorms_init, intr.msing) + odet = GeneralizedPerturbedEquilibrium.ForceFreeStates.OdeState(intr.numpert_total, 1, 1, intr.msing) GeneralizedPerturbedEquilibrium.ForceFreeStates.initialize_el_at_axis!(odet, ctrl, mats, equil.profiles, intr) # Default (bidirectional=false): all chunks should have direction=+1 @@ -253,7 +204,7 @@ end # on this small multi-surface fixture. function run_solovev(integrator) ctrl, intr, equil, metric, mats = _prepared_solovev(integrator) - odet, _, _, _ = GeneralizedPerturbedEquilibrium.ForceFreeStates.eulerlagrange_integration(ctrl, equil, mats, intr) + odet = first(GeneralizedPerturbedEquilibrium.ForceFreeStates.eulerlagrange_integration(ctrl, equil, mats, intr)) vac = GeneralizedPerturbedEquilibrium.ForceFreeStates.free_run(odet, ctrl, equil, mats, intr) return real(vac.et[1]), intr end @@ -298,35 +249,16 @@ end # physically meaningful. BVP Δ' regression is concentrated on the DIIID-like # fixture below (intrinsically stable, well-conditioned BVP Δ'). - @testset "Riccati leaves a sparse u_store in the Riccati basis" begin - # The Riccati path stores only chunk-endpoint snapshots, so u_store / psi_store have - # strictly fewer entries than the forward path's dense saved steps, and the stored - # state is not in the Euler-Lagrange axis basis. Downstream ξ consumers must therefore - # be fed by the forward integrator; this test pins that contract. - function run_and_capture(integrator) + @testset "both ODE integrators leave solution profiles" begin + # Either store materializes into dΞ_ψ/dψ and Ξ_s, which is what PerturbedEquilibrium consumes. + for integrator in ("Forward", "Riccati") ctrl, intr, equil, metric, mats = _prepared_solovev(integrator) - odet, _, _, _ = GeneralizedPerturbedEquilibrium.ForceFreeStates.eulerlagrange_integration(ctrl, equil, mats, intr) - # Derivatives are recomputed on demand; materialize so the stores can be compared. - GeneralizedPerturbedEquilibrium.ForceFreeStates.materialize_derivative_stores!(odet, equil, mats, intr) - return odet + odet = first(GeneralizedPerturbedEquilibrium.ForceFreeStates.eulerlagrange_integration(ctrl, equil, mats, intr)) + @test GeneralizedPerturbedEquilibrium.ForceFreeStates.materialize_derivative_stores!(odet, equil, mats, intr) + @test odet.du_store_populated + @test size(odet.du_store) == size(odet.xi_s_store) == (intr.numpert_total, intr.numpert_total, odet.step) + @test all(isfinite, odet.du_store) && all(isfinite, odet.xi_s_store) end - - odet_fwd = run_and_capture("Forward") - odet_ric = run_and_capture("Riccati") - - @test odet_ric.step < odet_fwd.step - @test length(odet_ric.psi_store) < length(odet_fwd.psi_store) - - # The forward path returns dense ξ in the axis basis, ready for PerturbedEquilibrium. - @test odet_fwd.u_store_el_basis - @test odet_fwd.du_store_populated - - # The Riccati solution is in the Riccati basis, so the derivative stores cannot be - # materialized from it and stay empty rather than holding unusable values. - @test !odet_ric.u_store_el_basis - @test !odet_ric.du_store_populated - @test isempty(odet_ric.du_store) - @test isempty(odet_ric.xi_s_store) end # Production DIII-D energies and Δ′ are tracked by the diiid_n1_riccati harness cases. diff --git a/test/runtests_resist_eval.jl b/test/runtests_resist_eval.jl index f82b5f6595..2f703a53f2 100644 --- a/test/runtests_resist_eval.jl +++ b/test/runtests_resist_eval.jl @@ -74,18 +74,10 @@ dq = deriv_view(equil.profiles.q_spline, 1) s1 = SingType(; psifac=0.3, rho=sqrt(0.3), m=[2], n=[1], q=2.0, q1=dq(0.3), - delta_prime=ComplexF64[], - delta_prime_col=zeros(ComplexF64, 0, 0), - ua_left=zeros(ComplexF64, 0, 0, 0), - ua_right=zeros(ComplexF64, 0, 0, 0), - psi_ua_left=0.0, psi_ua_right=0.0) + delta_prime=ComplexF64[]) s2 = SingType(; psifac=0.7, rho=sqrt(0.7), m=[3], n=[1], q=3.0, q1=dq(0.7), - delta_prime=ComplexF64[], - delta_prime_col=zeros(ComplexF64, 0, 0), - ua_left=zeros(ComplexF64, 0, 0, 0), - ua_right=zeros(ComplexF64, 0, 0, 0), - psi_ua_left=0.0, psi_ua_right=0.0) + delta_prime=ComplexF64[]) @test s1.restype === nothing @test s2.restype === nothing @@ -105,11 +97,7 @@ dq = deriv_view(equil.profiles.q_spline, 1) s1 = SingType(; psifac=0.3, rho=sqrt(0.3), m=[2], n=[1], q=2.0, q1=dq(0.3), - delta_prime=ComplexF64[], - delta_prime_col=zeros(ComplexF64, 0, 0), - ua_left=zeros(ComplexF64, 0, 0, 0), - ua_right=zeros(ComplexF64, 0, 0, 0), - psi_ua_left=0.0, psi_ua_right=0.0) + delta_prime=ComplexF64[]) intr = ForceFreeStates.ForceFreeStatesInternal(; sing=[s1], msing=1) ForceFreeStates.resist_eval_all!(intr, equil) @@ -138,11 +126,7 @@ @testset "ggj_parameters: errors when restype not populated" begin s_unpop = SingType(; psifac=0.5, rho=sqrt(0.5), m=[2], n=[1], q=2.0, q1=1.0, - delta_prime=ComplexF64[], - delta_prime_col=zeros(ComplexF64, 0, 0), - ua_left=zeros(ComplexF64, 0, 0, 0), - ua_right=zeros(ComplexF64, 0, 0, 0), - psi_ua_left=0.0, psi_ua_right=0.0) + delta_prime=ComplexF64[]) @test s_unpop.restype === nothing @test_throws ArgumentError ggj_parameters(s_unpop, equil; eta=1e-7, rho=1e-7) end @@ -151,11 +135,7 @@ dq = deriv_view(equil.profiles.q_spline, 1) s1 = SingType(; psifac=0.3, rho=sqrt(0.3), m=[2], n=[1], q=2.0, q1=dq(0.3), - delta_prime=ComplexF64[], - delta_prime_col=zeros(ComplexF64, 0, 0), - ua_left=zeros(ComplexF64, 0, 0, 0), - ua_right=zeros(ComplexF64, 0, 0, 0), - psi_ua_left=0.0, psi_ua_right=0.0) + delta_prime=ComplexF64[]) intr = ForceFreeStates.ForceFreeStatesInternal(; sing=[s1], msing=1) ForceFreeStates.resist_eval_all!(intr, equil) gs = [ggj_parameters(intr.sing[1], equil; eta=eta, rho=rho, ising=1)] diff --git a/test/runtests_result_struct.jl b/test/runtests_result_struct.jl index f38ada7aa3..4cf9ffbfc7 100644 --- a/test/runtests_result_struct.jl +++ b/test/runtests_result_struct.jl @@ -22,6 +22,7 @@ include("display_check.jl") merge!(inputs["ForceFreeStates"], ffs_overrides) new = get(inputs["ForceFreeStates"], "integrator", "Riccati") old == new || haskey(inputs, old) && (new == "Galerkin" ? delete!(inputs, old) : (inputs[new] = pop!(inputs, old))) + haskey(inputs, "Riccati") && delete!(inputs["Riccati"], "renorm_threshold") for (name, section) in extra_sections haskey(inputs, name) ? merge!(inputs[name], section) : (inputs[name] = section) end @@ -94,7 +95,7 @@ include("display_check.jl") end end - @testset "riccati run publishes Δ′ and diagnostics but no ξ solution" begin + @testset "riccati run publishes Δ′, diagnostics and a ξ solution" begin mktempdir() do dir _stage_deck(dir, template; ffs_overrides=Dict{String,Any}("integrator" => "Riccati")) ret = GeneralizedPerturbedEquilibrium.main([dir]) @@ -102,10 +103,9 @@ include("display_check.jl") ffs = ret.ffs riccati_result[] = ffs @test ffs.integrator === :riccati - # Chunk-endpoint Riccati states are not a ξ solution; only the raw state is carried. - @test ffs.solution === nothing + @test ffs.solution !== nothing && ffs.solution.basis === :el_axis + @test ffs.solution.step == length(ffs.solution.psi_store) == size(ffs.solution.du_store, 3) @test ffs.diagnostics !== nothing - @test !ffs.diagnostics.u_store_el_basis @test ffs.closure === :ideal @test iszero(ffs.bpen) @@ -130,7 +130,7 @@ include("display_check.jl") @test !(@test_logs (:warn,) FFS.require(ffs, :galerkin, "a calculation")) # `require_solution` is the ξ-specific presence gate. - @test !(@test_logs (:warn,) FFS.require_solution(ffs, "a calculation")) + @test @test_logs FFS.require_solution(ffs, "a calculation") end end @@ -153,11 +153,11 @@ include("display_check.jl") @test ret.pe !== nothing @test isempty(ret.pe.response.permeability) @test isempty(ret.pe.response.C_delta_prime) - @test ret.ffs.free_boundary !== nothing && ret.ffs.solution === nothing + @test ret.ffs.free_boundary !== nothing && ret.ffs.solution !== nothing @test ret.resp.compute_singular_coupling for result in (ret.resp, ret.pe) plain = _check_display(result) - @test occursin("needs stored ξ profiles", plain) + @test occursin("a Riccati result is not accepted yet", plain) @test !occursin("needs a free-boundary result", plain) @test !occursin("singular coupling off", plain) @test !occursin(" 0 J", plain) diff --git a/test/runtests_riccati.jl b/test/runtests_riccati.jl index b3904b9702..ec60ed0196 100644 --- a/test/runtests_riccati.jl +++ b/test/runtests_riccati.jl @@ -1,22 +1,34 @@ -using LinearAlgebra, Random, TOML +using LinearAlgebra, Logging, TOML const FFS = GeneralizedPerturbedEquilibrium.ForceFreeStates -# The Solovev fixture prepared as solve prepares it, up to the integration. Each integration mutates its +# A fixture (Solovev by default) prepared as solve prepares it, up to the integration. Each integration mutates its # intr (sing[s].delta_prime etc.), so every run takes a freshly prepared one. -function prepared_solovev() +function prepared_solovev(; fixture="regression_solovev_ideal_example", integrator="Riccati", singfac_min=nothing, nchunks=0, over=Dict{String,Any}()) GPEC = GeneralizedPerturbedEquilibrium - ex = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example") + ex = joinpath(@__DIR__, "test_data", fixture) inputs = TOML.parsefile(joinpath(ex, "gpec.toml")) - # The fixture's [Forward] settings drive both paths; the Riccati control also carries the chunking. - inputs["ForceFreeStates"]["integrator"] = "Riccati" - inputs["Riccati"] = pop!(inputs, "Forward") + # The fixture's [Forward] tolerance drives both integrators; the Riccati control also carries the chunking. + inputs["ForceFreeStates"]["integrator"] = integrator + if integrator == "Riccati" + inputs["Riccati"] = delete!(pop!(inputs, "Forward"), "renorm_threshold") + inputs["Riccati"]["nchunks"] = nchunks + end + singfac_min === nothing || (inputs["ForceFreeStates"]["singfac_min"] = singfac_min) + merge!(inputs["ForceFreeStates"], over) equil = GPEC.Equilibrium.PlasmaEquilibrium(GPEC.equilibrium_config(inputs, ex), GPEC.analytic_config(inputs)) prob = GPEC.ForceBalanceProblem(equil; GPEC.force_balance_keywords(inputs, ex)...) ctrl, intr, equil, metric, mats = GPEC._prepare_solve(prob, GPEC.force_balance_algorithm(inputs)) return ctrl, intr, equil, metric, mats end +# Free-boundary energy eigenvalues of a Solovev fixture, crossing each surface at `singfac_min`. +function solovev_energies(integrator, singfac_min; fixture="regression_solovev_ideal_example") + ctrl, intr, equil, _, mats = prepared_solovev(; fixture, integrator, singfac_min) + odet = first(FFS.eulerlagrange_integration(ctrl, equil, mats, intr)) + return FFS.free_run(odet, ctrl, equil, mats, intr).et +end + @testset "Riccati Integration Tests" begin # ── Pure matrix unit tests — no equilibrium needed ──────────────────────── @@ -59,69 +71,68 @@ end @test u[:, :, 1] ≈ S rtol = 1e-12 end - @testset "renormalize_riccati! (OdeState)" begin - N = 3 - rng = [1.0+0.5im 0.2im 0.1; - 0.0im 1.2+0.1im 0.0; - 0.1+0.1im 0.0 0.9+0.3im] - U1 = rng .+ 0.5 * I(N) - U2 = 0.2 * rng .+ I(N) - - odet = FFS.OdeState(N, 10, 5, 1) - odet.u[:, :, 1] .= U1 - odet.u[:, :, 2] .= U2 - - S_expected = U1 / U2 - intr = FFS.ForceFreeStatesInternal(; mpert=N, numpert_total=N) - - FFS.renormalize_riccati!(odet, intr) - - @test odet.u[:, :, 2] ≈ I(N) - @test odet.u[:, :, 1] ≈ S_expected rtol = 1e-12 - end - # ── Shared Solovev setup ────────────────────────────────────────────────── # # equil (Grad-Shafranov solve) and mats (metric matrices) are expensive and # immutable after construction — built ONCE and shared across all tests below. # intr is cheap to (re)initialize but is mutated by each integration run # (sing[s].delta_prime etc.), so a fresh copy is made for each integration. - # - # Integration runs: - # intr_ric / odet_ric — Riccati path (shared by most tests) - # intr_fwd / odet_fwd — Forward path (energy comparison only) ctrl, intr_ric, equil, metric, mats = prepared_solovev() N = intr_ric.numpert_total - # Riccati integration. The driver returns (odet, propagators, chunks, S_at_surface_left); - # only odet is used here. - odet_ric, _, _, _ = FFS.riccati_eulerlagrange_integration(ctrl, equil, mats, intr_ric) + odet_ric = first(FFS.riccati_eulerlagrange_integration(ctrl, equil, mats, intr_ric)) - # Save inline Δ' values before any test that calls compute_delta_prime_from_ca! - # (which overwrites intr_ric.sing[s].delta_prime) - delta_prime_inline = [copy(intr_ric.sing[s].delta_prime) for s in 1:intr_ric.msing] - - vac_ric = FFS.free_run(odet_ric, ctrl, equil, mats, intr_ric) - et_ric = real(vac_ric.et[1]) + # ───────────────────────────────────────────────────────────────────────── - # Forward integration (needed only for energy comparison). - _, intr_fwd, _, _, _ = prepared_solovev() - odet_fwd, _, _, _ = FFS.forward_eulerlagrange_integration(ctrl, equil, mats, intr_fwd) - vac_fwd = FFS.free_run(odet_fwd, ctrl, equil, mats, intr_fwd) - et_fwd = real(vac_fwd.et[1]) + @testset "energy does not depend on the crossing distance" begin + et_far = solovev_energies("Riccati", 1e-3) + et_near = solovev_energies("Riccati", 1e-5) + @test isapprox(et_far[1], et_near[1]; rtol=1e-6) + # Forward's trapezoid crossing reaches the same energy only as the distance shrinks. + et_forward = solovev_energies("Forward", 1e-6) + @test isapprox(real(et_far[1]), real(et_forward[1]); rtol=1e-3) + end - # ───────────────────────────────────────────────────────────────────────── + @testset "surface with two resonances" begin + fixture = "regression_solovev_ideal_example_multi_n" + et_far = solovev_energies("Riccati", 1e-3; fixture) + et_near = solovev_energies("Riccati", 1e-5; fixture) + et_forward = solovev_energies("Forward", 1e-6; fixture) + # The least-stable multi-n eigenvalue is near marginal, so its relative tolerances are looser. + @test isapprox(et_far[1], et_near[1]; rtol=1e-4) + @test isapprox(real(et_far[1]), real(et_forward[1]); rtol=5e-3) + end - @testset "Riccati integration matches forward ODE — Solovev example" begin - # The two formalisms solve the same system, so the leading energy eigenvalue must - # agree closely; the measured Solovev disagreement is ≈5e-5. A regression of the - # Riccati/renormalization algorithm to ~1 % error would fail here loudly. - @test isapprox(et_ric, et_fwd; rtol=1e-4) + @testset "stored solution" begin + # Every solver step of every chunk is stored, in increasing ψ, ending at the edge state. + @test odet_ric.step == length(odet_ric.psi_store) == size(odet_ric.u_store, 4) + @test all(diff(odet_ric.psi_store) .> 0) + @test odet_ric.u_store[:, :, :, end] ≈ odet_ric.u + + # Each stored column is a solution: integrated directly from one stored node it lands on the next. + psi = odet_ric.psi_store + for k in (1, odet_ric.step ÷ 4, odet_ric.step ÷ 2, odet_ric.step - 1) + any(s -> psi[k] < s.psifac < psi[k+1], intr_ric.sing) && continue + params = (ctrl, equil, mats, intr_ric, FFS.OdeState(N, 1, 1, 0), FFS.IntegrationChunk(psi[k], psi[k+1], false, 0, 1)) + direct = FFS.solve_el(FFS.ODEProblem(FFS.sing_der!, odet_ric.u_store[:, :, :, k], (psi[k], psi[k+1]), params), + ctrl.algorithm.eulerlagrange_tolerance).u[end] + @test isapprox(direct, odet_ric.u_store[:, :, :, k+1]; rtol=1e-5) + end - # Riccati stores chunk endpoints, the forward path stores every saved ODE step, - # so the Riccati store is always the sparser of the two. - @test odet_ric.step <= odet_fwd.step + # Across the surfaces they are the solutions Forward stores: scaled to a unit edge displacement, the + # conjugate momenta agree on each side of every surface. (ξ itself differs there by Forward's crossing + # condition, which leaves a big-solution term next to each surface.) + function edge_scaled_momentum(integrator) + ctrl, intr, equil, _, mats = prepared_solovev(; integrator, singfac_min=1e-6) + odet = first(FFS.eulerlagrange_integration(ctrl, equil, mats, intr)) + return odet, intr, [odet.u_store[:, :, 2, k] / odet.u_store[:, :, 1, odet.step] for k in 1:odet.step] + end + ric, intr, momentum_ric = edge_scaled_momentum("Riccati") + fwd, _, momentum_fwd = edge_scaled_momentum("Forward") + shared = [(findfirst(==(p), ric.psi_store), k) for (k, p) in enumerate(fwd.psi_store[1:fwd.step]) if p in ric.psi_store] + @test length(shared) >= 2 * intr.msing + @test maximum(norm(momentum_ric[i] - momentum_fwd[k]) / norm(momentum_fwd[k]) for (i, k) in shared) < 1e-2 end # Note: a Solovev per-surface Δ' regression testset previously lived here, @@ -132,79 +143,44 @@ end # STRIDE BVP Δ' matrix (see runtests_parallel_integration.jl). @testset "Riccati end state has U₂ ≈ I" begin - # The outer-region re-integration that closes riccati_eulerlagrange_integration - # leaves odet.u[:,:,2] as the identity (canonical Riccati convention after final renorm) + # The edge state is reduced to the canonical Riccati form (S, I). @test odet_ric.u[:, :, 2] ≈ I(N) rtol = 1e-10 end - @testset "riccati_der! formula — Glasser 2018 Eq. 19" begin - # Verify riccati_der! correctly evaluates dS/dψ = w†·F̄⁻¹·w − S·Ḡ·S, w = Q − K̄·S. - # - # Test states are Hermitian (physical constraint: the EL system preserves S†=S from - # the axis). Non-Hermitian states would give ~5% disagreement — not a bug, but a - # consequence of the derivation assuming the physical symmetry. - # - # See benchmarks/benchmark_riccati_der.jl for the extended version with output. - - # Use an initialized OdeState just for spline_hint and chunk bounds - odet_tmp = FFS.OdeState(N, ctrl.algorithm.numsteps_init, ctrl.algorithm.numunorms_init, intr_ric.msing) - FFS.initialize_el_at_axis!(odet_tmp, ctrl, mats, equil.profiles, intr_ric) - chunks = FFS.chunk_el_integration_bounds(odet_tmp, ctrl, intr_ric) - - # 30% into each chunk: away from singularities at psi_end - test_psis = [c.psi_start + 0.3 * (c.psi_end - c.psi_start) for c in chunks] - - rng = Random.MersenneTwister(42) - for psi in test_psis - # Hermitian S: physical Riccati matrix is Hermitian (preserved by EL symmetry) - A = randn(rng, ComplexF64, N, N) - S = (A + A') / 2 - - # Manual RHS: w†·F̄⁻¹·w − S·Ḡ·S - L = zeros(ComplexF64, N, N) - Kmat = zeros(ComplexF64, N, N) - Gmat = zeros(ComplexF64, N, N) - mats.ideal.F_spline_lower(vec(L), psi; hint=mats._hint) - mats.ideal.K_spline(vec(Kmat), psi; hint=mats._hint) - mats.ideal.G_spline(vec(Gmat), psi; hint=mats._hint) - q = equil.profiles.q_spline(psi) - singfac = vec(1.0 ./ ((intr_ric.mlow:intr_ric.mhigh) .- - q .* (intr_ric.nlow:intr_ric.nhigh)')) - w = -Kmat * S - for i in 1:N - w[i, i] += singfac[i] - end - v = copy(w) - ldiv!(LowerTriangular(L), v) - ldiv!(UpperTriangular(L'), v) - dS_manual = adjoint(w) * v - S * Gmat * S - - # riccati_der! RHS - u_ric = zeros(ComplexF64, N, N, 2) - du_ric = zeros(ComplexF64, N, N, 2) - u_ric[:, :, 1] .= S - u_ric[:, :, 2] .= Matrix{ComplexF64}(I, N, N) - dummy = FFS.IntegrationChunk(psi, psi, false, 0, 1) - params = (ctrl, equil, mats, intr_ric, odet_tmp, dummy) - FFS.riccati_der!(du_ric, u_ric, params, psi) - - rel_err = norm(du_ric[:, :, 1] - dS_manual) / max(norm(dS_manual), 1e-10) - @test rel_err < 1e-10 + @testset "Δ′ does not depend on the chunk layout" begin + # The TJ fixture has a well-conditioned Δ′; the smallest layout has two chunks per segment. + function delta_prime(nchunks) + ctrl, intr, equil, _, mats = prepared_solovev(; fixture="regression_tj_ideal_example", nchunks) + odet, riccati = FFS.riccati_eulerlagrange_integration(ctrl, equil, mats, intr) + wv = FFS.free_run(odet, ctrl, equil, mats, intr).wv + result = FFS.compute_delta_prime_matrix(intr, riccati.propagators, riccati.chunks, riccati.columns.axis_state; wv, psio=equil.psio) + return result.matrix, riccati.chunks end + dp_auto, chunks_auto = delta_prime(0) + dp_few, chunks_few = delta_prime(1) + @test length(chunks_few) < length(chunks_auto) + @test !isempty(dp_auto) && all(isfinite, dp_auto) + @test isapprox(dp_few, dp_auto; rtol=1e-4) end - @testset "compute_delta_prime_from_ca! matches inline Δ'" begin - # Verify the standalone Δ' formula matches the inline Riccati crossing computation. - # Both apply the identical diagonal formula to the same ca_l/ca_r arrays, so the - # result must be bit-for-bit identical (not just approximately equal). - # - # Note: this call overwrites intr_ric.sing[s].delta_prime; delta_prime_inline was - # saved before free_run above so it holds the original inline values. - # - FFS.compute_delta_prime_from_ca!(odet_ric, intr_ric, equil) - for s in 1:intr_ric.msing - @test intr_ric.sing[s].delta_prime == delta_prime_inline[s] - end + @testset "truncation at the edge energy peak" begin + # The chunks, their propagators and the Δ′ system follow the edge to the peak. + ctrl, intr, equil, _, mats = prepared_solovev(; over=Dict{String,Any}("truncate_at_dW_peak" => true)) + odet, riccati = FFS.riccati_eulerlagrange_integration(ctrl, equil, mats, intr) + (; propagators, chunks) = riccati + axis_state = riccati.columns.axis_state + @test intr.psilim == odet.psi_store[end] == chunks[end].psi_end + @test length(propagators) == length(chunks) && count(c -> c.needs_crossing, chunks) == intr.msing + @test odet.u_store[:, :, :, end] ≈ odet.u + free = FFS.free_run(odet, ctrl, equil, mats, intr) + @test isfinite(free.et[1]) + # The Δ′ system reproduces the plasma response of the axis-to-edge products, and says so when it does not. + delta_prime = @test_logs min_level = Logging.Warn FFS.compute_delta_prime_matrix(intr, propagators, chunks, axis_state; + wv=free.wv, wp=free.wp, psio=equil.psio) + @test size(delta_prime.matrix) == (intr.msing, intr.msing) && all(isfinite, delta_prime.matrix) + @test size(delta_prime.coil) == (2 * intr.msing, intr.numpert_total) + @test_logs (:warn, r"differs from the axis-to-edge one") FFS.compute_delta_prime_matrix(intr, propagators, chunks, axis_state; + wv=free.wv, wp=1.01 .* free.wp, psio=equil.psio) end end diff --git a/test/runtests_sing.jl b/test/runtests_sing.jl index ed477946f3..e338a8faf5 100644 --- a/test/runtests_sing.jl +++ b/test/runtests_sing.jl @@ -91,7 +91,7 @@ using FastInterpolations: cubic_interp, CubicFit, LinearBinarySearch, Series, Ex equil = load_equilibrium_from_gpec(joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example", "gpec.toml")) # default single toroidal mode used in test data; matches intr.nlow set below - ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; nn_low=1) + ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; nn_low=1, algorithm=GeneralizedPerturbedEquilibrium.Forward()) intr = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesInternal() intr.numpert_total = 32 # replacing mpert (we set equal to 32). This is the same as msol # set mode ranges so sing_der can form singfac_vec consistently diff --git a/test/runtests_slayer_inputs.jl b/test/runtests_slayer_inputs.jl index 89b27dfbdc..2d44a5877a 100644 --- a/test/runtests_slayer_inputs.jl +++ b/test/runtests_slayer_inputs.jl @@ -26,11 +26,7 @@ # Helper to build a minimal SingType without touching unused fields _mk_sing(; psi, q, q1, m, n, delta_prime=-10.0 + 0im) = SingType(; psifac=psi, rho=sqrt(psi), m=[m], n=[n], q=q, q1=q1, - delta_prime=ComplexF64[delta_prime], - delta_prime_col=zeros(ComplexF64, 0, 0), - ua_left=zeros(ComplexF64, 0, 0, 0), - ua_right=zeros(ComplexF64, 0, 0, 0), - psi_ua_left=0.0, psi_ua_right=0.0) + delta_prime=ComplexF64[delta_prime]) @testset "surface_minor_radius: continuity + outboard > 0" begin # Minor radius grows monotonically with ψ (outboard midplane). diff --git a/test/runtests_solve_api.jl b/test/runtests_solve_api.jl index 41883afea2..6fa6f17397 100644 --- a/test/runtests_solve_api.jl +++ b/test/runtests_solve_api.jl @@ -18,15 +18,16 @@ include("display_check.jl") ) wall = GPEC.Vacuum.parse_wall(get(deck, "Wall", nothing)) - # The deck's `[ForceFreeStates]` block as `solve` keywords, less `nn` and the wall; its `[Forward]` numerics - # serve both ODE integrators. + # The deck's `[ForceFreeStates]` block as `solve` keywords, less `nn` and the wall; its `[Forward]` tolerance + # serves both ODE integrators. ffs_kwargs = filter(p -> !(p.first in (:nn, :wall)), GPEC.force_balance_keywords(deck, template)) ode_kwargs = GPEC.section_keywords(deck["Forward"], "Forward", GPEC._struct_fields(Forward)) forward = Forward(; ode_kwargs...) - riccati = Riccati(; ode_kwargs...) + riccati_kwargs = GPEC.section_keywords(delete!(copy(deck["Forward"]), "renorm_threshold"), "Riccati", GPEC._struct_fields(Riccati)) + riccati = Riccati(; riccati_kwargs...) # Copy the fixture deck into `dir` and apply `[Equilibrium]` / `[ForceFreeStates]` overrides and extra sections. - # An `integrator` override moves the `[Forward]` numerics to `[Riccati]`, or drops them for Galerkin. + # An `integrator` override moves the `[Forward]` tolerance to `[Riccati]`, or drops the section for Galerkin. function _stage_deck(dir; eq_overrides=Dict{String,Any}(), ffs_overrides=Dict{String,Any}(), extra_sections=Dict{String,Any}()) for name in readdir(template) cp(joinpath(template, name), joinpath(dir, name)) @@ -36,7 +37,8 @@ include("display_check.jl") merge!(inputs["Equilibrium"], eq_overrides) merge!(inputs["ForceFreeStates"], ffs_overrides) ode = pop!(inputs, "Forward") - inputs["ForceFreeStates"]["integrator"] in ("Forward", "Riccati") && (inputs[inputs["ForceFreeStates"]["integrator"]] = ode) + inputs["ForceFreeStates"]["integrator"] == "Forward" && (inputs["Forward"] = ode) + inputs["ForceFreeStates"]["integrator"] == "Riccati" && (inputs["Riccati"] = delete!(ode, "renorm_threshold")) foreach(section -> delete!(inputs, section), ("LocalStability", "Ballooning")) merge!(inputs, extra_sections) open(io -> TOML.print(io, inputs), toml_path, "w") @@ -69,19 +71,20 @@ include("display_check.jl") mktempdir() do dir reference = GPEC.main([_stage_deck(dir; ffs_overrides=Dict{String,Any}("integrator" => "Riccati"))]).ffs prob = ForceBalanceProblem(equil; nn=1, wall=wall, ffs_kwargs...) - api = solve(prob, Riccati(; ode_kwargs..., nchunks=40)) + api = solve(prob, Riccati(; riccati_kwargs..., nchunks=40)) @test api.integrator === :riccati @test api.control.algorithm.nchunks == 40 - @test api.solution === nothing + @test api.solution !== nothing && api.solution.basis === :el_axis @test api.delta_prime !== nothing msing = length(api.surfaces) @test size(api.delta_prime.matrix) == (msing, msing) @test api.delta_prime.matrix ≈ FFS.pest3_decompose(api.delta_prime.raw).Δ - # Chunking is a decomposition of the same problem, so Δ′ tracks the TOML run. + # Chunking is a decomposition of the same problem, so Δ′ tracks the TOML run to the accuracy + # the propagator products keep at this fixture's singfac_min = 1e-6. @test size(reference.delta_prime.matrix) == size(api.delta_prime.matrix) - @test api.delta_prime.matrix ≈ reference.delta_prime.matrix rtol = 1e-6 + @test api.delta_prime.matrix ≈ reference.delta_prime.matrix rtol = 1e-4 # With the deck's own (auto) chunking the two runs are the same computation. same = solve(prob, riccati)