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

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

With grid_type = "auto" and mpsi = 0 (the defaults; "log_asymptotic" is a legacy alias), the radial grid is built by a two-pass measured-curvature refinement driven by the single accuracy knob psi_accuracy (τ):

Pass 1 forms the equilibrium on a coarse three-region layout (geometric in log(ψ) at the core, uniform in the middle, geometric in log(1−ψ) at the edge).
A knot density is derived from one sizing rule — the cubic-spline derivative error model err(f′) ≈ h³|f''''|/24 ≤ τ·|f| — applied to three sources of curvature |f''''|:
- Measured (mid-radius and edge): divided differences on the pass-1 nodes of the 1D profiles (F, P, dV/dψ, q), the 2D geometry channels (r², η-offset, ν, Jacobian) along sampled θ-lines (the bicubic ψ-axis shares the same knots), and the kinetic profiles (n, T, ω_E) when loaded, so steep pedestal gradients attract knots. Nodal values come from independent field-line integrations, so the estimate is immune to inter-knot spline ringing. Curvature is measured against ρ = √ψ, in which the equilibrium is regular at the magnetic axis — in ψ itself the geometry channels behave as ψ^(k/2) (R−R₀ ~ √ψ), so their ψ-curvature diverges under refinement while their ρ-curvature converges; the ρ-spacing maps back through dψ/dρ = 2√ψ, which by itself packs √ψ-tight toward the axis.
- Separatrix model (edge floor, ψ ≥ 0.9): the same rule on q ≈ −A·ln(1−ψ) gives geometric-in-log(1−ψ) packing with uniform relative q′ error, independent of A. Normally inactive — the pass-1 layout is already log-packed at the edge, so the measured density dominates — it only guards a pass-1 that under-sampled the divergence.
- Axis model (core, ψ ≤ 0.03): the same rule on the power-law axis form gives geometric-in-log(ψ) packing (constant-ratio spacing ∝ ψ down to psilow). Here the model replaces measurement: nodal data from the smallest flux surfaces is dominated by integration and axis-extrapolation noise, which grid refinement amplifies rather than resolves.
- Rational surfaces: local packing around every q = m/n in the requested n range (pinned as mandatory knots by the main driver), at any ψ including inside the modeled core region. This local resolution of the ψ-splined Euler-Lagrange coefficient matrices is what converges the Δ′ boundary-value problem.
The density integral is equidistributed into the knot vector, mandatory rational-surface knots are inserted with a minimum-spacing snap guard, and pass 2 re-forms the equilibrium on the refined grid from the in-memory input (no file re-read).

Every region's knot count scales as τ^(-1/3), so tightening psi_accuracy refines the core, pedestal, and edge proportionally. The legacy grid_type = "ldp" (sin²-spaced), "pow1", "uniform", and explicit mpsi > 0 (single-pass, fixed layout) are still supported. Library users calling setup_equilibrium directly with mpsi = 0 receive the coarse pass-1 grid; use refined_psi_grid and the override_psi_nodes keyword to apply the refinement manually.

The packing on the DIII-D-like example (n=1) compared to fixed ldp grids — note the coarse spacing across the smooth mid-radius, the spacing dips at each rational surface, and the core/pedestal/edge packing (benchmarks/plot_grid_knot_placement.jl regenerates this figure):

Radial knot packing: auto two-pass vs ldp

Decomposing the density by source on the same example shows the pedestal band (ψ_N ≈ 0.85–0.98) is driven by measured curvature, not the edge floor: the pressure, q, and dV/dψ profiles contribute comparably, and the rzphi geometry channels — the Grad-Shafranov response to the same pedestal p′ (Shafranov-shift / angle-offset steepening) — contribute the most at every node in the band. Because the density is measured from the formed solution, the packing follows the pedestal wherever it sits:

Knot-density decomposition by source

API Reference

GeneralizedPerturbedEquilibrium.Equilibrium.AnalyticEqSpec — Type

AnalyticEqSpec

Single-source-of-truth entry describing one analytic-equilibrium kind: the gpec.toml section that carries its parameters, the *Config type built from that section, and the run function that turns the config into solver input. tj_analytic and tj_analytic_direct share TJAnalyticConfig but bind different run functions, so run_fn is per-entry rather than derived from config_type.

Fields

section::String — gpec.toml section key (e.g. "SOL_INPUT")
config_type::DataType — *Config type, constructed via its (::Dict) / (::String) ctor
run_fn::Function — (EquilibriumConfig, config) -> eq_input solver entry point

field	meaning
`n_i`	main-ion density
`n_e`	electron density
`T_i`/`T_e`	ion / electron temperature
`omega_E`	ExB rotation
`omega_tor`	toroidal rotation (optional)
`chi_e`	perpendicular heat diffusivity χ⊥
`chi_phi`	toroidal momentum diffusivity χ_φ
`species_densities`	named per-species densities (e.g. `"n_D"`, `"n_T"`) for explicit multi-ion input; `nothing` for ASCII / no extras
`provenance`	short string recording the source file/format