NIMROD Conferences and Presentations

We report progress for the NIMROD Project for the period October-December 1999. Progress during this period includes:

1. Second order finite element basis functions have demonstrated quadratic convergence and greatly improved accuracy for smaller numbers of mesh points;

2. Pressure equilibration by parallel thermal conduction on magnetic surfaces associated with islands has been demonstrated in realistic toroidal geometry. Heat flux anisotopy approaching 10⁸ has been demonstrated;

3. Two-fluid simulations of the RFP are continuing;

4. A campaign to validate NIMROD on linear and nonlinear cases at realistic values of S has demonstrated accurate scaling of the linear 2/1 resistive tearing mode growth rate over the range ;

5. Work on equilibrium flows is continuing;

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

7. Progress has been made in implementing neo-classical closures into the physics model;

8. NIMROD has been applied to problems involving spheromak formation and sustainment in gun geometry;

9. NIMROD is being applied to the nonlinear tilting mode in FRCs;

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

The NIMROD project continues to suffer from lack of manpower. During the present reporting period, progress in implementing both the resistive wall and the vacuum region was halted because the responsible scientists were fully involved in other important projects. (Tom Gianakon: modeling of the neoclassical tearing mode; Scott Kruger: validating improved reduced MHD model for use at GA.) The team's energy was also diverted toward preparing the PSACI proposal during the month of December.

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;

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.

Pressure equilibration on magnetic surfaces associated with islands has been successfully demonstrated in both slab geometry and non-circular toroidal geometry for values of

exceeding 10⁸. This was achieved by using increased grid resolution (as well as finding a bug!) This has allowed further progress in neoclassical tearing modes (see below). [Tom Gianakon]

Two variants of the flow and current form of the neoclassical closure have been added to NIMROD. The first drops the flow effects in the evaluation of the pressure anisotropy. The second closure directly relates the current density used in the pressure anisotropy computation to the perturbed pressure gradient. With these additions, four neoclassical closures are now functional in NIMROD: flow and current, current only, perturbed pressure, and hole.

The motivation for the current-only form was the observation that initial parameters for the flow and current closure result in a condition where the flow term is comparable to the current term in the evaluation of the pressure anisotropy. In SI units this occurs when

, where n is the number density and e is the charge of an electron. It may be possible that this term plays a significant role in neoclassical threshold dynamics, so this could be a rather exciting observation. The perturbed current form puts the closure into a form that looks like a bootstrap current and is similar to forms previously used in reduced MHD codes.

A neoclassical closure for the momentum equations has also been added to NIMROD and is now being tested. We have considered three small amplitude initial conditions for the plasma flow: 1) a pure toroidal flow, 2) a pure parallel flow, and 3) a pure poloidal flow, where all tests have n=0 only. They include parallel thermal diffusion, are based on the resistive MHD equations, and have no kinematic viscosity. The pressure anisotropy is of the form

Since the initial flow perturbation is small, we observe that the current remains small, so the second term on the right can be ignored. Additionally, since the simulation has n=0 only, toroidal flows do not affect the pressure anisotropy. This has been verified in the first test, as the flows do not change with time. The second case is important, because the parallel projection of the momentum equation leads to a term that always produces damping in the parallel direction. The test result reproduces this damping, as shown in Figure 1, and the damping is at the expected rate.

The third test, where the initial condition is a pure poloidal flow, is more complicated. The kinetic energy is observed to increase, as illustrated in the Figure 2, and the flow is dominated by the generation of a toroidal flow.

To understand this result the lowest order momentum equation is approximated by the rather simple equation

where N is a vector with projections in all three directions. This implies that the presence of poloidal flows will cause some effect in the toroidal direction, producing a flow that is neither damped nor driven, and the lack of a continuity equation may lead to the kinetic energy error. We suspect that after the poloidal flows have damped away, a steady state toroidal flow will remain. This test case is being continued to determine if such a situation develops.

The simulation reported at the APS-DPP meeting in Seattle, which uses the flow and current closure model, is not an NTM, and recent simulations with the same equilibrium confirm that it is simply a resistive tearing mode. We have analyzed this case with new diagnostics that have been added to NIMPLOT. These diagnostics include:

Using these diagnostics we have determined that the perpendicular velocities are almost exactly

, and the parallel electric field is dominated by the resistivity. The contribution from the bootstrap current does appear and looks to have the right sign for an NTM, but it is about a factor 20 times smaller than other contributions in Ohm's law. The large poloidal current observed in this case can be attributed to a Grad-Shafranov effect. The pressure flattening in the vicinity of the island produces a net decrease in pressure between the island resonant surface and the magnetic axis, implying

near the island. (B₀ is mainly toroidal, and

is mainly radial.) This matches well the observed poloidal current. A new series of NTM investigations has begun from a new equilibrium that should be stable to tearing modes. [Tom Gianakon]

Previously, we discovered that Hall and diamagnetic terms can induce large shear flows in RFP mean fields. These large flows are caused by the nonlinear evolution of the m=1 tearing modes, which are also responsible for the RFP dynamo effect. The flow arises quasilinearly in the Hall MHD model, and its magnitude is significantly larger than the diamagnetic drift speed.

At present we do not have a self-consistent model for the two-fluid effects, however, since gyroviscosity has not been included. We have formed a study group at Los Alamos in order to get a better understanding of what these terms are, how they interact with the Hall and diamagnetic terms, and how to best incorporate them numerically in NIMROD.

We have also discovered a "spurious " mode in the two fluid simulations. It comes in at about 1/10 the wavelength of the standard tearing mode. Although we do not understand its origin at present, we have been able to eliminate several possibilities. It is not a numerical instability, a two-fluid remnant of the g-mode, or a spurious mode caused by the formulation of the semi-implicit operator. We are presently looking at whether it might be caused by

problems. [Rick Nebel]

Most of the physics kernel of NIMROD has been converted to use standard Lagrange-type basis functions of arbitrary degree. More work needs to be done on the conjugate-gradient preconditioners, regularity conditions, blocks of triangular elements, and on time step limitations. Basis function definitions have been coded for polynomials up to order 2. We can therefore run NIMROD on test cases to examine spatial convergence with biquadratic elements. This has been done on a tearing mode in linear, circular cross-section geometry. In short, the convergence properties are very encouraging. Accurate answers are achieved with much smaller meshes when biquadratic elements are used. Furthermore, quadratic convergence of the

error is borne out.

The biquadratic finite elements have greatly improved the convergence properties of the NIMROD physics kernel. These are illustrated in Figures 3 and 4. Figure 3 shows the convergence of the growth rate of resistive tearing mode as the number of cells is varies, and Figure 4 shows the similar convergence of the error in

. We expect these algorithmic improvements to be incorporated in the public release of NIMROD in the next few months. [Carl Sovinec]

Without any imposed equilibrium poloidal rotation, a 1/1 resistive linear kink mode case was used to benchmark NIMROD with the reduced MHD version of FAR . A strong consistency was found in the growth rate relative to S (magnetic Reynolds number) between these codes. These growth rates shown in Figure 5 also agree very well with the Full FAR code results and with theoretical predictions by Hastie et al., for this mode.

We then used reduced FAR to explore the effects of poloidal rotation on this mode. Preliminary results suggest that poloidal rotation profile effects can fall into three classes of solutions depending on where the maximum in that rotation is relative to the mode rational surface. The growth rate can easily increase from rotation if the profile peak is inside the q=1 surface. In this case, the rotation is likely just supplying centripetal energy which is effectively weakening the poloidal magnetic field. The result is an effective beta poloidal:

. Since the growth rate has been shown to increase as beta increases, we are likely seeing this effect. Here poloidal rotation would cause an effective total beta:

. Poloidal rotation should therefore be much more effective at creating effective pressure than toroidal rotation. The growth rate increase shown below is consistent with this picture. Another case is if the poloidal rotation peak is near the plasma edge or beyond. This case approximates rigid rotation and the poloidal rotation effectively drops out of the equations. Rigid rotation has therefore no effect on mode growth as expected. The way to decrease the growth rate while staying below the poloidal unity Mach number is to have the rotation peak very near the very core so that its shear effect can overcome the centripetal energy effect. The shear should start to stabilize the mode at

. Here we are concerned with shear in the poloidal angular rotation at the q=1 surface. Again the reduced FAR simulations are roughly consistent with this picture.

As we did in the above figure, an important milestone for NIMROD will be to reproduce the interesting physics described in Figure 6. [Jean-Noel Leboeuf, Michael Kissick]

The NIMROD code has been 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. NIMROD is being run first to compare the linear growth rates with GATO, with the goal of running the case nonlinearly at the beta threshold.

This case, based on Shot 87009 is an extremely challenging due to its large Shafranov shift, extreme triangularity, and large shear in the exterior region. Our primary problems in running this case have been in the preprocessing of the equilibrium data (the EQDSK files from EFIT) which resulted in inaccurate calculations of the magnetic field and skewed rectangular blocks which are extremely prone to numerical error. A number of improvements in FLUXGRID have been implemented to overcome these problems and grids can now be generated (see Figure 7) which can run successfully for the highest, most critical beta values (

Using this grid, we have been able to begin studies of the stability of this equilibrium as a function of b_N. The growth rates as a function of b_N obtained by both NIMROD and GATO are show in Figure 8. None of the NIMROD results are converged, and should be considered preliminary. The unstable modes for b_N = 5 and b_N = 6 are ideal interchange modes. The unstable eigenfunction for the case b_N = 6 is shown in Figure 9. GATO (which determines ideal MHD stability only) shows stability for b_N < 5. However, NIMROD shows instability to a slowly growing 2/1 resistive interchange mode for these lower values of b_N. [Alfonso Tarditi].

A campaign has begun to determine the range of Lundquist number S accessible by the NIMROD physics kernel. The goal of the study is to push practical calculations to high enough values of S to be relevant to present experiments. Initial studies have concentrated on the 2/1 tearing mode in a doubly periodic cylindrical tokamak [Holmes, et al., Phys. Fluids 26, 2569 (1983)]. The converged linear growth rate and eigenfunctions at for the range

have been determined. The results are shown in Figure 10. The growth rate is found to scale as S^-0.67, in good agreement with analytic theory. In all cases, grid packing around the low order rational surfaces has been found to be crucial for obtaining accurate results. In the future these studies will be extended to shaped cross-sections in toroidal geometry. [Dalton Schnack]

With some very slight modifications to boundary conditions, NIMROD has been used to simulate spheromak sustainment from a magnetized coaxial plasma gun. Results of toroidally averaged poloidal magnetic flux and RB_f are shown below for one simulation; the inner electrode extends from the bottom of the domain to Z=1.0. Figure 11 displays the amplification of flux embedded in the gun electrodes resulting from MHD activity, but it does not indicate the stochastic nature of the 3D field.

The MHD activity is qualitatively the same as that which occurs in simple can configurations with electrodes at the end plate surfaces. A summary of this physics study on spheromak sustainment was presented at the APS-DPP meeting in Seattle. [Carl Sovinec]

Our medium-term goal is to study how energetic particles (beam, RF and fusion-alpha) effect stability and nonlinear saturation of the internal kink mode. Our short-term goal is to follow gyrokinetic particles in NIMROD fields.

We now have all of NIMROD running on our local DEC Alpha workstation. We are currently working on setting up an internal kink test case working from the GATO-Ideal MHD comparisons that was suggested by Tom Gianakon. For now, we are using cylindrical geometry. We have formulated particle interpolation, deposition and tracking with quadrilaterals. We have already coded a particle tracking algorithm, but the particle motion is incorrect (work in progress). We have also formulated a df method which will allow straight-forward evolution of the particle weights given the equilibrium temperature and density of the energetic component on the NIMROD grid. [Scott Parker, Charlson Kim]

A new web server was put into operation in the past quarter to fix some of the performance problems which had plagued the web site since it moved from LANL to SAIC. At the same time the new web server went into operation, the domain name changed from http://nimrod.saic.com to http://nimrodteam.org to better reflect the multi-institutional nature of the team, and to allow for easy transfer of web site responsibilities in the future. The web site itself was improved with errors and outdated links corrected. The biggest improvement has been in the addition of an automated mailing list manager for the NIMROD User's Group. People can subscribe and change the mailing list options all from the web.

In concert with many of the improvements in XDraw, the XDraw web page has been greatly improved and now includes documentation written by Nina Popova. Although this documentation satisfies one of the most requested needs, documentation remains a very underdeveloped area. Current plans are to finish several chapters in the User's Manual, and then put it on the developer's web page to allow for a more interactive and distributed development of the documentation. [Scott Kruger]

Prepared by the NIMROD Team