GitLab

The Particle-Mesh (PM) solver is one of the oldest improvements over the PP solver. Still one of the best references is the book by R.W. Hockney & J.W. Eastwood [hockney]. The PM solver introduces a discretization of space. The rectangular computation domain \Omega:=[-L_x,L_x]\times[-L_y,L_y]\times[-L_t,L_t], just big enough to include all particles, is segmented into a regular mesh of M=M_x\times M_y\times M_t grid points. For the discussion below we assume N=M_x=M_y=M_t.

The solution of Poisson’s equation is an essential component of any self-consistent electrostatic beam dynamics code that models the transport of intense charged particle beams in accelerators. If the bunch is small compared to the transverse size of the beam pipe, the conducting walls are usually neglected. In such cases the Hockney method may be employed [hockney, eastwoodandbrownrigg,hockneyandeastwood]. In that method, rather than computing N_p^2 point-to-point interactions (where N_p is the number of macro-particles), the potential is instead calculated on a grid of size (2 N)^d, where N is the number of grid points in each dimension of the physical mesh containing the charge, and where d is the dimension of the problem. Using the Hockney method, the calculation is performed using Fast Fourier Transform (FFT) techniques, with the computational effort scaling as (2N)^d (log_2 2N)^d.

When the beam bunch fills a substantial portion of the beam pipe transversely, or when the bunch length is long compared with the pipe transverse size, the conducting boundaries cannot be ignored. Poisson solvers have been developed previously to treat a bunch of charge in an open-ended pipe with various geometries [qiangandryne,qiangandgluckstern].

The solution of the Poisson equation,

\nabla^2\phi=-\rho/\epsilon_0,

for the scalar potential, \phi, due to a charge density, \rho, and appropriate boundary conditions, can be expressed as,

\phi(x,y,z)=\int\int\int{\mathrm{d}x' \,\mathrm{d}y' \,\mathrm{d}z'}\rho(x',y',z') G(x,x',y,y',z,z'),

where G(x,x',y,y',z,z') is the Green function, subject to the appropriate boundary conditions, describing the contribution of a source charge at location (x',y',z') to the potential at an observation location (x,y,z).

For an isolated distribution of charge this reduces to

\phi(x,y,z)=\int\int\int{\mathrm{d}x' \,\mathrm{d}y' \,\mathrm{d}z'}\rho(x',y',z') G(x-x',y-y',z-z'),

where

G(u,v,w)={\frac{1}{\sqrt{u^2+v^2+w^2}}}.

A simple discretization of Equation [convolutionsolution] on a Cartesian grid with cell size (h_x,h_y,h_z) leads to,

\phi_{i,j,k}=h_x h_y h_z \sum_{i'=1}^{M_x}\sum_{j'=1}^{M_y}\sum_{k'=1}^{M_t}  \rho_{i',j',k'}G_{i-i',j-j',k-k'},

where \rho_{i,j,k} and G_{i-i',j-j',k-k'} denote the values of the charge density and the Green function, respectively, defined on the grid M.

FFT-based Convolutions and Zero Padding

\bar{\phi}_j=\sum_{k=0}^{K-1}\bar{\rho}_k \bar{G}_{j-k}\quad,
\begin{array}{l}
j=0,\ldots,J-1 \\
k=0,\ldots,K-1 \\
j-k=-(K-1),\ldots,J-1 \\
\end{array}

\rho_k=\left\{
\begin{array}{l l}
\bar{\rho}_k & \quad \text{if }k=0,\ldots,K-1 \\
0 & \quad \text{if }k=K,\ldots,N-1. \\
\end{array}\right.

G_m=\left\{
\begin{array}{l l}
\bar{G}_m & \quad \text{if }m=-(K-1),\ldots,J-1 \\
0 & \quad \text{if }m=J,\ldots,N-K, \\
G_{m+iN}=G_{m} & \quad \text{for } i \text{ integer }.
\end{array}\right.

{\phi}_j=\frac{1}{N}\sum_{k=0}^{N-1} W^{-jk}
                    \left(\sum_{n=0}^{N-1} \rho_n W^{nk}\right)
                    \left(\sum_{m=0}^{N-1} G_m W^{mk}\right),
~~~~~~0 \le j \le N-1,

{\phi}_j=
          \sum_{n=0}^{K-1}~
          \sum_{m=0}^{N-1} \bar{\rho}_n G_m
\frac{1}{N}\sum_{k=0}^{N-1} W^{(m+n-j)k}
~~~~~~0 \le j \le N-1.

\sum_{k=0}^{N-1} W^{(m+n-j)k}= N \delta_{m+n-j,iN} ~ ~ ~ ~ ~(i~\mathrm{an~integer}).

{\phi}_j=\sum_{n=0}^{K-1}~\bar{\rho}_n G_{j-n+iN}
~~~~~~0 \le j \le N-1.

{\phi}_j=\sum_{n=0}^{K-1}~\bar{\rho}_n G_{j-n}
~~~~~~0 \le j \le N-1.

\bar{\phi}_j=\sum_{k=0}^{K-1}\bar{\rho}_k \bar{G}_{j+k}\quad,
\begin{array}{l}
j=0,\ldots,J-1 \\
k=0,\ldots,K-1 \\
j-k=-(K-1),\ldots,J-1 \\
\end{array}

Algorithm used in OPAL

\phi_{i,j,k}=h_x h_y h_z \text{FFT}^{-1} \{ ( \text{FFT}\{\rho_{i,j,k}\}) ( \text{FFT}\{G_{i,j,k}\}) \}

Interpolation Schemes

This is a scalable parallel solver for the Poisson equation within a Particle-In-Cell (PIC) code for the simulation of electron beams in particle accelerators of irregular shape. The problem is discretized by Finite Differences. Depending on the treatment of the Dirichlet boundary the resulting system of equations is symmetric or `mildly' non-symmetric positive definite. In all cases, the system is solved by the preconditioned conjugate gradient algorithm with smoothed aggregation (SA) based algebraic multigrid (AMG) preconditioning. More details are given in [Adelmann:2009p543].

The beam out of a cathode or in a plasma wake field accelerator can have a large energy spread. In this case, the static approximation using one Lorentz frame might not be sufficient. Multiple Lorentz frames can be used so that within each Lorentz frame the energy spread is small and hence the electrostatic approximation is valid. More details will be given in Version 2.0.0.

See Table [fieldsolvercmd] for a summary of the Fieldsolver command.

At present only a FFT based solver is available. Future solvers will include Finite Element solvers and a Multigrid solver with Shortley-Weller boundary conditions for irregular domains.

The dimensions in x, y and z can be parallel (TRUE) or serial FALSE. The default settings are: parallel in z and serial in x and y.

Number of grid points in x, y and z for a rectangular grid.

Two boundary conditions can be selected independently among x, y namely: OPEN and for z OPEN & PERIODIC. In the case you select for z periodic you are about to model a DC-beam.

Two Greens functions can be selected: INTEGRATED, STANDARD. The integrated Green’s function is described in [qiang2005, qiang2006-1, qiang2006-2]. Default setting is INTEGRATED.

The bounding box defines a minimal rectangular domain including all particles. With BBOXINCR the bounding box can be enlarged by a factor given in percent of the minimal rectangular domain.

The list of geometries defining the beam line boundary. For further details see Chapter [geometry].

The iterative solver for solving the preconditioned system: CG, BiCGSTAB or GMRES.

The interpolation method for grid points near the boundary: CONSTANT, LINEAR or QUADRATIC.

The tolerance for the iterative solver used by the MG solver.

The maximal number of iterations the iterative solver performs.

The behavior of the preconditioner can be: STD, HIERARCHY or REUSE. This argument is only relevant when using the MG solver and should only be set if the consequences to simulation and solver are evident. A short description is given in the Table below.

Suppose \mathrm{d} E the energy spread in the particle bunch is to large, the electrostatic approximation is no longer valid. One solution to that problem is to introduce k energy bins and perform k separate field solves in which \mathrm{d} E is again small and hence the electrostatic approximation valid. In case of a cyclotron see Section [cyclotron] the number of energy bins must be at minimum the number of neighboring bunches (NNEIGHBB) i.e. \mathrm{ENBINS} \le \mathrm{NNEIGHBB}.

The variable MINSTEPFORREBIN defines the number of integration step that have to pass until all energy bins are merged into one.

The option AMR=TRUE enables further commands to be used in the Fieldsolver command. In order to run a proper AMR simulation one requires following ingredients:

• Upper bound of refinement, specified by AMR_MAXLEVEL

• Maximum allowable grid size, specified by AMR_MAXGRID. This size holds for all levels in all directions.

• Mesh refinement strategy, specified by AMR_TAGGING. Depending on the tagging scheme, further keywords can be used. A summary is given in the Table below.