**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]