NIMROD Progress Report

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

1. Work on generalizing the finite element basis functions is underway;

2. Implementation of the AZTEC linear solver into a version of the code;

3. Completed simulations of the destabilization of ballooning modes by growth of the internal kink in D3D-shaped geometry;

4. Pressure equilibration by parallel thermal conduction on magnetic surfaces associated with islands has been demonstrated in slab geometry, but it remains problematic in realistic toroidal geometry;

5. Two-fluid simulations of RFP configurations in a limited poloidal wavenumber range demonstrate large shear flows generated by nonlinear forces;

6. A campaign to validate NIMROD on linear and nonlinear cases at realistic values of S has begun;

7. An interface from FLUXGRID (the interface from experimental equilibria to NIMROD) to PEST3 has been completed;

8. Work was begun to improve code performance in the presence of equilibrium flows;

9. New gridding options have allowed high-b D3D equilibria to be studied;

10. An improved algorithm for the generation of regions with triangular grids is being developed;

11. Grid packing (non-uniform radial grid) has been implemented for both circular cross section cylindrical geometry and flux surfaces based grids;

12. The present model for a vacuum region has been shown to converge to the correct linear growth rate for external modes, but with the penalty of an increased iteration count;

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

The NIMROD physics kernel has matured to be able to handle a number of new and interesting physics problems. During the past three months the NIMROD Team has been using the kernel to expand the application areas of the code. However, pushing the envelope always leads to the identification, and hopefully the solution, of new problems in the code system itself. (In fact, this is a primary purpose of such an exercise.) For example:

Other difficulties are being experienced in the areas of pre-processing and post-processing. These include:

Successful resolution of these issues will require a significant increase in team resources in specialized computer science areas. Without this support the team will likely be unable to deliver the level of user-friendliness originally envisioned.

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 (eg,, 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. An interface to the AZTEC parallel linear solver package is being implemented;

6. The implementation of kinetic and energetic particle effects by particle pushing using the gyro-kinetic model has begun;

1. Linear stability of 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, and the linear stability of the SSPX spheromak;

5. Two-fluid effects in RFPs and FRCs;

6. High-S campaign.

Simulations of the destabilization of ballooning modes by growth of the internal kink in D3D-shaped geometry have been completed. Shaping effects appear to suppress the disruptive nature of the ballooning modes, but further work is required. An example of the saturated pressure contours is shown in Figure 1. [Tom Gianakon]

Pressure equilibration on magnetic surfaces associated with islands has been demonstrated in slab geometry, but it remains problematic in realistic toroidal geometry. A cylindrical tearing mode is being investigated to consider whether the problem stems from the toroidal geometry. The pressure equilibration is the primary stumbling block at this point to realistic simulations of neoclassical tearing modes. [Tom Gianakon]

Two-fluid simulations of RFP configurations and a limited poloidal-wavenumber range demonstrate large shear flows generated by nonlinear forces. A series of simulations has found temporal convergence when Dtg_linear~0.01 for magnetic Prandtl numbers of ~40. A high axial-wavenumber mode develops with magnetic Prandtl number at unity, and it saturates at an amplitude of about 1/10% with respect to the average magnetic field. This mode is not present in simulations using only the resistive MHD Ohm's law. It may be a numerical artifact, and much more investigation is warranted. Studies with a broader poloidal-wavenumber range are beginning. [Rick Nebel]

A number of structural changes to FLUXGRID (the interface from experimental equilibria to NIMROD) have been made to increase its usefulness to the fusion community at large. The new structures make it simpler to implement it as a preprocessor for other codes. By customer request, work on an interface from FLUXGRID to PEST3 has been completed to allow additional validation comparisons for NIMROD. This interface will soon be available in a released version of the NIMROD package and will be useful for D3D simulations of neoclassical tearing modes. [Tom Gianakon and Scott Kruger]

Nonlinear cases generating strong flows and linear cases with large equilibrium flows (³ 0.1 v_A) lead to numerical instabilities with NIMROD, unless substantial viscosity is used. In an attempt to address this problem, two separate development branches have been created. One applies the recommendations in "Stability of Algorithms for Waves with Large Flows," by R. Lionello, Z. Mikic, and J. Linker, submitted to JCP. Those recommendations are:

The second branch modifies the time step with a splitting similar to what is used in the ALE (Arbitrary Lagrangian Eulerian) scheme. Simplified numerical analyses of both schemes demonstrate stability with strong flow in 1D; however, neither seems to solve the problem in 3D NIMROD simulations. Literature on incompressible fluid computations with finite elements is being studied, and we have learned that our treatment of

may not be adequate. This is motivating a change in basis functions described below. In the meantime, correction 1) of the Lionello scheme has been incorporated in the recent NIMROD release. [Carl Sovinec]

Work on generalizing the finite element basis functions is underway. Besides addressing the difficulties with flow, the modifications are intended to improve spatial convergence, particularly in the poloidal direction and to reduce or eliminate difficulties associated with nonzero

. So far, the preprocessor NIMSET and the postprocessor NIMPLOT have been converted to the generalized basis functions. A test of interpolated B distributions (not results of finite element computations) have matched the theoretical prediction of second-order convergence in the

error for biquadratic elements. [n^th order derivatives are predicted to converge n-orders more slowly than the function itself.] A significant number of the modifications required for the kernel have been completed while converting the postprocessor, but much more development is required to finish the conversion. [Carl Sovinec]

Internal kinks are expected to determine the MHD stability of the Electric Tokamak (ET) under completion at UCLA because of its high aspect ratio. We have embarked on a linear study of 1/1 resistive internal kink modes for ET profiles and parameters using the latest version, Version 2.3.3, of NIMROD. Typical ET start-up plasma parameters used are: central density of 6.3x1018 m-3, peak beta value of 0.8%, magnetic field of .25 Tesla, aspect ratio of 7.1. The profiles used are shown in Figure 2. We have also applied the full MHD version of the FAR code to this study. Use of the TOQ-FAR interface written by Scott Kruger enables FAR and NIMROD to start from an identical TOQ equilibrium. We are delighted to report that the latest version of NIMROD has made it possible to access much higher and more realistic S values (S ~ 106) than had ever been possible before (S ~ 104) for this equilibrium and where overlap with FAR occurs. The results from an S scan using both NIMROD and FAR are displayed in Figure 3. The agreement in the linear growth rates is excellent, and the scaling of S^-1/3 agrees with the theoretical and numerical calculations reported in Hastie et al., Phys. Fluids 30, 1756 (1987). The slight discrepancy between NIMROD and FAR can be attributed to the fact that the FAR calculations were performed without any viscosity whereas the NIMROD calculations were carried out with a Prandl number of 1 for all S values.

Most importantly, we are testing NIMROD for poloidal flow. Since nearly all other tests in flow are toroidal, we are filling a gap in the testing. So far NIMROD is having trouble running without numerical instabilities for the reasonably strong flows which we expect on ET. We are using the same 1/1 mode low beta ET case as with the benchmark to FAR. The ideal situation would be if NIMROD could simulate this mode with a realistic ET poloidal rotation profile and a magnitude of about 104 to 105 m/s that is in between the poloidal sound and poloidal shear Alfven speed. This speed would correspond to an H-mode type bifurcation. We need this rotation study done at the S = 105 range where the mode is believable and matches FAR. So far the fastest poloidal rotation which we can do in this case is only <100 m/s (a couple orders of magnitude too low). Work is quickly progressing though, and with the success of this goal, NIMROD will be ready to handle many ET relevant problems where no other code could. We expect this to be very helpful for interpreting stability issues and data in the experiment.

The latest version of NIMROD is currently being compiled and tested on an SP cluster platform (AIX, RS/6000 processors) available at UCLA. Because AIX F90 compilers are very rigid, this exercise refines the many parts of the code that could be written more precisely. For Version 2.3.3, we found only two small places that needed very minor improvements. Since NERSC will be moving towards SP systems soon, this exercise is also helping the team with the proper makefile constructions and any changes for immediate use on the NERSC SP machines. We have applied for and may receive permission to be "early users" on the new NERSC SP machine and NIMROD would be one of the codes tested there. [Michael Kissick and Jean-Noel Leboeuf]

The team has a fall-99 milestone to extend NIMROD's linear solver options by integrating in the parallel iterative sparse-matrix library AZTEC, a DOE-sponsored library developed at Sandia National Labs. The biggest effort required for the integration is to perform a parallel remapping of the Ax = b equation from the data layout used by NIMROD to the layout needed by AZTEC. We have completed about 75% of the coding for this task, including options for real and complex, scalar and vector (multiple unknowns/grid point), and regular and irregular NIMROD grids. Currently we are awaiting a new version of AZTEC that will provide full support for complex pre-conditioners, partially motivated by NIMROD needs. We hope to complete the coding and testing of the new capability in time to present results at the November APS Plasma Physics meeting in Seattle. [Steve Plimpton]

The NIMROD code has been applied to the simulation of high-beta D3D discharges to investigate the growth of ideal MHD modes driven through their instability threshold. The D3D 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 4) which can run successfully for the highest, most critical beta values (

). The study is now progressing to complete the scan of the different values of the normalized beta and to explore different more options for the choices of the numerical parameters that control the grid generation. [Alfonso Tarditi and Scott Kruger]

As of September 17th, Charlson Kim (Physics PhD graduate student) and Scott Parker at CU-Boulder have began working on adding gyrokinetic particles to NIMROD. The present status of this project is as follows:

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 eigenfunction at S = 10⁵ have been determined. Linear studies at S = 10⁶ are almost complete. Nonlinear simulations at S = 10⁵ are being performed. 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]

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. One model is increase the resistive 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 and viscosity. The use of the shape variable, along with the advection of marker particles, allows for easy generalization to nonlinear regimes where the boundary can move. To validate this model, a simple cylindrical equilibrium which has analytic solutions [Shafranov, Soviet Physics - Technical Physics, 15, 175 (1970)] is studied. The equilibrium is pressureless and has a step-function current that produces the q-profile that is flat within the core and increases monotonically in the external vacuum region.

This equilibrium is unstable to an external kink mode. In Figure 5 we plot the ratio of the growth rate of this mode as determined by NIMROD to the analytic growth rate, as a function of the ratio of the vacuum resistivity to the core plasma resistivity. The growth rate asymptotes to the analytic value as the vacuum resistivity is increased.

This solution does not come without a price. In the vacuum region at high resistivity, the Ohm's Law becomes the elliptic problem

, while in the core region the equation is predominantly hyperbolic. This discontinuity of equation properties is a classic "stiff" problem. Computationally, it manifests itself as an increase in the number of iterations required by the linear solver. This is illustrated in Figure 6. This increase affects the overall performance of the code. Based on our experience with D3D equilibria, the time that would be required to run DI3D equilibria with a vacuum region is unacceptable. The problem is then not ``how can we converge to the proper solution", but ``how can we converge to the proper solution efficiently".

Two ways around the problem of increased iterations have been formulated. One is to multiply the Ohm's Law by the conductivity; which results in better-conditioned matrices. The other is to find an alternate method for de-coupling the velocity from the magnetic field in the vacuum region. The idea here is to find a way to affect the velocity in the vacuum region without either increasing the resistivity in the vacuum region or affecting the flows in the core plasma. Simply adding a drag force to the momentum equation in the vacuum region retards all vacuum displacements and makes the vacuum act as an effective solid boundary. This stabilizes the external modes. We are formulating an approach in the shape variable is used to define a direction perpendicular to the plasma-vacuum interface and a drag force will be implemented which acts on the velocity field tangential to plasma-vacuum interface in the vacuum region. This should provide the necessary de-coupling without eliminating the de-stabilizing vacuum displacements. Both approaches will be tested in the next few months. [Scott Kruger]

A number of improvements have been made in the past 3 months to the fluxgrid pre-processor. Of primary importance are the improvements made to handle highly-shaped equilibria which has been motivated by the D3D case which Alfonso Tarditi is working on. This case is extremely triangular, has a large Shafranov shift, and has strong magnetic shear within the last 5% of the plasma flux. To handle this case, we

This case has many of the same characteristics as the low aspect ratio tokamak equilibria which are being proposed as the target equilibria. The work done on this pre-processing should help for those cases in the future.

In the next 3 months, some major changes for FLUXGRID and NIMSET are being proposed to increase its flexibility for generating triangle meshes, and for improving the meshes in the vacuum region. [Scott Kruger]

A new User's Manual is now on the nimrod web site in draft form. Further improvements are being made slowly, with a more concerted effort being planned once the current kernel restructuring current underway is completed.

The nimrod web site has had several performance and network problems. To alleviate these problems, a new computer was ordered to serve as the nimrod web server. It should be online and tested by the APS meeting.

Most of the recent work on the nimrod web site has been aimed at improving the communication among the team. A method allowing for the automatic uploading of documents, web access to the CVS repository, and experimental mailing list server have all been implemented. Once the new web server is implemented, further improvements along these lines are being planned. [Scott Kruger]

Prepared by the NIMROD Team