EQUILIPY — Python Grad–Shafranov CutFEM Solver

EQUILIPY is a Python finite-element solver for the Grad–Shafranov (GS) equation, the nonlinear elliptic PDE that governs the magnetohydrodynamic (MHD) equilibrium of an axisymmetric magnetically confined plasma (e.g. a tokamak). Solving it yields the poloidal magnetic-flux surface configuration ψ(R, Z) at equilibrium.

The solver is built on the Cut Finite Element Method (CutFEM): the plasma boundary is defined implicitly by a level-set function and is free to cut through a fixed background mesh, so the plasma region can deform and evolve towards equilibrium without remeshing. A conforming FEM mode (body-fitted meshes) is also provided for verification.

---

Table of contents

• Goal & background
• Key features
• Repository structure
• Installation
• Quick start
• Solver description
• Meshes & element families
• Building a simulation (step by step)
• Tests, validation & convergence

---

Goal & background

The GS equation couples the plasma pressure and current profiles with the confining magnetic field. Because the current carried by the plasma depends on its own cross-section shape — which in turn is shaped by the self-induced field — the physical problem is intrinsically free-boundary. EQUILIPY solves two complementary formulations:

• FIXED-boundary — an artificial benchmark where the plasma shape (and thus its boundary) is known a priori. With the right current model and boundary condition, analytical solutions exist, making this mode ideal for verification and convergence studies.
• FREE-boundary — the physical problem, where the plasma shape is unknown. The external-magnet confinement is projected onto the computational-domain boundary via a Green's-function formalism, and the plasma domain evolves iteratively towards equilibrium.

The unfitted CutFEM discretisation is well suited to the large boundary deformations of the free-boundary problem: interfaces are tracked with level-set functions rather than by conforming the mesh.

---

Key features

• Two formulations: fixed-boundary (verification) and free-boundary (physical, iterative).
• Two discretisations: CutFEM (unfitted, level-set interface) and conforming FEM (body-fitted).
• Weak boundary/interface conditions via Nitsche's method, stabilised with ghost penalty.
• Element families up to cubic order — triangles (TRI03/06/10) and quadrilaterals (QUA04/09/16), i.e. P1–P3 / Q1–Q3, all verified against a finite-difference basis audit.
• Several plasma-current models (LINEAR, ZHENG, NONLINEAR, JARDIN, pressure/β-constrained, APEC, …) and analytical solutions for validation.
• External-magnet models: filament coils, rectangular multi-coils, quadrilateral and shaped coils.
• Critical-point detection (magnetic axis / saddle) via nonlinear optimization for normalisation.
• Reproducible test suite with automated convergence and (β, ζ) parameter-optimization scripts.

---

Repository structure

EQUILIPY/
├── src/            Solver source code
├── MESHES/         Pre-built meshes (TRI / QUA families, orders 1–3)
├── TESTs/          Test scripts (.py), notebooks (.ipynb), convergence/optimization runners
├── figures/        Figures referenced in this README
├── RESULTS/        Simulation output (created on demand; one pickle folder per run)
├── requirements.txt
└── README.md       ← you are here

Core source modules (src/)

File	Role
GradShafranovSolver.py	Main solver: assembly, linear solve, error computation, driver EQUILIPY()
Mesh.py	Mesh loading, element classification, level-set handling, plasma-boundary path, plotting
Element.py	Element integrals (domain, Nitsche, ghost), tessellation, reference geometry
FELagrangeanbasis.py	Lagrange shape functions & derivatives (orders 1–3); reference→physical mapping
GaussQuadrature.py	Gauss quadrature tables (TRI orders 1–10, QUA orders 1–5)
PlasmaCurrent.py	Toroidal-current models and analytical solutions
InitialPlasmaBoundary.py / InitialPSIGuess.py	Level-set / ψ initial models
Tokamak.py / Magnet.py	Tokamak device: first wall and external-coil models
Greens.py	Green's-function boundary projection (free boundary)
_critical.py, _update.py, _L2error.py, _Bfield.py, _plot.py, _output.py	Critical points, level-set update, error norms, magnetic field, plotting, I/O

---

Installation

Clone the repository:

git clone https://github.com/Elmanyer/EQUILIPY.git
cd EQUILIPY

Install the required Python packages (ideally in a virtual environment):

pip install -r requirements.txt

The core dependencies are numpy, scipy, and matplotlib; Jupyter is included to run the example notebooks.

---

Quick start

Ready-to-run examples live in TESTs/ (both .py scripts and .ipynb notebooks):

cd TESTs
python TestFIXED-CutFEM-LINEAR.py       # fixed-boundary CutFEM, analytical LINEAR case

Results are written to a new RESULTS/… folder. To reproduce the convergence study across all element families and solvers:

python run_optimization.py    # optimal (beta, zeta) for structured CutFEM -> optimization_results.txt
python run_convergence.py     # FEM + CutFEM convergence rates             -> convergence_results.txt

---

Solver description

The Grad–Shafranov problem

Δ*ψ = μ₀ R j_φ(ψ)      in the plasma region Ω⁻
Δ*ψ = 0                 in the vacuum region Ω⁺   (free boundary only)

with Δ* = R ∂/∂R (1/R ∂/∂R) + ∂²/∂Z². The plasma boundary Γ is the level-set zero contour; ψ = ψ_D is enforced weakly on Γ (Nitsche). In free boundary the outer-wall condition ψ = ψ_B comes from the Green's-function projection of the coil currents.

CutFEM discretisation

• The background mesh does not conform to Γ. Elements are classified as plasma interior (Dom = −1), cut (Dom = 0), vacuum (Dom = +1), or computational boundary (Dom = +2).
• Cut elements are tessellated into sub-elements with adapted quadratures.
• The interface condition is imposed with Nitsche's method (penalty beta); near-degenerate cut configurations are conditioned by ghost-penalty stabilization (penalty zeta).
• In fixed boundary the PDE is integrated only over the active mesh (plasma ∪ cut elements); the vacuum is decoupled.

Free-boundary workflow (double loop)

1. Outer loop — converge the boundary values ψ_B projected onto the computational-domain boundary via the Green's-function formalism.
2. Inner loop — with those boundary values, solve the GS boundary-value problem (updating the plasma level-set each iteration) until the poloidal flux converges.

Both loops take tolerances and maximum-iteration limits as inputs. The fixed-boundary problem uses a single outer iteration.

---

Meshes & element families

Meshes are named <FAMILY>_<KIND>_<h>, e.g. QUA16_REC_0.1, TRI06_FEM_0.05. Refinement levels available: 1.0, 0.5, 0.1, 0.05, 0.02.

Family	Element	Order	KIND (mesh kinds available)
TRI03 / TRI06 / TRI10	triangle	1 / 2 / 3	FEM (body-fitted), REC_STRUC (structured cut), REC_UNSTR (unstructured cut)
QUA04 / QUA09 / QUA16	quadrilateral	1 / 2 / 3	FEM (body-fitted), REC (structured cut)
TRI03_ITER_FIRSTWALL	—	—	Tokamak first-wall / computational-boundary mesh

• FEM meshes are body-fitted (no cut elements) — used by the conforming FEM solver.
• REC meshes are rectangular backgrounds cut by the level-set — used by CutFEM.
• Use QuadratureOrder2D = 5 for QUA families (orders 1–5 implemented) and 8 for TRI.

Meshes were generated with GiD.

---

Building a simulation (step by step)

The TESTs/ files show complete examples; the essential steps are:

1. Create the solver and set parameters

Equilibrium = GradShafranovSolver()
Equilibrium.FIXED_BOUNDARY     = True        # True: fixed boundary; False: free boundary
Equilibrium.GhostStabilization = True        # ghost penalty on/off
Equilibrium.QuadratureOrder2D  = 8           # 5 for QUA families
Equilibrium.QuadratureOrder1D  = 6
Equilibrium.beta = 100.0                     # Nitsche penalty
Equilibrium.zeta = 0.0                       # ghost penalty
Equilibrium.ext_maxiter = 1                  # outer loop (free boundary uses several)
Equilibrium.int_maxiter = 50; Equilibrium.int_tol = 1e-10
# critical-point solver initial guesses: R0_axis/Z0_axis, R0_saddle/Z0_saddle, opti_maxiter, opti_tol
Equilibrium.InitialiseParameters()

2. Computational-domain mesh

Equilibrium.MESH = Mesh('TRI06_REC_STRUC_0.1')

3. Tokamak device (first wall always recommended; magnets mandatory for free boundary)

wall = Mesh('TRI03_ITER_FIRSTWALL')
coil = QuadrilateralCoil(name='PF1', Itotal=5.73e6, Xcenter=np.array([3.94, 7.57]), Area=0.25)
Equilibrium.TOKAMAK = Tokamak(WALL_MESH=wall, MAGNETS=[coil])   # MAGNETS optional for fixed boundary

4. Initial plasma boundary (level-set; fixed for fixed-boundary, an initial guess otherwise)

Equilibrium.initialPHI = InitialPlasmaBoundary(EQUILIBRIUM=Equilibrium, GEOMETRY='LINEAR',
                                               R0=6.0, epsilon=0.32, kappa=1.7, delta=0.33)
Equilibrium.DomainDiscretisation(INITIALISATION=True)

If the initial plasma boundary coincides with the mesh boundary (a body-fitted FEM mesh), the whole mesh is the plasma domain and the problem reduces to standard FEM.

5. Initial ψ guess

Equilibrium.initialPSI = InitialGuess(EQUILIBRIUM=Equilibrium, PSI_GUESS='LINEAR',
                                      R0=6.0, epsilon=0.32, kappa=1.7, delta=0.33)
Equilibrium.InitialisePSI()

6. Plasma-current model

Equilibrium.PlasmaCurrent = CurrentModel(EQUILIBRIUM=Equilibrium, MODEL='LINEAR',
                                         R0=6.0, epsilon=0.32, kappa=1.7, delta=0.33)

7. Launch

Equilibrium.EQUILIPY("TS-CASE")     # writes results to RESULTS/TS-CASE-<mesh>/

Available GEOMETRY/PSI_GUESS models include LINEAR, ZHENG, CUBICHAM, …; current MODELs include LINEAR, ZHENG, NONLINEAR, JARDIN, APEC, and constrained variants (see the class docstrings in src/).

---

Tests, validation & convergence

The TESTs/ folder contains the verification suite:

File(s)	Purpose
TestFIXED-FEM-LINEAR.*, TestFIXED-CutFEM-LINEAR.*	Fixed-boundary, analytical LINEAR case (FEM / CutFEM)
TestFIXED-*-NONLINEAR.*	Nonlinear current model
TestFREE-ITER-*.ipynb	Free-boundary ITER-like cases (coils / shaped coils)
run_optimization.py	(β, ζ) grid search for structured CutFEM → optimization_results.txt
run_convergence.py	Convergence rates for FEM + CutFEM (structured & unstructured), all orders → convergence_results.txt
_runner.py	Shared execution core (run_case) used by the two runners

Validated convergence (LINEAR case, fixed boundary): all solvers reach the optimal O(h^(p+1)) rate — O(h²) for linear, O(h³) for quadratic, O(h⁴) for cubic elements (the finest cubic meshes saturate near solver precision), summarised in convergence_results.txt:

Relative L2 error versus 1/h for the fixed-boundary LINEAR case. Markers denote the solver/mesh type (○ FEM, △ CutFEM structured, □ CutFEM unstructured) and colours the element family; dashed grey lines are the reference h²/h³/h⁴ slopes for linear, quadratic and cubic elements.

---

EQUILIPY — developed at the Barcelona Supercomputing Center (BSC), Computer Applications in Science and Engineering (CASE), Nuclear Fusion group. Author: Pau Manyer Fuertes (pau.manyer@bsc.es).