Nonlinear Fusion MagnetoHydrodynamics with Finite
Elements
C. R. Sovinec
Los Alamos National Laboratory
and
The NIMROD Team
presented at the
International Sherwood Fusion Theory Conference
March 2729, 2000
Los Angeles, California
OBJECTIVE
The objective of this presentation is to discuss
the advantages and pitfalls of using finite elements for
timedependent magnetohydrodynamic simulations. Computations of
resistive tearing with the NIMROD code show convergence rates
that are near theoretical predictions.
OUTLINE
 Introduction

 Project summary
 Code features
 Finite element discretization

 Relation to incompressible NavierStokes
 Review of INS Analysis
 Formulation
 Divergencestability
 Convergence
 Timedependent problems
 Error diffusion
 Test results
 Semiimplicit advance
 Summary
PRESENT TEAM MEMBERS AND ADVISORS
Ahmet Aydemir IFS
James Callen UWI
Ming Chu GA
John Finn LANL
Tom Gianakon LANL
Charlson Kim CUBoulder
Scott Kruger SAIC
JeanNoel Leboeuf UCLA
Richard Nebel LANL
Scott Parker CUBoulder
Steve Plimpton SNL
Nina Popova MSU
Dalton Schnack SAIC
Carl Sovinec LANL
Alfonso Tarditi SAIC
EQUATIONS
NIMROD solves Maxwell's equations plus a set of
fluidmoment relations:
where quasineutrality, , is assumed,
w_{p} is the
plasma frequency, and r is the mass density.

 Density evolution will be added.

 T. Gianakon has implemented neoclassical closures
(1D33).

 C. Kim is implementing particlebased closures (1D07).
FINITE ELEMENT DISCRETIZATION

 The finite element approach starts with the solution space
and not the differential operators.

 After choosing a family of discrete spaces, a field is
represented by a sum over basis functions:
,
where is a local basis of the piecewise polynomial
representation.

 Substituting the discrete representation for each field
into a weak form of the set of partial differential equations,
or into its variational form, leads to either a linear or a
nonlinear matrix equations for the set of all
coefficients.

 A finite element representation is said to be
conforming if the chosen discrete space is a subset of
the admissible space for the variational form of the original
set of PDEs.
We note that our timeadvance for B has
a lot in common with the timeindependent incompressible
NavierStokes equations.
Magnetic field equation: [Simplifying to
resistive MHD]
Incompressible NavierStokes
The application of finite element methods to
incompressible NavierStokes 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., (SpringerVerlag, 1988).
C. Cuvelier, A. Segal, and A. A. van Steenhoven,
Finite Element Methods and Navier Stokes Equations, (D.
Reidel Publishing Co., 1986).
Incompressible NavierStokes Analysis
PDE
Weak Form
for all [Subscript implies b. c. satisfied.]
for all [Subscript implies 0 mean.]
Beside the nonlinear term, the divergencefree condition
makes INS a much more challenging application than ordinary
structural mechanics problems.

 Using discrete spaces for the primitive variables is the
most common and flexible approach, but the divergence issue
must be dealtwith directly.

 The "integrated method" solves the weak finite element
form with the divergence constraint equation.
 The "penalty method" modifies the continuity equation
to
.
This introduces finite compressibility, and solutions
should be checked in the limit of .


 Divergencefree elements are also possible (see R.
Gruber and J. Rappaz, Finite Element Methods in Linear
Ideal Magnetohydrodynamics (SpringerVerlag, 1985) but
they may be the most restrictive.
 A potential representation may be used to satisfy the
divergence condition analytically. [The unstructured M3D code
does thissee H. R. Strauss, et. al., "MHD Simulations on an
Unstructured Mesh," in Proceedings of the Workshop on Nonlinear
MHD and Extended MHD, March 2526, 1998, UWCPTC 981.]

 In general, this leads to lowerorder spatial
convergence on the physical fields, however.
 Applying boundary conditions may be more complicated
than a primitivevariable formulation.
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
divergencefree vector field approaches a continuous
divergencefree 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 {V^{h}, S^{h}}.

 Acceptable 2D element pairs include:


 continuous, 6pt quadratic V^{h} and
continuous, linear S^{h} on triangles
(TaylorHood)
 continuous, 9pt biquadratic V^{h} and
continuous, bilinear S^{h} on quadrilaterals
(TaylorHood)
 continuous, 6pt quadratic V^{h} and
piecewise constant S^{h} on triangles
(CrouzeixRaviart)
 continuous, 7pt quadratic V^{h} and
piecewise linear S^{h} on triangles
(CrouzeixRaviart)
 continuous, 9pt biquadratic V^{h} and
piecewise linear S^{h} on quadrilaterals
(CrouzeixRaviart)
 Unfortunately, simpler element pairs are not
acceptable:

 continuous, linear V^{h} and
piecewiseconstant S^{h} on triangles [This gives
an overdetermined system.]
 continuous, bilinear V^{h} and
piecewiseconstant S^{h} on quadrilaterals [A
checkerboard pressure mode allows any velocity field to
satisfy the divergence constraint.]
 Other element pairs fail in a more subtle way; see the
Gunzburger reference.
If divergencestability is satisfied and
a unique solution exists (conditions are not turbulent), then
convergence is assured.

 The difference between the finite element solution and the
best piecewise polynomial approximation for the solution is
independent of mesh spacing. The rate of spatial convergence is
then tied to the polynomial degree used for the basis
functions.

 This applies to regular or irregular meshes, hence the
attractiveness of the finite element method for domains with
complicated borders.
 If the solution is sufficiently smooth,
.
Timedependent Problems
A discrete time advance allows alternative
implementations of the divergence constraint. On a fundamental
level, all share characteristics with the steadystate
treatments.

 NIMROD uses an errordiffusion technique* for controlling .
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.§
 Advance to B* with physical terms.
 Solve
^{§} 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 NavierStokes.
The penalty method needs to satisfy
divergencestability.

 If p^{h} is eliminated numerically after
choosing an acceptable {V^{h}, S^{h}}
pair, the condition is still satisfied.

 If p is eliminated analytically, acceptability
depends on the form of .
Is the space represented by acceptable?

 This condition is not satisfied for linear
V^{h} on triangles, or for bilinear
V^{h} on quadrilaterals.

 For 7pt. quadratic V^{h} on triangles and
9pt. biquadratic V^{h} on quadrilaterals, the
piecewise linear scalar field is a subset of all possible
, so the
analytic elimination penalty method or error diffusion method
are acceptable. (Effectively producing CrouzeixRaviart
elements)
The NIMROD finite element representation has been
generalized to allow arbitrary degree polynomial bases of the
Lagrange type (H^{1}).

 By selecting second order (and higher?)
degree polynomials, error diffusion is divergencestable with
arbitrarily large values of k_{b.}


 This statement does not imply that every
computation with an acceptable solution space has zero
magnetic divergence,

 nor does it imply that magnetic divergence has no
effect on the rest of the solution.

 It implies that spatial convergence with an acceptable
solution space removes magnetic divergence and its
effects.

 We have not yet verified that the Fourier representation of
the periodic direction does not change acceptability for
toroidal components where .

 We have successfully used bilinear and linear elements for
B in many computations. However, some cases fail
catastrophically. Why does it work with unacceptable elements
in some cases?

 A clue is that k_{b must be ~} h.
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=10^{4} and P_{m}=1.
Equilibrium B
Linear Eigenfunction ~
We have varied k_{b} over orders of magnitude
to compare bilinear (32x32 mesh) and biquadratic solutions (16x16
mesh).
Linear Growth Rate vs.
log(k_{b})
Divergence Error, , vs. log(k_{b})
In a similar, but doublyperiodic slab, using
k_{b}=50h was
catastrophic for bilinear elements:

 A slow numerical mode arose after 0.15 diffusion
times.

 There is a
mode in the envelope of the eigenfunction.
Repeating the computation either with
k_{b}=h and
bilinear elements or with k_{b}=50h
and biquadratic elements avoided the
vertextovertex 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. [k_{b}=2000h]
Convergence of the solution and error are compared
for the cylindrical problem, varying only the number of cells in
the azimuthal direction. [k_{b}=2h]
Linear Growth Rate vs. Mesh Spacing
ln(Divergence Error) vs. ln(Mesh Spacing)
SEMIIMPLICIT ADVANCE
A semiimplicit 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 semiimplicit scheme replaces mass density in
the momentum equation with a new spatial operator
After normalizing pressure by
P_{0} , velocity by the sound speed, and assuming
an e^{ikx} dependence, the modified leapfrog advance has
the form,
The complex eigenvalue representing the linear operations of
the time step,
satisfies
where
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 semiimplicit operator has a strong influence
on accuracy at large Dt (large compared with propagation times of normal
modes).
The NIMROD semiimplicit scheme for the MHD advance
includes the following equation for velocity:

 After integration by parts in the finite element procedure,
the linear part of the operator has the form of the linear
potential energy of the ideal MHD system.

 This is very similar to the scheme by K. Lerbinger and J.
F. Luciani, J. Comput. Phys. 97, 444 (1991), but the
0subscripted fields here contain the n=0 part of the
solution, in addition to the equilibrium, and the coefficient
for the isotropic operator is updated with the magnitude of
nonsymmetric perturbations to avoid CFL constraints imposed by
nonlinear waves.

 This is computationally efficient for fusion problems where
the magnetic field is dominantly symmetric. There is no
toroidal coupling in the matrices, permitting parallel
decomposition over n.

 Hall terms, when used, are timesplit from the MHD advance
of B, and require a separate semiimplicit operator.
[Harned, J. Comput. Phys. 83, 1]
For nonlinear tokamak simulations, the coefficient
for the isotropic part of the operator may be much smaller than
that for the anisotropic part.

 A simulation of a (3,2) tearing mode in DIIID serves as an
example.

 Coefficients are proportional to the square of the
respective CFL number.

 See the presentation by Schnack, et al. (1C24) for more
details on this computation.
SUMMARY

 Numerical analysis of finite elements applied to the
incompressible NavierStokes equations has proven extremely
valuable for understanding and improving the treatment of
magnetic field in the NIMROD code.

 Our original formulation with linear and bilinear elements
did not satisfy the divergencestability condition.

 The error diffusion technique with k_{b}~h
works with loworder elements in cases where
the generation of is easily
controllable.

 The divergencestability condition is satisfied with
quadratic elements (and higher order elements?) and the error
diffusion technique.

 The rate of convergence increases with polynomial
degree, and results are less sensitive to the value of
k_{b}.

 A semiimplicit advance is computationally efficient for
fusion applications.
This poster will be available from
http://nimrodteam.org , along with many other presentations on
the NIMROD code development project.