Galerkin Δ′ Solver (RDCON outer region)

Construct an empty GalCell for mpert poloidal modes and Hermite order np.

GalCell

One grid cell of the Galerkin discretization. Port of Fortran cell_type (gal.f).

Fields

• etype::GalCellType — none / resonant / extension(1,2)
• extra::GalSide — none / left / right (side of the neighbouring singular surface)
• emap::Int — global index of the extra (small resonant) DOF; 0 if none
• map::Matrix{Int} — (mpert, np+1) local→global DOF map (Fortran cell%map(mpert, 0:np))
• x::NTuple{2,Float64} — [x_left, x_right] of the cell
• x_lsode::Float64 — resonant-integration truncation point (res cells only)
• erhs::ComplexF64, ediag::ComplexF64 — resonant RHS and diagonal contributions
• rhs::Matrix{ComplexF64} — (mpert, np+1) cell RHS (driving term from the big solution)
• emat::Matrix{ComplexF64} — (mpert, np+1) coupling of the extra DOF to the Hermite DOFs
• mat::Array{ComplexF64,4} — (mpert, mpert, np+1, np+1) cell stiffness block

Cell type flanking a singular surface (Fortran cell%etype, gal.f).

GalInterval

One interval between consecutive singular surfaces (or between a surface and a domain edge). Port of Fortran interval_type (gal.f). x/dx are length nx+1 (Fortran 0:nx).

GalMatchResult

Coil-driven RPEC matched solution from the outer↔inner asymptotic matching (Fortran rmatch match_rpec). Populated by gal_match_rpec (GalerkinMatch.jl) when gal_match_flag.

• cout::Matrix{ComplexF64} — (2·msing, mcoil) outer-region plasma-solution coefficients.
• cin::Matrix{ComplexF64} — (2·msing, mcoil) inner-region coefficients.
• xi::Array{ComplexF64,3}, xi_deriv::Array{ComplexF64,3} — (mpert, ngrid, mcoil) matched ξ(ψ) and analytic ξ′(ψ) on the gal grid, one column per coil drive (identity-at-edge basis).
• deltar::Matrix{ComplexF64} — (msing, 2) inner-layer matching data (Δ₁, Δ₂) per surface.
• bpen::Matrix{ComplexF64} — (msing, mcoil) inner-layer penetrated (reconnected) resonant field at each rational surface, one column per coil drive. Read off the GGJ inner solution at the layer center (X=0) exactly as Fortran match_output_solution builds intotsol_b (match.f) — cusp-free, fit-free. Zero for the ideal branch (gal_ideal_flag), where the inner layer is skipped.
• inner_psi::Vector{Vector{Float64}} — per surface, the inner-layer ψ grid ψ_s ± X·x0/v1 (left wing reversed then right wing, so ψ ascends through ψ_s). Empty in the ideal branch.
• inner_xi::Vector{Matrix{ComplexF64}} — per surface, the inner-layer displacement ξ_ψ(ψ) on inner_psi, (length(inner_psi[s]), mcoil), one column per coil drive. This is Fortran match_solution's intotsol (deltac component 2): resc·(±Ξ₁·cin[2s] + Ξ₂·cin[2s-1]), odd parity Ξ₁ antisymmetric across ψ_s, even parity Ξ₂ symmetric, resc=(v1/x0)^(1/2+p1). The raw inner-layer contribution only (no regular outer background added), so it shares the singular asymptote with the outer near ψ_s.
• rpec_eig::Vector{ComplexF64} — forced eigenvalues γ_s = 2πi·n·f_s per surface.
• residual::Float64 — relative linear-solve residual ‖mat·cof − rmat‖/‖rmat‖.

Which side of the adjacent singular surface a cell lies on (Fortran cell%extra).

GalSingAsymp

Two-sided singular asymptotics for the Galerkin solver: separate right (sig=+1) and left (sig=−1) power series, mirroring Fortran sing.f vmatr/vmatl. Each is evaluated at a positive distance |ψ−ψ_s| (real √, real pfac, no z<0 renormalization), so the directly-extracted Δ′ matches gal.f.

GalWorkspace

Global workspace for the Galerkin solve. Port of Fortran gal_type (gal.f).

intvl has length msing+1 (Fortran intvl(0:msing), accessed as intvl[ising+1]). kl = ku = mpert*(np+1). The banded mat storage layout depends on solver (gal.f):

• "LU" → ldab = 2*kl + ku + 1, full band (LAPACK gbtrf!/gbtrs!)
• "cholesky" → ldab = kl + 1, lower band only (LAPACK pbtrf!/pbtrs!, 'L')

GalerkinResult

Outputs of the outer-region Galerkin solve.

Fields

• delta::Matrix{ComplexF64} — Δ′ matrix, (nsol, 2*msing) with nsol = 2*msing; the small resonant coefficients (Fortran delta, gal.f).
• Ap, Bp, Gammap, Deltap::Matrix{ComplexF64} — PEST-3 matching blocks, each (msing, msing) (Fortran gal_write_pest3_data, gal.f).
• msing::Int — number of resonant singular surfaces included.
• sing_psi, sing_q::Vector{Float64}, sing_m, sing_n::Vector{Int} — per-surface identifiers.
• di::Vector{Float64}, alpha::Vector{ComplexF64} — Mercier index and exponent per surface.
• solution::Union{Nothing,GalerkinSolution} — reconstructed radial ξ(ψ) and analytic ξ′(ψ) on the gal-native grid; nothing if no resonant surfaces.

GalerkinSolution

Reconstructed outer-region radial solution on the gal-native grid (ξ and analytic ξ′ consumed by PerturbedEquilibrium). Populated by gal_output_solution (GalerkinSolution.jl).

• psi::Vector{Float64}, q::Vector{Float64} — radial grid (inner→edge) and its safety factor.
• issing::Vector{Bool} — grid points sitting on a singular surface (skipped; left zero).
• xi::Array{ComplexF64,3}, xi_deriv::Array{ComplexF64,3} — (mpert, ngrid, nsol) ξ and dξ/dψ.

gal_make_grid!(intvl::GalInterval, ising::Int, msing::Int, sings::Vector{SingType},
               nn::Int, psilow::Float64, psihigh::Float64, ctrl::ForceFreeStatesControl)

Set the cell types (res/ext/ext1/ext2/none), sides, resonant integration limits, and packed node positions for interval ising (0-based, 0:msing). Port of Fortran gal_make_grid (gal.f), dx1dx2_flag=true path. intvl.x/intvl.dx have length nx+1, intvl.cells length nx (preallocated by the caller). sings is the filtered list of resonant surfaces (1-based, sings[ising] = Fortran sing(ising)).

Note (faithful to Fortran): the lower-surface width scale uses |nn·q1| (gal.f) but the upper-surface scale uses |q1| without nn (gal.f). Identical for nn = 1.

gal_make_map!(ws::GalWorkspace, mpert::Int)

Build the local→global DOF numbering with the 2*mpert C¹ overlap between adjacent Hermite cells and the interleaved extra (small resonant) DOFs (emap). Port of Fortran gal_make_map (gal.f). Sets ws.ndim. The skip/emap bookkeeping reserves a global slot for each resonant coefficient and shares it between the res cell and its neighbouring ext cell.

gal_pack(nx::Int, pfac::Float64, side::String) -> Vector{Float64}

Packed logical grid on [-1, 1]. Port of Fortran gal_pack (gal.f). For side="both" returns 2*nx+1 nodes (Fortran -nx:nx) clustered toward both ends when pfac < 1 (strong packing near the bracketing singular surfaces). pfac > 1 packs toward the centre; pfac == 1 is uniform.

gal_resonant_surfaces(intr, equil) -> (sings, psilow, psihigh)

Resonant surfaces inside the integration domain, in increasing ψ (Fortran sing(1:msing)), with the domain bounds. psilow is raised above the axis by sing_min! when qlow > qmin (excludes the q<qlow core); it falls back to the equilibrium axis bound if sing_min! was not run (intr.psilow still 0).

gal_assemble_mat!(ws, mpert)

Scatter every cell's mat/emat/ediag into the banded global matrix ws.mat. Port of gal_assemble_mat (gal.f). Supports both "LU" (full band, offset = kl+ku+1) and "cholesky" (lower band only, offset = 1).

gal_assemble_rhs!(ws, mpert)

Scatter every cell's rhs/erhs into the global RHS ws.rhs. Port of gal_assemble_rhs (gal.f). A new RHS column (isol) begins at each big-solution driver: an ext2-left cell or a res-right cell.

gal_extension!(cell, ising, ffit, profiles, intr, asymps, sings, nn, nodes, weights)

Build cell.emat, cell.ediag, cell.rhs, cell.erhs for an extension cell (ext/ext1/ext2). Port of gal_extension (gal.f). The big solution (column 1 of the gal asymptotic) drives the RHS; for ext cells the small solution (column 2) also forms the resonant coupling emat/ediag.

gal_gauss_quad!(cell, ffit, profiles, intr, nodes, weights, swap_edge)

Accumulate the nonresonant Hermite stiffness block cell.mat by Gauss-Lobatto quadrature. Port of gal_gauss_quad (gal.f). When swap_edge (the final cell of the last interval), the right-node value/slope DOFs (Fortran pb(2)↔pb(3)) are swapped so the edge-value DOF lands at index 3, matching the free-boundary BC (gal.f).

gal_get_fkg(ffit, intr, x, q) -> (F, K, G)

Evaluate the mpert×mpert matrices F = Q F̄ Qᴴ, K = Q K̄, G = Ḡ at flux x with safety factor q, where Q = diag(singfac) and singfac = m - n q (direct). Port of gal_get_fkg (gal.f). Uses the un-factored reduced ffit.fmats_gal (F̄), ffit.kmats (K̄), ffit.gmats (Ḡ). F/K/G are the ideal-MHD Euler–Lagrange coefficient matrices of the outer-region weak form (Glasser 2016, PoP 23, 112506).

gal_hermite(x, x0, x1) -> (pb, qb)

Hermite cubic basis values pb and derivatives qb on [x0, x1] (length 4, Fortran 0:3). Port of gal_hermite (gal.f).

gal_resonant!(cell, ising, ffit, profiles, intr, asymps, sings, nn, gal_tol, gal_gnstep, verbose)

Compute the resonant-cell contributions cell.erhs, cell.ediag, cell.rhs, cell.emat by adaptive quadrature (QuadGK) of the residual-operator integrand. Port of gal_lsode_int/gal_lsode_der (gal.f): the Fortran LSODE accumulation is replaced by quadgk over [x_bdy, x_lsode].

gal_set_boundary!(ws, mpert, wv_edge)

Apply the axis (u(0)=0 identity) and edge boundary conditions to the boundary cells' mat, plus the rpec coil-source RHS injection, before global assembly. Port of gal_set_boundary (gal.f), including its full three-way edge branch (rpec count ncoil = ws.nsol - 2*msing):

• rpec (ncoil > 0): identity edge and a unit source at the edge value DOF of coil column 2*msing+ipert for each poloidal mode ipert (gal.f). ws.rhs is already allocated and gal_assemble_rhs! only fills columns 1:2*msing, so the injection is not clobbered.
• free boundary (wv_edge !== nothing): add the vacuum block wvac·psio² at the edge value DOF (gal.f).
• fixed boundary (else): identity edge (gal.f).

wv_edge is nothing for the rpec and fixed cases.

gal_get_solution(ws, asymps, sings, intr, x, iintvl, icell, isol; want_d=true)
    -> (sol, dsol, iintvl, icell)

Reconstruct the displacement sol (length mpert) — and, when want_d, its analytic radial derivative dsol = d(sol)/dψ — at flux x for solution column isol, from the solved Galerkin coefficients ws.sol. Port of gal_get_solution (gal.f, non-cut path). iintvl/icell are a monotone cell cursor advanced and returned for the next call. dsol is nothing if !want_d.

gal_output_solution(ws, asymps, sings, intr, profiles, psihigh) -> GalerkinSolution

Build the gal-native packed radial grid (inner→edge; GAL_INTERP_NP_RES points per resonant/extension cell, GAL_INTERP_NP elsewhere, plus the edge point ψ=psihigh) and evaluate ξ AND the analytic ξ′ over it for every solution column. Port of gal_output_solution (gal.f); the binary/ASCII writes and the b_flag ξ→b^ψ conversion (off by default) are not ported — we keep ξ itself.

gal_match_rpec(ctrl, equil, intr, gal_result) -> GalMatchResult

Solve the coil-driven RPEC matching from the outer Δ′ (gal_result) and the per-surface inner-layer Δ. Requires gal_result.solution (reconstructed outer ξ/ξ′) and the rpec coil block gal_result.delta_coil.

The resistive path uses per-surface inputs ctrl.gal_eta/gal_rho/gal_rotation (length msing) + ctrl.gal_gamma. When ctrl.gal_ideal_flag, the inner layer is skipped and cout=0, so the matched solution is the bare ideal coil column (Fortran rmatch coil%ideal_flag) — the DCON/EL reference.

gal_matched_odestate(gal_result, ffit, intr) -> OdeState

Pack the RPEC-matched outer solution into an OdeState shaped exactly like the shooting integrator's, so PerturbedEquilibrium consumes it unchanged. Mirrors Fortran idcon_build's gal branch (idcon.f) and globalsol.bin (the on-surface issing points are dropped, match.f):

• u_store[:,:,1,ip] = ξ_ψ (matched fundamental matrix, mode×coil-drive, identity-at-edge basis)
• ud_store[:,:,1,ip] = dξ_ψ/dψ (analytic)
• ud_store[:,:,2,ip] = ξ_s = −A⁻¹(B·ξ′ + C·ξ) via ffit (same outer ideal-MHD relation as sing_der!)
• u_store[:,:,2] = 0 (PE never reads it; matches Fortran's unused u2 in the gal path)

The grid is the gal-native grid (inner→edge); step indexes the edge so build_flux_matrix derives the edge BC from u_store[:,:,1,step].

Empty GalerkinResult for a domain with no resonant surfaces.

gal_make_arrays!(ws, ctrl, equil, ffit, intr, asymps, sings, nn, wv_edge)

Assemble the global banded matrix and RHS. For each cell: Gauss-Lobatto Hermite stiffness (gal_gauss_quad!), then the resonant (gal_resonant!) or extension (gal_extension!) contributions; then the boundary conditions and the scatter into ws.mat/ws.rhs. Port of gal_make_arrays.

gal_pest3_blocks(delta, msing) -> (Ap, Bp, Gammap, Deltap)

PEST-3 matching blocks (each msing×msing) as ± combinations of the left/right Δ′ entries. Port of gal_write_pest3_data. Row index 2i-1/2i = left/right of surface i.

galerkin_solve(ctrl::ForceFreeStatesControl, equil, ffit::FourFitVars,
               intr::ForceFreeStatesInternal; vac_data=nothing) -> GalerkinResult

Compute the outer-region Δ′ matching matrix by the singular Galerkin method. Port of gal_solve (gal.f). Single toroidal mode only (intr.npert == 1). Returns a GalerkinResult; if there are no resonant surfaces in the domain it returns an empty result.

vac_data (a VacuumData from free_run!) supplies the free-boundary edge term wv_edge = vac_data.wv · psio²; pass nothing for a fixed-boundary edge.

write_galerkin!(out_h5, result::GalerkinResult)

Write the Galerkin Δ′ outputs into the open HDF5 file under the galerkin/ group. Replaces the Fortran delta_gw/pest3_data ASCII/binary outputs.