Comparison of Preconditioning Strategies in Energy Conserving Implicit Particle in Cell Methods

Particle in cell plasma simulations [1, 13] are amongst the most successful in obtaining new physics and in running at top performance on supercomputers [2]. Recently, this success has attracted several new developments. Among them, we focus the attention here on implicit and semi-implicit PIC methods. Given their simplicity and flexibility, these PIC formulations have gained wide use as effective computational tools for the simulation of multiscale plasma phenomena in a variety of fields from space [19], astrophysics [30], laser-plasma interaction [29], magnetic fusion [6] and space propulsion [7]. The interest on (semi-)implicit PIC has been intensified by the presentation of exactly energy conserving schemes [27, 10]. Semi-implicit methods had been previously developed along two lines of research (e.g. [3, 20] for a review): moment implicit [4] and direct implicit [17]. In these methods, the coupling between particles and fields is represented via a linear response: the direct implicit method uses a sensitivity matrix while the implicit moment method uses a Taylor expansion similar to the Sonine polynomial expansion of the Chapman-Enskog moment method [8]. Both methods found very successful applications either in commercial codes (www.orbitalatk.com/lsp/) or in research codes used by the plasma physics community: Venus[4], Celeste[33, 18], Parsek2D [25] and iPic3D [26].

Very recently, a new variant of semi-implicit method, ECsim [21] has been proposed. The new method is based on constructing a so-called mass matrix that represents how the particles in the system produce a plasma current in response to the presence of magnetic and electric fields. This response is linear and the ECsim method belongs to the class of semi-implicit methods. However, its key difference compared with the previously mentioned semi-implicit methods is that its mass matrix formation leads to exact energy conservation. A theorem proves the energy conservation to be valid at any finite time step[21], a conclusion confirmed to machine precision in practical multidimensional implementations[22].

Fully implicit methods had been deemed impractical in the past decades, but recent developments in handling non linear coupled systems made them practical. The original particle mover and field solver of the implicit moment method were already exactly energy conserving[27]. This property is however broken by the Taylor expansion. If instead the discretised field and particle equations are solved directly as a coupled non-linear system, energy is conserved exactly [27]. Several variants can include charge conservation besides energy conservation[10], for electromagnetic relativistic systems [23] and for the so-called Darwin pre-Maxwell approximation [9].

Even in our times, non-linear systems still are a formidable numerical task, especially for a large number of unknowns. For an efficient implementation, the numerical method used for solving the system needs to be carefully designed. So far, the research has focused on the Jacobian-Free Newton Krylov method (JFNK) [14] and on Picard iterations [32]. The performance of the non linear solver can be improved dramatically if good preconditioners are designed [16]. Preconditioners are operators that invert part of the problem, increasing the speed of convergence of the solution of the Jacobian problem needed for each Newton iteration.

In the context of implicit PIC, the attention has focused on the so-called particle hiding (or particle enslavement) approach [15, 27], especially suitable to hybrid architectures [11], and the use of fluid preconditioning for the fields [12].

In the present paper we revisit the particle hiding preconditioning and compare it with a number of new preconditioning methods able to significantly accelerate the iterative computation.

The first new approach proposed has been named field hiding after its similarities with the particle hiding method, with the main difference being its acting on the field solver instead of on the particle mover.

Second, we show that this method can be further optimized in PIC methods if the JFNK is completely replaced by a direct Newton solver where a Schwarz decomposed Jacobian is computed and inverted analytically, a method we refer to as Direct Newton-Schwarz (DNS).

Finally, all the implementations of the fully implicit method described above are compared with the new semi-implicit method ECsim.

All methods considered (including ECsim) conserve energy exactly, to machine precision and are virtually identical in the accuracy of the results produced for the cases tested. They only differ in the computing performance. We focus then on comparing the different implementations in terms of CPU time used. In particular, the analysis is conducted in the form of four tests to measure the computing performance while: 1) varying the time step, 2) varying the tolerance of the NK iterations, 3) varying the spatial grid resolution and 4) varying the number of particles.

We structured the paper as follows. Section 2 gives an introduction to the formulation of the energy conserving fully implicit method. Section 3 describes the preconditioning methods mentioned earlier together with a skeleton of their implementation in MATLAB. The results of the analysis are given in Section 4. Finally, Section 5 provides the discussion of the results and some important conclusion.

For 1D electrostatic problems, the ECFIM has been previously presented in terms of the Vlasov-Ampère system [27, 10], but it can be equivalently expressed in terms of the Vlasov-Poisson system [10], which is the approach followed in this work. We illustrate here the application in 1D but for all methods presented there is no barrier to the extension to multiple dimension, other than coding complexity.

The particle-in-cell method is based on introducing a new entity called computational particle, or super-particle, which is intended to represent a particular cluster of physical particles ($w_{p}$), having a charge $q_{p}=w_{p}q_{s}$, where $q_{s}$ is the charge of the physical particles of species $s$. Similarly, the mass is $m_{p}=w_{p}m_{s}$. The charge to mass ratio is physically correct and the (super)-particles still move with the same Newton equations of the physical particles. The core idea of the ECFIM is to use an implicit time centered particle mover:

where quantities under bar are averaged between the time step $n$ and $n+1$ (e.g. $\bar{{v}}_{p}=({v}_{p}^{n}+{v}_{p}^{n+1})/2$). The new velocity at the advanced time is then simply:

The electric field $\bar{E}_{p}$ is the electric field acting of the computational particle $p$ with position ${x}_{p}$ and velocity ${v}_{p}$ at the average time $\bar{{E}}_{p}=({E}_{p}^{n}+{E}_{p}^{n+1})/2$. Unlike the Vlasov-Ampère method, we rely here on the Poisson’s equation to compute the electric field:

where $\rho$ is the charge density. However, to reach an exactly energy conserving scheme, we reformulate the equation using explicitly the equation of charge conservation:

We take the temporal partial derivative of the Poisson’s equation (3) and substitute eq. (4)

Assuming now a time centered discretization, we express the variation of the scalar potential in one time step, $\delta\phi=\frac{\partial\phi}{\partial t}\Delta t$ as

and compute the new electric field as:

This equation can then be further discretized in space. We use here a staggered 1D discretization where the electric field and current are computed in the nodes of a uniform grid and the potential in the centers. In 1D this leads to $x_{c}=\{(i+1/2)*\Delta x\}$ for $i\in[0,N_{g}]$ and $x_{v}=\{i*\Delta x\}$ for $i\in[0,N_{g}+1]$ where $N_{g}$ is the number of cells. The fully discretized field equation is then:

where $v$ is the node to the left of the center $c$. The time dependent Poisson formulation is Galilean invariant and does not suffer from the presence of any curl component in the current [10]. In electrostatic systems the current cannot develop a curl because such curl would develop a corresponding curl of the electric field, and in consequence electromagnetic effects. The formulation above prevents that occurrence since it is based on the divergence of the current and any curl component of J is eliminated. This is not an issue in 1D but it would be in higher dimensions: the Vlasov-Poisson system can be more directly applied for electrostatic problems.

The peculiarity of PIC methods is the use of interpolation functions $W\left({x}_{g}-\bar{x}_{p}\right)$ ($g$ the generic grid point, center or vertex) to describe the coupling between particles and fields. For example, in the case of a the Cloud-in-Cell approach, the interpolation function reads  [13]

The interpolation functions can then be used to compute the electric field acting on each computational particles as

where ${E}_{g}$ is the electric at the vertices of the grid. Likewise, the current is computed as

The set of equations for particles (1) and fields (8) is non-linear and coupled. The coupling comes from the need to know the electric field at the advanced time before the particles can be advanced and in the need to know the new particles position and velocity before the current can be computed and the field advanced. The set of equations is also non-linear due to the non-linear dependence of the currents and fields on the particle positions. A non-linear solver is then needed to address the problem.

The solution provided describes an exactly energy conserving scheme [27, 10]: energy will be conserved for any choice of the time step $\Delta t$ at the level of accuracy chosen in the interactive solution of the system equations.

The goal of the study is to compare under identical straightforward style of coding different strategies for conducting the non-linear iteration required by the fully implicit method. To this end, we rely on a MATLAB implementation that we include in the appendices. Naturally, this is a limited scope: more sophisticated implementations in other languages for more complex problems would likely be needed for real world applications. Our aim here is to consider the promise shown in the simplest cases by each of the following strategies for the ECFIM:

1. 1.

   Jacobian Free Newton Krylov (JFNK) method without any preconditioning.

2. 2.

   JFNK with particle hiding, also called particle enslavement in the literature [27, 10].

3. 3.

   JFNK with field hiding

4. 4.

   Direct Newton approach with field hiding and using a partial Jacobian computed in a Schwarz particle-by-particle decomposition. We refer to this new method as Direct Newton-Schwartz (DNS).

These implementations are additionally compared with the ECSIM method [21, 22], which does not require any non linear iteration and is used here as comparative benchmark.

In conclusion, we would like to point out that the results obtained here in a simple case might be overturned in more complex cases using more advanced implementations. Nevertheless, we believe this study may be considered as a hitchhiker guide for more complex design choices in production codes. Below we report the details of the implementation of each of the ECFIM schemes being compared.

To evaluate the performance of the implementations discussed above, we present a parametric study organized in four tests:

• •

  the first test measures the performance with a different time step amplitude,

• •

  the second test measures the performance with a different tolerance of the Newton method,

• •

  the third test measures the performance with a different resolution (i.e. cell size),

• •

  the fourth test measures the performance with a different number of particles.

In particular, we set the analysis around an electron two stream instability in a neutral plasma featuring cold and fixed ions, with an electrons drift velocity $V_{e}=0.2\,c$. The study of the two-stream instability is one of the first initial approaches normally taken to test a new plasma-related code. Explanations on how this instability works can be found in plasma-physics books, such as in [1]. Particles are distributed in a 1D domain with size $L=20\pi\,c/{\omega}_{pe}$. Furthermore, periodic boundary conditions are assumed and the system is let evolve for 200 cycles. To make a clearer illustration, since each parametric study performs the change of a single parameter, it is convenient to establish the parameters of the reference set-up: the tolerance of the Newton and GMRES algorithms are equal to $10^{-10}$ (where present), ${\omega}_{pe}\Delta t=0.1$, $N_{p}=10^{5}$ and $N_{g}=64$. Anywhere not mentioned otherwise the corresponding parameter equals the one in the reference set-up. All tests include the comparison of 5 different methods:

1. 1.

   ECFIM with no preconditioner,

2. 2.

   ECFIM with particle-hiding preconditioner,

3. 3.

   ECFIM with field-hiding preconditioner,

4. 4.

   the Direct Newton Schwarz (DNS) method,

5. 5.

   the semi-implicit code ECSIM.

An example of the output obtained from running the unpreconditioned version of the code is given in Fig. 1, where signatures of the two stream instability are clearly visible in the particles phase-space (upper panel), together with the energy evolution and its relative error respect to the initial value (lower panels). This test featured a non-linear iteration tolerance of $10^{-14}$, $10000$ particles and $256$ cells.

In this paper we presented a set of different strategies to improve the computing performance of fully implicit Particle-in-Cell (PIC) codes. Unlike semi-implicit PICs, which make use of a linear approximation to improve performance, we consider here no approximations of the original formulation and focus only on advanced non-linear solvers that still produce the same exact result. Due to the heavy CPU and memory consumption requested by the iteration of the entire non-linear system, the fully implicit approach has long been considered expendable over the semi-implicit counterpart. This concern, however, is largely eliminated when preconditioners are used. In this work we proved that the computation with fully implicit PICs can be significantly improved whether a careful application of ad-hoc preconditioning methods be applied. The semi-implicit method remains the fastest, but the fully implicit methods come now close second.

The first preconditioner method considered here is found in previous literature [27] and named particle-hiding after its approach of segregating the particle mover from the field solver in the non-linear JFNK iteration. Particles are still continuously updated within the iterative cycle using the fields information, but now the particles velocity is no longer considered as part of the residual convergence, consequently leading to a lower number of iterations and a strong reduction of the memory demanded due to the non storage of the particles information. The downside is that for each residual evaluation, all particles still need to be moved, so that this approach still relies on the full mutual coupling between particles and fields, which conducts to further computational costs.

We have then proposed a further improvement based upon the fact that are actually the fields being the coupling source in the system rather than the particles. Particles do not directly interact with each other, their interaction is provided by the field. Mathematically, the fields are coupled with each particles while the particles just need the updated field information to compute their new position and velocity. The approach is based on a complementary variant of the previous method, and is named field-hiding after its application to the field information instead of particles’. It is now the field solver to be segregated out of the residual computation, so that the iteration is only performed over the particles velocity. Having removed all direct couplings from the JFNK calculation and having relegated all field couplings to the preconditioner, this approach is shown to perform faster then the previous one. However, the memory consumption downside is yet present, as the method requires a memory allocation as large as the number of particles considered in the system, which may become very heavy for real-system simulations.

In order to further relax this issue, a different preconditioner based on the field-hiding approach has been presented. Given the independence of the particles from each other, the Jacobian matrix appears to have all the cross-over derivatives equal to zero. We can then consider to eliminate the use of the Krylov space and resolve the linear part of the Newton iteration using a Schwarz decomposition. The JFNK is then no longer needed, each particles is solved independently via a direct newton iteration based on the Schwarz decomposed Jacobian computed and inverted analytically particle by particle. In doing so, we combine the two main advantages of the two previous preconditioners: 1) not having to deal with a large dimensional array and 2) not dealing with the main couplings directly inside the non-linear iterative loop.

This latter approach is shown to lead to a lower number of iterations required in the non-linear solver and to the fastest CPU performance. From the tests presented, we can in fact confirm that the Direct-Newton-Schwarz (DNS) method is the fastest to converge. The DNS also eliminates the need to reserve memory for the particle iterative solution as each particle is solved independently and directly. DNS comes second only to the semi-implicit ECSIM method, which yet relies on a completely different approximation approach.

In conclusion, we observe how the field-hiding preconditioner approach is reported here as the most promising, especially when coupled with a DNS formulation, avoiding Krylov subspace solvers altogether. In fact, all the tests with different number of particles, different mesh, different Newton iteration tolerance and different time step have shown that the field-hiding approach, and especially its DNS implementation, is the best fully non-linear system. While the ECsim semi-implicit approach is still the fastest, we note that its extension to relativistic system does not appear to be straight-forward, while that of the DNS fully implicit method is.