Equilibrium Module

The Equilibrium module provides tools to read, construct, and analyze magnetohydrodynamic (MHD) plasma equilibria. It wraps a variety of equilibrium file formats (EFIT, CHEASE, SOL, LAR, and others), runs the appropriate direct or inverse solver, and returns a processed PlasmaEquilibrium object ready for downstream calculations.

Overview

Key responsibilities of the module:

Read equilibrium input files and TOML configuration (see EquilibriumConfig).
Provide convenient constructors for analytic / model equilibria (Large

Aspect Ratio, Solev'ev).

Build spline representations used throughout the code (1D cubic and

2D bicubic splines).

Run the direct or inverse equilibrium solver and post-process results

(global parameters, q-profile finding, separatrix finding, GSE checks).

The module exposes a small public API that covers setup, configuration, and common analyses used by other GPEC components (e.g. ForceFreeStates, vacuum interfaces).

EFIT solver strategies

When eq_type is one of the EFIT-based options, three solver strategies are available:

`eq_type`	Method	Best for
`"efit"`	Geometric-angle field-line ODE	Default; fast and robust for standard cases
`"efit_arclength"`	Arc-length field-line ODE	Near-separatrix surfaces where the geometric-angle ODE becomes singular
`"efit_by_inversion"`	Contour.jl marching-squares → inverse solver	Highest geometric accuracy; avoids ODE singularities entirely

efit integrates field lines using the geometric angle η (polar angle from the magnetic axis, 0→2π) as the independent variable. Position at each step is computed as R = R₀ + rfac·cos(η), so the RHS denominator is Bz·cos(η) − Br·sin(η) — the projection of B_pol onto the radial direction. Near an X-point where Bp → 0 this denominator passes through zero, causing a coordinate singularity: the ODE steps become extremely small or the solver fails to converge.

efit_arclength uses arc length along the flux surface as the independent variable instead. The tangent direction (dR/ds, dZ/ds) = ∇ψ⊥/|∇ψ| has unit magnitude by construction, so the position equations have no denominator and remain well-behaved all the way to the separatrix. The 1/Bp factors in the accumulated integrals still diverge near X-points, but their solver tolerances are set large so they do not restrict the step size — only the position tracking drives adaptivity. The two methods produce identical results when both succeed; efit_arclength extends the reliable range of psihigh closer to 1.

efit_by_inversion traces flux surface level sets directly from ψ(R,Z) using marching squares, resamples each closed curve to a uniform geometric-angle grid, and feeds the result into the inverse equilibrium solver (the same path used by CHEASE input). The Cartesian evaluation grid is clipped to the separatrix bounding box and its resolution is set adaptively from a bilinear interpolation error bound, so no manual tuning is needed.

Radial grid packing

The default grid_type = "log_asymptotic" uses a three-region grid that respects the asymptotic behavior of q near both the magnetic axis and the separatrix:

Core (ψ < 0.15): geometric spacing in log(ψ) — handles the axis where profiles behave as ψⁿ
Middle (0.15 ≤ ψ ≤ 0.95): uniform spacing — preserves resolution through the pedestal region
Edge (ψ > 0.95): geometric spacing in log(1−ψ) — tracks the logarithmic divergence q ~ −A·ln(1−ψ) near a diverted separatrix

With mpsi = 0 (the default), the number of radial knots is chosen automatically from the psi_accuracy parameter (target absolute error in q). Two probe field-line integrations near psihigh estimate the local log-slope A, and the knot count is set so the cubic spline error stays below psi_accuracy throughout the domain. The legacy grid_type = "ldp" (sin²-spaced) and explicit mpsi are still supported.

API Reference

GeneralizedPerturbedEquilibrium.Equilibrium.DirectBField — Type

DirectBField

Internal mutable struct to hold B-field components and their derivatives at a point. It is used as a temporary workspace to avoid allocations in tight loops.