Next Article in Journal
Quantum Configuration and Phase Spaces: Finsler and Hamilton Geometries
Next Article in Special Issue
The Impact of Radio Frequency Waves on the Plasma Density in the Tokamak Edge
Previous Article in Journal
Acknowledgment to the Reviewers of Physics in 2022
Previous Article in Special Issue
Probability Distribution Functions of Solar and Stellar Flares
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Advances in the Implementation of the Exactly Energy Conserving Semi-Implicit (ECsim) Particle-in-Cell Method

University of Leuven Center for mathematical Plasma Astrophysics (CmPA), Department of Mathematics, Katholieke Universiteit Leuven (KU Leuven), 3000 Leuven, Belgium
Physics 2023, 5(1), 72-89; https://doi.org/10.3390/physics5010007
Submission received: 6 November 2022 / Revised: 14 December 2022 / Accepted: 19 December 2022 / Published: 18 January 2023

Abstract

:
The energy-conserving semi-implicit (ECsim) method presented by the author in 2017, is a particle-in-cell (PIC) algorithm for the simulation of plasmas. Energy conservation is achieved within a semi-implicit formulation that does not require any non-linear solver. A mass matrix is introduced to linearly express the particle-field coupling. With the mass matrix, the algorithm preserves energy conservation to machine precision. The construction of the mass matrix is the central nature of the method and also the main cost of the computational cycle. Here, three methods that modify the construction of the mass matrix are analyzed. First, the paper considers how the sub-cycling of the particle motion modifies the mass matrix. Second, a form of smoothing that reduces the noise while retaining exact energy conservation is introduced. Finally, an approximation of the mass matrix is discussed that transforms the ECsim scheme to the implicit moment method.

1. Introduction

The energy-conserving semi-implicit (ECSim) method is an algorithm for plasma simulation based on the particle-in-cell (PIC) approach [1]. PIC methods can be explicit, semi-implicit or fully implicit. Plasmas are governed by two sets of equations: the equations for the motion of the particles and the equations for the evolution of the fields. The two sets of equations are coupled because the field equations need the sources (current and charge) from the particles, and the particles need the fields to compute the force. As a particle moves, the fields are modified, and as the fields change, the forces on the particles are modified. This link is central to the physics of plasmas: plasmas are collective sets of particles interacting via the fields. The coupling between particles and fields is non-linear, and to represent it in discretized equations in its fullness, one needs fully implicit methods. In fully implicit methods, the particle equations of motion and the field equations are solved together within a non-linear solver, such as the Newton–Krylov approach [2,3]. In explicit methods, conversely, the coupling between particles and fields is suspended for a small time step [4]. In that small time interval, one assumes that the known fields can be used unchanged for moving the particles, and the particle information can be used unchanged to evolve the fields. This has three major consequences.
First, the explicit method is straightforward with no iteration is needed, and explicit PIC can be implemented as one of the most efficient algorithms known in computer science, consistently being a top performance achiever on any new computer architecture introduced. For example, PIC was one of the first applications to reach petascale performance [5]. Implicit PIC is much more complex in its implementation. Especially on massively parallel computers, reaching high efficiency is a challenge.
Second, in explicit PIC, the time step becomes limited by numerical stability considerations, requiring the use of high resolution. The peculiarity of PIC is that the resolution needs to not just be refined in time to resolve the electron plasma frequency: ω p e Δ t < 2 [6], where ω p e denotes the electron plasma frequency and Δ t is the time step, but also in space to avoid the so-called finite-grid instability [4]. This limitation is removed by the implicit approach that allows one to select grid spacing and a time step based on the accuracy needed and not on the stability of the numerical algorithm [7].
Third, in explicit PIC, energy is not conserved. Using a good resolution, energy is acceptably maintained. There is always a secular trend of energy increase [4], but as the resolution is relaxed closer to the stability limit of the finite grid instability, the energy increase becomes more severe, until at the instability limit, it starts to grow exponentially. This effect cannot be avoided, but it can be improved by using smoothing and higher-order interpolation techniques. Recent structure-preserving geometric particle-in-cell methods use symplectic integrators to ensure local energy conservation at small time steps [8]. The implicit PIC method, instead, conserves energy exactly, whatever resolution is used [2,3]. This feature is physically important and practically impactful. Physically, of course, confidence comes from using an algorithm that preserves one of the most established properties in physics: conservation of energy. If energy starts to spontaneously increase, confidence in the results is shaken. Practically, lack of energy conservation requires a tedious and careful tuning of the parameters to make sure the simulation does not increase its energy excessively, leading in some situations to excessively resolved models that need to use much more resolution than the processes of interest require.
The semi-implicit PIC method tries to make a compromise and retain some of the advantages of both approaches. In semi-implicit methods, the particles and the fields are still advanced together, and an iteration is needed, but the coupling is linearized, and the iteration uses linear solvers. Different methods are used to linearize the coupling. The implicit moment method formulates the particle response to changes in the fields using the moment closure method [9]. The direct implicit method uses a formal linear expansion of the coupling operator [10]. In all these approaches, the stability properties of the semi-implicit method are superior to explicit methods and allow a good compromise in the resolution needed [7]. However, energy is not conserved. Unlike explicit PIC, energy can either increase or decrease depending on the implementation because dissipation terms are included in the algorithm to suppress energy growth.
The ECsim approach is the first semi-implicit PIC to retain exact energy conservation as in the fully implicit PIC. ECsim uses a mathematical construct called a mass matrix to express the coupling between particles and fields. With the mass matrix, the coupling is linear, but energy is conserved exactly. Below, it is reviewed how this was achieved [1].
Since its recent introduction, ECsim has found application in a number of applications in space [11,12,13,14,15] and in fusion [16,17]. An important improvement has removed the lack of charge conservation in the original scheme [18,19]. The current paper reports some extensions of the method that can widen its practical applicability.
First, a method to introduce smoothing to reduce noise while retaining exact energy conservation is described. Smoothing is an affordable way to effectively introduce a higher-order interpolation scheme. It removes the high frequency part of the spectrum. Without special attention, smoothing will tend to break energy conservation by just removing high-frequency fluctuations from the system. Instead, a method that conserves energy is presented.
Second, sub-cycling is a convenient approach in plasma simulations to address the faster scales seen by the particles. In some applications, the particle response (or that of a subset of the particle population) is faster than the evolution of the fields, and it is beneficial to move the particles several times without advancing the fields. In explicit PIC, the operation is straightforward since fields and particles are not advanced together. In implicit and semi-implicit PIC, instead, sub-cycling also needs to be done in the coupling, requiring the modification of the algorithm to compute the current used in the field solution.
Finally, we illustrate how the mass matrix formulation opens up the opportunity to approximate the mass matrix in certain limit cases to reduce the cost of the simulation.
The paper is organized as follows. Section 2 recaps the key properties of the ECsim needed for the discussion. Section 3 introduces how smoothing can be implemented to reduce noise while retaining the property of energy conservation. Section 4 is dedicated to the algorithm for sub-cycling the particle motion within the ECsim scheme. Section 5 derives a limit case of the mass matrix formulation that transforms the ECsim algorithm into the standard implicit moment method [9]. Results are presented in Section 6, and final conclusions are drawn in Section 7.

2. Summary of the ECsim Method

The ECsim method [1] is based on a formulation of the mover similar to the classic leap-frog scheme:
x p n + 1 / 2 = x p n 1 / 2 + Δ t v p n , v p n + 1 = v p n + q p Δ t m p E n + θ ( x p n + 1 / 2 ) + v ¯ p × B n ( x p n + 1 / 2 ) ,
where x p denotes the particle (labeled with p index) space position, v p the particle velocity, q p is particle charge, m p is the particle mass, n denotes the number of a particle, and E and B are electrical and magnetic fields, respectively, v ¯ p = ( v p n + 1 + v p n ) / 2 is the averaged speed, and E n + θ = θ E n + 1 + ( 1 θ ) E n .
This mover differs from both the explicit leap-frog mover [20] and the implicit θ -mover [21]: it combines the first equation from the explicit mover with the second equation from the implicit mover, but with an important difference: the electric and magnetic fields are computed at the known position x p n + 1 / 2 rather than at the unknown position x ¯ p . The important consequence is that the particle equations can be solved directly without any iteration needed among themselves. Instead, in the standard θ -mover, a predictor-corrector iteration is required. Yet, the mover is still implicit because the new fields are not known until the field equations are solved. The ECsim method retains the coupling between advanced fields and advanced particles, requiring the solution of a linear coupled system. However, the mover itself does not require any iteration, a substantial simplification.
In the ECsim scheme, the electric and magnetic fields are computed at the known position x p n + 1 / 2 . In the θ -scheme, instead, they are computed at the unknown position x ¯ p . These two positions are conceptually the same; they express the particle position at the mid-time between the old and new evaluations of the velocity. However, one is computed explicitly, in the leap-frog sense, while the other is computed as part of a predictor-corrector iteration [22,23]. Both methods are second-order accurate, but the ECsim scheme is simpler to compute. This simplicity is not just a virtue in itself but leads to an important consequence: the simplicity allows us to formulate the coupling with the fields in a way that insures exact energy conservation without requiring non-linear iterations. The θ -scheme can be made energy conserving but at the cost of a fully non-linear iteration requiring a non-linear solver [2,3].The ECsim scheme allows exact energy conservation without requiring any non-linear iteration.
The properties of stability are determined by the field-particle coupling, and in this sense, the method is still implicit. For this reason, not requiring any non-linear iteration but still requiring a liner solver to deal with the field-particle coupling, the method is semi-implicit. This nomenclature is to distinguish it from the fully implicit method that requires the non-linear iteration.
Note that the force term is expressed using the magnetic field, B n ( x p n + 1 / 2 ) , at the initial time level, but the electric field is considered at the advanced intermediate level: E n + θ ( x p n + 1 / 2 ) . The reason for this choice is simplicity and the fact that the magnetic field does no work and using the old time level does not introduce any loss of energy conservation.
The coupling of particles and fields requires us to interpolate the fields to the particle positions:
E p n + θ = E n + θ ( x p n + 1 / 2 ) = g E g n + θ W ( x p n + 1 / 2 x g ) ,
B p n = B n + θ ( x p n + 1 / 2 ) = g B g n W ( x p n + 1 / 2 x g ) .
A generic index g is used for the grid. In the specific implementation within the iPic3D code [24,25], the electric field and magnetic field are not colocated, and g labels either centers (for B ) or vertices (for E ). Here, the notation is simplified: B p n B n ( x p n + 1 / 2 ) and E p n + θ E n + θ ( x p n + 1 / 2 ) . In the implementation here, the interpolation function, W, is a b-spline of order = 1 [26]:
W ( x p x g ) = b ( x p x g ) b ( y p y g ) b ( z p z g ) .
This expression reduces trivially in 1D (one-dimensional) for the examples reported below.
For Maxwell’s equation, the standard θ -scheme [24] is used:
g × E n + θ + 1 c B g n + 1 B g n Δ t = 0 , g × B n + θ 1 c E g n + 1 E g n Δ t = 4 π c J ¯ g .
The spatial operators in Equation (5) are discretized on the grid labeled by g introduced above.
The coupling of the field equations with the particles is expressed by the current for each species:
J ¯ s g = 1 V g p s q p v ¯ p W ( x p n + 1 / 2 x g ) ,
where the summation is over the particles of the same species, labeled by s, and V g is the cell volume.
As with the θ -mover, the velocity equation can be rewritten in the equivalent form [27],
v ¯ p = v ^ p + β s E ^ p ,
with
v ^ p = α p n v p n , E ^ p = α p n , E p n + θ
and the rotation matrix α p n given by
α p n = 1 1 + ( β s B p n ) 2 I β s I × B p n + β s 2 B p n B p n ,
where I is the dyadic tensor (a matrix with a diagonal of 1) and β s = q p Δ t / 2 m p (independent of the particle weight and unique to a given species). The elements of the rotation matrix are indicated as α p i j , n , with labels i and j referring to the 3 components of the vector space (x, y, z).
Substituting Equation (7) into Equation (6), one obtains without any approximation or linearization,
J ¯ s g = 1 V g p q p v ^ p W p g + β s V g p q p E ^ p n + θ W p g ,
where the notation is shortened, W p g = W ( x p n + 1 / 2 x g ) , and the summation is intended over all particles of species s.
Using Equation (8), the expression for the current becomes
J ¯ s g = J ^ s g + β s V g p q p α p n E p n + θ W p g ,
where
J ^ s g = 1 V g p q p v ^ p W p g .
Computing then the electric field on the particles by interpolation form the grid as in Equation (3), it follows that
J ¯ s g = J ^ s g + β s V g p g q p α p n E g n + θ W p g W p g .
Exchanging the order of summation, one obtains:
J ¯ s g = J ^ s g + β s V g g M s , g g E g n + θ ,
where the actor in the leading role of the ECsim scheme—the mass matrix [28] is defined:
M s , g g i j = p q p α p i j , n W p g W p g .
There are 3 v (where v is the number of velocity directions) mass matrices, and in matrix notation, they can be written as M g g , that is, without the indices i , j for the vector directions.
The mass matrices M s , g g are the most important aspect of the ECsim method and are also the most expensive part of the computation [25]. A number of symmetries can be used to reduce the cost. Speeding up the construction can be achieved using offloading to accelerator processors (e.g., graphical processing units, GPU) [29]. The mass matrices, as shown in Equation (14), provide an explicit linear link between the advanced current at the mid-point of the time step and the electric field at the advanced time. This linear relationship can be substituted into the discretized Maxwell’s Equation (5) to form a linear set of equations:
g × E n + θ + 1 c B n + 1 B n Δ t = 0 , g × B n + θ 1 c E n + 1 E n Δ t = 4 π c J ^ g + g M g g E g n + θ ,
where the total current is J ^ g = s J ^ s g , and the species summed mass matrices that are written by elements are:
M g g i j = s β s V g M s , g g i j .
The direct link provided by the mass matrix is analytically exact for the original set of discretized equations. Unlike the implicit moment method, where the equations have to be approximated by Taylor series expansion [27], here, the link is still exactly the same as in the original set of discretized equations. Having eliminated the need for any approximation or Taylor series expansion is the reason why ECsim conserves energy exactly.
Let us consider now how energy conservation can be shown in the case θ = 1 / 2 . This is important for the derivations below to consider what key steps enable energy conservation.
The inner product of the velocity from Equation (1) with the average speed, v ¯ p , gives, by summing over all particles,
1 2 p m p ( v p n + 1 ) 2 ( v p n ) 2 = Δ t p q p g v ¯ p · E ¯ g W p g ,
where the electric field is computed as an average consistent with the choice θ = 1 / 2 , and the magnetic field drops out as obvious from the properties of the cross-product.
Exchanging the summation over particles and cells leads to:
1 2 p m p ( v p n + 1 ) 2 ( v p n ) 2 = Δ t g J ¯ g · E ¯ g ,
where it is recognized that J ¯ g = p q p v ¯ p W p g .
Multiplying the first Equation (5) by B ¯ g and the second by E ¯ g and summing them leads to:
( B g n + 1 ) 2 ( B g n ) 2 2 c + ( E g n + 1 ) 2 ( E g n ) 2 2 c = Δ t 4 π c J ¯ g · E ¯ g + E g · g × B B g · g × E .
Assuming a mimetic grid discretization that preserves the continuum properties of the operators and summing over all grid points gives:
g ( B g n + 1 ) 2 ( B g n ) 2 4 π + g ( E g n + 1 ) 2 ( E g n ) 2 4 π = Δ t g J ¯ g · E ¯ g + c Δ t 4 π g g · ( E g × B g ) .
This conservation law states that the variation of the magnetic and electric energy, as measured on the grid, equals the amount exchanged with the particles and carried by the grid-discretized divergence of the Poynting flux. For energy to be conserved in the system, the energy exchange term on the particle equations (Equation (19)) needs to be identical to that on the field equations, (Equation (21)). This term is indeed identically equal to g J ¯ g · E ¯ g in both equations. Energy conservation is enforced exactly to round off.
In addition to guaranteeing physical conservation of energy, a cornerstone in any physical model, the existence of this conservation constraint also guarantees a form of non-linear stability of the discretized equations [1], expanding the stability of the semi-implicit method compared with the moment implicit scheme [11,16].

3. Smoothing with the Mass Matrix Formulation

Smoothing can be designed to be compatible with the energy conserving properties of the mass matrix. Only the electric field is chosen to be smoothed, since the magnetic field tends to be much smoother in PIC simulations and smoothing is not needed.
From the proof of energy conservation recapped above, it is clear that for energy conservation, the smoothing of the current must be done in the same way as that of the electric field in the mover. Starting from the mover and calling S g g the smoothing operator, a smoothed electric field on the grid is defined as:
E g S M = g S g g E g n + θ .
From Equation (22), the smoothed electric field acting on a particle can be computed as:
E p S M = g E g S M W ( x p n + 1 / 2 x g ) .
The second equation of motion, Equation (1), then uses the smoothed electric field as:
v p n + 1 = v p n + q p Δ t m p E p S M + v ¯ p × B p n ,
where B p n is still computed as above.
From Equation (24), one can again compute the current as
J ¯ s g = J ^ s g + β s V g g M s , g g E g S M .
The energy exchange term for the particles then becomes
1 2 p m p ( v p n + 1 ) 2 ( v p n ) 2 = Δ t s g J ^ s g + β s V g g M s , g g E g S M · E g S M .
Applying smoothing to the source term of the second of the Maxwell Equation (16), one obtains:
g × B n + θ 1 c E g n + 1 E g n Δ t = 4 π c s g S g g J ^ s g + β s V g g M s , g g E g S M ,
where the last term is expressed from Equation (25).
For the fields, the energy integral then becomes:
g ( B g n + 1 ) 2 ( B g n ) 2 4 π + g ( E g n + 1 ) 2 ( E g n ) 2 4 π c Δ t 4 π g g · ( E g × B g ) = Δ t s g g S g g J ^ s g + β s V g g M s , g g E g S M · E ¯ g .
For the two energy integrals to be the same, the exchange term seen by the particles must be equal to that seen by the fields. The right-hand sides of Equation (26) must then equal that of Equation (28):
g g S g g J ^ s g + β s V g g M s , g g E g S M · E ¯ g = g J ^ s g + β s V g g M s , g g E g S M · g S g g E ¯ g .
Switching g with g (just names) in the right-hand side, the equivalence above holds when the smoothing operator is symmetric (i.e., the matrix repenting it is symmetric), a common property shared by many smoothing operators [3].
Note that the electric field is smoothed but not the magnetic field, which tends to be less noisy by its nature.

4. Sub-Cycling with the Mass Matrix Formulation

A mass matrix can also be defined in the presence of sub-cycling or orbit-averaging movers. The velocity update of ECsim can be reformulated for sub-cycling as
v p ν + 1 = v p ν + q p Δ t ν m p E n + θ ( x p ν ) + v p ν + 1 + v p ν 2 × B n ( x p ν ) .
It is assumed assume that the time step Δ t between field updates is subdivided into N ν , not necessarily equal sub-steps Δ t ν .
The positions for the field evaluations, x p ν , during the sub-cycle can be computed in different ways. The simplest is to assume a straight orbit within Δ t , similar to the leap-frog approach:
x p ν = x p n 1 / 2 + v p n ν = 0 ν = ν Δ t ν ,
which can all be computed at once since the same velocity is used for all points along the trajectory. The first step starts from the old position x p ν = 0 = x x n and old velocity v p ν = 0 = v x n , and the last step leads to the final position x p ν = N ν = x p n + 1 and final velocity v p ν = N ν = v p n + 1 . The fields are assumed to be those computed at the time level θ within the field update time step Δ t .
Another promising approach is to recall that most often in plasma physics, particles are not moving in straight lines, but rather they are frozen into the field lines, moving in cyclotron orbits with drifts due to the inhomogeneity of the fields. In the spirit of gyro-averaging, certain applications of the implicit method might need to step over the gyration time scale, and the positions of the particles used in Equation (30) would then be chosen to achieve accurate gyro-averaging, for example taking N ν positions along the gyro-orbit of a particle [30] to compute an average force on the particle’s center of gyration.
The example of the two strategies above for computing the intermediate positions can be made in a single explicit step that generates all positions at once: in this case, each substep contribution to the mass matrix and the moments can be computed in parallel, greatly improving the parallel performance. In practice, the N ν operations required by the substepping algorithm can all be done in parallel in an embarrassingly parallel approach that scales ideally on supercomputers: no communication between the particles and between the substeps is needed.
Equation (30) can be inverted with the same vector manipulations used for Equation (1) to obtain:
v p ν + 1 + v p ν 2 = v ^ p ν + β s E ^ p ν ,
where hatted quantities have been rotated by the magnetic field computed at the location x p ν :
v ^ p ν = α p ν v p ν , E ^ p ν = α p ν E n + θ ( x p ν ) ,
via a rotation matrix, α p n , defined as in the case of a singe step, but with the magnetic field computed at the last substep position:
α p ν = 1 1 + ( β s B n ( x p ν ) ) 2 I β s I × B n ( x p ν ) + β s 2 B n ( x p ν ) B n ( x p ν ) .
From Equation (32), one can obtain directly from its definition (6) the mean current over all sub-steps without any further approximation or linearization:
J ¯ s g = 1 V g p q p ν Δ t ν Δ t v ^ p ν + β s E ^ p ν W p g ν ,
where W p g ν = W ( x p ν x g ) . Using now the definitions of the hatted quantities, Equation (33), one can cast the mean current in the same form as in the single step formulation:
J ¯ s g = J ^ s g + β s V g g M s , g g E g n + θ ,
but with the new definition for
J ^ s g = 1 V g p q p ν Δ t ν Δ t v ^ p ν W p g ν .
and the mass matrix for a sub-cycled trajectory defined by
M s , g g i j = p q p ν Δ t ν Δ t α p i j , ν W p g ν W p g ν .
Note that the definition of the sub-cycled mass matrix (Equation (38)) and the sub-cycled hatted current (Equation (37)) treats each sub-interval for each particle as if they were independent: the equations see each sub-interval as a particle. It is as if the system has N p N ν particles made by each particle for each sub-interval Δ ν . This feature lends itself to a simpler computing implementation where each particle is spawned into N ν treated as independent particles in the interpolation, mass matrix computation and current gathering step. This is a valuable approach in parallel and vectorized computer architectures, such as GPUs.
Regardless of how the positions for the particles during the sybcycling are chosen, if the current defined in Equation (36) is computed with the mass matrices defined in Equation (38), energy is conserved. In fact, the term g J ¯ g · E ¯ g is again identical when computed from the equations for the particles and for the fields.

5. Simplification of the Mass Matrix: The Limit of the Implicit Moment Method

There is a specific case where the mass matrix formulation takes a much simplified form: in the case of the interpolation of the nearest grid point (NGP). In that case, the interpolation function W p g is simple: 1 for the nearest grid point, which is labeled here as g p , and 0 for everywhere else:
W p g = δ g p g ,
where δ a b is the Kronecker delta.
In this case, the mass matrix becomes
M s , g g i j = p q p α p i j , n δ g p g δ g p g .
When substituted in the expression for the current, this leads to
J ¯ s g = J ^ s g + β s V g p g q p α p n E g n + θ δ g p g δ g p g .
Given the properties of Kronecker delta, the summation over g can be done first:
g E g n + θ δ g p g = E g p n + θ ,
and substituting:
J ¯ s g = J ^ s g + β s V g p q p α g p n E g p n + θ δ g p g ,
where the point that the α s are computed using the magnetic fields of the nearest grid point is used, consistent with the NGP interpolation. Using again the properties of Kronecker delta, the summation over the particles can be done directly:
1 V g p q p = ρ s g ,
to obtain
J ¯ s g = J ^ s g + β s ρ s g α g n E g n + θ .
This is the same expression that links the electric field and the current in the implicit moment method [31] used in Venus2D [9], in Celeste3D [23] and iPic3D [24]. It is then worth considering this limit expression for the mass matrix formulation. Most often, the NGP cannot be used in practice for its well known excessive noise, but the expression (45) still can be used even in the presence of other interpolation schemes.
Equation (45) is exact only for the NGP scheme where it still leads to exact energy conservation. When, instead, other orders of interpolation are used, energy conservation is lost, but it still provides a meaningful link between the current and electric field: it is, in fact, the same used for decades by the implicit moment method. If Equation (45) is used, ECsim becomes more similar to the implicit moment method, but it still differs in one key aspect: the interpolations are computed using the position provided by the leap-frog algorithm for the particle position. This is known explicitly and does not require any iteration. In the implicit moment method, the mover requires to iterate between the velocity update and position update using a predictor-corrector scheme. This is not required in ECsim, where the particle position is known explicitly from the previous time step.

6. Results

6.1. Effects of Sub-Cycling

To check how the innovations described above perform in practice, let us first consider the effect of sub-cycling. The goal is twofold. First, one needsto verify that the energy is indeed conserved and sub-cycling does not break the energy conservation. Second, one tests on a specific case how far one can push the sub-cycling before the correctness of the results deteriorates. The first task above is just a confirmation of the rigorously exact calculations above, and it is merely a verification test for the code implemented here. The second task is only an illustration because the usefulness of sub-cycling is highly problem-dependent.
A two-stream instability was initiated by considering a domain of L / d e = 2 π divided into 64 cells with 10,000 particles. A mini-app was useful to readily study implementations in different computer architectures: it was intended merely as a test, not as a real world problem. Space was normalized in electron skin depth, d e = c / ω p e , and time was normalized in electron plasma frequency, ω p e t . The ions were kept immobile, while the electrons had a thermal speed of v the / c = 0.02 (with c denoting the speed of light), and the two beams had a net speed of V 0 / c = ± 0.1 . The reference case without sub-cycling ( N ν = 1 ) has a CFL (Courant–Friedrichs–Lewy) condition Courant number, C = V 0 Δ t / Δ x = 0.1 . We then proceeded to carry out 10 runs with N ν = 1 –10, and Δ t was increased by the same factor. That is, the particle Δ t was kept fixed and the Δ t of the field was increased. As the number of sub-cycles increased, the Δ t increased, but the total time remains the same: ω p e t = 50 . A perturbation of the particle velocity was added with mode number, m = 5 :
u p = u p + V 1 sin 2 π m x p L ;
with an amplitude V 1 = V 0 / 10 .
Here, a MATLAB implementation was used that was tested on a MAC OSX 10.14.6 with an Intel Core i7 2,6GHz processor complemented with 16GB DDR3 1600 MHz memory.
The results of the study are summarized in Table 1 and in Figure 1 and Figure 2.
Table 1 reports the reduction of the cost of the simulation as the number of sub-cycles increased. The MATLAB implementation is well vectorized and efficient in handling the particle projections, so much so that as the number of sub-cycles is increased, the cost per Δ t hardly increases at all. If one compares the cost of a time step without sub-cyclig (0.08395 s) with that with 10 sub-cycles (0.0958 s), the increase is minimal even though particles have been moved and projected to the moments and the construction of the mass matrix 10 times more. While this result is specific to the vectorization done via MATLAB, other forms of vectorization are also possible on GPUs using OpenACC [32] and CUDA [33] programs opening up a similar opportunity for this type of optminization in other implementations.
As required by the theoretical derivations above, energy should be conserved in all cases. This is indeed the case in Table 1, a confirmation that the MATLAB implementation used is bug-free.
However, if sub-cycling is beneficial in reducing the cost of the simulation, it introduces a degradation in the physics sense. Particles are moved with the same time step in each sub-cycle regardless of N ν , but as N ν increases, the overall time step Δ t for the field recalculation is increased. In a number of problems, it is still beneficial to accept this compromise. Figure 1 and Figure 2 show the degradation of the results with N ν . Unfortunately, the instability modeled is a complex non-linear process, and the linear phase is not important: an easy quantitative metric of accuracy is not available. The study was already started from a significant perturbation. If not added then several modes develop at the same time, and still, the linear phase is muddied by the interaction of many modes. There is no single metric that can easily summarize the quality of the evolution. One has to rely on qualitative visual comparison. For this reason, let us look carefully at the velocity space distribution (Figure 1) and phase space (Figure 2). As N ν is increased, the correct physics description is progressively lost.
The most characteristic feature of the two-stream instability is the formation of a flat top distribution, a distribution where, for a range of velocities, the distribution remains flat. This feature is prominent without sub-cycling, and it is still reasonably well represented at higher sub-cycling numbers. However, the tail of higher-energy particles is reduced when sub-cycling is increased, a reflection of the fact that the electric field responsible for particle acceleration [34] is less accurately computed. These are non-physical artifacts of excessive sub-cycling.
The same conclusion is reached by analysing phase space. In this case, the most important feature is the formation of electron holes, regions of phase space depleted of electrons and in fact completely void of them. These features, which are also observed in experimental measurements, are distorted and expelled to the edge of the velocity range as N ν is increased.
Sub-cycling is compatible with the mass matrix formulation, retaining the property of exact energy conservation. It can be implemented effectively enough, only minimally increasing the cost of the computational cycle despite the increase in the number of particle operations while reducing the number of time steps for the same total time. However, it must be used with care because not updating the fields after the particles are moved in a sub-cycle introduces physical errors even though energy is still exactly conserved.

6.2. Effects of Smoothing

To test the effects of smoothing, let us consider another type of streaming instability: the transverse electromagnetic instability driven by counter streaming beams [35]. Again, a 1D plasma (x is the only spatial variable) is considered, but now, all three components of the particle velocity are taken into account. The two counter-streaming beams are directed along y . Using the speed of light for normalization, a case with the beams having v the / c = 0.01 counter-streaming with speed v 0 y / c = 0.2 to be considered. No initial perturbation is added and the natural noise is let to initiate the instability. The transverse electromagnetic streaming instability is common is often used in the context of understanding how magnetic fields can be generated in the universe [36] a form of dynamo. A setup similar to that reported by Innocenti et al. [37] is used. A characteristic of this instability is that it segregates in phase space particles with opposite signs of v y , initially residing in the two different beams.
Figure 3 shows the evolution of the phase space cross-section ( x , v x ) : the red particles have v y < 0 and the blue particles have v y > 0 . The evolution initially retains their separation, but in time, phase-space mixing takes over. More details of the physics of this simulation can be found by Innocenti et al. [37]. Let us focus here on comparing the normal ECsim case without smoothing with one where the smoothing kernel, S = [ 1 / 4 , 1 / 2 , 1 / 4 ] , is applied by convolution 3 times. Figure 3 compares the smoothed and non-smoothed case. Smoothing does not alter the evolution in any profound ways, but it affects it. Since the simulation is started from its natural noise, no two simulations will be identical, and of course, the smoothed simulation will have less noise. The locations of the islands formed in phase space are not the same, but the overall features are similar.
The energy exchange is also similar in the smoothed and non-smoothed runs, as shown in Figure 4. In particular, the kinetic energy is lost primarily to produce magnetic energy. The transverse counter-streaming instability is a form of magnetic dynamo that spontaneously creates a magnetic field by using the kinetic energy of the counter-streaming beams. The total energy is conserved to machine precision in both runs. Smoothing, as shown in Section 3 theoretically, indeed conserves energy exactly.
The main effect of smoothing is to eliminate the high-frequency part of the spectrum. Figure 5 and Figure 6 compare the k ω spectrum of selected fields. The noise at high values of k is reduced. The spectrum identifies the low k part of the spectrum as dominant. The characteristic arch of the electromagnetic waves (light waves) is also prominent.

7. Discussion

The ECsim method allows to make a PIC simulation within a semi-implicit approach that conserves energy exactly. This feature, besides having its own intrinsic value, also leads to improved stability, allowing, for example, to consider realistic plasma conditions in the heliosphere, including the colder solar wind where the Debye length is very small compared with the scales of interest.
This paper reports on three new developments.
First, a method for smoothing the electric field in an algorithm that preserves energy conservation is introduced. When it can be avoided, smoothing should be avoided, but when there is a need for it, the algorithm presented here achieves smoothing without breaking energy conservation. This can be, for example, the case when excessive noise alters the correct transport properties of a plasma, an issue of great importance, for example, in fusion energy studies [38].
Second, a method for computing the mass matrix in the presence of particle sub-cycling is introduced. Sub-cycling can be useful in a number of situations. Most often, if it can be avoided, it should be avoided because particles and fields move together, with the frozen-in condition being a cardinal property of plasmas. However, frozen-in properties are valid at the large scales typical of MHD. There are many situations where the fields evolve slowly while particles move quickly. For example, in gyro motion. Sub-cycling can be used together with gyroaveraging. The method presented here allows us to use sub-cycling while constructing a mass matrix that continues to preserve the exact energy conservation.
Finally, a limit case when the mass matrix calculation becomes especially simple is discussed and it is shown that, in this limit, the plasma particle response in the ECsim method becomes identical to that of the implicit moment method (IMM). This has two implications. First, it allows us to use the same code either in the full energy-conserving ECsim mode or in the less expensive IMM mode. Comparisons can then be made more readily. Second, the derivation clarifies the theoretical links between ECsim and IMM, revealing in what limit the two become identical.
These three new steps have both a theoretical significance and practical value, bringing understanding of the properties of ECsim and broadening its range of applications.

Funding

The work reported received funding from the Katholieke Universiteit Leuven (KULeuven) Bijzonder Onderzoeksfonds (BOF) under the C1 project TRACESpace and from the European Union project DEEP-SEA (grant agreement 955606).

Data Availability Statement

The results presented in this article were obtained using the codes available at http://www.ecsimpic.org (accessed on 5 November 2022).

Conflicts of Interest

The author declares no conflict of interest.

References

  1. Lapenta, G. Exactly energy conserving semi-implicit particle in cell formulation. J. Comput. Phys. 2017, 334, 349–366. [Google Scholar] [CrossRef] [Green Version]
  2. Markidis, S.; Lapenta, G. The energy conserving particle-in-cell method. J. Comput. Phys. 2011, 230, 7037–7052. [Google Scholar] [CrossRef] [Green Version]
  3. Chen, G.; Chacón, L.; Barnes, D.C. An energy-and charge-conserving, implicit, electrostatic particle-in-cell algorithm. J. Comput. Phys. 2011, 230, 7018–7036. [Google Scholar] [CrossRef] [Green Version]
  4. Birdsall, C.K.; Langdon, A.B. Plasma Physics via Computer Simulation; CRC Press/Taylor & Francis Group: Boca Raton, FL, USA, 1991. [Google Scholar] [CrossRef]
  5. Bowers, K.J.; Albright, B.J.; Yin, L.; Daughton, W.; Roytershteyn, V.; Bergen, B.; Kwan, T.J.T. Advances in petascale kinetic plasma simulation with VPIC and Roadrunner. J. Phys. Conf. Ser. 2009, 180, 012055. [Google Scholar] [CrossRef]
  6. Hockney, R.W.; Eastwood, J.W. Computer Simulation Using Particles; CRC Press/Taylor & Francis Group: Boca Raton, FL, USA, 1988. [Google Scholar] [CrossRef]
  7. Lapenta, G. Particle simulations of space weather. J. Comput. Phys. 2012, 231, 795–821. [Google Scholar] [CrossRef]
  8. Xiao, J.; Qin, H.; Liu, J. Structure-preserving geometric particle-in-cell methods for Vlasov-Maxwell systems. Plasma Sci. Technol. 2018, 20, 110501. [Google Scholar] [CrossRef] [Green Version]
  9. Brackbill, J.U.; Forslund, D.W. An implicit method for electromagnetic plasma simulation in two dimensions. J. Computat. Phys. 1982, 46, 271–310. [Google Scholar] [CrossRef]
  10. Langdon, A.B.; Cohen, B.I.; Friedman, A. Direct implicit large time-step particle simulation of plasmas. J. Computat. Phys. 1983, 51, 107–138. [Google Scholar] [CrossRef]
  11. Lapenta, G.; Gonzalez-Herrero, D.; Boella, E. Multiple-scale kinetic simulations with the energy conserving semi-implicit particle in cell method. J. Plasma Phys. 2017, 83, 705830205. [Google Scholar] [CrossRef] [Green Version]
  12. Walker, R.J.; Lapenta, G.; Berchem, J.; El-Alaoui, M.; Schriver, D. Embedding particle-in-cell simulations in global magnetohydrodynamic simulations of the magnetosphere. J. Plasma Phys. 2019, 85, 905850109. [Google Scholar] [CrossRef] [Green Version]
  13. Zhou, H.; Tóth, G.; Jia, X.; Chen, Y.; Markidis, S. Embedded kinetic simulation of Ganymede’s magnetosphere: Improvements and inferences. J. Geophys. Res. Space Phys. 2019, 124, 5441–5460. [Google Scholar] [CrossRef]
  14. Lapenta, G.; El Alaoui, M.; Berchem, J.; Walker, R. Multiscale MHD-kinetic PIC study of energy fluxes caused by reconnection. J. Geophys. Res. Space Phys. 2020, 125, e2019JA027276. [Google Scholar] [CrossRef]
  15. Lapenta, G.; Schriver, D.; Walker, R.J.; Berchem, J.; Echterling, N.F.; El Alaoui, M.; Travnicek, P. Do we need to consider electrons’ kinetic effects to properly model a planetary magnetosphere: The case of Mercury. J. Geophys. Res. Space Phys. 2022, 127, e2021JA030241. [Google Scholar] [CrossRef]
  16. Gonzalez-Herrero, D.; Micera, A.; Boella, E.; Park, J.; Lapenta, G. ECsim-CYL: Energy Conserving Semi-Implicit particle in cell simulation in axially symmetric cylindrical coordinates. Comput. Phys. Commun. 2019, 236, 153–163. [Google Scholar] [CrossRef] [Green Version]
  17. Park, J.; Lapenta, G.; Gonzalez-Herrero, D.; Krall, N.A. Discovery of an electron gyroradius scale current layer: Its relevance to magnetic fusion energy, Earth’s magnetosphere, and sunspots. Front. Astron. Space Sci. 2019, 6, 74. [Google Scholar] [CrossRef] [Green Version]
  18. Chen, Y.; Tóth, G. Gauss’s law satisfying energy-conserving semi-implicit particle-in-cell method. J. Comput. Phys. 2019, 386, 632–652. [Google Scholar] [CrossRef] [Green Version]
  19. Campos Pinto, M.; Pagès, V. A semi-implicit electromagnetic FEM-PIC scheme with exact energy and charge conservation. J. Comput. Phys. 2022, 453, 110912. [Google Scholar] [CrossRef]
  20. Boris, J.P. Relativistic plasma simulation-optimization of a hybrid code. In Proceedings of the Fourth Conference on Numerical Simulation of Plasmas, Washington, DC, USA, 2–3 November 1970; pp. 3–67. [Google Scholar]
  21. Brackbill, J.U.; Cohen, B.I. (Eds.) Multiple Time Scales; Academic Press, Inc.: Orlando, FL, USA, 1985. [Google Scholar] [CrossRef]
  22. Vu, H.X.; Brackbill, J.U. Accurate numerical solution of charged particle motion in a magnetic field. J. Comput. Phys. 1995, 116, 384–387. [Google Scholar] [CrossRef]
  23. Lapenta, G.; Brackbill, J.U.; Ricci, P. Kinetic approach to microscopic-macroscopic coupling in space and laboratory plasmas. Phys. Plasmas 2006, 13, 055904. [Google Scholar] [CrossRef] [Green Version]
  24. Markidis, S.; Lapenta, G.; Rizwan-uddin. Multi-scale simulations of plasma with iPIC3D. Math. Comput. Simul. 2010, 80, 1509–1519. [Google Scholar] [CrossRef] [Green Version]
  25. Gonzalez-Herrero, D.; Boella, E.; Lapenta, G. Performance analysis and implementation details of the Energy Conserving Semi-Implicit Method code (ECsim). Comput. Phys. Commun. 2018, 229, 162–169. [Google Scholar] [CrossRef] [Green Version]
  26. de Boor, C. A Practical Guide to Splines; Springer: New York, NY, USA, 1978. [Google Scholar]
  27. Vu, H.X.; Brackbill, J.U. CELEST1D: An implicit, fully kinetic model for low-frequency, electromagnetic plasma simulation. Comput. Phys. Comm. 1992, 69, 253–276. [Google Scholar] [CrossRef]
  28. Burgess, D.; Sulsky, D.; Brackbill, J.U. Mass matrix formulation of the FLIP particle-in-cell method. J. Comput. Phys. 1992, 103, 1–15. [Google Scholar] [CrossRef]
  29. Boella, E.; Innocenti, M.E.; Bettencourt, M.; Chimeh, M.; Lapenta, G.; Parodi, P.; Shukla, N.; Spiga, F. On Adding GPU Support to the Particle-in-Cell code ECsim. OpenACC and Hackathons Summit 2022, 2 August 2022, Virtual. Available online: https://www.youtube.com/watch?v=9FK8zty8BgA (accessed on 15 December 2022).
  30. Lee, W.W. Gyrokinetic particle simulation model. J. Comput. Phys. 1987, 72, 243–269. [Google Scholar] [CrossRef] [Green Version]
  31. Ricci, P.; Lapenta, G.; Brackbill, J.U. A simplified implicit Maxwell solver. J. Computat. Phys. 2002, 183, 117–141. [Google Scholar] [CrossRef] [Green Version]
  32. OpenACC. Available online: https://www.openacc.org/ (accessed on 15 December 2022).
  33. CUDA Zone. Available online: https://developer.nvidia.com/cuda-zone (accessed on 15 December 2022).
  34. Lapenta, G.; Markidis, S.; Marocchino, A.; Kaniadakis, G. Relaxation of relativistic plasmas under the effect of wave-particle interactions. Astrophys. J. 2007, 666, 949–954. [Google Scholar] [CrossRef]
  35. Fried, B.D. Mechanism for instability of transverse plasma waves. Phys. Fluids 1959, 2, 337. [Google Scholar] [CrossRef]
  36. Lazar, M.; Schlickeiser, R.; Wielebinski, R.; Poedts, S. Cosmological effects of Weibel-type instabilities. Astrophys. J. 2009, 693, 1133–1141. [Google Scholar] [CrossRef] [Green Version]
  37. Innocenti, M.E.; Lazar, M.; Markidis, S.; Lapenta, G.; Poedts, S. Electron streams formation and secondary two stream instability onset in the post-saturation regime of the classical Weibel instability. Phys. Plasmas 2011, 18, 052104. [Google Scholar] [CrossRef] [Green Version]
  38. Nevins, W.M.; Hammett, G.W.; Dimits, A.M.; Dorland, W.; Shumaker, D.E. Discrete particle noise in particle-in-cell simulations of plasma microturbulence. Phys. Plasmas 2005, 12, 122305. [Google Scholar] [CrossRef]
Figure 1. Particle velocity distribution at the end of 10 different simulations of same total duration varying the number of sub-cycles, N ν = 1 (top), 2, 3 to 10 (bottom), and correspondingly increasing time step, Δ t , by the same factor. See text for details.
Figure 1. Particle velocity distribution at the end of 10 different simulations of same total duration varying the number of sub-cycles, N ν = 1 (top), 2, 3 to 10 (bottom), and correspondingly increasing time step, Δ t , by the same factor. See text for details.
Physics 05 00007 g001
Figure 2. Scatter plot of the phase space for the particles at the end of 10 different simulations of same total duration varying the number of sub-cycles, N ν = 1 (top), 2, 3 to 10 (bottom), and correspondingly increasing Δ t by the same factor. See text for details.
Figure 2. Scatter plot of the phase space for the particles at the end of 10 different simulations of same total duration varying the number of sub-cycles, N ν = 1 (top), 2, 3 to 10 (bottom), and correspondingly increasing Δ t by the same factor. See text for details.
Physics 05 00007 g002
Figure 3. Evolution of the phase space section ( x , v x ) in the transverse counter-streaming instability. The red particles have v y < 0 , and the blue particles have v y > 0 . Three times are shown, and the left (right) panels report the non-smoothed (smoothed) run: ω p e t = 6.35 (a,b), 25 (c,d), 43.75 (e,f), and 62.5 (g,h). See text for details.
Figure 3. Evolution of the phase space section ( x , v x ) in the transverse counter-streaming instability. The red particles have v y < 0 , and the blue particles have v y > 0 . Three times are shown, and the left (right) panels report the non-smoothed (smoothed) run: ω p e t = 6.35 (a,b), 25 (c,d), 43.75 (e,f), and 62.5 (g,h). See text for details.
Physics 05 00007 g003
Figure 4. Evolution of the energy in the transverse counter-streaming instability without smoothing (left) and with smoothing (right) with the exchange between kinetic and magnetic energy, with a minority contribution going to the electric field energy (top) and the total energy conservation (bottom).
Figure 4. Evolution of the energy in the transverse counter-streaming instability without smoothing (left) and with smoothing (right) with the exchange between kinetic and magnetic energy, with a minority contribution going to the electric field energy (top) and the total energy conservation (bottom).
Physics 05 00007 g004
Figure 5. Transverse counter-streaming instability in absence of smoothing for E x (left,) E z (middle), and B y (right) in the spatiotemporal plane ( x , t ) (top) and in the spectral plane ( k , ω ) (bottom). See text for details.
Figure 5. Transverse counter-streaming instability in absence of smoothing for E x (left,) E z (middle), and B y (right) in the spatiotemporal plane ( x , t ) (top) and in the spectral plane ( k , ω ) (bottom). See text for details.
Physics 05 00007 g005
Figure 6. Transverse counter-streaming instability in presence of smoothing. Three fields are shown from left to right: E x , E z and B y . The top row shows the spatiotemporal plane ( x , t ) , and the bottom row shows the spectral plane ( k , ω ) . E x (left), E z (middle), and B y (right) in the spatiotemporal plane ( x , t ) (top) and in the spectral plane ( k , ω ) (bottom).
Figure 6. Transverse counter-streaming instability in presence of smoothing. Three fields are shown from left to right: E x , E z and B y . The top row shows the spatiotemporal plane ( x , t ) , and the bottom row shows the spectral plane ( k , ω ) . E x (left), E z (middle), and B y (right) in the spatiotemporal plane ( x , t ) (top) and in the spectral plane ( k , ω ) (bottom).
Physics 05 00007 g006
Table 1. Effect of sub-cycling on a two-stream instability run. Values attained while varying the number N ν of sub-cycles, from 1 to 10 and correspondingly increasing the duration of the simulation in seconds, the seconds spent per the time step, Δ t , and the energy error, Δ E , by the same factor Δ t are reported.
Table 1. Effect of sub-cycling on a two-stream instability run. Values attained while varying the number N ν of sub-cycles, from 1 to 10 and correspondingly increasing the duration of the simulation in seconds, the seconds spent per the time step, Δ t , and the energy error, Δ E , by the same factor Δ t are reported.
N ν Time (s)Seconds/ Δ t Δ E
1.042.730.08395 1.75 × 10 16
2.020.940.08244 1.1 × 10 16
3.014.630.08657 1.015 × 10 16
4.010.880.08567 7.388 × 10 17
5.08.970.08881 6.764 × 10 17
6.07.360.08762 8.342 × 10 17
7.06.460.08972 7.988 × 10 17
8.05.740.09111 8.492 × 10 17
9.05.270.09411 8.14 × 10 17
10.04.790.0958 8.61 × 10 17
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Lapenta, G. Advances in the Implementation of the Exactly Energy Conserving Semi-Implicit (ECsim) Particle-in-Cell Method. Physics 2023, 5, 72-89. https://doi.org/10.3390/physics5010007

AMA Style

Lapenta G. Advances in the Implementation of the Exactly Energy Conserving Semi-Implicit (ECsim) Particle-in-Cell Method. Physics. 2023; 5(1):72-89. https://doi.org/10.3390/physics5010007

Chicago/Turabian Style

Lapenta, Giovanni. 2023. "Advances in the Implementation of the Exactly Energy Conserving Semi-Implicit (ECsim) Particle-in-Cell Method" Physics 5, no. 1: 72-89. https://doi.org/10.3390/physics5010007

APA Style

Lapenta, G. (2023). Advances in the Implementation of the Exactly Energy Conserving Semi-Implicit (ECsim) Particle-in-Cell Method. Physics, 5(1), 72-89. https://doi.org/10.3390/physics5010007

Article Metrics

Back to TopTop