An Introduction to the NIMROD Fusion Magnetohydrodynamics Simulation Project

1
An Introduction to the NIMROD Fusion
Magnetohydrodynamics Simulation Project
Prof. Carl Sovinec Department of Engineering
Physics University of Wisconsin-Madison presented
at Argonne National Laboratory November 21, 2002
CEMM
2
Goals for NIMROD(Non-Ideal Magnetohydrodynamics
with Rotation, an Open Discussion Project)

• Develop a simulation code package for studying
  three-dimensional, nonlinear electromagnetic
  activity in laboratory fusion experiments.
• Allow flexibility in the geometry and physics
  models used in simulations.
• Allow efficient computation on a wide range of
  platforms from PCs to massively parallel
  supercomputers.
• Provide user-friendly features, such as a
  graphical interface and documentation, and make
  the code publicly available. http//nimrodteam.o
  rg
• Apply techniques such as integrated product
  development and quality function deployment in
  design and development.

3
The project has been a multi-institutional effort
since 1996.
Curt Bolton, OFES Dan Barnes, LANL Dylan Brennan,
GA James Callen, Univ. of WI Ming Chu, GA Tom
Gianakon, LANL Alan Glasser, LANL Chris Hegna,
Univ. of WI Eric Heldstudents, Utah
State Charlson Kim, CU-Boulder Michael Kissick,
Univ. of WI
Scott Kruger, SAIC-San Diego Jean-Noel Leboeuf,
UCLA Rick Nebel, LANL Scott Parker,
CU-Boulder Steve Plimpton, SNL Nina Popova,
MSU Dalton Schnack, SAIC-San Diego Carl
Sovinecstudents, Univ. of WI Alfonso
Tarditi, NASA-JSC
Presently, there is non-team-member use of the
NIMROD code at LLNL, IFS, Univ. of WA, UCLA, and
AIST-Japan.
4
Characteristics of Magnetically Confined Fusion
Plasmas

• Reactor-grade conditions require a mix of
  deuterium and tritium at ntE?1020 m-3s and T?10
  KeV. Both species are fully ionized at these
  conditions.

• The Lorentz force, qV?B, confines the
  perpendicular motion of charged particles in a
  magnetic field.
• A toroidal container can have lines of B
  completely enclosed, but the field must be
  twisted in order to avoid rapid perpendicular
  particle drifts.

5
Ideally, Magnetic Configurations Consist of
Nested Toroidal Flux Surfaces

• Each charged particle will tend to remain
  confined on a magnetic surface.
• Surfaces provide insulation for hot and dense
  conditions in the center.

Equilibrium Flux Surfaces and Pressure from the
General Atomics DIII-D Tokamak
6
Pressure and current density gradients can drive
asymmetric collective modes unstable, changing
the magnetic topology.
Puncture-plots show simulated magnetic topology
changing gradually over 300 wave transit times
around the DIII-D tokamak. Initial conditions
are taken from experimental measurements.
7
Integrated modeling is the new horizon.
Simulations of the Pegasus tokamak at the Univ.
of WI are suggestive of what is possible.

• Plasma current and separatrix evolve
  self-consistently with applied loop voltage and
  vertical-field ramp.
• Transport has a strong influence on dynamics.
• High-order spatial accuracy is essential for
  distinguishing closed flux and open flux through
  modeled transport effects.

8
This study integrates MHD and transport effects
with realistic geometry and experimental
parameters.
loop voltsg
EMF (V)
Current (kA)
fplasma current
Time (ms)
Pegasus data courtesy of A. Sontag.
Axisymmetric simulation results.

• Study has focused on 2D evolution, but 3D
  tearing-mode simulation is a straightforward
  extension for NIMROD.
• Results emphasize interaction between MHD,
  transport effects, and overall performance.

9
Physical models for these macroscopic dynamics
are based on fluid-like magnetohydrodynamic (MHD)
descriptions.

• Density and magnetic-divergence diffusion are
  for numerical purposes.

10

• Conditions of interest possess two properties
  that pose great challenges to numerical
  approachesanisotropy and stiffness.
• Anisotropy produces subtle balances of large
  forces, nearly singular behavior at rational
  surfaces, and vastly different parallel and
  perpendicular transport properties.
• Stiffness reflects the vast range of time-scales
  in the system, and targeted physics is slow
  (transport scale).

11
The NIMROD code has a unique combination of
advanced numerical methods for solving systems of
PDEs that describe high-temperature plasmas

• High-order finite element representation of the
  poloidal plane
• accuracy for MHD and transport anisotropy at
  realistic parameters Sgt106, c/cperpgt109
• flexible spatial representation
• Temporal advance with semi-implicit and implicit
  methods
• multiple time-scale physics from ideal MHD (ms)
  to transport (10-100 ms)
• Coding modularity for physics model development
• Large-scale parallel computing capability

12
The finite element method provides an approach to
spatial discretization that has the needed
flexibility and accuracy.
NIMROD uses 2D finite elements (that are general
with respect to the degree of polynomials used
for basis functions) for the poloidal plane and
finite Fourier series for the periodic direction,
which may be toroidal, azimuthal, or a periodic
linear coordinate.
13
The semi-implicit advance is derived through the
differential approximation for an implicit time
advance for ideal linear MHD with arbitrary time
centering, q.
Using the alternative differential approximation
to the resulting wave equation leads to
where L is the ideal MHD force operator. We may
drop the Dt-term on the rhs to avoid numerical
dissipation and arrive at a semi-implicit advance
stable for all Dt where V is leap-frogged with B
and p.
14
A linear resistive tearing study in a periodic
cylinder shows that nearly singular behavior can
be reproduced with packed finite elements and a
large time-step.

• This S106 computation has a 32x32 mesh of
  bicubic elements and Dt100tA (1.8x105 times
  explicit). g is within 2.

15
Accuracy while varying the mesh and degree of
polynomial basis functions meets expectations for
biquadratic and bicubic elements.

• Divergence errors are too large with bilinear
  elements for these S106 conditions and the
  numerical parameters.

16
A nonlinear simulation of a classical
tearing-mode demonstrates full application
inlow-dissipation conditions.

• NIMROD Simulation
• StR/tA106
• PmtR/tn0.1
• tA1ms
• b ltlt1 to avoid GGJ stabilization

• DIII-D L-mode Startup Plasma
• R. LaHaye, Snowmass Report
• S1.6x106
• Pm4.5
• tA0.34ms
• tE0.03s

17
Small D (linear gtA5x10-4) leads to nonlinear
evolution over the energy confinement time-scale.
Saturation of Coupled Island Chains
Magnetic Energy vs. Time

• 5th-order accurate biquartic finite elements
  resolve anisotropies.
• 20,000 semi-implicit time-steps evolve solution
  for times gt tE.
• Explicit computation is impossibleg2x108
  time-steps.

18
Thermal conduction also exercises spatial
accuracy for realistic ratios of thermal
conductivity coefficients (109).

• Adaptive meshing alone cannot provide the needed
  accuracy in nonlinear 3D simulations magnetic
  topology changes across islands and stochastic
  regions.
• High-order finite elements provide a solution.

A simple but revealing quantitative test is a
box, 1m on a side, with source functions to drive
the lowest eigenmode, cos(px) cos(py), in T(x,y)
and Jz (x,y). Mass density is large to keep V
negligible.

• Analytic solution is independent of c,

• Computed T-1(0,0) measure effective cross-field
  conductivity.
• Any simple rectangular mesh has poor alignment.

19
Convergence of the steady state solution shows
that even bicubic elements are sufficiently
accurate for realistic parameters.

• Bilinear elements have severe difficulties with
  the test by conductivity-ratio values of 106.

20
Simulations of realistic configurations bring
together the MHD influence on magnetic topology
and rapid transport along field lines to show the
net effect on confinement.
SWINDLE these plots were handy but the
computation ran the MHD first, then thermal
conduction.
21
Tests of anisotropic thermal conduction at
various times during the nonlinear classical
tearing evolution reproduce an analytic wd-4
scaling. Fitzpatrick, PoP 2, 825 (1995)

• Conductivity ratio is scaled until an inflection
  in T within (2,1) island is achieved.
• Power-law fit is c/cperp3.0x103 (wd /a)-4.2.
• Result is for toroidal geometry.

• High-order spatial convergence is required for
  realistic anisotropy.
• Implicit thermal conduction is required for
  stiffness.

22
Solving ill-conditioned matrices is often the
most performance-limiting aspect of the algorithm.

• The condition number of the velocity-advance
  matrix can be estimated as

which can be gt 1011 in some computations.

• We have been using a home-grown conjugate
  gradient method solver with a parallel
  line-Jacobi preconditioner.
• It has been running out of wind on some of the
  more recent applications, forcing a reduction of
  time-step.
• We are presently implementing calls to Sandias
  AZTEC library, but we are interested in other
  possibilities (PETSc), too.

23
Conclusions

• Test results and past and present physics
  applications show the effectiveness of combining
  the semi-implicit method with a variational
  approach to spatial representation.
• Improved performance is expected from algorithm
  refinements.
• Iterative solution methods
• Adaptive meshing
• Advection (not discussed here)

24
Directions for the Project

• Hall and other two-fluid terms are in the NIMROD
  code, but the implementation requires small
  time-steps for accuracy.
• We are working on improved formulations.
• The ability to solve nonsymmetric matrices is
  important for this.
• Kinetic physics
• Parallel electron streaming effects E. Held,
  USU
• Gyrokinetic hot ion effects C. Kim and S.
  Parker, CU
• Resistive wall and external vacuum fields T.
  Gianakon, S. Kruger, and D. Schnack