ForceFreeStates.Riccati - Reconstruct ξ profiles from the Δ′ BVP solution #671

Summary

An idea for discussion, not a committed plan. The Riccati integrator in ForceFreeStates (FFS) currently publishes no ξ profiles, so a Riccati result cannot feed PerturbedEquilibrium (PE), and a Riccati-fed MatchProblem gives bpen and the matched wp but no matched ξ or composite inner-layer profiles. The Δ′ boundary-value problem (BVP) that Riccati already solves contains enough information to reconstruct those profiles.

The idea

Each BVP solve returns a coefficient vector: the axis coefficients, and the coefficients of the asymptotic (Frobenius) solutions on the left and right of every rational surface. A solution column over the whole plasma is those coefficients expanded on the solutions the BVP was built from:

• Axis region. Back-substitute the axis coefficients through the per-chunk QR factors of the conditioned axis propagator, then apply the chunk propagators to get the state inside each chunk.
• Between surfaces. The shooting solutions started from the asymptotic basis on each side of a surface, multiplied by that side's coefficients, out to the midpoint where the two sides meet.
• Inside the crossing standoff. The Frobenius series itself, multiplied by the same coefficients. The on-surface point is flagged and excluded.
• ξ′. From the Euler-Lagrange relation at each node, using the conjugate momentum that the reconstruction already carries.

This gives two families of columns, in the same layout the Galerkin solver publishes:

• 2·msing big-solution columns (unit big solution on one side of one surface), from the fixed-boundary Δ′ solves;
• N error-field (EF) columns (unit edge displacement in each mode, no big solution), from the Δ_EF solve. These are the ideal-closure solution.

A cut solution (the leading-order resonant big and small terms removed near each surface) can be evaluated alongside, as the background for MatchProblem's composite inner-layer profiles.

What it would enable

• plasma_response and PE on a Riccati result, using the EF columns as the solution.
• A matched ξ and composite inner profiles from a Riccati-fed MatchProblem, as Galerkin gives today.
• ξ and dξ/dψ in gpec.h5 for Riccati runs.

Design points to settle

• Naming. The result should not be typed or labelled as Galerkin (GalerkinSolution, :gal_native, the GalerkinIntegration HDF5 group). A formalism-neutral type and basis name would let consumers stay unaware of which solver ran.
• One home. Galerkin's basis lives in ffs.galerkin.solution. Publishing both formalisms' bases in one field would avoid consumers branching on the integrator.
• Opt-in. Riccati is the default integrator, and the reconstruction needs dense solver output and stores several N × ngrid × (2·msing + N) arrays. It should sit behind a Riccati(; …) keyword and deck key with a deck/API parity case, as cut_solution does for Galerkin.
• Separation. Solving the BVP and reconstructing profiles are separate jobs and should be separate functions.
• HDF5. Where the Riccati profiles and basis are written, following docs/development/hdf5-conventions.md.

Open questions

• wp and the free-boundary energies come from the Riccati sweep, whose surface crossing differs from the BVP's. Are the reconstructed EF columns consistent with that wp to an acceptable tolerance, or should the edge momentum of the EF columns be used instead?
• How closely should the PE singular-coupling outputs (for example C_delta_prime) agree with a Forward solve and with Galerkin's EF columns, and how does that depend on the node spacing near the surfaces and on singfac_min?
• What is the time and memory cost on a DIII-D case, and how dense does the stored grid need to be?
• Which regression-harness cases should track the new output?

Status

A local prototype on feature/riccati-xi-profiles implements the reconstruction and passes ξ through to MatchProblem and PE. It is parked and not proposed for merge as it stands; it is always-on and reuses the Galerkin types and HDF5 group.