🔬 physics

A MATLAB-Based Graphical User Interface for Finite Element Solution of the Grad--Shafranov Equation

This paper introduces \textsc{MHDEquilSolver}, an interactive MATLAB-based graphical user interface that employs the finite element method to solve the Grad--Shafranov equation for axisymmetric tokamak plasmas, enabling real-time parameter adjustment, sensitivity analysis, and educational exploration of magnetohydrodynamic equilibria.

Original authors: Yahya Sadeghi, Samaneh Yarmahmoodi

Published 2026-09-28

📖 4 min read☕ Coffee break read

CC BY 4.0

Original authors: Yahya Sadeghi, Samaneh Yarmahmoodi

Original paper licensed under CC BY 4.0 (https://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Inside the heart of a fusion reactor, a superheated gas called plasma swirls in a donut-shaped chamber, held in place not by physical walls, but by an invisible cage of magnetic fields. Keeping this plasma stable is one of the most difficult challenges in modern physics. If the magnetic cage wavers, the plasma cools instantly, and the energy production stops. To design a reactor that works, scientists must first understand exactly how the plasma settles into a balanced state, a condition known as equilibrium. This balance is a tug-of-war between the pressure of the hot gas trying to expand outward and the magnetic forces squeezing it inward. Because the shape of the plasma and the strength of the magnetic fields change constantly, calculating this balance requires solving a complex mathematical puzzle that has no simple, direct answer.

For decades, researchers have relied on powerful computer programs to solve this puzzle, but these tools were often built for experts, written in older languages, and difficult to modify or understand. They acted like black boxes: you put numbers in, and a result came out, but the process in between was hidden. A new study by Yahya Sadeghi and Samaneh Yarmahmoodi from the Nuclear Science and Technology Research Institute in Iran aims to open that box. They have created a new, interactive software tool called MHDEquilSolver. Unlike previous programs, this tool is designed to be transparent and adaptable, allowing scientists and students to see exactly how changing a single setting ripples through the entire system to alter the shape and stability of the plasma.

The researchers built their tool using a popular programming environment known for its ability to handle complex math and create visual displays. At its core, the software solves the fundamental equation that governs the plasma's equilibrium. To do this, it breaks the donut-shaped space into thousands of tiny triangles, a technique called the finite element method. This allows the computer to approximate the solution by calculating the balance of forces in each small piece and then stitching them together. The software then uses a step-by-step mathematical process to refine this approximation until the solution is precise enough to be useful. The team ensured the program is robust, meaning it can handle difficult shapes or extreme settings without crashing, often by automatically simplifying the calculation if the initial attempt gets too complicated.

What makes this tool unique is its interactive interface, which puts the user in direct control of the plasma's shape and behavior. The screen displays a central image of the plasma cross-section, surrounded by six different graphs that update instantly as the user turns virtual knobs. These knobs control the physical properties of the reactor, such as the size of the donut, how stretched out the plasma is vertically, and how much it is indented on one side to form a "D" shape. Other controls adjust the internal pressure of the gas and the strength of the magnetic fields. As a user changes these values, the software immediately redraws the magnetic field lines, the pressure distribution, and the electric currents flowing within the plasma. This real-time feedback allows a user to see, for example, that increasing the vertical stretch of the plasma makes the magnetic cage stronger, or that changing the pressure profile shifts the center of the plasma slightly outward.

The authors tested the software on three distinct plasma shapes: a simple circle, an elongated oval, and the complex "D" shape used in modern reactors. In every case, the tool successfully calculated the equilibrium state and produced results that matched known theoretical solutions. The software is fast enough to run on a standard office computer, typically taking between ten and fifteen seconds to generate a full solution, even when the calculation involves thousands of tiny triangles. The researchers also performed a sensitivity analysis, a systematic check to see how much each input setting influences the final result. They found that small changes in the pressure settings can significantly shift the position of the magnetic center, while changes in the magnetic field settings alter the safety margin that keeps the plasma stable.

This work does not claim to solve the ultimate problems of fusion energy or to replace the massive, supercomputer-based codes used for final reactor design. Instead, it fills a specific gap by providing a clear, educational, and rapid prototyping environment. It serves as a bridge between abstract theory and practical application, allowing students to visualize the invisible forces at play and helping researchers quickly test new ideas before committing to more complex simulations. By making the mathematics of plasma equilibrium visible and adjustable, the tool offers a new way to explore the delicate balance required to harness the power of the stars.

Technical Summary: A MATLAB-Based Graphical User Interface for Finite Element Solution of the Grad–Shafranov Equation

Problem Statement
The Grad–Shafranov equation is the fundamental nonlinear elliptic partial differential equation governing ideal magnetohydrodynamic (MHD) equilibrium in axisymmetric toroidal plasmas, such as those in tokamaks. While established equilibrium codes like EFIT and CHEASE exist, they often rely on legacy Fortran/C++ codebases, require significant expertise to operate, lack real-time interactivity, and offer limited extensibility for educational purposes or rapid parameter scanning. There is a need for a tool that bridges the gap between complex experimental operations and academic training by providing an interactive, transparent, and modular environment for solving the Grad–Shafranov equation.

Methodology
The paper introduces MHDEquilSolver, a MATLAB-based Graphical User Interface (GUI) designed to solve the Grad–Shafranov equation using the Finite Element Method (FEM) with mesh refinement.

Mathematical Formulation: The solver addresses the equilibrium condition $J \times B = \nabla p$ $J \times B = \nabla p$ , coupling the poloidal flux function $\psi(R, Z)$ $ψ (R, Z)$ with the plasma pressure $p(\psi)$ $p (ψ)$ and the toroidal field function $F(\psi) = RB_\phi$ $F (ψ) = R B_{ϕ}$ . The code implements quadratic polynomial parameterizations for these profiles:
- $p(\psi) = p_0 + p_1\psi + p_2\psi^2$
- $F(\psi) = f_0 + f_1\psi + f_2\psi^2$
  This approach balances computational efficiency with the flexibility to model experimentally relevant profiles (e.g., peaked H-mode or broad L-mode pressure profiles).
Geometry and Meshing: The GUI supports three axisymmetric plasma boundary shapes: circular, elliptical (with elongation $\kappa$ ), and D-shaped (with triangularity $\delta$ ). The geometry is decomposed into minimal regions to satisfy MATLAB's PDE Toolbox requirements. An initial triangular mesh is generated with a maximum edge length constraint ( $H_{max} = a/5$ ); for typical parameters ( $R_0/a = 3, \kappa = 1.5$ ), this yields approximately 1,350 elements. The mesh is subsequently refined using a longest-edge bisection algorithm.
Numerical Solution: The weak form of the Grad–Shafranov equation is solved using a damped Newton method with an analytically computed Jacobian to ensure quadratic convergence near the solution. The solver imposes homogeneous Dirichlet boundary conditions on the Last Closed Flux Surface (LCFS).
Robustness and Error Handling: The implementation includes fallback strategies for mesh generation, such as detecting self-intersecting boundaries, resolving acute angles, and automatically coarsening the mesh if the primary generation fails. This ensures the tool degrades gracefully rather than crashing under extreme parameter combinations.
Interface Architecture: The GUI features a control panel for real-time adjustment of geometric parameters ( $R_0, a, \kappa, \delta$ ) and profile coefficients ( $f_0, f_1, f_2, p_1, p_2$ ). It dynamically updates six subplots upon parameter modification, displaying poloidal flux contours, mesh visualization, magnetic field components, pressure profiles, and toroidal current density.

Key Contributions

Algorithmic Transparency: The paper provides a detailed description of the numerical framework, including specific implementation insights into geometry generation, mesh adaptation, and post-processing routines, facilitating reproducibility.
Systematic Sensitivity Analysis: The authors present a quantitative analysis of how input parameters influence the equilibrium solution. The study derives partial derivatives showing that pressure coefficients ( $p_1, p_2$ ) strictly increase plasma pressure and toroidal current, thereby increasing the Shafranov shift. Similarly, it analyzes the impact of toroidal field coefficients on the safety factor ( $q$ -profile) and geometric parameters on domain scaling and stability limits.
Educational and Prototyping Utility: The tool serves as a dual-purpose platform: an educational instrument for visualizing MHD equilibrium principles and a rapid prototyping environment for screening operational scenarios and profile effects without the overhead of complex legacy codes.

Results

Verification: The numerical implementation was verified against the analytical Solov'ev equilibrium solution. Convergence studies demonstrated an $O(H^2)$ error reduction rate for linear finite elements, consistent with theoretical expectations.
Performance: On a standard workstation (Intel Core i3 12th Gen, 16 GB RAM), the solver converges within 7–10 Newton iterations. For a typical mesh of approximately 5,000 elements, total execution times range from 10 to 15 seconds depending on mesh size and refinement steps.
Case Studies: The solver successfully reproduced physical behaviors for circular, elliptical, and D-shaped configurations. Notably, it captured the outward Shafranov shift in circular cross-sections, the volume increase due to elongation, and the edge stability improvements associated with positive triangularity in D-shaped plasmas.

Significance and Claims
The paper positions MHDEquilSolver not as a replacement for high-fidelity design codes used in final reactor engineering, but as a specialized tool for education and rapid prototyping. Its significance lies in:

Accessibility: Lowering the barrier to entry for students and researchers new to MHD equilibrium by providing an intuitive, interactive interface with built-in visualization.
Interactivity: Enabling real-time exploration of parameter dependencies, which is difficult with traditional batch-processing codes.
Modularity: Offering a codebase (approx. 900 lines) that separates the control interface, computation engine, and visualization, making it extensible for future enhancements.

The authors acknowledge limitations, including the assumption of fixed plasma boundaries, axisymmetry, and the restriction to quadratic profile representations. They explicitly state that the model does not incorporate kinetic effects or non-ideal MHD physics, limiting its use for detailed design studies or capturing certain macroscopic instabilities. Future work is envisioned to include general profile functions, coupling with transport solvers, and the addition of free-boundary capabilities.

Technical Summary: A MATLAB-Based Graphical User Interface for Finite Element Solution of the Grad–Shafranov Equation

More like this