**NIMROD
Progress Report**

**2 ^{nd}
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:

Full two-fluid Ohm's law, including ideal and resistive MHD, Hall and electron pressure gradient terms, electron inertia, and electron advection;

Anisotropic thermal conduction;

A choice of 3 neo-classical closures for the parallel viscous stress tensor;

Full MPI implementation for massively parallel computer architectures.

3. The following grids have been implemented:

Flux coordinates based on axisymmetric nested flux within a separatrix;

Doubly periodic circular cylinder;

Singly periodic circular cylinder ("soup can") geometry, with

*R*= 0 boundary (regularity) conditions;For all of the above, grid packing (non-uniform grid) in the radial (

*y*) coordinate direction;Two dimensional cartesian slab geometry.

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 (*m*_{e})
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.7x10^{6} and
Pr = 10^{4}. The two black triangles are at
*S* = 2.7x10^{5} and Pr = 10^{3}.
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.7x10^{6} and Pr = 10^{4}
and those in black are at *S* = 2.7x10^{5} and
Pr=10^{3} for a TFTR-like equilibrium based on the perturbed
pressure variant of the flow damping viscous stress tensor that uses
. Here, *m*_{e}/*n*_{e}
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.7x10^{6} and Pr=10^{4} for a TFTR-like
equilibrium based on the perturbed pressure variant of the flow
damping viscous stress tensor that uses . Here, *m*_{e}/*n*_{e}
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* = 1x10^{5}
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/W_{ci} 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 *n*^{2}. 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 *T _{jn}* can be
formed by multiplying the appropriate partial differential equation
for

. (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 {*T _{jn}*} 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 {*T _{jn}*}.

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.23x10^{8},
and D*t*=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 *b*_{N} ~ 2.
Since NIMROD presently requires a conducting wall to be at the
separatrix, the instability threshold is in the higher range of
4.5 < *b*_{N} < 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 *b*_{N}
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 *b*_{N}
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 *b*_{N} = 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 *b*_{N} = 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 *b*_{N} = 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 *b*_{N} = 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__)