NIMROD Progress Report
3^{nd} Quarter FY 2000
July 22, 2000
1. Overview
We report progress for the NIMROD Project for the period April-June, 2000.. Progress during this period includes:
1. A thin resistive wall boundary condition bounded by a vacuum region has been implemented in the bilinear element version of NIMROD for periodic slab and cylinder cases;
2. Changes for the higher-order elements have been merged into the main line NIMROD code repository, and the combined code has been made available to the Team for testing and evaluation;
3. Work has started on the M3D-NIMROD comparisons for two-fluid simulations;
4. The FLUXGRID pre-processor has been rewritten to improve this accuracy as well as allowing for more flexible grids, including packing in the theta direction;
5. The NIMROD web site has been modified to improve its functionality and appearance;
6. NIMROD has successfully computed the nonlinear resistive MHD evolution of D3D discharge 86144 at 2250 msec. with S = 10^{7}, albeit with relatively large viscosity (Pr = 10^{3}).
Progress in some areas continues to be slow due to lack of manpower.
2. Status
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:
1. The physics kernel contains the following physical models:
2. The following grids have been implemented:
3. Interfaces with several commonly used equilibrium file formats are available.
1. Beta-testing of the high-order finite elements;
2. An algorithm for a vacuum region and a moving separatrix is being implemented;
3. The implementation of kinetic and energetic particle effects by particle pushing using the gyro-kinetic model is continuing;
4. 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.
3. Highlights
Resistive Wall and Outer Vacuum
A thin resistive wall boundary condition bounded by a vacuum region has been implemented in the bilinear element version of NIMROD for periodic slab and cylinder cases. The implementation extends to the toroidal case once an appropriate response matrix becomes available through Morrell Chance's vacuum code [M. S. Chance Phys. Plasmas 4 (1997), 2161]. The thin resistive wall approximation presumes that the normal component of magnetic field is constant through the wall. This is used to evaluate the vacuum solution outside the wall and provides the normal component of the normal gradient of the magnetic field outside the resistive wall. The plasma response provides a similar gradient on the inside of the wall, and the difference is used to advance the normal component in the wall. The updated normal component is then applied as a boundary condition to the NIMROD solution and as such is explicitly applied in the time-step.
Our initial set of test simulations have focused on the periodic slab geometry with a tearing unstable equilibrium [J. E. Finn, Phys. Plasmas 5 (1998), 461]. The equilibrium current is allowed to flow in the periodic direction or in the finite element plane so that convergence properties of the finite element representation and the Fourier representation can be studied. In zero-flow cases that compare results from a resistive wall versus a conducting wall, the island widths agree with the reference when the mode is oriented in the Fourier direction. When the mode is oriented in the finite-element direction, the island width tends to be larger and distorted in shape. The cause may be a lack of poloidal resolution in our test.
Fig. 1. Magnetic flux surfaces for the saturated tearing mode in a sheared slab with a resistive wall on the right side (no flow).
Fig. 2. Magnetic flux surfaces for the saturated tearing mode in a sheared slab with a conducting wall on the right side (no flow).
The introduction of flow across the island in the finite element direction has proven problematic with a possible cause identified as the synergistic /flow instability. Therefore, our first efforts with flow have focused on the Fourier representation, which should exhibit a weaker problem. In this orientation with slow flow velocities, locking is observed. At high flow velocities, the time-step had to be reduced to approximately 1/100 of an Alfvén time (t_{A}) such that the flow CFL is around 0.03. In one case the equilibrium flow is given a magnitude such that . The simulation exhibits a slow decay of the total flow as radial magnetic field slowly diffuses through the wall, and translation continues for three resistive times (S=10^{5}). The net flow eventually reaches a value of , at which point rapid locking occurs, and the island grows to the zero-flow resistive wall result. The diagnostic drops at least an order of magnitude between the unlocked and locked states, suggesting that the synergistic instability may be playing a role prior to locking. [Tom Gianakon]
Fig. 3. Evolution of the radial component of magnetic field at a position near the non-ideal wall as locking forces slow the net flow, and the magnetic island becomes stationary in the laboratory frame.
Higher-order Finite Elements
During this quarter, changes for the higher-order elements have been merged into the main line NIMROD code repository, and the combined code has been made available to the Team for testing and evaluation. However, prior to this limited distribution, two important changes were made, and nonlinear tests were performed.
The basis function definitions in quadrilaterals now use general formulas for computing the value of a basis function at an arbitrary point in its domain. In the 1D domain [0,1], the n-th degree Lagrange polynomial has nodes at
.
Since each Lagrange basis function is nonzero at only one of the nodes in the domain, they must have the form
.
Normalizing so that gives
.
It is extremely unlikely that this is a new result, but Team members had not encountered it in the literature. For 2D quadrilaterals, the n^{2} basis functions are
.
These general definitions have been implemented in NIMROD, along with general Gaussian quadrature integration for the required level of accuracy, so that the polynomial degree in quadrilaterals is now unrestricted. In principle spatial convergence may be achieved through grid resolution, polynomial degree resolution, or both.
The second change is the elimination of interior-centered data prior to calling the iterative matrix equation solver. Coupling to the interior nodes is always local across a grid cell for Lagrange elements. Therefore, a direct elimination of this data only affects local connections that are nonzero in the original matrix. The amount of data that is eliminated scales as in grids of quadrilaterals, where n is the polynomial degree. As n is increased, the iterative solver does less work, helping overall efficiency.
Final tests before the limited distribution included nonlinear simulations. An RFP computation was performed with bicubic elements, and convergence on a spheromak result was obtained by running a bilinear computation with a 40x40 grid, then a biquadratic computation with a 30x30 grid, and finally a bicubic computation with a 30x30 grid. [Carl Sovinec]
Benchmarks for Computational Initiative
Work has started on the M3D-NIMROD comparisons for two-fluid simulations. We have received a test case from Linda Sugiyama. It is an ideal m=1 kink mode for a tokamak equilibrium in a periodic cylinder. The M3D profiles have been incorporated in NIMROD, and we have started the simulation comparisons. To date we have only compared MHD results to serve as a guide before attempting two-fluid simulations.
Initial results are very encouraging. NIMROD achieved growth rates within a few percent of the M3D results almost immediately. Figure 4 shows typical linear eigenfunctions for the radial flow. Note that the mode has the ideal parity with B_{r} vanishing at the m=1 resonance (located at r/a~.4), and V_{r} doesn't change sign at the resonance. This is consistent with an ideal mode and is as expected. Convergence studies with viscosity, time-step, and mesh size are now underway. [Rick Nebel]
Fig. 4a. B_{r} vs. r as a function of axial position
Fig. 4b. V_{r} vs. r as a function of axial position
Pre-processor
It had been observed that our pre-processor fluxgrid had difficulty maintaining accuracy of the Grad-Shafranov equation, with Grad-Shafranov equation residuals routinely exceeding unity near the edge of the plasma. The pre-processor has been rewritten to improve this accuracy as well as allowing for more flexible grids, including packing in the theta direction. Various minor improvements, such as making it compatible with the new higher-order element version of nimrod, have been done also. [Scott Kruger]
Web Site and Documentation
A summer intern has been hired at SAIC to maintain the web site and improve it's appearance. Various improvements have already been done to improve it's organization and readability. Any suggestions for improvement are greatly welcome and may be sent webmaster@nimrodteam.org. [Scott Kruger]
High-S Campaign
The High-S Campaign has previously demonstrated accurate scaling of linear tearing mode growth rates at values of the Lundquist number S up to 10^{9}. Recently the High-S Campaign has been concentrating on the much more difficult problem of the nonlinear evolution of an equilibrium corresponding to D3D shot 86144 at 2250 msec. This case has been previously modeled at low-S [Ref: Gianakon, et al, http://www.nimrodteam.org/presentations/Sherwood99/sherwood99_gianakon.ps] Shot 86144 displays continuous n = 1 sawtooth oscillations. At 2250 msec. a 3/2 magnetic island appears after a sawtooth crash. The size of this island is amplified by successive sawtooth crashes until it reaches a saturated level. This perturbation has been interpreted as a neo-classical tearing mode triggered by a sawtooth crash. The discharge eventually terminates when a large 2./1 island appears and locks to the wall. Successful modeling of this shot is extremely challenging, and requires:
The present campaign concentrates on the first point, above, and aims at computing the nonlinear resistive MHD evolution of this discharge at high S.
We begin with an equilibrium determined from an efit reconstruction of the discharge at 2250 msec. A flux-based grid computed from this equilibrium is shown in Fig. 5. The equilibrium has q(0) = 0.98, and q_{sep} = 4.1. The grid has 128 radial (y) points and 64 angular (q) points. Three Fourier modes (n = 0, 1, and 2) are included in the simulation. The grid is packed non-uniformly about the q = 1, 3/2,and 2 surfaces.
Fig. 5. Poloidal grid (128 ¥ 64) used for initial high-S studies of D3D shot 86144 at 2250 msec. The grid is packed about the q = 1, 3/2, and 2 surfaces. Different colored regions indicate grid decomposition for parallel computing (here 6 blocks).
The equilibrium is perturbed at low amplitude (dB ~ 10^{-9} Tesla) and allowed to evolve in time. Calculations have been made with increasing S, but with fixed viscosity. For S = 10^{5}, 10^{6}, and 10^{7}, this corresponds to magnetic Prandtl numbers of Pr = 10, 10^{2}, and 10^{3}, respectively.
The kinetic energy for the three Fourier modes [n = 0 (yellow), n = 1 (blue), and n = 2 (magenta)] as a function of time for the case S = 10^{7}, Pr = 10^{3} is shown in Fig.6. Both the n = 1 and n = 2 modes are linearly unstable (g_{n=1} = 9.883 ¥ 10^{3} sec^{-1}, g_{n=2} = 6.155 sec^{-1}). The linear eigenmode structure is shown in Figures 7 and 8, where we display vector plots of the poloidal components of the velocity for each mode.. The n = 1 is dominantly 1/1, with a small 2/1 component, and the n = 2 is dominantly 2/2, with a small 3/2 component.
Fig. 6. Kinetic energy vs. time for three Fourier modes [n = 0 (yellow), n = 1 (blue), and n = 2 (magenta)] for the case S = 10^{7}, Pr = 10^{3}.
Fig. 7. Poloidal velocity for the n = 1 linear eigenfunction in the core of the discharge. Note the dominant 1/1 structure, and the small 2/1 "halo".
Fig. 8. Poloidal velocity for the n = 2 linear eigenfunction in the core of the discharge. Note the dominant 2/1 structure, and the small 3/2 "halo".
Fig. 9. Magnetic field topology in saturated state for D3D shot 86144.2250 with S = 10^{7}, Pr = 10^{3}. Note the 1/1, 2/2, and 2/1 structures, and the absence of 3/2 magnetic islands.
At t ~ 2 msec. the stable n = 0 mode becomes nonlinearly driven by the n = 1 mode. At t ~ 2.2 msec. the n = 2 mode makes a similar transition. The driven character of these modes is indicated by their growth rates becoming twice that of the n = 1. The growth of all modes saturates at t ~ 3.1 msec.
A Poincaré plot showing the magnetic field topology at t = 4 msec. is shown in Fig. 9. The 1/1 and 2/2 distortions in the center and the 2/1 magnetic island are clearly visible. Note the absence of a 3/2 magnetic island for this case, although a moderately sized 3/2 structure was seen with the less realistic parameters S = 10^{5} and Pr = 10.
NIMROD has successfully computed the nonlinear resistive MHD evolution of D3D discharge 86144 at 2250 msec. with S = 10^{7}, albeit with relatively large viscosity (Pr = 10^{3}). The driven 3/2 magnetic island that was seen in the experiment (and which presumably serves at the seed island for a neo-classical tearing mode) is not seen in this simulation, although it does appear at lower values of S. Calculations are underway at S = 10^{7} with lower viscosity (Pr = 10). Preliminary indications are that a substantially finer grid will be required for these cases. [Dalton Schnack]