The following is published in Plasma Phys.
Control. Fusion **41** (1999) A747-A755. http://www.iop.org/EJ/welcome

The NIMROD code: a new approach to numerical plasma physics

A. H. Glasser, C. R. Sovinec, and R. A. Nebel

Los Alamos National Laboratory, MS K717, Los Alamos, NM 87545, USA

T. A. Gianakon

University of Wisconsin-Madison, 1500 Engineering Drive, Madison, WI 53706, USA

S. J. Plimpton

Sandia National Laboratory, MS 1111, Albuquerque, NM 87185, USA

M. S. Chu

General Atomics Corporation, P. O. Box 8909, San Diego, CA 92186, USA

D. D. Schnack

Science Applications International Corporation, 10260 Campus Point Drive, San Diego, CA 92121, USA

and the NIMROD Team

**Abstract.** NIMROD is a code development project
designed to study long-wavelength, low-frequency, nonlinear
phenomena in toroidal plasmas with realistic geometry and
dynamics. The numerical challenges of solving the fluid
equations for a fusion plasma are discussed, and our
discretization scheme is presented. Simulations of a resistive
tearing mode show that time steps much greater than the Alfven
time are possible without loss of accuracy. Validation tests of
a resistive interchange mode in a shaped equilibrium, a
ballooning mode, and nonlinear activity in reversed-field
pinches are described.

**PACS number:** 52.65.-y

1. Introduction

Among the most important phenomena in present-day fusion experiments are various forms of nonlinear global magnetic activity. In tokamaks, saturated magnetohydrodynamic (MHD) modes can lock to walls and fields errors, often leading to plasma disruption [1,2]. Nonlinear neoclassical effects produce magnetic islands, inducing transport that precludes operation near ideal MHD stability limits; hence, overall performance is diminished [3,4]. In alternative concepts such as the reversed-field pinch (RFP) and spheromak, nonlinear global magnetic activity determines the very nature of the configuration. Thus, numerical simulation and theoretical understanding of these phenomena are crucial for advancing magnetically confined fusion. The NIMROD [Non-Ideal MHD with Rotation, Open Discussion] code development project seeks to address these issues.

There are many pitfalls in modeling this activity numerically. The fluid model of a magnetized plasma forms a system of equations that is stiff and very anisotropic. The nonlinear behavior of interest occurs on resistive tearing to resistive diffusion time scales, while waves propagate many orders of magnitude faster. Furthermore, the equations need to be solved in realistic geometric configurations-including open field lines beyond a separatrix if appropriate for a particular experiment-to capture the complete global nature of the activity. Equally important is the need to accurately represent magnetic resonant layers that conform to equilibrium magnetic surfaces in the plasma core. An additional challenge stems from the fact that the fastest present-day computers have massively parallel architectures. NIMROD addresses these difficulties by using a finite element spatial discretization of the poloidal plane, an accurate semi-implicit time-advance, and three-dimensional domain decomposition linked by message passing communication. In this paper, we briefly describe essential elements of the equations and numerical algorithm, and we present some validation simulations and initial research efforts.

2. Physics model

The fluid equations we solve are based on velocity moments of the Boltzmann equation, where closure terms may be selected for particular conditions. We presently consider two species; however, it is numerically advantageous to form a center-of-mass velocity equation and a generalized Ohm's law to permit a semi-implicit time advance [5,6]. The system of equations is then Maxwell's equations plus

where quasineutrality, , is assumed, *w _{p}* is the plasma frequency, and

The stress tensor terms can be important, particularly in non-collisional plasmas. While most of our simulations have not included them, initial efforts at simulating neoclassical magnetic islands use a Chew-Goldberger-Low form

where the magnetic pumping coefficients are computed from an analytic model for the Braginskii and banana regimes [8]. We have ordered the contributions in the perturbed fields, and only first-order terms are retained at present. Other analytic closures will be tested, and some type of particle closure may eventually be used [9,10].

3. Numerical algorithm

The discrete form of the equations in NIMROD uses a finite element representation of the poloidal plane and Fourier series for the periodic direction (either toroidal or periodic linear). The two-dimensional finite elements allow arbitrary poloidal cross sections, but the Fourier series requires geometric uniformity in the periodic direction. We employ two types of finite elements to achieve flexibility without compromising accuracy: 1) quadrilateral elements where geometric and equilibrium quantities are bicubic but dependent variables are bilinear, and 2) triangular elements where all quantities are linear. A complete grid may be assembled from blocks of each element type in the fashion shown in Figure 1. A grid generator aligns the core grid lines to equilibrium flux surfaces, and packing may be used to provide extra resolution at resonant surfaces. A triangulation routine is used to discretize the rest of the domain up to a realistic experiment wall. The grid blocks may be allocated to different processors on parallel computers, with message passing communication at common borders. We plan to implement a means of adapting the grid to changing flux surfaces.

The Fourier series representation for the periodic direction provides accuracy for derivatives in that direction, and they ease several numerical tasks. In our semi-implicit time advance, the matrix equations do not couple different Fourier components. This considerably reduces the burden of inverting matrices, and it permits parallel decomposition of the Fourier direction (in addition to the poloidal decomposition). Different Fourier components couple through explicit terms, which are treated in a pseudospectral fashion [11]. An efficient parallel fast-Fourier-transform (FFT) is used to transform fields to functions of toroidal angle, where nonlinear terms are easily computed, then another FFT transforms the product back to Fourier series.

The stiffness of the equations is another critical issue, since the phenomena we wish to simulate are many orders of magnitude slower than wave propagation times. [We have chosen to use the fluid equations without aspect ratio orderings or assumptions on divergence of flow, so that the code may be applied to a variety of situations without adding special terms to capture the relevant physics.] Optimally, the equations are advanced on time steps that are only one to two orders of magnitude below the physics time scale (based on a Taylor expansion of the evolution). Many implicit and semi-implicit methods permit numerically stable simulations at large time step. However, retaining accuracy is challenging. Spatial truncation errors may couple different normal modes, permitting accuracy only at small time step. This is less catastrophic than 'spectral pollution' in eigenvalue computations [12], although the source of error is similar. In resistive MHD, for example, the shearing behavior near a singular layer is nearly force free, since . If truncation error proportional to Dt coupled this shearing behavior to a compressional response, small time steps would be required to accurately reproduce the mode.

We have addressed this issue with a time-split, semi-implicit method [5]. For each time step, velocity is advanced with an operator that is based on ideal MHD. Prior to spatial discretization and integration by parts, the numerical equation has the form,

(1)

where the subscripted terms are the toroidal average fields,
and viscous terms have been omitted for clarity. The differential
operators on the left side of the equation form the semi-implicit
operator for MHD physics. When spatial discretization and
integration by parts are applied, the terms with numerical
coefficient *f _{a}* represent the linear potential
energy of the ideal MHD system. The operator with coefficient

The magnetic field is updated in two segments, with Hall terms advanced separately from MHD terms. The Hall terms require another semi-implicit operator to avoid numerical instabilities due to whistler waves [6]. Our present approach solves the following twice for the Hall segment,

(2)

where , and
*f _{h}* is a numerical parameter. The predictor step
advances half of the total time step, and the term uses the predicted
values in the corrector step. Linear and nonlinear tests of this
approach are being conducted to determine its effectiveness.

A disadvantage of using low order polynomials for primitive
physical fields is that this representation of **B** cannot
have zero divergence everywhere, except in special cases. To
prevent an accumulation of divergence over many time steps, we
have employed an error diffusion technique [16]. The term
, which is
diffusive for the longitudinal part of the field, is either added
to the MHD segment or a separate parabolic diffusion equation is
solved to 'clean' the field. Minimum acceptable values for
*k* are problem dependent, but we have found that
*k* must not significantly exceed the electrical
diffusivity to avoid unphysical behavior. Truncation error from
the term can affect
the physical, solenoidal part of **B**, and this error is
proportional to *k*.

We have investigated the numerical properties of NIMROD with
dispersion relations of the discrete equations and many
validation tests. The numerical dispersion relations-derived in
uniform conditions to lowest order in -demonstrate that the choice of primary
dependent variables to be expanded as bilinear elements is
critical [17]. Formulations with electromagnetic and
electrostatic potentials as the primary variables couple
different normal modes through truncation errors, although the
divergence of **B** problem is avoided. Formulations with
magnetic field as the primary variable factor shear and
compressional waves to all orders of Dt
and Dx. This property has been borne out
by computations of shear Alfven waves where ; they converge to accurate frequencies at
time steps that are large compared to compressional time scales.
Closer to simulations of research value are linear tearing modes
for a 0-*b* paramagnetic equilibrium in a
periodic cylinder. Figure 2 shows temporal convergence of the
growth rate at time steps above the propagation time of a
compressional wave across the entire minor radius. This
demonstrates that our time step is not tied to Alfvenic time
scales for numerical stability or accuracy of low frequency
phenomena, making it possible to simulate tearing behavior at
large Lundquist number (, where *t _{r}* is the
diffusion time and

The matrices resulting from the semi-implicit operators in equations (1) and (2) are symmetric, positive definite, but they are very ill-conditioned. A large burden is therefore placed on the numerical linear algebra used to solve the matrix equations in order that computational gains be realized at large time step. We have used preconditioned conjugate gradients with an additive Schwarz technique [18] applied to the task of preconditioning across grid blocks that may be allocated to different processors. Within each grid block, we precondition with point Jacobi, line Jacobi, incomplete Cholesky, or a direct solve. All of these methods deteriorate as the number of vertices and grid blocks increases. To make them more scalable, we are investigating the use of a coarse-grid component in the preconditioner.

4. Testing and validation

An essential task for a new code like NIMROD is to verify that
physics results it produces are consistent with previous
analytical and computational results. A recent linear result of
interest is an unstable resistive interchange mode in the D-IIID
tokamak [19,20]. The profile of the safety factor *q* in
this equilibrium has negative central shear, with a central value
, a minimum value
, and an edge
value of . There
are singular
surfaces at and
0.573, with *y* the normalized poloidal flux. The
mode is unstable because of a positive value of the
pressure-curvature term *D _{R}* at the inner
singular surface.
Figure 3 shows the structure of the eigenfunction from an MHD
simulation with .
It is strongly localized to the neighborhood of the inner
surface, as
expected from analytical theory. The poloidal mode structure is
as predicted by the analytical theory, with dominant poloidal
mode number on each
surface. The growth rate and eigenfunction are in good agreement
with previous results [19,20]. Further work is needed to verify
that this agreement persists at higher values of

Another linear MHD test was performed on a high-*n*
ballooning mode in a toroidal equilibrium with major radius 2.6
m, minor radius 0.8 m, aspect ratio 3.26, and peak toroidal
*b* ranging from 0.5% to 3%, with parallel
current profile and pressure profile shape held fixed and safety
factor on axis .
The infinite-*n* ideal ballooning criterion predicts
instability for peak *b* above 2.19%. NIMROD
simulations were performed for with 128 radial cells and 128 poloidal cells; the
time step is 1/10 of the Alfven time, and . Figure 4 shows a contour plot of the
perturbed pressure, illustrating a typical ballooning mode
structure. Below the predicted ideal stability criterion, there
is a sharp reduction in the growth rate due to a transition from
an ideally unstable ballooning mode to a slower-growing resistive
ballooning mode.

To demonstrate the code's nonlinear capabilities, we present
some preliminary work on toroidal simulations of the MHD dynamo
in RFPs. To our knowledge, toroidal geometry effects have not
been previously studied with nonlinear simulations. The effects
may have an impact on the dynamo and magnetic fluctuation level
at small aspect ratio, where poloidal coupling is strong. The
simulation presented here has and ,
with Fourier components ; pressureless MHD conditions are used. The simulation
starts from an unstable paramagnetic pinch, where is the most unstable mode.
Figure 5 shows the evolution of the reversal parameter and the
kinetic energies in the different components. The field reverses
only after nonlinear interactions broaden the spectrum. The RFP
simulations are being extended to large *S* values with more
Fourier components to achieve more realistic conditions.
Simulations in toroidal and periodic linear geometry will be
compared for a range of aspect ratios to determine the degree to
which geometry affects the MHD activity in RFPs.

Tests of NIMROD's two-fluid capabilities are being conducted
in large Hall parameter conditions, similar what was produced in
the ZT-40 RFP [21]. The geometry is a periodic linear cylinder,
equilibria have central *b* values of 0-20%, and
. There is an
ideally unstable kink mode that rotates at a frequency slightly above the
diamagnetic frequency when Hall terms are included in Ohm's law;
we have observed that the frequency increases with *b*. Quasilinear simulations with and Fourier components show much stronger interaction
between the components when the Hall terms are included. This
work is being extended to investigate two-fluid effects on the
RFP dynamo.

5. Summary

The NIMROD code is being developed as a unique computational tool for fusion plasma applications. It combines geometric flexibility with an accurate representation of resonant layers in the core plasma and a time advance algorithm that can use large time steps for studying slow phenomena. Two-fluid equations are solved in non-reduced form to provide comprehensive information of relevant physical phenomena.

Future code development efforts will focus on the following. A
coarse-grid preconditioner will be added to the existing Schwarz
domain decomposition method to improve convergence of the linear
solver. A vacuum region and resistive wall will be added.
Boundary conditions for will be implemented to allow treatment of spheromaks and
non-fusion applications. Grid packing and adaptation will be
added to improve resolution while economizing on grid size.
Anisotropic thermal conductivity and analytical closure schemes
will be added. Improved 3D visualization techniques will be
developed. Resonant particles will be added through a
*df* method to study kinetic
effects.

Acknowledgment

This work is supported by the U. S. Department of Energy.

References

[1] La Haye R J, Fitzpatrick R, Hender T C, Morris A W,
Scoville J T and Todd T N 1992 *Phys. Fluids* **B4**
2098

[2] Strait E J, Taylor T S, Turnbull A D, Ferron J R, Lao L
L, Rice B, Sauter O, Thompson S J and Wroblewski D 1995
*Phys. Rev. Lett.* **74** 2483

[3] La Haye R J *et al* 1996 *Proc. 16th Int. Conf.
on Plasma Phys. and Controlled Nucl. Fusion* (IAEA, Vienna,
1997) **1** 747

[4] Sauter O *et al* 1997 *Phys. Plasmas*
**4**, 1654

[5] Harned D S and Kerner W 1985 *J. Comput. Phys.*
**60** 62; Harned D S and Schnack D D 1986 *J. Comput.
Phys.* **65** 57; and Schnack D D, Barnes D C, Mikic Z,
Harned D S and Caramana E J 1987 *J. Comput. Phys.*
**70** 330

[6] Harned D S and Mikic Z 1989 *J. Comput. Phys.*
**83** 1

[7] Krall N A and Trivelpiece A W 1986 *Principles of
Plasma Physics* (San Francisco Press) pp 88-92

[8] Hirschman S P and Sigmar D J 1981 *Nucl. Fusion*
**9** 1079

[9] Barnes D C, Nebel R A, Sovinec C R and Nystrom W D 1996
"Three-dimensional Calculations Using the Quiet Implicit PIC
Method," presented at the *1996 Int. Sherwood Fusion Theory
Conf.* Philadelphia, Pennsylvania

[10] Park W *et al* 1992 *Phys. Fluids* **B4**
2033

[11] Schnack D D, Baxter D C and Caramana E J 1984 *J.
Comput. Phys.* **55** 485

[12] Gruber R and Rapaz 1985 *Finite Element Methods in
Linear Ideal Magnetohydrodynamics* (Berlin:
Springer-Verlag)

[13] Lerbinger K and Luciani J F 1991 *J. Comput.
Phys.* **9**7 444

[14] Lionello R, Mikic Z and Linker J A 1998 "Stability of
Algorithms for Waves with Large Flows," submitted to *J.
Comput. Phys.*

[15] Courant R, Friedrichs K O and Lewy H 1928
*Mathematische Annalen* **100** 32

[16] Marder B 1987 *J. Comput. Phys.* **68** 48

[17] Glasser A H and Sovinec C R 1997 "Numerical Analysis of
the NIMROD Formulation," presented at the *1997 Int. Sherwood
Fusion Theory Conf.* Madison, Wisconsin

[18] Smith B, Bjorstad P and Gropp W 1996 *Domain
Decomposition, Parallel Multilevel Methods for Elliptic Partial
Differential Equations* (Cambridge Univ. Press)

[19] Glasser A H, Greene J M and Johnson J L 1975 *Phys.
Fluids* **18** 875

[20] Chu M, Greene J M, Lao L L, Miller R L, Bondeson A,
Sauter O, Rice B W, Rice E J, Taylor T S and Turnbull A D 1996
*Phys. Rev. Lett.* **27** 2710

[21] Baker DA *et al* 1984 *Proc. 10th Int. Conf. on
Plasma Phys. and Controlled Nucl. Fusion* (IAEA, Vienna,
1985) 439

Figures

Figure 1. A computational grid of the D-IIID poloidal plane for NIMROD, illustrating the flux-surface conforming grid in the core and the unstructured triangulation in the edge. Discontinuities of color indicate grid block boundaries. The triangulation is coarse for purpose of illustration only; a simulation would require more refinement.

**Figure 2.** Temporal convergence of an MHD tearing mode with
kinematic viscosity equivalent to electrical diffusivity in a
periodic cylinder with close-fitting conducting shell. The grid
has 61 radial cells and 31 azimuthal cells. The vertical axis is
, and the
horizontal axis is , where *t _{A}* is
based on the minor radius.

**Figure 3.** Contour plot of the D-IIID resistive
interchange eigenfunction; the flux-surface normal velocity is
shown.

**Figure 4.** Contour plot of the ballooning mode eigenfunction; perturbed
pressure is shown.

**Figure 5.** Temporal evolution of the reversal parameter,
, and kinetic
energies in the Fourier components of the nonlinear RFP
simulation. The different components are represented from
to as cyan through magenta in
hue. Parameters are chosen such that the resistive diffusion time
is 1s.