NIMROD Progress Report

We report progress for the NIMROD Project for the period April-June, 2000.. Progress during this period includes:

1. A thin resistive wall boundary condition bounded by a vacuum region has been implemented in the bilinear element version of NIMROD for periodic slab and cylinder cases;

2. Changes for the higher-order elements have been merged into the main line NIMROD code repository, and the combined code has been made available to the Team for testing and evaluation;

3. Work has started on the M3D-NIMROD comparisons for two-fluid simulations;

4. The FLUXGRID pre-processor has been rewritten to improve this accuracy as well as allowing for more flexible grids, including packing in the theta direction;

5. The NIMROD web site has been modified to improve its functionality and appearance;

6. NIMROD has successfully computed the nonlinear resistive MHD evolution of D3D discharge 86144 at 2250 msec. with S = 10⁷, albeit with relatively large viscosity (Pr = 10³).

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

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.

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.

1. Beta-testing of the high-order finite elements;

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

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

4. Resistive wall boundary condition is being implemented

1. Linear stability studies for the UCLA Electric Tokamak (ET);

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

3. High-b disruption in D3D;

4. Relaxation in spheromaks;

5. Two-fluid effects in RFPs and FRCs;

6. High-S campaign.

A thin resistive wall boundary condition bounded by a vacuum region has been implemented in the bilinear element version of NIMROD for periodic slab and cylinder cases. The implementation extends to the toroidal case once an appropriate response matrix becomes available through Morrell Chance's vacuum code [M. S. Chance Phys. Plasmas 4 (1997), 2161]. The thin resistive wall approximation presumes that the normal component of magnetic field is constant through the wall. This is used to evaluate the vacuum solution outside the wall and provides the normal component of the normal gradient of the magnetic field outside the resistive wall. The plasma response provides a similar gradient on the inside of the wall, and the difference is used to advance the normal component in the wall. The updated normal component is then applied as a boundary condition to the NIMROD solution and as such is explicitly applied in the time-step.

Our initial set of test simulations have focused on the periodic slab geometry with a tearing unstable equilibrium [J. E. Finn, Phys. Plasmas 5 (1998), 461]. The equilibrium current is allowed to flow in the periodic direction or in the finite element plane so that convergence properties of the finite element representation and the Fourier representation can be studied. In zero-flow cases that compare results from a resistive wall versus a conducting wall, the island widths agree with the reference when the mode is oriented in the Fourier direction. When the mode is oriented in the finite-element direction, the island width tends to be larger and distorted in shape. The cause may be a lack of poloidal resolution in our test.

Fig. 2. Magnetic flux surfaces for the saturated tearing mode in a sheared slab with a conducting wall on the right side (no flow).

The introduction of flow across the island in the finite element direction has proven problematic with a possible cause identified as the synergistic

/flow instability. Therefore, our first efforts with flow have focused on the Fourier representation, which should exhibit a weaker

problem. In this orientation with slow flow velocities, locking is observed. At high flow velocities, the time-step had to be reduced to approximately 1/100 of an Alfvén time (t_A) such that the flow CFL is around 0.03. In one case the equilibrium flow is given a magnitude such that

. The simulation exhibits a slow decay of the total flow as radial magnetic field slowly diffuses through the wall, and translation continues for three resistive times (S=10⁵). The net flow eventually reaches a value of

, at which point rapid locking occurs, and the island grows to the zero-flow resistive wall result. The

diagnostic drops at least an order of magnitude between the unlocked and locked states, suggesting that the synergistic instability may be playing a role prior to locking. [Tom Gianakon]

Fig. 3. Evolution of the radial component of magnetic field at a position near the non-ideal wall as locking forces slow the net flow, and the magnetic island becomes stationary in the laboratory frame.

During this quarter, changes for the higher-order elements have been merged into the main line NIMROD code repository, and the combined code has been made available to the Team for testing and evaluation. However, prior to this limited distribution, two important changes were made, and nonlinear tests were performed.

The basis function definitions in quadrilaterals now use general formulas for computing the value of a basis function at an arbitrary point in its domain. In the 1D domain [0,1], the n-th degree Lagrange polynomial has nodes at

Since each Lagrange basis function is nonzero at only one of the nodes in the domain, they must have the form

It is extremely unlikely that this is a new result, but Team members had not encountered it in the literature. For 2D quadrilaterals, the n² basis functions are

These general definitions have been implemented in NIMROD, along with general Gaussian quadrature integration for the required level of accuracy, so that the polynomial degree in quadrilaterals is now unrestricted. In principle spatial convergence may be achieved through grid resolution, polynomial degree resolution, or both.

The second change is the elimination of interior-centered data prior to calling the iterative matrix equation solver. Coupling to the interior nodes is always local across a grid cell for Lagrange elements. Therefore, a direct elimination of this data only affects local connections that are nonzero in the original matrix. The amount of data that is eliminated scales as

in grids of quadrilaterals, where n is the polynomial degree. As n is increased, the iterative solver does less work, helping overall efficiency.

Final tests before the limited distribution included nonlinear simulations. An RFP computation was performed with bicubic elements, and convergence on a spheromak result was obtained by running a bilinear computation with a 40x40 grid, then a biquadratic computation with a 30x30 grid, and finally a bicubic computation with a 30x30 grid. [Carl Sovinec]

Work has started on the M3D-NIMROD comparisons for two-fluid simulations. We have received a test case from Linda Sugiyama. It is an ideal m=1 kink mode for a tokamak equilibrium in a periodic cylinder. The M3D profiles have been incorporated in NIMROD, and we have started the simulation comparisons. To date we have only compared MHD results to serve as a guide before attempting two-fluid simulations.

Initial results are very encouraging. NIMROD achieved growth rates within a few percent of the M3D results almost immediately. Figure 4 shows typical linear eigenfunctions for the radial flow. Note that the mode has the ideal parity with B_r vanishing at the m=1 resonance (located at r/a~.4), and V_r doesn't change sign at the resonance. This is consistent with an ideal mode and is as expected. Convergence studies with viscosity, time-step, and mesh size are now underway. [Rick Nebel]

It had been observed that our pre-processor fluxgrid had difficulty maintaining accuracy of the Grad-Shafranov equation, with Grad-Shafranov equation residuals routinely exceeding unity near the edge of the plasma. The pre-processor has been rewritten to improve this accuracy as well as allowing for more flexible grids, including packing in the theta direction. Various minor improvements, such as making it compatible with the new higher-order element version of nimrod, have been done also. [Scott Kruger]

A summer intern has been hired at SAIC to maintain the web site and improve it's appearance. Various improvements have already been done to improve it's organization and readability. Any suggestions for improvement are greatly welcome and may be sent webmaster@nimrodteam.org. [Scott Kruger]

The High-S Campaign has previously demonstrated accurate scaling of linear tearing mode growth rates at values of the Lundquist number S up to 10⁹. Recently the High-S Campaign has been concentrating on the much more difficult problem of the nonlinear evolution of an equilibrium corresponding to D3D shot 86144 at 2250 msec. This case has been previously modeled at low-S [Ref: Gianakon, et al, http://www.nimrodteam.org/presentations/Sherwood99/sherwood99_gianakon.ps] Shot 86144 displays continuous n = 1 sawtooth oscillations. At 2250 msec. a 3/2 magnetic island appears after a sawtooth crash. The size of this island is amplified by successive sawtooth crashes until it reaches a saturated level. This perturbation has been interpreted as a neo-classical tearing mode triggered by a sawtooth crash. The discharge eventually terminates when a large 2./1 island appears and locks to the wall. Successful modeling of this shot is extremely challenging, and requires:

The present campaign concentrates on the first point, above, and aims at computing the nonlinear resistive MHD evolution of this discharge at high S.

We begin with an equilibrium determined from an efit reconstruction of the discharge at 2250 msec. A flux-based grid computed from this equilibrium is shown in Fig. 5. The equilibrium has q(0) = 0.98, and q_sep = 4.1. The grid has 128 radial (y) points and 64 angular (q) points. Three Fourier modes (n = 0, 1, and 2) are included in the simulation. The grid is packed non-uniformly about the q = 1, 3/2,and 2 surfaces.

Fig. 5. Poloidal grid (128 ¥ 64) used for initial high-S studies of D3D shot 86144 at 2250 msec. The grid is packed about the q = 1, 3/2, and 2 surfaces. Different colored regions indicate grid decomposition for parallel computing (here 6 blocks).

The equilibrium is perturbed at low amplitude (dB ~ 10^-9 Tesla) and allowed to evolve in time. Calculations have been made with increasing S, but with fixed viscosity. For S = 10⁵, 10⁶, and 10⁷, this corresponds to magnetic Prandtl numbers of Pr = 10, 10², and 10³, respectively.

The kinetic energy for the three Fourier modes [n = 0 (yellow), n = 1 (blue), and n = 2 (magenta)] as a function of time for the case S = 10⁷, Pr = 10³ is shown in Fig.6. Both the n = 1 and n = 2 modes are linearly unstable (g_n=1 = 9.883 ¥ 10³ sec^-1, g_n=2 = 6.155 sec^-1). The linear eigenmode structure is shown in Figures 7 and 8, where we display vector plots of the poloidal components of the velocity for each mode.. The n = 1 is dominantly 1/1, with a small 2/1 component, and the n = 2 is dominantly 2/2, with a small 3/2 component.

Fig. 6. Kinetic energy vs. time for three Fourier modes [n = 0 (yellow), n = 1 (blue), and n = 2 (magenta)] for the case S = 10⁷, Pr = 10³.

Fig. 7. Poloidal velocity for the n = 1 linear eigenfunction in the core of the discharge. Note the dominant 1/1 structure, and the small 2/1 "halo".

Fig. 8. Poloidal velocity for the n = 2 linear eigenfunction in the core of the discharge. Note the dominant 2/1 structure, and the small 3/2 "halo".

Fig. 9. Magnetic field topology in saturated state for D3D shot 86144.2250 with S = 10⁷, Pr = 10³. Note the 1/1, 2/2, and 2/1 structures, and the absence of 3/2 magnetic islands.

At t ~ 2 msec. the stable n = 0 mode becomes nonlinearly driven by the n = 1 mode. At t ~ 2.2 msec. the n = 2 mode makes a similar transition. The driven character of these modes is indicated by their growth rates becoming twice that of the n = 1. The growth of all modes saturates at t ~ 3.1 msec.

A Poincaré plot showing the magnetic field topology at t = 4 msec. is shown in Fig. 9. The 1/1 and 2/2 distortions in the center and the 2/1 magnetic island are clearly visible. Note the absence of a 3/2 magnetic island for this case, although a moderately sized 3/2 structure was seen with the less realistic parameters S = 10⁵ and Pr = 10.

NIMROD has successfully computed the nonlinear resistive MHD evolution of D3D discharge 86144 at 2250 msec. with S = 10⁷, albeit with relatively large viscosity (Pr = 10³). The driven 3/2 magnetic island that was seen in the experiment (and which presumably serves at the seed island for a neo-classical tearing mode) is not seen in this simulation, although it does appear at lower values of S. Calculations are underway at S = 10⁷ with lower viscosity (Pr = 10). Preliminary indications are that a substantially finer grid will be required for these cases. [Dalton Schnack]

Prepared by the NIMROD Team