---
title: 'M3D-C1: High-Order Extended-MHD Code'
url: https://www.emergentmind.com/topics/m3d-c1
type: topic
---

# M3D-C1: High-Order Extended-MHD Code

M3D-C\(^1\) is a high-order finite-element extended-magnetohydrodynamics code for global, nonlinear simulation of magnetically confined plasmas in realistic toroidal geometry. Across the cited literature it appears as a fully 3D-capable code used in axisymmetric and non-axisymmetric settings, with fully implicit or split-implicit time advance, high-order \(C^1\)-continuous finite elements, and applications spanning vertical displacement events, energetic-particle–MHD interaction, bootstrap-current closure, and internal-kink and sawtooth dynamics in SPARC-like plasmas [1908.02387] [2107.13663] [2507.05166] [2604.02172].

## 1. Definition, scope, and code family

M3D-C\(^1\) is described as a high-fidelity extended-MHD code designed for global, nonlinear simulations of magnetically confined plasmas on transport as well as MHD timescales. In the 2025 bootstrap-current study it is characterized technically as solving single-fluid, and more general extended-MHD, equations on a high-order finite-element grid with \(C^1\) continuity, with fully implicit or split-implicit time advance so that the timestep is not limited by the Alfvén Courant condition and can reach transport timescales [2507.05166]. In the 2019 vertical-displacement-event benchmark it is likewise described as a high-order finite element code that solves the nonlinear time-dependent extended MHD equations and uses a split-implicit time advance to enable simulations over transport time scales [1908.02387].

The code family is explicitly broader than a single equilibrium class. One source states that M3D-C\(^1\) was originally developed for axisymmetric tokamaks and has been generalized to fully 3D domains so it can evolve strongly non-axisymmetric stellarators [2507.05166]. Another shows the code run in strictly axisymmetric mode for an \(n=0\) benchmark even though it is fully 3D capable [1908.02387]. A further extension, M3D-C1-K, is developed on top of the M3D-C1 finite-element code to add a kinetic energetic-particle module while retaining the MHD solver [2107.13663].

A concise view of representative uses is given below.

| Study | Problem class | Role of M3D-C\(^1\) |
|---|---|---|
| [1908.02387] | Axisymmetric VDE benchmark | One of three large-scale nonlinear MHD codes |
| [2107.13663] | EP–MHD hybrid simulation | Base code for M3D-C1-K |
| [2507.05166] | Bootstrap-current closure | Self-consistent Sauter/Redl implementation |
| [2604.02172] | SPARC sawtooth modeling | High-fidelity 3D extended-MHD solver |

This suggests that M3D-C\(^1\) is best understood not as a single-purpose disruption code, but as a general extended-MHD platform onto which additional closures and kinetic couplings can be grafted.

## 2. Numerical architecture and discretization

The numerical core is consistently described as finite-element based. In the VDE benchmark, M3D-C\(^1\) uses triangular wedge finite elements; in axisymmetry these reduce to high-order triangular finite elements in the poloidal \(R\)–\(Z\) plane [1908.02387]. In the bootstrap-current and hybrid-kinetic descriptions, the code is associated with high-order \(C^1\)-continuous finite elements and realistic toroidal geometry, including unstructured meshes and implicit time integration [2507.05166] [2107.13663].

The time-advance strategy is central to its operating regime. The split-implicit formulation is used to treat stiff MHD components implicitly, permitting long-time evolution across resistive and transport timescales rather than only Alfvénic timescales [1908.02387]. In the hybrid code, the MHD system remains implicit or semi-implicit, while energetic particles are advanced explicitly in subcycles between MHD updates [2107.13663]. In SPARC sawtooth simulations, implicit time stepping is highlighted as enabling long multi-timescale runs from Alfvénic to resistive dynamics [2604.02172].

The mesh and domain structure depend on application. For axisymmetric VDEs, a single unified mesh spans a central plasma region, a finite-thickness resistive wall region, and an outer vacuum region, so the wall is modeled as a true volumetric conductor rather than a thin shell [1908.02387]. For SPARC baseline-like scenarios, two meshes are used: a whole-device mesh including plasma, first wall, conducting vacuum vessel, and coils, and a cheaper plasma-only mesh with the boundary close to the last closed flux surface [2604.02172]. In the QA stellarator nonlinear runs reported in the bootstrap-current paper, the 3D grid uses 36 toroidal planes and about \(1.88\times 10^5\) 3D elements [2507.05166].

A plausible implication is that M3D-C\(^1\) derives much of its versatility from the combination of global geometry, high-order continuity, and implicit advancement rather than from any single plasma model.

## 3. Governing equations and physical closures

The code is used with several related MHD models. In the VDE study, the plasma region employs single-fluid, resistive MHD with mass continuity, a momentum equation including viscosity, an induction equation via Ohm’s law with Spitzer resistivity, an electron temperature equation with anisotropic heat conduction and Ohmic heating, and particle diffusion [1908.02387]. In the SPARC study, the governing model is presented as a single-fluid extended MHD with two temperatures, \(T_i\) and \(T_e\), including ion continuity, momentum balance, induction, and separate electron and ion energy equations [2604.02172].

The SPARC formulation writes
\[
\frac{\partial \mathbf{B}}{\partial t}
= \nabla\times (\mathbf{V}\times\mathbf{B} - \eta \mathbf{J}),\qquad
\mu_0 \mathbf{J} = \nabla\times\mathbf{B},
\]
with Ohm’s law
\[
\mathbf{E} = \eta \mathbf{J} - \mathbf{V}\times\mathbf{B},
\]
and total pressure \(P=P_i+P_e=n_i(T_i+Z_iT_e)\) under quasineutrality \(Z_i n_i=n_e\) [2604.02172]. The same study states explicitly that Hall, electron-inertia, and FLR terms are not used there, so the operative model is effectively single-fluid resistive or ideal MHD [2604.02172].

Geometry-dependent field representations are also specified. In the SPARC work the magnetic field and velocity are written in Clebsch-like form,
\[
\mathbf{B} = \nabla \psi_p \times \nabla\Phi - \nabla_\perp f' + I(\psi)\,\nabla\Phi,
\]
\[
\mathbf{V} = R^2 \nabla U \times \nabla\Phi + \omega R^2 \nabla\Phi + R^{-2}\,\boldsymbol{\chi},
\]
with primary unknowns \(n_i,\; T_i,\; T_e,\; U,\; \omega,\; \chi,\; \psi,\; f\) [2604.02172]. In the VDE benchmark, the wall region evolves according to
\[
\partial_t \mathbf{B} = - \nabla \times (\eta_{\textnormal{wall}} \mathbf{j}),\qquad
\mathbf{j} = \frac{\nabla \times \mathbf{B}}{\mu_0},
\]
while the vacuum region is treated with \(\mathbf{j}=0\) [1908.02387].

Boundary conditions are problem-specific and physically consequential. In the VDE benchmark, no-slip velocity is imposed at the plasma–wall interface, temperature and density are fixed there, and no boundary condition is imposed on the magnetic field at the resistive wall within the mesh, allowing halo currents to flow into and out of the wall [1908.02387]. In SPARC-like plasma-only runs, perturbations vanish at a boundary placed near the last closed flux surface [2604.02172].

## 4. Bootstrap-current modeling and quasisymmetric stellarators

A major 2025 extension adds self-consistent bootstrap-current models to M3D-C\(^1\). Two analytical closures are implemented: a generalized Sauter model and a revised Sauter-like model due to Redl et al. [2507.05166]. In both cases the quantity supplied by neoclassical theory is the flux-surface-averaged parallel current density \(\langle j_\parallel B\rangle\).

For the original Sauter model, the implementation uses
\[
\langle j_{\parallel} B\rangle =
- F(\psi)\Big( p\,L_{31}\,\frac{\partial \ln p}{\partial\psi}
+ p_e\,L_{32}\,\frac{\partial \ln T_e}{\partial\psi}
+ p_i\,\alpha L_{34}\,\frac{\partial \ln T_i}{\partial\psi} \Big),
\]
while the Redl model is implemented as
\[
\langle j_{\parallel}B\rangle =
- F(\psi)\Big( p\,L_{31}\frac{\partial \ln n}{\partial\psi}
+ p_e (L_{31}+L_{32})\frac{\partial \ln T_e}{\partial\psi}
+ p_i (L_{31}+\alpha L_{34})\frac{\partial \ln T_i}{\partial\psi} \Big)
\]
[2507.05166]. The same source notes that Sauter’s model works well at low collisionality but becomes inaccurate for \(\nu_e^* \gtrsim 1\), which motivates inclusion of Redl’s revision [2507.05166].

The bootstrap-current source enters the MHD system through a modified Ohm’s law,
\[
\mathbf{E} = -\,\mathbf{v}\times\mathbf{B} + \eta(\mathbf{J} - \mathbf{J}_x),
\]
so the induction equation becomes
\[
\frac{\partial \mathbf{B}}{\partial t} = \nabla\times(\mathbf{v}\times\mathbf{B}) - \nabla\times\big[\eta(\mathbf{J} - \mathbf{J}_x)\big].
\]
The model current is taken as purely parallel and divergence-free, with
\[
J_\parallel(\psi_t) = \frac{\langle \mathbf{J}_x\cdot\mathbf{B}\rangle}{\langle B^2 \rangle},\quad
\mathbf{J}_x = J_\parallel(\psi_t)\,\frac{\mathbf{B}}{B}
\]
[2507.05166].

The same study generalizes the closure to quasisymmetric stellarators using the Landreman isomorphism, which maps neoclassical transport in a QS stellarator to an effective axisymmetric problem [2507.05166]. In this setting M3D-C\(^1\) uses effective tokamak-like quantities such as \(qR\), \(\epsilon\), and \(f_t\), derived from stellarator geometry, and applies a Landreman-modified Sauter-Redl formula to compute bootstrap current in QA configurations [2507.05166]. The implementation is benchmarked against NEO, XGCa, and SFINCS, with percent-level agreement reported for the benchmark cases, and nonlinear QA simulations show that enabling the bootstrap model largely maintains the toroidal current density whereas disabling it allows strong decay [2507.05166].

This suggests that M3D-C\(^1\) has evolved from an extended-MHD solver with prescribed transport coefficients into a platform that can incorporate neoclassical current closures in both axisymmetric and quasisymmetric geometries.

## 5. Hybrid energetic-particle extension: M3D-C1-K

M3D-C1-K is a hybrid PIC–MHD extension built on M3D-C\(^1\) for simulations of energetic particles interacting with global plasma modes [2107.13663]. The bulk plasma remains fluid, while energetic particles are advanced kinetically. The kinetic module uses a slow-manifold “classical Pauli particle” Boris pusher, a \(\delta f\) formulation for low-noise moment calculation, and two coupling schemes, pressure coupling and current coupling [2107.13663].

The slow-manifold algorithm is designed to preserve long-term orbit fidelity. Instead of explicit full-orbit integration or conventional guiding-center RK4 alone, it advances a fictitious particle whose orbit lies on the slow manifold of the Lorentz dynamics by using the Boris update with an effective field \(\mathbf{E}^\dagger=\mathbf{E}-\mu\nabla B\) [2107.13663]. The hybrid paper reports that errors in toroidal canonical momentum and energy remain bounded under the slow-manifold Boris method, whereas RK4 guiding-center integration exhibits secular drift, especially for passing particles with large \(v_\parallel\) [2107.13663].

Energetic-particle moments are coupled back into the MHD solver either through \(-\nabla\cdot\mathbf{P}_h\) in pressure coupling or through \((\mathbf{J}-\mathbf{J}_h)\times\mathbf{B}\) in current coupling [2107.13663]. The code also supports FLR effects through 4-point gyro-averaging of fields and deposition [2107.13663]. GPU acceleration is implemented for particle pushing and weight evolution with OpenACC. On the cited Summit benchmark, 4 million particles pushed for 50 time steps gave 134.22 s for RK4 guiding center on CPUs, 12.65 s for RK4 guiding center on GPUs, and 10.32 s for slow-manifold Boris on GPUs [2107.13663].

Validation is performed on linear fishbone, toroidal Alfvén eigenmode, and reversed-shear Alfvén eigenmode problems, with good agreement reported against M3D-K, NIMROD, and other eigenvalue, kinetic, and hybrid codes [2107.13663]. The paper states that differences between pressure and current coupling are negligible in the benchmarked linear problems [2107.13663].

## 6. Benchmarks, applications, and scientific role

A defining use of M3D-C\(^1\) is benchmark-oriented cross-code comparison. In the 2019 VDE study it is one of three 3D nonlinear MHD codes, together with JOREK and NIMROD, used for axisymmetric simulations of vertical displacement events in a vertically unstable NSTX-like equilibrium [1908.02387]. That benchmark reports good agreement in linear growth rates between linear NIMROD and M3D-C\(^1\), and excellent agreement in full nonlinear VDE evolution regarding plasma location, plasma currents, and eddy and halo currents in the wall [1908.02387]. The same work identifies several physical regimes: for sufficiently small wall resistivity and large edge resistivity, the linear VDE growth rate scales as \(\gamma \propto \eta_w\), whereas at higher edge temperature or larger wall resistivity the growth rate becomes limited by the open-field-line plasma response rather than the wall alone [1908.02387].

The SPARC study uses M3D-C\(^1\) as a high-fidelity 3D extended-MHD code to investigate low-\(n\) instabilities and sawtooth crashes in relaxed baseline cases based on the SPARC Primary Reference Discharge [2604.02172]. The linear simulations identify a dominant \(n=1\) internal kink at the \(q=1\) surface, with growth rates that are strongly sensitive to the temperature profile, plasma \(\beta\), and \(q_0\) near unity [2604.02172]. Nonlinear runs show two qualitatively distinct regimes: low-\(\beta\) current-driven sawtooth-like oscillations with moderate pressure flattening, and a baseline high-\(\beta\) case in which combined current and pressure drive yields a strong sawtooth with magnetic reconnection and a hollowed pressure profile [2604.02172]. The paper interprets the baseline event as a mixture of Kadomtsev and Wesson mechanisms rather than a purely one-model crash [2604.02172].

Taken together, these studies place M3D-C\(^1\) in a particular niche. It is repeatedly used where global geometry, current redistribution, wall coupling, or profile-driven mode structure matter enough that reduced or local models are inadequate. This suggests that its principal scientific role is as a predictive extended-MHD framework for problems in which equilibrium evolution, stability, and nonlinear topology change must be handled self-consistently.

## 7. Limits, assumptions, and nomenclature

The cited literature is explicit about modeling limits. In the VDE benchmark all calculations are axisymmetric, so inherently non-axisymmetric stages of VDEs are not simulated there, and the thermal quench is imposed artificially through transport coefficients [1908.02387]. In the bootstrap-current work, Sauter and Redl are local closures with finite validity ranges, and the Landreman extension is expected to fail once strong non-QS structure or magnetic stochasticity develops; in such cases the paper advocates post-processing M3D-C\(^1\) equilibria with SFINCS or NEO 3D to quantify error [2507.05166]. In the SPARC sawtooth study, external heating, energetic particles, Hall physics, electron inertia, and FLR effects are omitted, so the reported crash thresholds and recovery times are intentionally an MHD-only baseline rather than a full burning-plasma prediction [2604.02172]. In the hybrid energetic-particle extension, the published benchmarks focus on linear physics, with nonlinear EP–MHD interaction identified as a target for future work [2107.13663].

A recurring source of ambiguity is the string “M3D” itself. Several unrelated arXiv works use “M3D” for topics such as 3D medical image analysis, dataset condensation, methylation-profile testing, and monolithic 3D DRAM, and those papers explicitly do not define “M3D-C1” as a concept within their own nomenclature [2404.00578] [2312.15927] [1410.6677] [2008.11367]. In the plasma-physics literature, by contrast, M3D-C\(^1\) consistently denotes the finite-element MHD code and its extensions [1908.02387] [2107.13663] [2507.05166] [2604.02172].

In this sense, M3D-C\(^1\) is not a generic “M3D” variant but a specific computational framework in fusion-plasma simulation: a high-order, implicit, geometry-resolving extended-MHD code that has been extended toward neoclassical bootstrap-current closure and hybrid energetic-particle kinetics while retaining its role as a benchmark-grade nonlinear MHD solver.

Source: https://www.emergentmind.com/topics/m3d-c1