Recent Algorithmic and Computational Efficiency Improvements in the NIMROD Code

S. J. Plimpton, Sandia National Laboratory

C. R. Sovinec and T. A. Gianakon, Los Alamos National Laboratory

S. E. Parker, University of Colorado-Boulder

and The NIMROD Team

presented at the

41st Annual Meeting of the American Physical Society, Division of Plasma Physics, November 15-19, 1999

OUTLINE

1. Introduction
2. 1. Code and Team descriptions
   2. Algorithm features
   3. Equations
   4. Challenges
3. Matrix Solution
4. 1. Global preconditioning for CG
   2. Linking AZTEC
5. Anisotropic thermal conduction
6. Generalized basis elements
7. Strong flow
8. 1. Time advance
   2. Divergence condition
9. Modeling kinetic effects

The NIMROD Code:

•  

• Developed to simulate nonlinear electromagnetic phenomena that fundamentally alter the magnetic confinement of fusion-grade plasmas.

•  

• Designed to handle the stiff mathematical equations that describe hot magnetized plasmas, where time and spatial scales vary by many orders of magnitude.

•  

• Implemented in a parallel algorithm to run on the fastest computers available today and in the future.

•  

• NIMROD is publicly available from http://nimrodteam.org

The NIMROD Team:

•  

• NIMROD is a multi-institutional project commissioned by the Department of Energy's Office of Fusion Energy Science. The product is publicly available to provide the greatest possible benefit to the fusion community.

•  

• The NIMROD Team presently includes members from SAIC-San Diego, SNL, UT-Austin, General Atomics, UCLA, CU-Boulder, MSU, and LANL.

PRESENT TEAM MEMBERS AND ADVISORS

NIMROD FEATURES

The NIMROD spatial representation is finite element for one plane and Fourier series for the third (periodic) direction.

•  

• Permits simulating detailed configurations for comparison with experiment and for simulating non-fusion applications.

•  

• Allows refinement in regions of interest.

•  

• Present and future developments will take advantage of the flexibility of finite elements to change basis functions and to dynamically adapt the mesh to solutions.

•  

Our time-split, semi-implicit advance of the coupled Maxwell/velocity-moment (electrons and ions) set of equations permits efficient parallel computation of slow, nonlinear electromagnetic phenomena.

•  

• The following magnetic field Poincare plots show magnetic topology changing gradually over ~300 wave transit times around the toroidal configuration of the DIII-D tokamak at General Atomics in San Diego. Initial conditions are taken from experimental measurements.

EQUATIONS

NIMROD presently solves Maxwell's equations plus the following set of fluid moment relations:

where quasineutrality, , is assumed, w_p is the plasma frequency, and r is the mass density.

CHALLENGES

•  

• Our semi-implicit time advance solves the full ideal MHD operator for accuracy at large time steps. The matrix is symmetric positive definite but is often very ill conditioned. Subdomain-based preconditioners lose efficiency in large grids. Below we describe a global parallel preconditioning scheme that scales well for regions of structured grid. We are investigating the AZTEC package for help with regions of unstructured grid.

•  

• Equilibration of temperature along field lines in realistic conditions stresses spatial truncation errors with extreme anisotropy. A semi-implicit form of anisotropic thermal conduction has been implemented and is now being used in simulations of neoclassical tearing modes.

•  

• Linear and bilinear basis elements are convenient, but convergence properties have been disappointing, especially with respect to the required poloidal resolution. We are generalizing the NIMROD formulation to use polynomials of larger degree.

•  

• Our time split semi-implicit scheme seems susceptible to numerical instabilities when fluid Reynolds numbers are significant. Efforts to improve the time and spatial discretization are expected to help.

•  

• We have begun code development to incorporate kinetic effects via gyrokinetic particles.

SOLVING ILL CONDITIONED MATRICES

When applied to high temperature magnetized plasma, the NIMROD time advance requires solution of an ill conditioned matrix.

•  

• The condition number of our semi-implicit matrix may be estimated as

, from the fastest normal mode represented in the numerical system on a given grid. [The smallest eigenvalue is approximately unity.] When simulating slow resistive behavior, we may need to take time steps that are much larger than the Alfven time for propagation across the entire domain. Thus, the condition number may be very large.

•  

• A further complication stems from parallel decomposition of the poloidal plane. Conjugate gradient iterations require only local communication of decomposed data. However, preconditioning needs some form of global communication to relax errors in the lowest wavenumbers.

•  

• A set of block-based preconditioners have been developed for NIMROD, but all become ineffective as the problem sizes are pushed to solve realistic conditions.

We have recently found that global line Jacobi preconditioning in the structured grid region scales well with problem size and time step. [See below.]

As a refresher, a Jacobi iteration scheme is the simplest matrix splitting technique for solving a matrix equation. If we split the matrix A into L+D+U, where D is the diagonal and L and U are the strictly lower and upper triangles, respectively, Jacobi iterations find the solution of

by iterating

or

The scheme converges provided the magnitude of the largest eigenvalue of D-1(L+U) is less than unity. Furthermore, a Jacobi iteration step may be used to precondition the conjugate gradient scheme.

The definition of D may be generalized to improve the effectiveness of each iteration step.

We have found that using the 1D matrices formed by matrix elements along the radial grid lines and along the azimuthal grid lines in the structured grid region produces an effective preconditioner when the semi-implicit matrix is ill conditioned. [This forms two sets of 1D matrices, and the solutions from the two sets are averaged.]

This line-Jacobi preconditioning is most effective when the 1D solves are extended as far as possible across the grid. Data for each line may be distributed on different processors. We use non-blocking MPI point-to-point communication to make the exchanges illustrated in the following schematic:

Decomposition swaps reorganize from the normal block pattern (above) to the line pattern (below) for global preconditioning, then swap back to the block pattern.

SCALING RESULTS

Fixed problem size scaling for global line Jacobi on the T3E for a 100x100 polar grid case. The reported times are the total time step loop ("loop"), the line swapping time ("line") including data copies and parallel communication, and the seam communication time ("seam"). Times are per time step, averaged over 5 steps. The efficiency is (1-block loop time)/(# processors * loop time). The number of blocks equals the number of processors in each case, and the nxbl / nybl decomposition is balanced.

Comparison of preconditioner effectiveness on the 100x100 case. Block-based line Jacobi with periodic blocks (nybl=1, nxbl=nprocs), and balanced decomposition are compared with the global scheme.

Comparison of time step rates (per CPU time per processor). The rate is the inverse of the loop time per node per time step.

Increasing problem size scaling for global line Jacobi on the T3E. The nxbl / nybl decomposition is balanced.

Comparison of preconditioner effectiveness with increasing problem size.

Comparison of time per time step on the T3E with increasing problem size.

NIMROD & AZTEC

Steve Plimpton

Sandia National Labs

Seattle APS Meeting - Nov 99

NIMROD support: Carl Sovinec

AZTEC support: Mike Heroux

Good News: Fully implemented in 2.3.2

Bad News: Doesn't compete with NIMROD rblock line-solve methods

What is AZTEC?

• Parallel linear solver library

• Solves Ax = b

• Variety of pre-conditioners -- Jacobi, ILU, ICC, ILUT, ...

• Variety of solvers -- CG, GMRES, ...

• Point or vector unknowns -- MSR, VBR matrix formats

• Real and complex

• Designed for unstructured grids, asymmetric matrices

• Highly parallel

• Developed at Sandia National Labs

• Similar in scope to Argonne's PETSc

Why do This

• Rapid trial of different solvers & preconditioners on NIMROD physics.
• NIMROD preconditioners work well on rblocks, but only have Jacobi for tblocks.
• AZTEC has automated multilevel solvers for general unstructured grids (still in beta test).

Key point:Ax = b is different in NIMROD vs AZTEC

•  

• NIMROD:

• •  

  • duplicated unknowns at processor boundaries

  •  

  • matrix rows/elements not fully formed at boundaries

  •  

  • F90 data structures with local indexing

   

• AZTEC:

• •  

  • unknowns assigned uniquely to a processor

  •  

  • processor owns all matrix elements (entire row of matrix)

  •  

  • C data structures with global indexing

AZTEC into NIMROD: One-time Setup

•  

• Determine matrix structure/decomposition in AZTEC format

•  

• Communication parameters to map NIMROD data -> AZTEC

•  

• Messy setup enables run-time remap to be fast.

•  

• Similar flavor to NIMROD inter-block "seaming" but uglier.

AZTEC into NIMROD: Matrix Creation

•  

• Done only when a NIMROD matrix is formed, factored.

•  

• Allows for multiple RHS solves on subsequent timesteps.

AZTEC into NIMROD: Solve

•  

• Map x,b vectors from NIMROD -> AZTEC format.

•  

• Set AZTEC options.

•  

• Call AZTEC solver.

•  

• Remap x from AZTEC -> NIMROD format via seam communication

Results

•  

• 3 test problems, bad condition numbers, convergence criterion = 1.0E-6

•  

• Timings in CPU secs on 1 proc of T3E

•  

• Complex is for hardest NIMROD solve = VMHD equation

Test #1: 64x100 circular grid w/ center point, 6465 unknowns

Test #2: 2 external tblocks + 48x32 circular grid w/ center point, 4494 unknowns

Test #3: same as #2 w/ 10x smaller timestep, 4494 unknowns

Comments on Results

•  

• Also tried AZTEC CG/BILU and CG/ILUT to no avail.

•  

• NIMROD line-solve pre-conditioners work best even when tblocks exist.

•  

• Only AZTEC method faster than NIMROD equivalent is real Jacobi.

•  

• GMRES is about 4x slower per iteration than CG.

•  

• AZTEC scales well in parallel, but so does NIMROD line-solve, so no win at this point.

What is hard for AZTEC ?

•  

• This is a (mostly) structured grid and a symmetric (HPD) matrix.

•  

• NIMROD line-solve preconditioner leverages azimuthal nature of solution.

•  

• AZTEC complex is stored as 2 real/imag matrices:

What to try with AZTEC?

•  

• Multilevel

•  

• should do well on problems where Jacobi works but is slow

•  

• AZTEC developers have seen significant speed-ups in some cases

•  

• Dual preconditioning

•  

• let NIMROD precondition rblocks, AZTEC the tblocks

•  

• let AZTEC do the solve

•  

• Coupled equations

•  

• formulate all NIMROD equations as one large matrix

•  

• motivated by other NIMROD issues

•  

• AZTEC may do better on this harder problem

ANISOTROPIC THERMAL CONDUCTION

See "Simulations of Neoclassical Tearing Modes with NIMROD," by T. A. Gianakon

GENERALIZED BASIS ELEMENTS

The truncation error of the finite element method is intimately tied to interpolation errors of the chosen discrete space.

For piecewise polynomials of order , the interpolation of a sufficiently smooth function u satisfies:

and

where h is a measure of the grid spacing, and and are the norm and semi-norm, respectively, of the Sobolev space over the domain W with the first subscript indicating order of differentiation. [Jiang, The Least-squares Finite Element Method, Springer-Verlag 1998, p. 64.]

Hence, order of convergence depends on the choice of polynomials used.

Until now, the NIMROD algorithm has been restricted to using only bilinear basis functions in regions of quadrilateral elements and linear basis functions in regions of triangular elements.

This restriction has led to at least two difficulties:

1.  

2. Spatial convergence in the azimuthal direction takes a large number of elements, often comparable to the resolution required in the radial direction. Besides the obvious costs associated with using more elements, this also increases matrix condition numbers and tightens flow CFL restrictions.

3.  

4. Our magnetic field representation has nonzero divergence, and controlling that error without being too invasive is often tricky.

We expect that using quadratic basis elements will improve both issues significantly.

•  

• Quadratic elements will better represent poloidal curvature in magnetic fusion configurations.

•  

• Quadratic elements will allow us to satisfy the Ladyzhenskaya-Babuska-Brezzi condition for the magnetic field, which is critical for the velocity in finite element computations of incompressible flows. [Gunzburger, "Mathematical Aspects of Finite Element Methods for Incompressible Flows," in Finite Elements Theory and Application, eds. Dwoyer, Hussaini, and Voigt, ICASE/NASA LaRC Series, Springer-Verlag, 1986.]

NIMROD is presently being modified to use any Lagrange type element, where only first-order derivatives are square-integrable.

The generalized version of NIMROD is partially functional, and we have been able to compare biquadratic and bilinear elements on two test problems.

1) Sheared slab with a linear n=0 tearing instability, where then number of cells in the "radial" direction is twice the number in the "azimuthal" direction.

•  

• The curves labeled "0" result with bilinear elements and the curves labeled "1" result with biquadratic elements.

•  

• The quantity k_divb is a wavenumber associated with error in magnetic field divergence. It is defined as

2) Resonant n=2 instability in an annulus. Azimuthal resolution is the same as radial resolution in these cases.

Comments on these test cases:

•  

• Biquadratic elements use twice as much data per direction as bilinear elements, so the CPU times for a biquadratic simulation scales like CPU time for a bilinear simulation with twice the number of mesh points per direction.

•  

• Only the divergence error for the sheared slab case shows clear asymptotic scaling. For bilinear elements the error converges linearly and for biquadratic elements it converges quadratically. This is expected from the convergence estimate for first-order derivatives.

•  

• Neither test is difficult with bilinear elements. Harder tests will provide more conclusive evidence regarding the efficacy of higher-order elements.

PROBLEMS WITH SIGNIFICANT FLOWS

•  

• NIMROD simulations with significant poloidal flow, even linear simulations, frequently result in numerical instability.

•  

• We have tested two schemes that are proven stable in 1D analyses.

•  

• The semi-implicit scheme described by Lionello, Mikic, and Linker [J. Comput. Phys. 152, 346 (1999)] has a leap-frog character. They recommend: 1) including the wave generating terms in the advection predictor step but with a coefficient of 1/2, 2) using the semi-implicit operator in both predictor and corrector steps, and 3) increasing the semi-implicit coefficient from the minimum value required without flow.

•  

• A time-split scheme related to an Arbitrary Lagrangian-Eulerian approach, where all variables are first advanced without advection, then advection is applied to each.

•  

• Neither approach solved the problem.

•  

• Present suspicions are that the Ladyzhenskaya-Babuska-Brezzi condition needs to be applied to the magnetic field representation in finite element MHD computations. The development of quadratic basis functions will make this possible.

KINETIC EFFECTS THROUGH SIMULATION PARTICLES

See "Kinetic MHD: Making the Micro-Macro Connection," by Scott Parker

SUMMARY

•  

• For cases without significant regions of unstructured grid, NIMROD now has an effective, scalable linear solver for conditions where the semi-implicit operator is very ill conditioned.

•  

• We have linked NIMROD with the AZTEC package to find better performance for regions of unstructured grid and to open the possibility of solving nonsymmetric matrices.

•  

• Anisotropic thermal conduction has been implemented with a semi-implicit approach, which seems to be effective even with widely disparate parallel and perpendicular diffusivities.

•  

• The finite element basis functions in NIMROD are being generalized to reduce resolution requirements and to alleviate problems associated with nonzero divergence of magnetic field.

•  

• Simulating strong flows remains a concern. We are expecting to make progress through the generalized basis functions.

•  

• Code development for modeling kinetic effects is underway.

Please visit us at http://nimrodteam.org or http://nimrod.saic.com, where this poster will be made available.