Nonlinear Fusion Magneto-Hydrodynamics with Finite Elements


C. R. Sovinec

Los Alamos National Laboratory







presented at the


International Sherwood Fusion Theory Conference

March 27-29, 2000

Los Angeles, California





The objective of this presentation is to discuss the advantages and pitfalls of using finite elements for time-dependent magnetohydrodynamic simulations. Computations of resistive tearing with the NIMROD code show convergence rates that are near theoretical predictions.










Ahmet Aydemir IFS

James Callen U-WI

Ming Chu GA

John Finn LANL

Tom Gianakon LANL

Charlson Kim CU-Boulder

Scott Kruger SAIC

Jean-Noel Leboeuf UCLA

Richard Nebel LANL

Scott Parker CU-Boulder

Steve Plimpton SNL

Nina Popova MSU

Dalton Schnack SAIC

Carl Sovinec LANL

Alfonso Tarditi SAIC





NIMROD solves Maxwell's equations plus a set of fluid-moment relations:







where quasineutrality, , is assumed, wp is the plasma frequency, and r is the mass density.










where is a local basis of the piecewise polynomial representation.






We note that our time-advance for B has a lot in common with the time-independent incompressible Navier-Stokes equations.


Magnetic field equation: [Simplifying to resistive MHD]


Incompressible Navier-Stokes


The application of finite element methods to incompressible Navier-Stokes equations has received considerable mathematical analysis over the last three decades. And, we can take advantage of that work!


See, for example:

M. D. Gunzburger, "Mathematical aspects of finite element methods for incompressible viscous flows," in Finite Element Theory and Application,D. L. Dwoyer, M. Y. Hussaini, and R. G. Voigt, eds., (Springer-Verlag, 1988).

C. Cuvelier, A. Segal, and A. A. van Steenhoven, Finite Element Methods and Navier Stokes Equations, (D. Reidel Publishing Co., 1986).




Incompressible Navier-Stokes Analysis





Weak Form

for all [Subscript implies b. c. satisfied.]

for all [Subscript implies 0 mean.]




Beside the nonlinear term, the divergence-free condition makes INS a much more challenging application than ordinary structural mechanics problems.


This introduces finite compressibility, and solutions should be checked in the limit of .




A necessary condition for convergence with the integrated and penalty function methods is that the solution space satisfies the Ladyzhenskaya - Babuska - Brezzi or divergence - stability condition:


There exists a g>0, independent of the mesh spacing, h, such that

for the chosen family of solution spaces,



This condition ensures that the discrete divergence-free vector field approaches a continuous divergence-free field in the limit as





Z is defined as the set of subspaces satisfying



The end result is a limit on the acceptable sets of {Vh, Sh}.









If divergence-stability is satisfied and a unique solution exists (conditions are not turbulent), then convergence is assured.










Time-dependent Problems

A discrete time advance allows alternative implementations of the divergence constraint. On a fundamental level, all share characteristics with the steady-state treatments.

Faraday's law is modified to:

* B. Marder, "A Method for Incorporating Gauss' Law into Electromagnetic PIC Codes," J. Comput. Phys. 68, 48 (1987).


This is related to the "projection method" of Brackbill and Barnes.§

§ J. U. Brackbill and D. C. Barnes, "The Effect of a Nonzero on the Numerical Solution of Magnetohydrodynamic Equations," J. Comput. Phys. 35, 426 (1980).




Error diffusion is also closely related to the penalty method for incompressible Navier-Stokes.




The penalty method needs to satisfy divergence-stability.




Is the space represented by acceptable?






The NIMROD finite element representation has been generalized to allow arbitrary degree polynomial bases of the Lagrange type (H1).







A linear tearing mode serves as a test for divergence error and convergence.

The geometry is a circular cross section, periodic cylinder. The polar grid is composed of quadrilaterals. S=104 and Pm=1.


Equilibrium B


Linear Eigenfunction ~




We have varied kb over orders of magnitude to compare bilinear (32x32 mesh) and biquadratic solutions (16x16 mesh).

Linear Growth Rate vs. log(kb)

Divergence Error, , vs. log(kb)




In a similar, but doubly-periodic slab, using kb=50h was catastrophic for bilinear elements:


Repeating the computation either with kb=h and bilinear elements or with kb=50h and biquadratic elements avoided the vertex-to-vertex oscillations.




Convergence of the solution and error improve with polynomial degree in the cylindrical problem, as predicted.

Linear Growth Rate vs. Mesh Spacing

ln(Divergence Error) vs. ln(Mesh Spacing)

Radial and azimuthal resolution are varied together. [kb=2000h]




Convergence of the solution and error are compared for the cylindrical problem, varying only the number of cells in the azimuthal direction. [kb=2h]

Linear Growth Rate vs. Mesh Spacing

ln(Divergence Error) vs. ln(Mesh Spacing)





A semi-implicit algorithm [Schnack, J. Comput. Phys. 70, 330, for example] is used to advance the solution from initial conditions.

A simple example with a 1D linear sound wave provides a review:

A semi-implicit scheme replaces mass density in the momentum equation with a new spatial operator

After normalizing pressure by P0 , velocity by the sound speed, and assuming an eikx dependence, the modified leap-frog advance has the form,

The complex eigenvalue representing the linear operations of the time step,



The magnitude of l is unity (advance is numerically stable without dissipation) if

Setting the coefficient

ensures stability for all Dt.


A coefficient exceeding the linear numerical stability requirement stabilizes nonlinear terms, which may then appear in explicit form.




The choice of semi-implicit operator has a strong influence on accuracy at large Dt (large compared with propagation times of normal modes).

The NIMROD semi-implicit scheme for the MHD advance includes the following equation for velocity:







For nonlinear tokamak simulations, the coefficient for the isotropic part of the operator may be much smaller than that for the anisotropic part.










This poster will be available from , along with many other presentations on the NIMROD code development project.