This paper discusses a novel fully implicit formulation for a one-dimensional electrostatic particle-in-cell (PIC) plasma simulation approach. Unlike earlier implicit electrostatic PIC approaches (which are based on a linearized Vlasov-Poisson formulation), ours is based on a nonlinearly converged Vlasov-Ampère (VA) model. By iterating particles and fields to a tight nonlinear convergence tolerance, the approach features superior stability and accuracy properties, avoiding most of the accuracy pitfalls in earlier implicit PIC implementations. In particular, the formulation is stable against temporal (Courant-Friedrichs-Lewy) and spatial (aliasing) instabilities. It is charge-and energy-conserving to numerical round-off for arbitrary implicit time steps (unlike the earlier "energy-conserving" explicit PIC formulation, which only conserves energy in the limit of arbitrarily small time steps). While momentum is not exactly conserved, errors are kept small by an adaptive particle sub-stepping orbit integrator, which is instrumental to prevent particle tunneling (a deleterious effect for long-term accuracy). The VA model is orbit-averaged along particle orbits to enforce an energy conservation theorem with particle sub-stepping. As a result, very large time steps, constrained only by the dynamical time scale of interest, are possible without accuracy loss. Algorithmically, the approach features a Jacobian-free Newton-Krylov solver. A main development in this study is the nonlinear elimination of the new-time particle variables (positions and velocities). Such nonlinear elimination, which we term particle enslavement, results in a nonlinear formulation with memory requirements comparable to those of a fluid computation, and affords us substantial freedom in regards to the particle orbit integrator. Numerical examples are presented that demonstrate the advertised properties of the scheme. In particular, long-time ion acoustic wave simulations show that numerical accuracy does not degrade even with very large implicit time steps, and that significant CPU gains are possible.

No abstract

No abstract

This paper presents an initial study that is intended to explore the development of a scalable fully-implicit stabilized unstructured finite element (FE) capability for low-Mach-number resistive MHD. The discussion considers the development of the stabilized FE formulation and the underlying fully-coupled preconditioned Newton-Krylov nonlinear iterative solver. To enable robust, scalable and efficient solution of the large-scale sparse linear systems generated by the Newton linearization, fully-coupled algebraic multilevel preconditioners are employed. Verification results demonstrate the expected order-of-acuracy for the stabilized FE discretization of a 2D vector potential form for the steady and transient solution of the resistive MHD system. In addition, this study puts forth a set of challenging prototype problems that include the solution of an MHD Faraday conduction pump, a hydromagnetic Rayleigh-Bernard linear stability calculation, and a magnetic island coalescence problem. Initial results that explore the scaling of the solution methods are presented on up to 4096 processors for problems with up to 64M unknowns on a CrayXT3/4. Additionally, a large-scale proof-of-capability calculation for 1 billion unknowns for the MHD Faraday pump problem on 24,000 cores is presented.

The conceptual development and implementation of a scalable nonlinear solver for the time-dependent three-dimensional (3D) compressible resistive magnetohydrodynamics model (MHD) is discussed. The approach is based on Jacobian-free Newton–Krylov technology, preconditioned with multigrid methods for algorithmic scalability. The key to the approach is the reformulation of the hyperbolic MHD system into a parabolic one, which is amenable to multigrid techniques. Such reformulation (parabolization) aims to render the modified system block diagonally dominant (unlike the original MHD system, which is diagonally submissive for implicit time steps Δt larger than the explicit Courant–Friedrichs–Lewy time step). The algorithm has been tested on a variety of 2D and 3D configurations which demonstrate its excellent algorithmic scalability properties. In particular, it is shown that, serially, the CPU time for a given simulation scales linearly with the number of unknowns, and that large CPU gains (∼30) are attainable with respect to an explicit approach even for moderate grids (256×256). In parallel, the algorithm features excellent scalability properties up to thousands of processors (4096) and millions of unknowns (∼1.3×108).