Abstract
Based on the previously developed Energy Conserving Semi Implicit Method (ECsim) code, we present its cylindrical implementation, called ECsim-CYL, to be used for axially symmetric problems. The main motivation for the development of the cylindrical version is to greatly improve the computational speed by utilizing cylindrical symmetry. The ECsim-CYL discretizes the field equations in two-dimensional cylindrical coordinates using the finite volume method . For the particle mover, it uses a modification of ECsim’s mover for cylindrical coordinates by keeping track of all three components of velocity vectors, while only keeping radial and axial coordinates of particle positions. In this paper, we describe the details of the algorithm used in the ECsim-CYL and present a series of tests to validate the accuracy of the code including a wave spectrum in a homogeneous plasmas inside a cylindrical waveguide and free expansion of a spherical plasma ball in vacuum. The ECsim-CYL retains the stability properties of ECsim and conserves the energy within machine precision, while accurately describing the plasma behavior in the test cases.
1. Introduction
Particle In Cell (PIC) algorithms are widely used to study natural and laboratory plasmas using fullykinetic description [1, 2]. In a typical PIC algorithm, particles and fields are advanced alternatively in a cycle.Particles are moved with the fields obtained by interpolation from a grid and the particle information is thenused to project the sources of the fields, i.e. current and charge, onto the grid advancing the field values intime. This alternation of particle movement and field solution is the typical feature of explicit PIC codes.Explicit PIC codes are straight-forward to design and lead to extremely efficient parallel implementations.As such, explicit PIC codes are always among the most successful applications on new supercomputers.
However, explicit PIC codes suffer from a scale limitation [1, 2]: the alternative field and particle solverapproach of the explicit cycle leads to two limitations: the time step needs to resolve the electron plasmafrequency and the grid spacing need to resolve a multiple of the Debye length (the multiple depends on theinterpolations scheme). Plasma problems are typically multi scale, with ions responding on much longer andslower scales than electrons. The need for explicit PIC codes to resolve all electron scales, led to developmentof implicit PIC methods to overcome this difficulty.
Implicit methods retain the coupling of particle motion and field solution at each time step. They areno longer done separately in sequence but rather become part of an iteration. This coupling is non-linearand it requires the use of non-linear solvers [3, 4].
However, when the coupling of particles and fields is linearized, the iteration becomes linear and typicallybased on Krylov solvers: in this case the methods are considered semi-implicit. This semi-implicitness should not be confused with the methods that solve implicitly only the field equations but still remain explicit inalternating particle mover and field solver [5]. This class of algorithms are designed to overcome the speedof light limitation in the electromagnetic waves and are not considered here.
All implicit and semi-implicit PIC codes enable the user to use larger time steps and grid spacing, anadvantage when dealing with multi scale problems where the electron physics need not be represented indetail but still needs to be captured in macroscopic terms [6]. For these algorithms, the electron responseto large scales is captured correctly but the small scale fluctuations (e.g. Langmuir waves) are damped[7, 8]. While fully implicit methods require an advanced arsenal of non-linear solvers, supported by efficientpreconditioning [9] to reduce the number of iterations, semi-implicit methods can use more standard linearsolvers.
Semi-implicit methods were the first to be deployed in large practical applications. The literature reportstwo lines of semi-implicit codes: the direct implicit method [10] and the moment implicit methods [7].The most modern implementations are, respectively, LSP [11] and iPic3D [12]. Among the semi-implicitalgorithms, we have recently developed a new algorithm called ECsim [13, 14, 15] that adds one moredesirable property to the previous semi-implicit schemes: energy conservation.
ECsim uses a new tool to represent the linear response of the plasma to the fields: the mass matrix[16, 13]. The use of a new mover allows the introduction of this linear formulation without requiring anylinearization: the result is that ECsim is exactly energy conserving. This has two aspects. First, thereis a formal analytical energy integral that has been derived to show energy conservation by the scheme.Second, when coded into a computer implementation, energy is conserved to machine precision, a powerfuldebugging feature as well as a desirable physical property [13].
ECsim has the extra computational complexity of requiring the computation of the mass matrix byinterpolation from the particles. The extra cost of this operation was evaluated in Ref. [15]. Comparedwith other PIC codes, ECsim requires more memory to store the mass matrix and more computational coststo compute it. When compared with specifically with the moment implicit methods, such as iPic3D [12],ECsim increases the number of operations needed to project particle information to the grid, as required tocompute the mass matrix. However, this increase is in part compensated by the reduced complexity of themover. The cost per cycle increases by a factor of 2 to 4 but it is more than compensated by the ability touse larger time steps allowed by energy conservation [15]. The increase in memory is not significant as it doesnot increase the particle memory, by far the overwhelming majority of the memory used: by comparison,the additional memory required to store the mass matrix is not significant.
In many applications, there is an easy way to reduce the cost: use the intrinsic symmetry of the problem.Many experiments and natural processes have a natural axis of symmetry and cylindrical symmetry allowsone to treat them as two-dimensional problems. Symmetry reduces the computational effort by reducing thedimensionality, enabling higher spatial resolutions or reducing the computational cost. We have implementeda new version of ECsim in cylindrical coordinates (hereinafter ECsim-CYL) for these axially symmetric cases.This new version, analogously to the original ECsim code, conserves the total energy exactly and retainsthe same stability properties.
Cylindrical symmetry has been used in several previous explicit and even semi-implicit codes, based onsimpler algorithms that do not conserve energy exactly. In the explicit approach, cylindrical symmetry isdescribed in textbooks [1] and has been applied in production codes [17]. Also for the moment implicitapproach a cylindrical version was developed [18], based on the same formulation of the implicit momentmethod [8] in cartesian geometry, at the basis of Venus, Celeste and iPic3D: this method does not conserveenergy exactly.
However, ECsim-CYL is the first semi-implicit code which implements the Energy Conserving SemiImplicit algorithm in cylindrical coordinates and is the first semi-implicit code that conserves energy.
More recently, general curvilinear geometry codes have been developed both for explicit [19, 20] and forimplicit codes [21].
ECsim-CYL has been developed in C++ and uses MPI for parallel communications. It shares theinfrastructure (particle communication, particle-to-grid interpolations, etc.) of the fully 3D Cartesian ECsim[15] and has a similar scaling behavior and performance.
The paper is organized as follows. In section 2, we describe the algorithm used in the ECsim-CYL implementation regarding the discretization of the field equations using cylindrical coordinates and theimplementation of the particle mover. In section 3, we test the accuracy of both the field solver and particlemover. In section 4, we test the ECsim-CYL results against the analytical formula regarding the plasma wavein an homogeneous plasma inside a cylindrical waveguide. This is followed by the test results of ECsim-CYLin the simulation of homogeneous expansion of a plasma sphere in vacuum with different spatio-temporalresolution and ion mass ratio in section 5. In Section 6, we summarize the results and present the conclusions.
