ForceFreeStates - IMPROVEMENT - Compute solution derivatives on demand

src/ForceFreeStates/EulerLagrange.jl

@@ -174,7 +174,6 @@ function serial_eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil::E
 
     # Initialization
     odet = OdeState(intr.numpert_total, ctrl.numsteps_init, ctrl.numunorms_init, intr.msing)
-    odet.du_store_populated = true
     if ctrl.sing_start <= 0
         initialize_el_at_axis!(odet, ctrl, ffit, equil.profiles, intr)
     elseif ctrl.sing_start <= intr.msing
@@ -643,9 +642,8 @@ function cross_ideal_singular_surf!(
     # so the result is in a different convention. The canonical Δ' is the STRIDE BVP matrix
     # (compute_delta_prime_matrix!) populated by the parallel FM path.
 
-    sing_der!(du1, odet.u, params, odet.psifac)
-
     # Store values after crossing step and advance
+    odet.q = equil.profiles.q_spline(odet.psifac; hint=odet.spline_hint)
     store_ode_data!(odet, odet.psifac, odet.u)
 end
 
@@ -686,10 +684,8 @@ function cross_kinetic_singular_surf!(
     sing_der!(du2, odet.u, params, odet.psifac)
     odet.u .+= (du1 .+ du2) .* dpsi
 
-    # re-evaluate at the final post-crossing u so the stored derivatives match u_store
-    sing_der!(du1, odet.u, params, odet.psifac)
-
     # Store crossing step
+    odet.q = equil.profiles.q_spline(odet.psifac; hint=odet.spline_hint)
     store_ode_data!(odet, odet.psifac, odet.u)
 end
 
@@ -741,7 +737,6 @@ function integrate_el_region!(
     q_range = abs(q_end - q_start)
 
     steps_in_segment = Ref(0)
-    du_buffer = zeros(ComplexF64, intr.numpert_total, intr.numpert_total, 2)
 
     function segment_callback!(integrator)
         ctrl, _, _, intr, odet, chunk = integrator.p
@@ -760,7 +755,8 @@ function integrate_el_region!(
         in_edge_scan = ctrl.psiedge < intr.psilim && integrator.t >= ctrl.psiedge
 
         if near_start || near_end || (odet.total_steps % ctrl.save_interval == 0) || in_edge_scan
-            sing_der!(du_buffer, integrator.u, integrator.p, integrator.t)
+            # q at the accepted point, not the last internal Runge-Kutta stage
+            odet.q = equil.profiles.q_spline(integrator.t; hint=odet.spline_hint)
             store_ode_data!(odet, integrator.t, integrator.u)
         end
     end
@@ -773,7 +769,7 @@ function integrate_el_region!(
     # Guarantees the pre-crossing (or pre-edge) state is always stored in u_store,
     # regardless of where the last accepted step landed relative to the near_end band.
     if odet.step == 1 || odet.psi_store[odet.step-1] != sol.t[end]
-        sing_der!(du_buffer, sol.u[end], (ctrl, equil, ffit, intr, odet, chunk), sol.t[end])
+        odet.q = equil.profiles.q_spline(sol.t[end]; hint=odet.spline_hint)
         store_ode_data!(odet, sol.t[end], sol.u[end])
     end
 
@@ -1007,12 +1003,14 @@ function transform_u!(odet::OdeState, intr::ForceFreeStatesInternal)
             odet.u_store[:, :, 1, istep] .= gauss_buffer
             mul!(gauss_buffer, odet.u_store[:, :, 2, istep], transforms[:, :, ifix])
             odet.u_store[:, :, 2, istep] .= gauss_buffer
-            mul!(gauss_buffer, odet.du_store[:, :, 1, istep], transforms[:, :, ifix])
-            odet.du_store[:, :, 1, istep] .= gauss_buffer
-            mul!(gauss_buffer, odet.du_store[:, :, 2, istep], transforms[:, :, ifix])
-            odet.du_store[:, :, 2, istep] .= gauss_buffer
-            mul!(gauss_buffer, odet.xi_s_store[:, :, istep], transforms[:, :, ifix])
-            odet.xi_s_store[:, :, istep] .= gauss_buffer
+            # Derivative stores are empty unless a path filled them analytically (galerkin);
+            # materialized ones are computed after this transform and need no fixup.
+            if !isempty(odet.du_store)
+                mul!(gauss_buffer, odet.du_store[:, :, istep], transforms[:, :, ifix])
+                odet.du_store[:, :, istep] .= gauss_buffer
+                mul!(gauss_buffer, odet.xi_s_store[:, :, istep], transforms[:, :, ifix])
+                odet.xi_s_store[:, :, istep] .= gauss_buffer
+            end
         end
         jfix = kfix + 1
     end

src/ForceFreeStates/ForceFreeStatesStructs.jl

@@ -505,11 +505,22 @@ and a small set of temporary matrices and factors used to compute singular-layer
   - `u_store::Array{ComplexF64,4}` - Stored solution arrays at each saved step with shape
     `(numpert_total, numpert_total, 2, numsteps_init)` (complex solution state used by the solver).
 
-  - `du_store::Array{ComplexF64,4}` - du₁/dψ and du₂/dψ at the accepted point of each saved step, same shape as `u_store`.
+  - `du_store::Array{ComplexF64,3}` - dΞ_ψ/dψ (the u₁ block only) at each saved step, shape
+    `(numpert_total, numpert_total, step)`. Empty until `materialize_derivative_stores!` fills it,
+    except on the galerkin-matched path which supplies the analytic derivative at construction.
+    du₂/dψ is never stored densely — its only consumer evaluates it on demand at bracket nodes.
 
-  - `xi_s_store::Array{ComplexF64,3}` - Clebsch displacement Ξ_s at each saved step, eq. 18 of Glasser 2016, shape `(numpert_total, numpert_total, numsteps_init)`.
+  - `xi_s_store::Array{ComplexF64,3}` - Clebsch displacement Ξ_s at each saved step, eq. 18 of Glasser 2016,
+    shape `(numpert_total, numpert_total, step)`. Empty until materialized, same as `du_store`.
 
-  - `du_store_populated::Bool` - True once the serial EL path has filled `du_store`/`xi_s_store`; other paths leave it false.
+  - `u_store_el_basis::Bool` - True when `u_store` holds the Euler-Lagrange state `(u₁, u₂)`, so the
+    derivative kernel can be re-applied to it. False on the sparse parallel path, whose stored columns
+    are chunk-endpoint Riccati matrices; `materialize_derivative_stores!` refuses to run there.
+
+  - `du_store_populated::Bool` - True once `du_store`/`xi_s_store` hold valid data in the final
+    (post-transform, post-normalization) basis. Set by `materialize_derivative_stores!` or by the
+    galerkin-matched constructor; stays false where the stores cannot be materialized, e.g. the
+    sparse parallel path whose solution is in the Riccati basis.
 
   - `crit_store::Vector{Float64}` - Stored crit parameter values (smallest eigenvalue of W⁻ꜝ) (length `numsteps_init`).
 
@@ -527,10 +538,6 @@ and a small set of temporary matrices and factors used to compute singular-layer
 
   - `u::Array{ComplexF64,3}` - Current working solution arrays with shape `(numpert_total, numpert_total, 2)`.
 
-  - `du::Array{ComplexF64,3}` - du/dψ from the latest `sing_der!` call, shape `(numpert_total, numpert_total, 2)`.
-
-  - `xi_s::Matrix{ComplexF64}` - Ξ_s from the latest `sing_der!` call, shape `(numpert_total, numpert_total)`.
-
   - `ising_start::Int` - Index of the starting singular surface to be crossed during integration.
 
     # Initialization parameters
@@ -578,8 +585,9 @@ and a small set of temporary matrices and factors used to compute singular-layer
     psi_store::Vector{Float64} = Vector{Float64}(undef, numsteps_init)
     q_store::Vector{Float64} = Vector{Float64}(undef, numsteps_init)
     u_store::Array{ComplexF64,4} = Array{ComplexF64}(undef, numpert_total, numpert_total, 2, numsteps_init)
-    du_store::Array{ComplexF64,4} = Array{ComplexF64}(undef, numpert_total, numpert_total, 2, numsteps_init)
-    xi_s_store::Array{ComplexF64,3} = Array{ComplexF64}(undef, numpert_total, numpert_total, 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} = Array{ComplexF64}(undef, numpert_total, numpert_total, 2, msing)
@@ -592,8 +600,6 @@ and a small set of temporary matrices and factors used to compute singular-layer
     psifac::Float64 = 0.0
     q::Float64 = 0.0
     u::Array{ComplexF64,3} = zeros(ComplexF64, numpert_total, numpert_total, 2)
-    du::Array{ComplexF64,3} = zeros(ComplexF64, numpert_total, numpert_total, 2)
-    xi_s::Matrix{ComplexF64} = zeros(ComplexF64, numpert_total, numpert_total)
     ising_start::Int = 0
     psimax::Float64 = 0.0
     needs_crossing::Bool = false

src/ForceFreeStates/Free.jl

@@ -28,9 +28,10 @@ end
 """
     normalize_eigenfunctions!(odet::OdeState, wt::AbstractMatrix{ComplexF64}, psio::Float64) -> OdeState
 
-Rescale the stored EL solution vectors in `odet.u_store`, `du_store`, and `xi_s_store` so the
-edge displacement matches the free-boundary eigenvectors `wt` (scaled by `2π·psio·1e-3`).
-Modifies `odet` in place. Call after `free_run` when downstream code consumes the stored ξ profiles.
+Rescale the stored EL solution vectors in `odet.u_store` (and `du_store`/`xi_s_store` when a
+path has already filled them) so the edge displacement matches the free-boundary eigenvectors
+`wt` (scaled by `2π·psio·1e-3`). Modifies `odet` in place. Call after `free_run` when
+downstream code consumes the stored ξ profiles.
 """
 @with_pool pool function normalize_eigenfunctions!(odet::OdeState, wt::AbstractMatrix{ComplexF64}, psio::Float64)
     N = size(wt, 1)
@@ -41,12 +42,14 @@ Modifies `odet` in place. Call after `free_run` when downstream code consumes th
         odet.u_store[:, :, 1, istep] .= tmp_mat
         mul!(tmp_mat, odet.u_store[:, :, 2, istep], coeffs)
         odet.u_store[:, :, 2, istep] .= tmp_mat
-        mul!(tmp_mat, odet.du_store[:, :, 1, istep], coeffs)
-        odet.du_store[:, :, 1, istep] .= tmp_mat
-        mul!(tmp_mat, odet.du_store[:, :, 2, istep], coeffs)
-        odet.du_store[:, :, 2, istep] .= tmp_mat
-        mul!(tmp_mat, odet.xi_s_store[:, :, istep], coeffs)
-        odet.xi_s_store[:, :, istep] .= tmp_mat
+        # Only the galerkin path carries derivative stores this early; materialized ones
+        # are built from the normalized u_store afterwards.
+        if !isempty(odet.du_store)
+            mul!(tmp_mat, odet.du_store[:, :, istep], coeffs)
+            odet.du_store[:, :, istep] .= tmp_mat
+            mul!(tmp_mat, odet.xi_s_store[:, :, istep], coeffs)
+            odet.xi_s_store[:, :, istep] .= tmp_mat
+        end
     end
     return odet
 end

src/ForceFreeStates/Galerkin/GalerkinMatch.jl

@@ -264,33 +264,21 @@ function gal_matched_odestate(gal_result::GalerkinResult, ffit::FourFitVars, int
     ngrid_f = length(psi_f)
 
     u_store = zeros(ComplexF64, mpert, mpert, 2, ngrid_f)
-    du_store = zeros(ComplexF64, mpert, mpert, 2, ngrid_f)
+    du_store = zeros(ComplexF64, mpert, mpert, ngrid_f)
     xi_s_store = zeros(ComplexF64, mpert, mpert, ngrid_f)
 
-    amat = Matrix{ComplexF64}(undef, mpert, mpert)
-    bmat = Matrix{ComplexF64}(undef, mpert, mpert)
-    cmat = Matrix{ComplexF64}(undef, mpert, mpert)
-    tmp = Matrix{ComplexF64}(undef, mpert, mpert)
     hint = Ref(1)
     for ip in 1:ngrid_f
         ξ = @view xi_f[:, ip, :]
         ξ′ = @view dxi_f[:, ip, :]
         @views u_store[:, :, 1, ip] .= ξ
-        @views du_store[:, :, 1, ip] .= ξ′
-        # ξ_s = −A⁻¹(B·ξ′ + C·ξ), exactly as the ideal path of sing_der! (Sing.jl:1015-1049)
-        ffit.amats(vec(amat), psi_f[ip]; hint=hint)
-        ffit.bmats(vec(bmat), psi_f[ip]; hint=hint)
-        ffit.cmats(vec(cmat), psi_f[ip]; hint=hint)
-        LinearAlgebra.LAPACK.potrf!('U', amat)
-        LinearAlgebra.LAPACK.potrs!('U', amat, bmat)   # bmat ← A⁻¹ B
-        LinearAlgebra.LAPACK.potrs!('U', amat, cmat)   # cmat ← A⁻¹ C
-        xs = @view xi_s_store[:, :, ip]
-        mul!(tmp, bmat, ξ′)
-        xs .= .-tmp
-        mul!(tmp, cmat, ξ)
-        xs .-= tmp
+        @views du_store[:, :, ip] .= ξ′
+        # ξ_s = −A⁻¹(B·ξ′ + C·ξ), the same node quantity the Euler-Lagrange path computes.
+        @views compute_node_xi_s!(xi_s_store[:, :, ip], ξ′, ξ, ffit, psi_f[ip]; hint=hint)
     end
 
+    # Derivatives here are the analytic galerkin ξ′, not recomputable from the ODE kernel.
     return OdeState(; numpert_total=mpert, numunorms_init=1, msing=gal_result.msing, numsteps_init=ngrid_f,
-        step=ngrid_f, total_steps=ngrid_f, psi_store=psi_f, q_store=q_f, u_store=u_store, du_store=du_store, xi_s_store=xi_s_store)
+        step=ngrid_f, total_steps=ngrid_f, psi_store=psi_f, q_store=q_f, u_store=u_store, du_store=du_store,
+        xi_s_store=xi_s_store, du_store_populated=true)
 end

src/ForceFreeStates/Riccati.jl

@@ -1009,12 +1009,6 @@ See: Glasser (2018) Phys. Plasmas 25, 032507 — Eq. 19 (dual Riccati form)
     # dS = w†·v - S·Ḡ·S  [Glasser 2018 eq. 19, dual Riccati]
     mul!(dS, adjoint(w), v)   # dS = w†·v
 
-    # Store du1/dψ = Q·v as a diagnostic before v is reused
-    # Q·v = diag(singfac_vec)·v = Ξ'_Ψ (displacement gradient, with U₂ = I)
-    @. odet.du[:, :, 1] = singfac_vec * v
-    @view(odet.du[:, :, 2]) .= 0
-    odet.xi_s .= 0
-
     # Subtract S·Ḡ·S (reuse v and tmp to avoid extra allocation)
     mul!(tmp, gmat, S)        # tmp = Ḡ·S
     mul!(v, S, tmp)           # v   = S·Ḡ·S
@@ -1794,9 +1788,11 @@ function _assemble_propagators_serially!(odet::OdeState, propagators::Vector{Chu
             riccati_cross_ideal_singular_surf!(odet, ctrl, equil, ffit, intr, chunk.ising)
             last_crossing_step = odet.step - 1
         else
-            # Save non-crossing end-of-chunk state. du_store is not meaningful here — when
-            # ctrl.populate_dense_xi=true the entire odet is replaced by a serial-EL pass
-            # at the end of parallel_eulerlagrange_integration.
+            # Save non-crossing end-of-chunk state. These columns are FM/Riccati chunk
+            # endpoints, not the Euler-Lagrange state — when ctrl.populate_dense_xi=true the
+            # entire odet is replaced by a serial-EL pass at the end of
+            # parallel_eulerlagrange_integration.
+            odet.u_store_el_basis = false
             if odet.step >= size(odet.u_store, 4)
                 resize_storage!(odet)
             end