NIMROD Progress Report

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

1. High-order finite elements have been generalized, and extended to bi-cubic order;

2. Problems with the 2-fluid calculations have been traced to the form of the semi-implicit operator;

3. A campaign to validate NIMROD on linear and nonlinear cases at realistic values of S has computed linear and nonlinear tearing modes in D3D geometry at S = 10⁵;

4. Improvements to the NIMROD web site have led to increased communication between team members.

5. Progress has been made in implementing neo-classical closures into the physics model, and in applying them to neo-classical tearing modes (NTMs);

6. NIMROD is being applied to modeling high-b disruptions in D3D.

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. Generalization of the finite element representation to include higher order (e.g., quadratic) basis functions;

2. Algorithms to improve code performance in the presence of equilibrium flow are being developed;

3. Calculations on triangular grids are being validated;

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

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

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

Numerous neoclassical closures have been implemented and tested in the NIMROD code. They fall into three broad categories: 1) flow and current, where one assumes a Chew-Goldberger-Low form for the stress-tensor and implements a neoclassical form for the pressure anisotropy; 2) flow damping, where electron and ion flows in the poloidal direction are explicitly damped; and 3) perturbed pressure models where an ideal MHD approximation is made to write the poloidal current that appears in the first two models as a function of the pressure gradient. The first model is generally unable to reproduce poloidal flow damping of ions and generates unphysical toroidal angular momentum. Therefore, it has not been tested in the context of perturbed bootstrap currents. The second model successfully reproduces ion poloidal flow damping, but the generation of perturbed bootstrap currents has not been demonstrated and requires more extensive convergence studies. Preliminary results with the second model exhibit significant pressure profile degradation, attributed to insufficient spatial resolution of anisotropic heat flux when the toroidally symmetric part of the magnetic field loses alignment with the grid. The third closure form is used only for bootstrap current generation. A numerical instability has been identified with this form, which requires further investigation.

A mini-workshop on NTMs with scientists from the NIMROD Team, the University of Wisconsin, and from the Princeton Plasma Physics Laboratory was held during the 2000 International Sherwood Fusion Theory Conference. The issue of which vector projections of the damping force survive flux surface averaging was discussed, providing guidance on the suitability of the different forms of the neoclassical closure. [Tom Gianakon]

In an effort to better elucidate some of our two-fluid simulation problems, a 1D numerical stability analysis (taking k parallel to B₀) was performed on the two fluid algorithm including the semi-implicit operator. Shear Alfvén waves, ion cyclotron waves, whistler waves, and electron cyclotron waves were all considered. The analysis indicates that the algorithm is numerically stable for all of these waves. At large time step, stability for the whistler and electron cyclotron waves is achieved by raising the effective mass of the electrons, essentially following the plasma on the effective electron cyclotron motion time-scale. The resulting effective electron mass may be comparable to the ion mass.

RFP-like linear cases previously showing an unphysical high-wavenumber mode when run with the two-fluid Ohm's law and bilinear finite elements have been investigated with a temporal convergence study using biquadratic finite elements. At large time step (Dt~0.01/g), using quadratic elements does not remove the unphysical mode. At a greatly reduced time step (Dt~0.0001/g), a rotating mode similar to the resistive MHD mode is recovered. (See Fig. 1)

Since quadratic elements (with different levels of magnetic divergence diffusivity) do not prevent the unphysical mode, but a reduced time step does, we conclude that our semi-implicit Hall advance leads to the undesirable behavior. The existing semi-implicit operator for Hall physics has the form,

which acts on

. This operator has the same form as electron inertia, which appears as the operator

which also acts on ¶B/¶t. Our Hall semi-implicit operator therefore acts like electron inertia with an effective electron mass of eB₀Dt.

Fig. 1. Comparison of eigenfunctions for linear two-fluid and resistive MHD stability computations for an RFP-like equilibrium. The different colors represent different phases for the same eigenfunction. The small time step, 2-fluid case shows the influence of differential electron flow, unlike the MHD case.

Electron inertia gives rise to magnetic reconnection if the electron skin depth, c/w_pe, is comparable to the resistive tearing layer width. If this depth, computed with the effective electron mass, is large, we can expect an unphysical amount of reconnection. For the computation displayed above, the effective electron skin depth exceeds 10% of the minor radius for the large time step. This is larger than the expected tearing layer width for S=10⁴, which is used in this computation. When the time step is reduced two orders of magnitude, the effective electron skin depth decreases one order of magnitude, and a smooth eigenfunction is recovered. We note that other computations successfully use the 2-fluid advance at large time step, but those cases compute non-resonant ideal modes or have S = 10³, so they are less sensitive to the effective electron skin depth.

Although temporal convergence with the existing algorithm seems to produce accurate results, the convergence properties are not acceptable. As S is increased (by decreasing dissipation) the physical mode slows, and an optimal algorithm allows larger time steps. In contrast, the existing 2-fluid algorithm requires us to decrease the time step with S to keep the effective electron skin depth below the tearing layer width (until the tearing layer width is smaller than the actual electron skin depth). We will implement alternative semi-implicit operators, such as the one described by Harned and Mikic [JCP 83, 1 (1989)] to improve the 2-fluid convergence properties of NIMROD. [Rick Nebel and Carl Sovinec]

The local and global line Jacobi preconditioning schemes for regions of quadrilaterals have been updated to use data centered along element sides and within element interiors. This preconditioning now appears to be as effective for simulations with higher-order elements as it is for bilinear elements.

Definitions for cubic elements have been added for quadrilateral cells. The generalized version of NIMROD can therefore use bilinear, biquadratic, or bicubic finite elements for dependent quantities, and a user may specify the type of element with a single parameter in the input file. We have verified that bicubic elements further improve spatial convergence on a tearing mode, as shown in Figs. 2 and 3.

Fig. 2. Convergence on growth rate for a tearing mode using different types of quadrilateral finite elements. The inverse of the number of cells used in both the radial and azimuthal directions is indicated by the horizontal axis.

Fig. 3 Magnetic divergence error, , vs. resolution, plotted on a log-log scale. The slope of each line shows the convergence order for this error.

The generalized version of NIMROD will be fully functional for simulations not requiring triangular elements during the next quarter. It will be made available to the team for testing and evaluation. [Carl Sovinec]

The NIMROD code is being applied to the simulation of high-beta D3-D discharges to investigate the growth of ideal MHD modes driven through their instability threshold. The D3-D shot #87009 was chosen as basic reference and different beta scenarios where considered, from a stable region through marginal stability and beyond it. Results from the GATO code indicate the onset of unstable behavior (internal kink modes) for a normalized beta b_N of about 2. A comparison with NIMROD linear results was reported in the 1^st Quarter 2000 Progress Report.

We have begun to study this case nonlinearly. The initial conditions correspond to b_N < 2, so they are stable to the kink mode. Heating is then applied to slowly increase the pressure profile, thus driving the discharge through marginal stability. A sequence of pressure profiles at various times during the simulation is shown in Fig. 4.

Fig. 4. The n=0 pressure profile at various times during the nonlinear simulation of D3D shot #87009. The heating model increases the pressure profile. The lowest profile is just below the marginal stability point.

A mode with spatial structure similar to the linear eigenfunction reported in the 1^st Quarter 2000 Progress Report begins to grow. The magnetic energy in the n = 1 component as a function of time is shown in Fig. 5.

Fig. 5. Magnetic energy in the n = 1 Fourier mode as a function of time after the initiation of heating. The original profile is stable. The heating drives the discharge through the marginal stability point.

Since the profiles are evolving dynamically in response to the heating, the growth rate of the mode will not be steady, but rather will be a function of time. Analytic predictions [J. D. Callen, C. C. Hegna, B. W. Rice, E. J. Strait, and A. D. Turnbull. Physics of Plasmas, 6, 2963 (1999)] indicate that the mode should grow as exp(at^3/2), rather than as exp(gt), which implies that the instantaneous growth rate of the mode should evolve as g(t) ~ t^1/2. The instantaneous growth rate of the n=1 mode shown in Fig. 5 is plotted in Fig. 6 as a function of time. We have not yet determined the time dependence of this function.

The eigenfunction corresponding to this case shows signs of loss of spatial resolution. We therefore have good reason to believe that the above results are numerical rather than physical. We are continuing to study this case with finer grids. [Alfonso Tarditi]

The high-S campaign has progressed to studies of tearing modes in toroidal geometry with non-circular cross section. The case under consideration is D3D shot #86144 at 2250 msec. In the experiment a 3/2 magnetic island appears at about this time. This island eventually evolves into a neo-classical tearing mode.

Shot #86144 has q(0) < 1 and q(a) about 4.. Since we are using resistive MHD for the high-S campaign, this profile would be unstable to an n=1 mode. This is not seen in the experiment; presumably it is stabilized by non-MHD effects. In order to avoid this mode in the MHD study, we have artificially raised the q-profile such that q(0) = 1.05. The equilibrium flux surfaces are shown in Figure 7.

We have applied NIMROD to study the n = 2 resistive stability of this case, initially at S = 10⁵. Using grid packing in the outer regions of the plasma, we find an exponentially growing n = 2 eigenfunction. The spatial structure of the poloidal velocity, the toroidal velocity, the poloidal magnetic field, the toroidal magnetic field, and the pressure for the unstable mode are shown in Figures 8-12, respectively. A careful examination of these figures reveals structures indicating resonances at the 3/2, 4/2 and 5/2 surfaces.

This case was then run nonlinearly at S = 10⁵ with the n = 0,1, and 2 modes. The kinetic energy for each mode as a function of time is shown in Figure 13. The mode is seen to grow and saturate nonlinearly. The "blip" at about t = 4.2 ¥ 10^-3 sec. is the result of a restart with a faulty time step. It takes a few time steps to readjust and recover the previous growth rate. The increasing growth rate late in the run is due to time step reduction in order to satisfy the advective CFL stability condition. (The initial time step corresponded to a linear wave CFL number exceeding 50,000.), and may indicate the effect of the semi-implicit term on the growth rate. The mode saturates benignly with a 3/2 magnetic island that is too small to make an interesting surface of section. However, it is possible (but not shown) that it is of sufficient size to act as the seed island for a neo-classical tearing mode.

Fig. 13. Kinetic energy in the toroidal modes (yellow: n = 0; light blue: n = 1; magenta: n = 2) for nonlinear simulation of D3D shot #86144 at 2250 msec.

Since this case does not have robust nonlinear behavior, it is not a good candidate for testing NIMROD's nonlinear behavior at higher S. To obtain more interesting dynamics we have lowered q(0) to it's original value (~ 0.96) in order to destabilize the robust 1/1 resistive mode. These calculations are presently underway. [Dalton Schnack]

We are in the process of incorporating energetic particle effects into the NIMROD nonlinear MHD simulation. We have already implemented integrating the particle equations of motion in the NIMROD magnetic fields. This is substantial progress because we needed to formulate how to include particles in an unstructured-grid finite element representation of the fields. The energetic particles will couple through the pressure tensor in the MHD momentum equation.[Reference: W. Park, et al., Phys. Fluids B 4, 2033 (1992)] Our target physics problem is energetic particle stabilization of the internal kink mode in DIII-D.

NIMROD is an MHD simulation code that uses the finite element method to advance the field equations. The use of finite elements allows for quite general geometries, including divertor regions and complex vacuum chambers. However finite elements presents a challenge for the implementation of particles. The grid is highly nonuniform and/or unstructured and the particle gather and scatter becomes nontrivial. There is also the issue of is tracking each particle as it moves from one finite element to the next. We have extended the ideas for particles in a triangular finite element grid [Reference: F.Kazeminezhad, S. Zalesak, D. Spicer, ``A particle model on an unstructured mesh,'' Computer Physics Communications, 90, 267 (1995)] to the quadrilateral finite elements used in NIMROD. We have formulated the shape functions which determines the weighting of particles onto the nodes. The shape functions have also been suggested as a means of tracking the particles, cite Ref. 2. Kazeminezhad suggests the use of the shape functions to search for the particles. However, after some study, the tracking algorithm of Kazeminezhad was found to be computationally inefficient. The present tracking algorithm uses a uniform grid to deduce which finite element the particle is located in. Once implemented, it is our intention to study energetic particle effects on internal kink modes. [Scott Parker and Charlson Kim]

The goal of this work is to implement a model for the vacuum region outside the separatrix that:

The solution of equations in regions of disparate material properties is a classic computationally stiff problem, and the choices that are made in this work represent a balance of the three goals. For example, eigenvalue codes such as GATO and MARS use the Green's function approach to give a response matrix for the vacuum region, but this method is not easily generalized to nonlinear simulations where the boundaries can move.

The problem is to de-couple the fluid motion from the magnetic field in the vacuum region without paying an excessive computational penalty. The model that NIMROD uses is to increase the resistivity in the vacuum region as compared to the core plasma region. The vacuum region is distinguished from the core region by the use of a shape variable that controls the spatial distribution of the resistivity. The use of the shape variable, along with the advection of a "shape variable" which will map the plasma-vacuum interface.

The method has been shown to give the correct behavior in the linear regime. The advection of the concentration variable and the corresponding resistivity mapping has been implemented, and is being tested for robustness. The next step is to incorporate the three-dimensional resistivity into the Ohm's Law along with a suitable semi-implicit operator.

Much of the development on the pre-processors fluxgrid and nimset have been related to the developments of the vacuum region, specifically in preparation for having a GATO-NIMROD benchmark with a vacuum region.

NIMROD has two grid-types: structured rectangular grids and unstructured triangular grids. The structured grids give better performance so we want to maximize their use, but previously they were only defined inside the separatrix. The method that fluxgrid uses to extend the rblocks into the vacuum region is based on the same algorithms that GATO uses for the self-similar walls and conformal walls. This allows for easy comparison of NIMROD and GATO for these simple wall shapes. To run NIMROD with the full DIII-D wall geometry, both r-blocks and t-blocks must be used. Figure 14 shows a grid with a self-similar wall and r-blocks extended beyond the separatrix. Details of this grid near the separatrix are shown in Figure 15. A sample grid with t-locks used to capture the wall geometry is shown in Figure 16. [Scott Kruger]

users@nimrodteam.org.

This mailing list for the NIMROD User's Group is to aid in using the code and discuss technical issues related to NIMROD.

This mailing list for general NIMROD announcements for those interested in the code, but don't wish to hear the technical details found in the User's Group.

People are able to subscribe, unsubscribe, see past mailings, and view the subscription list from the web. The NIMROD User's group held its first public meeting at the Sherwood Theory Conference. The next User's Group meeting will be held at the APS meeting.

The NIMROD documentation web page has been reorganized and improved. Most of the improvements are on the backend to make it easier for all of the developer's to contribute to the documentation. [Scott Kruger]

Prepared by the NIMROD Team