NIMROD Progress Report

2nd Quarter FY 2001

Prepared by the NIMROD Team

May 11, 2001

1. Overview

We report progress for the NIMROD Project for the period January - March, 2001. Progress during this period includes:

1. A paper entitled "Limitations on the stabilization of resistive tearing modes", by Tom Gianakon, was submitted to Physics of Plasmas. The paper describes a current drive model implemented in a variant of the NIMROD code;

2. Comprehensive testing of neoclassical closures implemented in NIMROD are nearing completion;

3. Work has proceeded on the Two-fluid comparison case with M3D;

4. A procedure for creating and solving three-dimensional matrix equations, i.e. linear systems of equations with implicit coupling among different Fourier components, has been developed for the NIMROD formulation, and the required subroutines have been added to the physics kernel;

5. Anisotropic thermal conduction has been reformulated using the algorithm referenced in 4, above. Extremely accurate results have been obtained;

6. Calculations of a high-b disruption in DIII-D shot 87009 are continuing;

7. Applications of NIMROD to plasma propulsion problems are continuing.

Progress in some areas continues to be slow due to lack of manpower.

2. Status

The present version of the NIMROD code is 3.0.0. It is publicly available from http://nimrod.saic.com/download/home.html. The status of this version of the NIMROD code is as follows:

1. Arbitrary order finite elements are now available;

2. The physics kernel contains the following physical models:

3. The following grids have been implemented:

4. Interfaces with several commonly used equilibrium file formats are available.

The following development projects are underway:

1. An algorithm for a vacuum region and a moving separatrix is being implemented;

2. The implementation of kinetic and energetic particle effects by particle pushing using the gyro-kinetic model is continuing;

3. Resistive wall boundary condition is delayed due to lack of manpower

The following simulation studies are ongoing:

1. Nonlinear growth and saturation of the neo-classical tearing mode in D3D;

2. High-b disruption in D3D;

3. Two-fluid effects in RFPs and FRCs;

4. Simulation of plasma propulsion systems.

3. Highlights

Tearing Mode Stabilization Studies

Paper entitled: "Limitations on the stabilization of resistive tearing modes", by Tom Gianakon was submitted to Physics of Plasmas that describes a current drive model implemented in a variant of the NIMROD code. Stabilization is limited by the smaller of the RF current drive width or an anisotropic diffusion scale length.

Neoclassical Tearing Modes

Comprehensive testing of neoclassical closures implemented in NIMROD are nearing completion. The stability boundary can be mapped by considering multiple simulations at different initial island widths and also electron poloidal flow damping rates (me) and noting whether the simulation is unstable (solid pointing up triangle) or stable (solid pointing down triangle). Simulations have not been completed for those points that are double marked. The two prime candidate models for the electron viscous stress tensor that are based on the perturbed pressure gradient variant of electron poloidal flow damping have been completed and shown to agree reasonably well with the analytic threshold prediction based on anisotropic thermal diffusion. The two models differ in that one presumes either the flow damping acts in the grad theta direction (Figure 1) or in the direction (Figure 2). The red triangles in both figures are based on S = 2.7x106 and Pr = 104. The two black triangles are at S = 2.7x105 and Pr = 103. The lower S causes a significantly lower neoclassical threshold.

Figure 1. Neoclassical tearing mode simulations with NIMROD (triangles) reproduce the analytic anisotropic thermal diffusion threshold and saturation (solid line). The upward arrows are unstable points, the downward arrows are stable, and the double marked points are simulations underway. The simulations marked in red are at S = 2.7x106 and Pr = 104 and those in black are at S = 2.7x105 and Pr=103 for a TFTR-like equilibrium based on the perturbed pressure variant of the flow damping viscous stress tensor that uses . Here, me/ne is the ratio of the electron poloidal flow damping rate to the electron collision frequency.

Figure 2. Neoclassical tearing mode simulations with NIMROD (triangles) reproduce the analytic anisotropic thermal diffusion threshold and saturation (solid line). The upward green arrows are unstable points, the downward red arrows are stable, and the double marked points are simulations underway. The simulations are at S=2.7x106 and Pr=104 for a TFTR-like equilibrium based on the perturbed pressure variant of the flow damping viscous stress tensor that uses . Here, me/ne is the ratio of the electron poloidal flow damping rate to the electron collision frequency.

As illustrated in Figure 3, the flow damping form of the ion viscous stress-tensor has been used to demonstrate the stabilizing effect from the neoclassical enhancement of the polarization current (NEPC) on a standard resistive tearing mode. The NEPC effect adds to the effective inertia of modes with growth rates smaller or on the order of the ion flow damping. Neoclassical tearing mode threshold dynamics are presumably significantly modified by this effect. Other resistive instabilities such as ballooning modes may also be modified by this effect. [Tom Gianakon]

Figure 3. The neoclassical enhancement of the polarization current that enters through the ion viscous stress tensor acts to stabilize slow growing modes. The equilibrium is based on the 1a equilibrium used under the PSACI benchmark activity with S = 1x105 and Pr = 10-3. Here, parameterizes the magnitude of the poloidal flow damping. (The results are not fully converged.)

Two-fluid Benchmark

Work has proceeded on the Two-fluid comparison case with M3D. Some profile problems that were present in the results presented at Sherwood have been straightened out. The Figures 4 and 5 below summarize the present state of the results. Note that for the MHD case, the NIMROD growth rates are larger than M3D. The major difference between the codes for these cases is that M3D has pressure flattening in the energy equation from their treatment of effective thermal conductivity. Whether this accounts for the difference is unclear at this point.

The cases labeled H1 and H2 are two fluid results for two different Hall parameters (1/Wci in M3D units). Note that the real frequencies increase with the Hall parameter (as expected, with w ~ 2w*) and the growth rates increase with the Hall parameter as well (not expected). Convergence studies (in time-step) are presently underway to see if the semi-implicit operator is responsible for the increase in growth rate. [Rick Nebel]



Figure 4. Growth rate versus viscosity for the two-fluid PSACI benchmark case.



Figure 5. Real frequency versus Hall parameter for two different cases for the PSACI two-fluid benchmark case.



3D Matrix Equation Capabilities

A procedure for creating and solving three-dimensional matrix equations, i.e. linear systems of equations with implicit coupling among different Fourier components, has been developed for the NIMROD formulation, and the required subroutines have been added to the physics kernel. The procedure retains the full three-dimensional domain decomposition of the NIMROD algorithm by using a matrix-free approach, so parallel scaling is not affected adversely.

Like 2D matrix equations, a 3D equation in NIMROD requires a computation of the right hand side (with information from the old time-step only) for each Fourier component. Unlike the 2D equations, however, a full 3D matrix is not formed. In many cases, such a matrix would be dense in n (Fourier component index), requiring storage that increases like n2. Fortunately, Krylov-space iterative methods are based on what happens to a vector when the matrix is applied to it, and they do not require information about individual matrix elements. [Many preconditioning algorithms need information about some of the elements, however; see below.] These matrix-vector products can be formed directly from integrations over finite elements without storing individual matrix elements.

Consider the diffusivity operator for example. When the scalar field T is approximated with the NIMROD spatial representation, we have

(1)

where is the basis function that is a finite element in r and z and a Fourier basis function in j. A linear equation for Tjn can be formed by multiplying the appropriate partial differential equation for T by each of the set {}, integrating over the volume, integrating second-order terms by parts, then substituting the expansion (1) for T. Ignoring the surface integral, the diffusivity operator leads to the following term for each (k,m)

. (2)

When is independent of j, the only nonzero contribution is from n = m, so the operation is diagonal in Fourier index. When is a function of j, there are nonzero contributions from . In either case, the operator is linear as long as is not a considered a function of the unknown T. Since {Tjn} is just a set of coefficients, the order of summation and integration may be exchanged:

. (3)

Though (2) and (3) are equivalent, (2) is written as a finite element computation of a vector quantity, where the factor in parentheses is the expansion of T used in the volume integration (with a "generic_eval" of T in an "rhs" integrand). In contrast, (3) is written as a sum of the matrix elements appearing in the parentheses, with the elements of the vector {Tjn}.

For differential operators, the matrices are always sparse in k-j (coupling is local in the poloidal plane). Diagonal in n, i.e. 2D, matrices then require only a moderate amount of storage. For 2D matrix equations, it is more efficient to compute and store the matrix elements and perform matrix-vector products directly when solving the system of linear equations. With storage for 3D matrices being prohibitive as the number of Fourier components is increased, it is more practical to perform a finite element integration like (2) during each iteration of the linear solve. Since a large fraction of computation time is typically devoted to linear solves, the efficiency of performing finite element integrations is more important when using this matrix-free approach.

The new iter_3d_cg Fortran 90 module performs matrix-free conjugate gradient iterations using the same integration routines that are used to find the right-hand side vector for a linear system. When calling the 3D solver, the name of the integrand routine appropriate for matrix-vector computations like (2) is passed into the solver routine. Preconditioning is based on the diagonal in n part of the linear system. Since our preconditioning routines are algebraic, a 2D matrix that is equivalent to the diagonal in n part of the 3D matrix is required. These matrices are generated with the "matrix_create" routines used for 2D linear systems, and a separate integrand routine is required.

A disadvantage of the matrix-free approach is that it is not possible to eliminate cell-interior data (see the FY2000 3rd quarter report), as that requires an explicit computation of matrix elements. The preconditioning operations have something of a substitute for the interior elimination, but it only uses diagonal in n couplings, and it must be performed at each preconditioning step.

The first use of a 3D linear equation in NIMROD is for anisotropic diffusion, described below. Other uses arising immediately are resistive diffusion with 3D variations in the diffusivity, particularly with a distorting plasma-vacuum interface, and the velocity equation with an evolving 3D density distribution. [Carl Sovinec]

Anisotropic Diffusion

Our original formulation of anisotropic thermal conduction used a semi-implicit approach [Sovinec and Schnack, "Semi-implicit Thermal Conduction in Magnetohydrodynamic Simulations," Los Alamos report LAUR 96-4599], where toroidal variations in the diffusivity tensor appear through explicit terms only. This numerically stable approach only requires 2D matrices; however, our experience with NIMROD has shown that the temporal accuracy is unacceptable for tokamak parameters unless the time-step is many orders of magnitude smaller than that required for the rest of the algorithm. To address this shortcoming, implicit diffusion with a 3D matrix equation arising from the toroidally varying diffusivity, , has been incorporated as described in the previous section. [The operator is presently applied to pressure, since density is not being evolved. When the continuity equation is fully incorporated, the anisotropic diffusion will be applied to temperature.]

As a demonstration of the new algorithm, anisotropic diffusion is applied to the pressure profile resulting from a nonlinear simulation of the PSACI macroscopic benchmark case 1b (a tearing mode in a DIII-D-like equilibrium). The magnetic field profile is held fixed for the diffusion test, and the initial pressure profile is not constant along the flux surfaces of the magnetic island. The spatial representation is a radially packed 48x28 mesh of bicubic finite elements with 0?n?2.

With the 3D implicit operator, k?=0.423, k||=4.23x108, and Dt=1x10-4, the island structure appears immediately in the pressure profiles. Figure AD-1 shows the pressure profile and Poincaré surface of section for the magnetic field at j=0 after 15 time-steps. Qualitatively, the result is very good. Furthermore, when the test is continued for a perpendicular diffusion time, the configuration reaches a new steady state (Fig. AD-2) with the profile interior to the island reduced uniformly (Fig. AD-3). The long-time behavior indicates that artificial cross-field diffusion is not a problem, despite the large k||/k? ratio and the fact that the magnetic field is not aligned with the finite element mesh in the vicinity of the island. Overall, the numerical performance is very encouraging for tokamak simulation. Quantitative tests will also be carried out, however. [Carl Sovinec and Tom Gianakon]





Figure AD-1. Pressure profiles with a Poincaré surface of section for the magnetic field overlaid at j=0 after 15 time-steps of the anisotropic diffusion test.

Figure AD-2. "Probe" traces of pressure for n=0 (sans equilibrium, bottom trace), n=1 (middle trace), and n=2 (top trace) at (R=1.31, Z=-0.25).

Figure AD-3. Pressure profiles at j=0 along the horizontal line outward from the magnetic axis. Different traces show different times from initial conditions to the final steady state.

Testing the continuity equation: application to plasma propulsion problems

Tests on the 2D (n=0) continuity equation have been continued. A customized version of NIMROD is used featuring, in addition to the continuity equation, a provision for external loading of arbitrary plasma profiles in cylindrical (or rectangular) coordinates and a model for open-end boundary conditions.

At first, to check the basic hydrodynamic functionality, the case of a traveling plasma pulse without magnetic field was considered. A pulse with velocity along z-axis is generated by a plasma source located near z=0. In fig. 1 the transient evolution of the pulse with shock-like features is shown.



Fig 1. Density profile of plasma pulse generated by a plasma source near z=0 and propagating along the z-axis

Other runs have been performed also in the context of advanced propulsion applications, for the simulation of the plasma in the magnetic nozzle/exhaust region of the Variable Specific Impulse Magnetoplasma Rocket experiment (VASIMR [1]).

NIMROD will eventually provide 3D, MHD/two-fluid simulation capabilities to support the present analysis of the VASIMR nozzle/exhaust plasma and will be able to reproduce quantitatively the exhaust plasma profiles up to the region where the detachment will take place. .

In these simulations a plasma pulse is now drifting through a nozzle-shaped magnetic field. A 2D (r-z) section of a cylindrical domain is considered. Open-end boundary conditions are imposed along both the longitudinal and the radial direction. The source is injecting plasma in a low-density background environment and over time a plasma pulse is formed while propagating through the magnetic nozzle. In figure 2 the initial transient of the pulse formation and propagation is shown.



Fig 2. Density profile of plasma pulse building up and propagating through a VASIMR- like magnetic field (plasma source near z=0 is producing a flow with velocity along the z-axis)

Eventually, a regime is reached where the flow through the open-end boundary is balanced by the injection rate. A condition close to this steady state is illustrated in fig. 3: while the longitudinal flow maintains an essentially constant profile, some residual perturbation is still propagating outward in the radial direction. [Alfonso Tarditi]

High-b Disruption in DIII-D Shot 87009

DIII-D shot 87009 is a reversed shear discharge that exhibited a high-b disruption while undergoing neutral beam heating. The offending instability was observed to grow at a rate faster than exponential. One theoretical explanation [Callen, Hegna, Rice, Strait, and Turnbull, Phys. Plasmas 6, 2963 (1999)] is that the exponential growth rate changes with time as the plasma is driven through the ideal stability boundary, yielding a growth characterized by exp[(t/t)3/2], where t is a hybrid time scale containing the heating and incremental ideal MHD growth rates. We are attempting to verify this behavior with NIMROD, and to study the nonlinear nature of the high-b disruption. In shot 87009 instability occurs at bN ~ 2. Since NIMROD presently requires a conducting wall to be at the separatrix, the instability threshold is in the higher range of 4.5 < bN < 5.

We have verified the linear stability properties of this discharge with NIMROD. The grid used for these studies, which is based on the equilibrium flux surfaces, is shown in Figure D-1. This grid is packed around the low order rational surfaces. Bi-cubic finite elements have been used for all calculations.

Figure D-1. Computational grid based on equilibrium flux surfaces for DIII-D shot 87009. The grid is packed around the low order rational surfaces.

The linear stability properties of this discharge are illustrated in Figure D-2, where we plot the linear growth rate of the n = 1 mode as a function of bN as obtained from both NIMROD and GATO. The marginal ideal stability point as determined by DCON is also shown. Above this value, both NIMROD and GATO find an ideal mode, with NIMROD consistently exhibiting larger converged growth rates. A comparison of the NIMROD and GATO eigenfunctions is shown in Figure D-3. Below the ideal marginal point NIMROD finds a resistive interchange mode. The eigenfunctions for that mode are shown in Figure D-4.



Figure D-2. Linear growth rate versus bN for DIII-D shot 87009 with a conducting wall on the separatrix. Both NIMROD and GATO results are shown. NIUMROD finds a resistive interchange mode below the DCON ideal marginal point.

Figure D-3. Comparison of linear eigenfunctions for the n = 1 ideal MHD instability of DIII-D shot 87009 at bN = 5. Poloidal flow vectors are shown. The agreement is quite good.

Figure D-4. Poloidal flow vectors (left) and pressure perturbed pressure contours (right) for the n = 1 resistive interchange mode in DIII-D shot 87009 at bN = 4.

Our goal is to simulate the heating and destabilization of the discharge beginning below the ideal stability limit. As a preliminary calculation, we have computed the nonlinear evolution of the ideal n = 1 instability at bN = 5. The evolution of the pressure contours for this case is shown in Figure D-5. A true disruption is not observed, and the mode saturates in the presence of a conducting wall. However, the filamentary structure of the toroidal current (not shown) could lead to reconnection and enhanced thermal losses along wandering field lines.

t = 0 msec. t = 0.120 msec.

t = 0.131 msec. t = 0.135 msec.

Figure D-5. Nonlinear evolution of the pressure surfaces in DIII-D shot 87009 at bN = 5. The modes saturates with a conducting wall on the separatrix.

Calculations with heating beginning below the ideal marginal stability point are presently underway. [Dalton Schnack]



[1] VASIMR: http://spaceflight.nasa.gov/mars/technology/propulsion/aspl)

13