
1. Field Solver
 1.1. FFT Based ParticleMesh (PM) Solver
 1.2. Iterative Space Charge Solver
 1.3. Energy Binning
 1.4. The
FIELDSOLVER
Command  1.5. Define the Fieldsolver to be used
 1.6. Define Domain Decomposition
 1.7. Define Number of Grid Points
 1.8. Define Boundary Conditions
 1.9. Define Greens Function
 1.10. Define Bounding Box Enlargement
 1.11. Define Geometry
 1.12. Define Iterative Solver
 1.13. Define Interpolation for Boundary Points
 1.14. Define Tolerance
 1.15. Define Maximal Iterations
 1.16. Define Preconditioner Behavior
 1.17. Define the number of Energy Bins to use
 1.18. Define AMR Solver
1. Field Solver
Space charge effects are included in the simulation by specifying a field solver described in this chapter and attaching it to the track command as described in Chapter [track]. By default, the code does not assume any symmetry i.e. full 3D. In the near future it is planed to implement also a slice (2D) model. This will allow the use of less numbers of macroparticles in the simulation which reduces the computational time significantly.
The space charge forces are calculated by solving the 3D Poisson equation with open boundary conditions using a standard or integrated Green function method. The image charge effects of the conducting cathode are also included using a shifted Green function method. If more than one Lorentz frame is defined, the total space charge forces are then the summation of contributions from all Lorentz frames. More details will be given in Version 2.0.0.
1.1. FFT Based ParticleMesh (PM) Solver
The ParticleMesh (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
selfconsistent 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
pointtopoint
interactions (where N_p
is the number of macroparticles),
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 openended 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(xx',yy',zz'),
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_{ii',jj',kk'},
where \rho_{i,j,k}
and G_{ii',jj',kk'}
denote the values of the charge density and the Green function,
respectively, defined on the grid M
.
FFTbased Convolutions and Zero Padding
FFTs can be used to compute convolutions by appropriate zeropadding of the sequences. Discrete convolutions arise in solving the Poisson equation, and one is typically interested in the following,
\bar{\phi}_j=\sum_{k=0}^{K1}\bar{\rho}_k \bar{G}_{jk}\quad,
\begin{array}{l}
j=0,\ldots,J1 \\
k=0,\ldots,K1 \\
jk=(K1),\ldots,J1 \\
\end{array}
where \bar{G}
corresponds to the free space Green
function, \bar{\rho}
corresponds to the charge density,
and \bar{\phi}
corresponds to the scalar potential. The
sequence \{\bar{\phi}_j\}
has J
elements,
\{\bar{\rho}_k\}
has K
elements, and
\{\bar{G}_m\}
has M=J+K1
elements.
One can zeropad the sequences to a length N\ge M
and use
FFT’s to efficiently obtain the \{\bar{\phi}_j\}
in the
unpadded region. This defines a zeropadded charge density,
\rho
,
\rho_k=\left\{
\begin{array}{l l}
\bar{\rho}_k & \quad \text{if }k=0,\ldots,K1 \\
0 & \quad \text{if }k=K,\ldots,N1. \\
\end{array}\right.
Define a periodic Green function,
G_m
, as follows,
G_m=\left\{
\begin{array}{l l}
\bar{G}_m & \quad \text{if }m=(K1),\ldots,J1 \\
0 & \quad \text{if }m=J,\ldots,NK, \\
G_{m+iN}=G_{m} & \quad \text{for } i \text{ integer }.
\end{array}\right.
Now consider the sum
{\phi}_j=\frac{1}{N}\sum_{k=0}^{N1} W^{jk}
\left(\sum_{n=0}^{N1} \rho_n W^{nk}\right)
\left(\sum_{m=0}^{N1} G_m W^{mk}\right),
~~~~~~0 \le j \le N1,
where W=e^{2\pi i/N}
. This
is just the FFTbased convolution of \{\rho_k\}
with
\{G_m\}
. Then,
{\phi}_j=
\sum_{n=0}^{K1}~
\sum_{m=0}^{N1} \bar{\rho}_n G_m
\frac{1}{N}\sum_{k=0}^{N1} W^{(m+nj)k}
~~~~~~0 \le j \le N1.
Now use the relation
\sum_{k=0}^{N1} W^{(m+nj)k}= N \delta_{m+nj,iN} ~ ~ ~ ~ ~(i~\mathrm{an~integer}).
It follows that
{\phi}_j=\sum_{n=0}^{K1}~\bar{\rho}_n G_{jn+iN}
~~~~~~0 \le j \le N1.
But G
is periodic with period
N
. Hence,
{\phi}_j=\sum_{n=0}^{K1}~\bar{\rho}_n G_{jn}
~~~~~~0 \le j \le N1.
In the physical (unpadded) region,
j\in \left[0,J1\right]
, so the quantity jn
in Equation [finaleqn] satisfies (K1)\le jn \le J1
. In
other words the values of G_{jn}
are identical to
\bar{G}_{jn}
. Hence, in the physical region the FFTbased
convolution, Equation [fftconvolution], matches the convolution in
Equation [bruteforceconvolution].
As stated above, the zeropadded sequences need to have a length
N \ge M
, where M
is the number of elements
in the Green function sequence \left\{x_m\right\}
. In
particular, one can choose N=M
, in which case the Green
function sequence is not padded at all, and only the charge density
sequence, \left\{r_k\right\}
, is zeropadded, with
k=0,\ldots,K1
corresponding to the physical region and
k=K,\ldots,M1
corresponding to the zeropadded region.
The above FFTbased approach – zeropadding the charge density array,
and circularshifting the Green function in accordance with
Equation [periodicgreenfunction] – will work in general. In addition, if
the Green function is a symmetric function of its arguments, the value
at the end of the Green function array (at grid point J1
)
can be dropped, since it will be recovered implicitly through the
symmetry of Equation [periodicgreenfunction]. In that case the approach
is identical to the Hockney method [hockney,
eastwoodandbrownrigg,hockneyandeastwood]. Lastly, note that the above
proof that the convolution, Equation [fftconvolution], is identical to
Equation [bruteforceconvolution] in the unpadded region, works even when
W^{j k}
and W^{m k}
are replaced by
W^{j k}
and W^{m k}
, respectively, in
Equation [fftconvolution]. In other words, the FFTbased approach can be
used to compute
\bar{\phi}_j=\sum_{k=0}^{K1}\bar{\rho}_k \bar{G}_{j+k}\quad,
\begin{array}{l}
j=0,\ldots,J1 \\
k=0,\ldots,K1 \\
jk=(K1),\ldots,J1 \\
\end{array}
simply by changing the direction of the Fourier transform of the Green function and changing the direction of the final Fourier transform.
Algorithm used in OPAL
As a result, the solution of Equation [openbruteforceconvolution] is then given by
\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}\}) \}
where the notation has been introduced that
\text{FFT}\{ . \}
denotes a forward FFT in all 3
dimensions, and \text{FFT}^{1}\{ . \}
denotes a backward
FFT in all 3 dimensions.
1.2. Iterative Space Charge Solver
This is a scalable parallel solver for the Poisson equation within a ParticleInCell (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' nonsymmetric 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].
1.3. Energy Binning
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.
1.4. The FIELDSOLVER
Command
See Table [fieldsolvercmd] for a summary of the Fieldsolver command.
Command  Purpose 


Specify a fieldsolver 

Specify the type of field solver: 

If 

If 

If 

Number of grid points in 

Number of grid points in 

Number of grid points in 

Boundary condition in 

Boundary condition in 

Boundary condition in 

Defines the Greens function for the FFT Solver 

Enlargement of the bounding box in % 

Geometry to be used as domain boundary 

Type of iterative solver 

Interpolation used for boundary points 

Tolerance for iterative solver 

Maximum number of iterations of iterative solver 

Behavior of the preconditioner 

Maximum adaptive mesh refinement level (single level if


Maximum grid size. It has to be smaller than 

Meshrefinement strategy 
1.5. Define the Fieldsolver to be used
At present only a FFT based solver is available. Future solvers will include Finite Element solvers and a Multigrid solver with ShortleyWeller boundary conditions for irregular domains.
1.6. Define Domain Decomposition
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
.
1.8. Define Boundary Conditions
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 DCbeam.
1.9. Define Greens Function
Two Greens functions can be selected: INTEGRATED
, STANDARD
. The
integrated Green’s function is described in [qiang2005, qiang20061,
qiang20062]. Default setting is INTEGRATED
.
1.10. Define Bounding Box Enlargement
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.
1.11. Define Geometry
The list of geometries defining the beam line boundary. For further details see Chapter [geometry].
1.12. Define Iterative Solver
The iterative solver for solving the preconditioned system: CG
,
BiCGSTAB
or GMRES
.
1.13. Define Interpolation for Boundary Points
The interpolation method for grid points near the boundary: CONSTANT
,
LINEAR
or QUADRATIC
.
1.16. Define Preconditioner Behavior
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.
Value  Behavior 


The preconditioner is rebuilt in every time step (enabled by default) 

The hierarchy (tentative prolongator) is reused 

The preconditioner is reused 
1.17. Define the number of Energy Bins to use
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.
1.18. Define AMR Solver
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.
Value  Behavior 


Mark each cell if


Mark each cell if the electric field component of any
direction satisfies


It performs a loop over all particles of a level and
computes the dot product of the momenta. Every cell that contains a
particle with


If the charge density of a cell is greater or equal
to the value specified with 

Cells with equal or more particles are
refined. The bound is specified with 

Cells with equal or less particles are
refined. The bound is specified with 