The DOE Princeton Plasma Physics Laboratory (PPPL), part of the Center for Extended Magnetohydrodynamics Modeling (CEMM), the Scientific Computation Research Center (SCOREC, RPI), part of the Interoperable Tools for Advanced Petascale Simulations (ITAPS) center and the Argonne National Laboratory, part of the Toward Optimal Petscale Simulation (TOPS) project, are jointly working on developing tools and methods to provide a realistic assessment of the mechanisms leading to disruptive and other stability limits in present and next generation of fusion devices.

This joint work, which has to aim to provide new insights into some of the most critical and complex phenomena in plasma and fusion science, requires the development of:

  • 3D nonlinear magneto-fluid based models of hot, magnetized fusion plasma.
  • Tools and methods to support the adaptive discretization, partitioning and parallel management of computational domain that could include several billion degrees of freedom.
  • Efficient parallel linear algebraic solvers to be used on petascale machines

Description of ITAPS tools

The implementation of such complex modeling to be solved on petascale machines requires the team developing new capabilities that include:


  • Construction of 3D computational domains encountered in plasma fusion (Figure 1) supported by the development of GMI (iGeom implementation) which supports queries on the geometric model, FMDB (iMesh/iMeshP implementation) which supports the parallel management of the domain discretization and IPComMan [1] service which improves the efficiency of parallel communications.
  • Development and implementation of specific domain interfaces using iMesh/iMeshP/iGeom to support specific queries of the mesh during the construction of the local FEA matrices.
  • Development and implementation of numbering mechanisms to store the defined degrees of freedom to support the construction of local FEA matrices in a parallel environment.
  • Development and implementation of matrix/vector container interfaces to support the efficient construction and compact storage of local matrices.
  • Development of parallel assembly algorithms to support the construction and storage of distributed global matrices.
  • Development of AMSI software [3] that integrates (Figure 2) standalone components (parallel mesh, geometry, FEA analysis (CEMM), linear algebraic solver (TOPS)) within a efficient parallel analysis framework (M3D-C1) (Figures 3).
  • Construction of an adaptive loop (Figure 4) aiming at better capturing physical phenomenon using the following four additional software components:
    • ErrorEstim which supports the estimation of the discretization error based on the calculation of hessian matrix (second derivative).
    • AdapMapGen which supports the calculation of both isotropic and anisotropic size fields based on calculation of the discretization error
    • MeshAdapt service which supports the parallel anisotropic mesh adaptation of unstructured mesh
    • MeshCurving service which supports the parallel higher order mesh adaptation of unstructured mesh
computational domains

Figure 1. 3D Domain made of 2D planes deployed along the toroidal direction

M3D-C1 framework

Figure 2. Integration of standalone components into M3D-C1 framework.

equi-pressure surface

Figure 3. (a) Initial equi-pressure surface at the boundary and some field lines. (b) NSTX sawtooth, showing pressure contours and surface with some B-lines.

adaptive loop

Figure 4. Adaptive loop applied to fusion application


  • Development of a high-order triangular element with continuous first derivatives (C1 continuity) which allows the Galerkin method to be applied without introduction of new auxiliary variables and reach an order of accuracy of h4 [2].
  • The extension of the 2D C1 element used in (r, z) planes to a full 3D spatial configuration where 1D cubic Hermite basis functions are used along the toroidal direction (Figure 1).
  • Development of a resistive wall boundary condition to allow both two-dimensional and three-dimensional magnetic perturbations to penetrate the wall.


  • The extension of PETSc software package to support the efficient parallel resolution of block diagonal matrices to be formed within the resolution of 3D plasma fusion models.


Adaptive loop: Adaptive methods reduce the total number of degrees of freedom (Figures 5, 6 and 7) by refining the mesh only in domain of interest and modeling general curved reactor domains (Figure 8).

anisotropic mesh adaptation

Figure 5. Anisotropic mesh adaptation to better capture fusion phenomenon.

density of unstable tearing mode

Figure 6. (a) Linear perturbed current density of unstable tearing mode, (b) close-up of linear perturbed current density with mesh superimposed.

Tokamat peeling-ballooning eigenmode

Figure 7. Tokamak peeling-ballooning eigenmode.

Tokamak euilibrium flux surfaces

Figure 8. (a) Tokamak equilibrium flux surfaces, (b) flux surfaces with adaptive mesh superimposed.

Figure 9. M3D strong scaling study on Jaguar (2008)

M3D-C1 Solver scalability: has been improved by using SuperLU_dist to factor the large sparse matrix associated with the M3D-C1 implicit formulation for the 2F MHD equations every 5-20 time-steps instead of every time-step as in the original software. In the intervening time steps, the last factored matrix is used as a very effective preconditioner for a GMRES-based iterative solve in PETSc, which then requires only a few iterations to converge. In the initial application, this resulted in a net speed increase of about 70% for a 2D "approach to 2F equilibrium with flow in NSTX" application.

As can be seen from Figure 9, the solvers are scaling even better than ideal, but the data copy is bringing the entire scaling down, confirming the need for efficient parallel IO.

Additional funding

PPPL funding: Center for Extended Magnetohydrodynamics Modeling (CEMM)


[1] A. Ovcharenko, O. Sahni, K.E. Jansen, C.D. Carothers and M.S. Shephard, Neighborhood communication paradigm to increase scalability in large-scale dynamic scientific applications. Parallel Comput., under review, 2010.

[2] S.C. Jardin, J. Breslau, N. Ferraro, A high-order implicit finite element method for integrating the two-fluid magnetohydrodynamic equations in two dimensions, J. Comp. Phys. 226 (2) (2007) 2146–2174].

[3] F. Delalondre, C. Smith, M.S. Shephard, Collaborative software infrastructure to support adaptive multiple model simulation, Computer Methods in Applied Mechanics and Engineering, Computer Methods in Applied Mechanics and Engineering, Volume 199, Issues 21-22, 1 April 2010, Pages 1352-1370.


FonTier Front Tracking: Simulation of Formation and Implosion of Plasma Liners for Magneto-Inertial Fusion

ITAPS Personnel: Roman Samulyak (SBU / BNL), Lingling Wu (SBU)
Fusion Personnel: Scott Hsu (LANL), Paul Parks (General Atomics)
Project Status: Active

Web Accelerator

FronTier simulation of the implosion of plasma liner formed by the merger of 625 jets

ITAPS front tracking technologies have been used in computational studies of the plasma jet driven magneto-inertial fusion (PJMIF). In the PJMIF concept, a plasma liner, formed by merging of a large number of radial, highly supersonic plasma jets, implodes on the target in the form of two compact plasma toroids, and compresses it to conditions of the fusion ignition. By avoiding major difficulties associated with both the traditional laser driven inertial confinement fusion and solid liner driven MIF, the plasma liner driven magneto-inertial fusion potentially provides a low-cost and fast R&D path towards the demonstration of practical fusion energy. The goal of PJMIF simulations is to evaluate the method by estimating the fusion energy gain and to provide guidance for the plasma liner experiment being built at Los Alamos. In the first phase of research, we optimized the nuclear fusion gain in the liner target parameter space via spherically symmetric simulations. Single and double-layer deuterium and xenon liners have been investigated as well as liners to be used in the PLX experiment. By varying target and liner parameters, the implosion process was optimized for maximum fusion energy gain and compared with theoretical predictions and scaling laws. In the most optimal setup, fusion ignition and energy gain of 10 was achieved with energy release of 10 GJ. Current work focuses on large scale 3D simulations of the propagation and merger of high Mach number plasma jets, the formation and implosion of liners, and compression of targets. The merger of 125, 144 and 625 jets have been simulated and the uniformity and Mach number reduction of the corresponding liners have been investigated. During late stages of the implosion, the Mach number of 3D liners was about half of that of spherically symmetric liners. The uniformity of the liner is critical for the reduction of Rayleigh-Taylor instabilities in the target while maintaining a high Mach number in the liner is necessary to achieve high target compression rates.

Simulations support DOE program in High Energy Density Laboratory Plasmas, in particular the Plasma Liner Experiment at LANL.


R. Samulyak, P. Parks, L. Wu, Spherically symmetric simulation of plasma liner driven magnetoinertial fusion, 17 (2010), No 9

R. Samulyak, L. Wu, P. Parks, On the structure of plasma liners for plasma jet induced magneto inertial fusion, Physics of Plasmas, 2010. Submitted to Physics of Plasmas