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.

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 r is the mass density [7]. In its present form, NIMROD solves all but the continuity equation; however, the current advection terms in Ohm’s law and all nonadiabatic terms in the pressure equations have not been implemented, and the stress tensor term in the momentum equation is presently an isotropic kinematic viscosity. The Hall, electron inertia, and stress tensor terms may be omitted in a simulation to solve the resistive MHD equations.

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].

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 f_i is an isotropic operator used to numerically stabilize nonlinear interactions (f_i<<f_a~1) [13]. The advection term appearing here and in other equations is addressed with a predictor/corrector scheme, where all terms are used in both predictor and corrector steps [13,14]. This treatment of advection requires limiting time steps to the Courant-Friedrichs-Lewy condition [15]. Satisfying this condition is usually not too restrictive in fusion plasmas, where velocities are much less than the sound and Alfven speeds. However, some type of implicit advection may be required for strong flow conditions. The Lorentz and pressure forces on the right side of Equation (1) are evaluated at the old time for both prediction and correction.

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 t_A is the Alfven time).

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.

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 S, as the region of radial localization shrinks. The results shown here use a radial grid that is uniform in r. Future work will use grid packing in the neighborhood of the singular surfaces to improve resolution at higher values of S, and we will explore the nonlinear evolution of the mode.

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.

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.

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

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.