NIMROD Progress Report

1^st Quarter FY 2001

Prepared by the NIMROD Team

February 6, 2001

1. Overview

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

1. An improved algorithm for evaluating surface terms arising in the finite-element method has been implemented;

2. NIMROD is now ready for production use on the IBM SP3.

3. Benchmark calculations with M3D under the auspices of the Macroscopic Modeling Project are continuing;

4. Neo-classical tearing modes have been calculated with the "pressure and flow" closure. Threshold island widths have been observed;

5. Calculations of the feedback stabilization of neoclassical tearing modes by RF-induced current have begun;

6. Calculations of a high-b disruption in DIII-D shot 87009 are continuing;

7. Version 3.0.0 with arbitrary order finite elements has been released for production use.

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 being implemented

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;

3. Highlights

Surface Integrals

Neumann boundary conditions are typically treated as surface integrals in the functional of variational problems. The finite element method used in NIMROD is a variational approach, and surface integrals arise when there is a flux of something through the domain boundary. For example, a weak form of Faraday's law is

where g is a vector test function (of finite support in the finite element method). The second integral on the right represents a flux of magnetic flux due to a tangential electric field at the surface. In practice, surface integrals lead to local source densities in the discrete equations at the surface of the domain.

The numerical machinery for computing explicit surface integral contributions has been added to the physics kernel of NIMROD. The numerical integration routine is generic, and physical terms are specified in integrand routines--both are similar to what is used for volumetric integrals for quadrilateral and triangular cells. The first application is applied voltages in Faraday's law, as shown in the example above. Before the finite element basis functions were generalized, we were able to substitute a volumetric source that was localized near the wall. This approach led to qualitative errors with the generalized basis functions, where larger finite elements are treated more accurately. [Carl Sovinec]

Neoclassical Tearing Modes (NTM)

Neoclassical tearing mode simulations have continued based on the perturbed pressure gradient variants of the poloidal flow damping ansatz, e.g., the bootstrap current form. The effort has met limited success on reproducing the nonlinear threshold caused by anisotropic thermal diffusion. The caveat is that while a nonlinear threshold is observed, the threshold appears larger than what one expects from previous reduced MHD simulations and also analytic theory. The discrepancy is suspected to lie in a difference in direction for the perturbed bootstrap current between the simulations, which introduce the bootstrap current in the poloidal direction exclusively, and the analytic theory, which deals only with a parallel Ohm's law. The results are presented in the Figure 1.

Figure 1. Neoclassical tearing modes grow when sufficient perturbed bootstrap current as measured by the ratio of the electron poloidal flow damping rate (m_e) to the electron collisionality (n_e) is generated. A nonlinear threshold for the instability is implied.

An analysis of the neoclassical island evolution has required the development/enhancement of the existing field-line integration tool nimfl. The limitation of the old version of nimfl was that starting positions for field line integration were very restrictive especially when dealing with single helicity islands. A version of the nimfl tool has now been modified to use the same structure developed for resetting from old dump files, e.g. dump_reset.f. In principle, this also allows field line integration in grids that contain triangles. Flexibility in starting positions includes: 1) arbitrary (R,Z) positions, 2) a scheme motivated by Stuart Hudson's method for choosing positions along a resonant surface, and 3) logical coordinate grid positions. The second method involves computing a set of flux-surfaces from the total n = 0 magnetic field and determination of the q-profile on these new surfaces. Minimum, maximum, and zero values of the radial magnetic field for a given resonant surface correspond to the separatrix and o-point and are then used as starting positions. The tool has now also evolved to include the ability to compute a poloidal decomposition of the magnetic energy. This involves computing a PEST angle on the n = 0 surfaces and a Fourier decomposition of the appropriate projections of the magnetic field. This requires the Fourier decomposition of both covariant and contravariant projections of the magnetic field. An example decomposition from a reverse field pinch is presented in the Figure 2. The power is observed to be peaked in modes resonant with q = 1/5. The tool may be expected to also evolve into a mechanism for determination of an island width based on single helicity arguments. Such a tool is needed for the study of neoclassical tearing modes.

Figure 2. Poloidal decomposition of power from an RFP simulation demonstrates the functionality of a newly developed diagnostic capability for use with NIMROD.

Work is also progressing on the study of the feedback stabilization of neoclassical tearing modes by RF-induced current. The study is being conducted with a resistively unstable tearing mode in slab geometry. The desire in these simulations is to have an unstable tearing mode even at zero island width so that the effect of diffusion of the RF-induced current can be studied. In particular, the thermal diffusion effects associated with the NTM threshold that are expected to also be active set a minimum current channel width that limits the ability of the induced RF current to stabilize the island. Such an effect is illustrated in the accompanying Figure 3, where the width of the current source is varied and the minimum island width is determined. The sudden increase at the end of each line is the phase instability caused when the o-point and x-point of the island switch and the RF becomes destabilizing. [Tom Gianakon]

Figure 3. Stabilization of a tearing mode by auxiliary RF current drive indicates that the minimum island width is determined by a combination of the current channel width (Dx) and anisotropic thermal diffusion effects.

Energetic Particles

Progress on the implementation of kinetic particles in NIMROD continues. The capability of loading particles uniformly in logical coordinate space and evolving the particles in the NIMROD fields is now implemented. Current work includes further optimization of gathering of the NIMROD field quantities at the particle location and developing filtering techniques.

Although df will help reduce discrete particle noise, filtering is still expected to be needed to run with a manageable number of particles. To address this, we are developing a method for filtering the deposited quantities on the non-uniform grid. It is observed that there are dominant spatial non-uniformities on the order of the grid size in the deposited pressure. This is seen in Figure 4. As the grid size increases with the minor radius, so too does spatial extent of the noise. This is not unexpected. The deposition is done in a uniform logical space. The mapping from logical space to coordinate space stretches the non-uniformities.

Figure 4. Pressure non-uniformities in the deposited pressure from the unfiltered df-method.

To deal with this, the current implementation of the filter is performed in logical space. However, this filtering process is a function of position and is not applied uniformly in coordinate space. This may cause the introduction of an additional aliasing effect. The appropriateness of this is being considered and possibly better methods are being explored. This also raises a subtle issue concerning the particle shape, which is uniform in logical space but not in coordinate space. This means that as the particle orbits and enters different grids, the physical particle shape is changing. The implications of the effects of the nonuniform particle shape need to be investigated further.

In addition, we are implementing methods of loading arbitrary particle distribution functions onto the NIMROD grid. This will allow for the control of the driving profiles (density or temperature) and consistent implementation of the delta-f method. Loading arbitrary particle distributions will add to the flexibility of the kinetic particles in NIMROD and provide a step towards the eventual kinetic closure. [Charlson Kim and Scott Parker]

NIMROD on the NERSC IBM SP3

Drs. Hisaya Sugimoto and Akira Yaoita of the Japanese Electrotechnical Laboratory visited the LANL T-15 group during the month of December for the purpose of getting NIMROD running at the Tsukuba Advanced Computing Center. One of the two machines available there is an IBM SP3, similar to the gseaborg machine at NERSC. We have had problems creating an optimized version of NIMROD that will run accurately on an IBM SP3, but Drs. Sugimoto and Yaoita were able to consult colleagues in Japan. They recommended specifying the 32-bit executable option for any code using IBM's MPI library, and this led to major progress. We also broke some of the larger fortran 90 modules into smaller pieces, so that the optimizer would not take an excessive amount of CPU time. The entire code (Version 3_0) now compiles in an interactive log-in session and is ready for production use on IBM SP3s. [Carl Sovinec]

Two-Fluid Effects

At the last meeting of the Macroscopic Modeling Group, initial results on two fluid simulations were presented. Although the growth rated appeared to be similar to the MHD case, unphysical oscillations were observed in the instability growth. Further investigation indicated that these oscillations were growing from the boundary. Imposing a resistive layer at the boundary appears to have mitigated this problem.

However, a problem has emerged in the MHD comparisons. The growth rate results presented at the last meeting were too low by a factor of 2^1/2. This increases the discrepancy between our results and the M3D results. [Rick Nebel]

Vacuum Implementation

The numerical instability reported last time has been resolved by smoothing the initial distribution of the concentration variable. By eliminating the sharp discontinuity in the initial condition, the nonlinear advection of the concentration variable works well as can be seen in Figure 5, which shows the perturbed part of the concentration variable at various slices (recall that concentration variable > 0.5 implies vacuum). When the calculations are extended to have the three-dimensional resistivity be derived from the concentration variable, numerical problems occur when running on parallel processors. The cause of this problem is being investigated. [Scott Kruger]

Figure 5. The perturbed concentration variable at various poloidal slices shows the shifting of the surfaces as the concentration variable is advected.

Post-Processor Developments

The post-processor routines required to upload NIMROD data into the MDS+ data base have been ported over to the latest architecture. The data tree has been altered to reduce the amount of storage by taking advantage of advanced MDS+ features. Several suggestions were given at the Princeton meeting by other team members on better ways of storing the data, primarily on the best way of storing structured data versus unstructured data. These suggestions will be incorporated as greater familiarity with the MDS+ system occurs, and will be analogous to the M3D data structure to the extent possible. [Scott Kruger and Jeff Schachter]

Density Equation Studies

Tests on the continuity equation (so far implemented only for n=0) have been performed. At first, to check the basic hydrodynamic functionality, the case of a traveling plasma pulse without magnetic field was considered.

A customized version of NIMROD is used where, in addition to the continuity equation, a provision for external loading of an arbitrary plasma profile in cylindrical (or rectangular) coordinates was introduced.

Another important modification required for this test was the introduction of a model for open-end boundary conditions. The model essentially imposes a "sink term" in a small axial extent of the grid near one end of the cylindrical domain. The density, pressure and flow velocity profiles are then imposed in this region, rather than computed, to bring the fields smoothly to a background value that matches the Dirichlet boundary conditions imposed on the actual boundary of the domain.

The test case considers a cylindrical domain in 2D with axis along the z direction. For the initial condition a plasma pulse is loaded with gaussian shaped density and pressure profiles both in r and z. The pulse is initialized with a velocity component along z.

Fig. 6 shows the evolution of the density and pressure profiles. There is still a small effect of the open boundary on the plasma pulse that leaves the simulation domain, as can be seen in the ripple appearing on the profile of t=2.

Another test was performed in the context of developing a simulation of the plasma in the magnetic nozzle/exhaust region of the Variable Specific Impulse Magnetoplasma Rocket, VASIMR concept (for advanced propulsion applications, http://spaceflight.nasa.gov/mars/technology/propulsion/aspl).

Fig. 6. A non-magnetized, drifting plasma pulse along z in presence of open-end boundary conditions: the evolution of density and pressure profiles is shown at arbitrary times. The last 10% of the z-axis (form z=2.7 to z=3.0) is used to implement the open boundary condition.

Fig. 6 shows the time evolution of the plasma pulse as moves through the field lines: because the curved field lines and the density gradient the plasma acquires a significant radial velocity.

The simulation domain is shown in Fig. 7, where the actual VASIMR magnetic configuration is also depicted. A plasma pulse, similar to the previous test, is now drifting through a nozzle-shaped magnetic field, towards the open-end boundary

Fig. 7. VASIMR magnetic field lines and sketch of the 2D r-z simulation domain

Fig. 8. Density profile of plasma pulse traveling along the VASIMR magnetic field: because the curved field lines and the density gradient the plasma acquires a significant radial velocity. Limited radial resolution and a development of large density gradient close to the wall (still with ideal conductor boundary condition) are causing the large radial oscillation at the end of the run.

These simulations are of a very preliminary nature. The present effort is now directed towards developing a series of basic cases for quantitative theory checks.

Eventually, this model will be applied to study the exhaust plasma detachment and the optimization of the magnetic nozzle for maximum performances in realistic VASIMR configurations designed for space flight. [Alfonso Tarditi]