ForceFreeStates.Riccati - REFACTOR, FEATURE! - Overhaul Riccati integrator, removing extra serial integration, making Float64 BVP by default with Double64 opt-in, speedup ~x4, add riccati xi profile output

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

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"

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.
 

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).
 

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 χ⊥/χ_φ

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

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

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

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.

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]

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]

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]

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"

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