Equilibrium Solver

The equilibrium solver is the numerical core of SCPN-Fusion-Core. It solves the Grad-Shafranov equation for axisymmetric magnetostatic equilibria, the foundation upon which all transport, stability, and control computations are built.

The Grad-Shafranov Equation

In an axisymmetric toroidal plasma, force balance between the plasma pressure gradient, the magnetic pressure, and the magnetic tension reduces to a single nonlinear elliptic PDE known as the Grad-Shafranov (GS) equation (Grad & Rubin 1958, Shafranov 1966):

(1)\[\Delta^{*}\psi \;=\; -\mu_0 R^2\,p'(\psi) \;-\; F(\psi)\,F'(\psi)\]

where:

  • \(\psi(R,Z)\) is the poloidal magnetic flux function

  • \(\Delta^{*} = R \frac{\partial}{\partial R}\!\left(\frac{1}{R}\frac{\partial}{\partial R}\right) + \frac{\partial^2}{\partial Z^2}\) is the toroidal elliptic (Grad-Shafranov) operator

  • \(p(\psi)\) is the plasma pressure profile

  • \(F(\psi) = R B_\phi\) is the poloidal current function (diamagnetic function)

  • \(\mu_0\) is the permeability of free space

The magnetic field components are recovered from \(\psi\) as:

\[B_R = -\frac{1}{R}\frac{\partial\psi}{\partial Z}, \qquad B_Z = \frac{1}{R}\frac{\partial\psi}{\partial R}, \qquad B_\phi = \frac{F(\psi)}{R}\]

Numerical Method

SCPN-Fusion-Core discretises equation (1) on a uniform \((R, Z)\) grid using a five-point finite-difference stencil with the correct \(1/R\) toroidal geometry factor.

Picard Iteration

The nonlinear system is solved by Picard (fixed-point) iteration with under-relaxation:

\[\psi^{(k+1)} = (1 - \alpha)\,\psi^{(k)} + \alpha\,\psi_{\text{new}}^{(k)}\]

where \(\alpha\) is the relaxation factor (default 0.1) and \(\psi_{\text{new}}\) is obtained by solving the linearised GS equation at the current iterate.

The inner linear system is solved by one of two methods:

Red-Black SOR (default)

Successive Over-Relaxation with red-black ordering for cache-friendly memory access. This is the production default and is well-suited for grids up to 128x128.

Multigrid V-Cycle (selectable)

A geometric multigrid solver with standard restriction, prolongation, and Gauss-Seidel smoothing. Activated via:

kernel.set_solver_method("multigrid")

The multigrid path is asymptotically faster for large grids but carries slightly more implementation complexity.

Profile Models

Two pressure/current profile parametrisations are supported:

L-mode (default)

Linear profiles: \(p(\hat\psi) = p_0(1 - \hat\psi)\) and \(FF'(\hat\psi) = c(1 - \hat\psi)\) where \(\hat\psi\) is the normalised flux.

H-mode

Modified hyperbolic tangent (mtanh) pedestal profiles following the Hirshman-Whitson parametrisation:

\[T(\rho) = T_{\text{ped}} \cdot \frac{1}{2}\left[1 + \tanh\!\left(\frac{\rho_{\text{ped}} - \rho}{\Delta}\right)\right] + T_{\text{core}} \cdot (1 - \rho^2)\]

Boundary Conditions

The solver supports both fixed-boundary (Dirichlet) and free-boundary modes. In free-boundary mode, the vacuum field from external coils is computed and matched at the computational boundary.

Convergence Criteria

Convergence is declared when the relative \(L^2\) residual of the GS operator falls below a threshold (default \(10^{-6}\)):

||Delta* psi + mu_0 R^2 p' + F F'||_2 / ||psi||_2 < tol

Solver Tuning

The relaxation factor, grid resolution, and convergence threshold are the primary tuning parameters. Detailed guidance is provided in docs/SOLVER_TUNING_GUIDE.md.

Relaxation factor recommendations

Range

Behaviour

When to use

0.02–0.05

Very conservative, slow but stable

High-beta plasmas, strongly shaped equilibria

0.08–0.12

Default range (0.1 shipped)

Most ITER-class and medium-beta scenarios

0.15–0.25

Aggressive, faster convergence

Low-beta L-mode, small grids (33x33)

> 0.3

Likely to diverge

Not recommended

For compact high-field designs (SPARC-class, \(B > 10\,\text{T}\)), a relaxation factor of 0.08 is usually more stable than 0.1.

Configuration

The solver reads its parameters from a JSON configuration file:

{
    "solver": {
        "max_iterations": 500,
        "convergence_threshold": 1e-6,
        "relaxation_factor": 0.1,
        "grid_size": 65,
        "method": "sor"
    }
}

GEQDSK I/O

The eqdsk module reads and validates the standard G-EQDSK equilibrium file format used throughout the fusion community (Lao et al. 1985). Eight SPARC GEQDSK files from CFS are included for validation.

from scpn_fusion.core.eqdsk import read_geqdsk

eq = read_geqdsk("validation/reference_data/sparc/lmode_vv.geqdsk")

The reader extracts the full set of equilibrium data: poloidal flux grid \(\psi(R,Z)\), pressure and current profiles \(p(\psi)\) and \(FF'(\psi)\), boundary/limiter geometry, and scalar parameters (\(B_\text{T}\), \(I_p\), axis position, separatrix flux values).

Neural Equilibrium Solver

For real-time control applications where millisecond-scale equilibrium updates are required, the neural equilibrium solver (neural_equilibrium.py) provides a PCA + MLP surrogate that predicts flux-surface geometry from an 8-dimensional input vector:

\[(I_p,\; B_T,\; R_\text{axis},\; Z_\text{axis},\; p'_\text{scale},\; FF'_\text{scale},\; \psi_\text{axis},\; \psi_\text{bry}) \;\longrightarrow\; \hat{\psi}(R, Z)\]

The surrogate is trained on SPARC GEQDSK data with profile perturbations. Inference time is approximately 5 microseconds per point (Rust backend, synthetic weights).

Warning

No pre-trained neural weights are shipped with the package. Users must train the surrogate on their own simulation data using the provided train_on_sparc() convenience function.

Force Balance and 3D Geometry

The force_balance module provides iterative force-balance checking against the GS equation. The geometry_3d module generates 3D flux-surface meshes (OBJ format) from Fourier boundary parametrisations for visualisation and CAD integration.