\input{header} \chapter{Field Solver} \label{chp:fieldsolver} \index{Field Solver|(} \ifthenelse{\boolean{ShowDebug}}{ \TODO{AA will rewrite} }{} 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 \chpref{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 macro-particles 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. \section{FFT Based Particle-Mesh (PM) Solver} \label{sec:fieldsolver:fftbased} \index{Field Solver!FFT based} 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 \ref{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 \ref{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 \ref{qiangandryne,qiangandgluckstern}. The solution of the Poisson equation, $$%\label{ } \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'), \label{eq:convolutionsolution}$$ where $$G(u,v,w)={\frac{1}{\sqrt{u^2+v^2+w^2}}}. \label{eq:isolatedgreenfunction}$$ A simple discretization of Equation~\ref{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'}, \label{eq:openbruteforceconvolution}$$ 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$. \subsubsection{FFT-based Convolutions and Zero Padding} FFTs can be used to compute convolutions by appropriate zero-padding 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}^{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} \label{eq:bruteforceconvolution}$$ 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+K-1$ elements. One can zero-pad 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 zero-padded charge density, $\rho$, $$\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.$$ Define a periodic Green function, $G_m$, as follows, $$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. \label{eq:periodicgreenfunction}G$$ Now consider the sum $${\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, \label{eq:fftconvolution}$$ where $W=e^{-2\pi i/N}$. This is just the FFT-based convolution of $\{\rho_k\}$ with $\{G_m\}$. Then, $${\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.$$ Now use the relation $$\sum_{k=0}^{N-1} W^{(m+n-j)k}= N \delta_{m+n-j,iN}~~~~~(i~\rm an~integer).$$ It follows that $${\phi}_j=\sum_{n=0}^{K-1}~\bar{\rho}_n G_{j-n+iN} ~~~~~~0 \le j \le N-1.$$ But $G$ is periodic with period $N$. Hence, $${\phi}_j=\sum_{n=0}^{K-1}~\bar{\rho}_n G_{j-n} ~~~~~~0 \le j \le N-1. \label{eq:finaleqn}$$ In the physical (unpadded) region, $j\in \left[0,J-1\right]$, so the quantity $j-n$ in Equation~\ref{finaleqn} satisfies $-(K-1)\le j-n \le J-1$. In other words the values of $G_{j-n}$ are identical to $\bar{G}_{j-n}$. Hence, in the physical region the FFT-based convolution, Equation~\ref{fftconvolution}, matches the convolution in Equation~\ref{bruteforceconvolution}. As stated above, the zero-padded 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 zero-padded, with $k=0,\ldots,K-1$ corresponding to the physical region and $k=K,\ldots,M-1$ corresponding to the zero-padded region. %As just shown, FFT's can be used to compute the convolution shown in Equation~\ref{bruteforceconvolution}, with no zero-padding of the Green function, %as long as there are $J$ values of the potential, $K$ values of the charge density, and $J+K-1$ values of the Green function. %%%This is normally the case in particle-in-cell beam dynamics codes, in fact it is usually the case that J=K (the region containing the charge is also the %%%region within which the potential is computed), and the Green function is defined on... %These conditions are usually satisfied in particle-in-cell beam dynamics code, which have a a grid containing the charge density, %and which have a grid containing the Green function ({\it e.g.}, Equation~\ref{isolatedgreenfunction}, Equation~\ref{rgreenfunction}, or Equation~\ref{rintgreenfunction}) %for both positive and negative arguments. The above FFT-based approach -- zero-padding the charge density array, and circular-shifting the Green function in accordance with Equation~\ref{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 $J-1$ in this discussion) (at grid point $J-1$) can be dropped, since it will be recovered implicitly through the symmetry of Equation~\ref{periodicgreenfunction}. In that case the approach is identical to the Hockney method \ref{hockney, eastwoodandbrownrigg,hockneyandeastwood}. %Naively one might expect that, for a Green function defined at grid points $(0,\pm 1, \pm 2,\ldots)$, there would be an odd number of grid points. %But since the last point can be dropped, Green function arrays usually have an even number of grid points. %When the free space Green function is a symmetric function of $\mathbf{x}-\mathbf{x'}$, a circular shift of the Green function array, moving $\mathbf{x}-\mathbf{x'}=0$ from the %center of the grid to the corner corresponding to the origin of grid points, produces a periodic Green function that is identical to the %periodic Green function used in the Hockney method. Viewed this way, the only difference between an "ordinary" FFT-based convolution %and one described by Hockney is a circular shift of the Green function. The usual description of the Hockney approach involves making the Green function periodic and symmetric, %but that is because the Green function for the free space potential is symmetric; if the Hockney approach were used to directly compute the electric fields %by convolving the charge density with the electric field Green functions, then the Green functions would have to be anti-symmetrized. %But when viewed simply as a convolution in which the charge density is zero padded, no explicit period-ization is required (because it it %implicit in the FFT-based approach), so it will work for any Green function, regardless of symmetry. %In the area of signal processing, the values of charge density, $\left\{r_k\right\}$, correspond to the response function, %and the values of the Green function, $\left\{x_{j-k}\right\}$ or $\left\{x_m\right\}$, correspond to the signal. %Based on the above discussion it follows that, if two arbitrary sequences $\left\{r_0,\ldots,r_{K-1}\right\}$ and $\left\{x_0,\ldots,x_{M-1}\right\}$ are to be convolved %to produce a sequence $\left\{p_0,\ldots,p_{J-1}\right\}$, and if the number of points in $\left\{x_m\right\}$ is less than $J+K-1$, then, %in order to use FFT's, the sequence $\left\{p_j\right\}$ needs to be zero-padded to extend it to have $J+K-1$ points. Lastly, note that the above proof that the convolution, Equation~\ref{fftconvolution}, is identical to Equation~\ref{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~\ref{fftconvolution}. In other words, the FFT-based approach can be used to compute $$\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} \label{eq:bruteforcecorrelation}$$ simply by changing the direction of the Fourier transform of the Green function and changing the direction of the final Fourier transform. \subsubsection{Algorithm used in \textit{OPAL}} As a result, the solution of Equation~\ref{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}\}) \} \label{eq:oneterm}$$ 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. %In order to obtain ${\cal M}_2$ in Equation~\ref{splitOper1} we must solve Poisson's equation, %where $\rho$ stands for the %charge density and $\phi$ for the scalar electrostatic potential: %$$\label{eq:Poisson} %\nabla^2 \phi(\mathbf{q}) = - \frac{\rho(\mathbf{q})}{\epsilon_0} %$$ %subject to open boundary conditions in all spatial directions: $\phi (\mathbf{q}) \rightarrow 0$ as $|\mathbf{q}| \rightarrow \infty$ or %imposing periodic boundary conditions in longitudinal directions. The assumption of using an "isolated system" is physically %motivated by observing the ratio of the beam size to vacuum vessel domensions. It has the computational advantages that one can use cyclic convolution in Equation~\ref{FourierPoisson}. %The computational domain $\Omega \subset \R^3$ is simply connected and has a cylindrical %or rectilinear shape. %The corresponding integral equation reads: %$$%\phi(\mathbf{q}) = \dis\int\limits_\Omega \, G(\mathbf{q} - \mathbf{q}\primed) \,\rho (\mathbf{q'}) \;\mathrm{d} \mathbf{q}\primed , \;\;\Omega \subset \R^3 %$$ %where $G$ is the Green's function which gives the response to a unit source term. In 3D we have %$$%G(\mathbf{q}-\mathbf{q}\primed) = \dis\frac{1}{4 \pi \,|\mathbf{q} - \mathbf{q}\primed|} \,. %$$ %The electric field then follows from the electrostatic potential %$$\label{eq:Electric} %\mathbf{E} = - \nabla \phi. %$$ %The charges are assigned from the particle positions in continuum, onto the grid using %one of two available interpolation schemes: cloud in cell (CIC) or nearest grid point (NGP). %Then the Poisson equation is solved on the mesh and the electric field at the particle positions is obtained by %interpolating back from the mesh. The use of the convolution theorem to solve the discretized Poisson %equation \eqref{Poisson} on the grid can dramatically improve performance. %Let $\Omega^D$ be spanned by a mesh of $l \times n \times m$ with $l= 1 \dots M_x$, $n= 1\dots M_x$ and $m= 1 \dots M_t$. %The solution of the discretized Poisson equation with $\mathbf{k}=(l,n,m,)$ %$$\label{eq:DiscretizedPoisson} %\nabla^{2} \phi^D(\mathbf{k}) = - \frac{\rho^D(\mathbf{k})}{\epsilon_0}, \mathbf{k} \in \Omega^D. %$$ %$\phi^D$ then is given by convolution with the appropriate discretized Green's function $G_D$: %$$%\phi^D = \rho^D * G^D. %$$ %In Fourier space (hats) the convolution becomes a simple multiplication, with %$%G(\mathbf{q})=\frac{1}{4\pi}\frac{1}{|\mathbf{q}|}\longrightarrow \widehat{G}(\mathbf{k})=-\frac{1}{4\pi}\frac{1}{|\mathbf{k}|^2} %$ %and we get: %$$\label{eq:FourierPoisson} %\widehat{\phi}^D = \widehat{\rho}^D \cdot \widehat{G}^D. %$$ %Thus, the convolution sum is converted to a single multiplication at the cost of a Fourier transform. Fortunately, Fast Fourier Transform (FFT) on the mesh is a very fast and accurate method %of transforming mesh-defined quantities to Fourier space. It needs $\mathcal{O}$($M\log M$) computational effort, so that, together with the interpolations, an overall scaling of $\mathcal{O}$($N+M\log M$) is achieved. %In order to have good spatial resolution, small grid sizes are often necessary, which again require a %large number of particles. Therefore, both the grid size $M$ and the particle number $N$ %are limiting factors. %The PM Solver Algorithm is summarized in the following algorithm: %\begin{tabbing} %{\bf PM Solver Algorithm}\\ %\quad $\triangleright$ Assign particle charges $q_i$ to nearby mesh points to obtain $\rho^D$ \\ %\quad $\triangleright$ Use FFT on $\rho^D$ and $G^D$ to obtain $\widehat{\rho}^D$ and $\widehat{G}^D$ \\ %\quad $\triangleright$ Determine $\widehat{\phi}^D$ on the grid using $\widehat{\phi}^D = \widehat{\rho}^D \cdot \widehat{G}^D$ \\ %\quad $\triangleright$ Use inverse FFT on $\widehat{\phi }^D$ to obtain $\phi^D$ \\ %\quad $\triangleright$ Compute $\mathbf{E}^D= -\nabla \phi^D$\footnote{using a second order finite difference method} \\ %\quad $\triangleright$ Interpolate $\mathbf{E(\mathbf{q})}$ at particle positions $\mathbf{q}$ from $\mathbf{E}^D$ \\ %\end{tabbing} %\subsubsection{Open and Periodic Boundary Conditions} %In order to meet open boundary conditions and to remove the intrinsic periodicity of the FFT, the grid size needs to be doubled in all spatial dimensions and the charge distribution is located at only one octant. The charge distribution is set equal to zero elsewhere. If the potential is then calculated in the entire enlarged region, the correct potential for an isolated system is obtained in the physical' octant. This is referred to as the Hockney Trick' \ref{hockney}. For periodic boundary conditions in the longitudinal direction, the Hockney trick is applied to the transverse directions only. %The main drawback of this method is its high storage requirement. However, using symmetries one is %able to bound the required storage to $2N_g$ where $N_g$ is the grid size used in the physical %region of the calculation (see \ref{hockney} on p. 213). \subsubsection{Interpolation Schemes} \index{Field Solver!Interpolation Schemes} More details will be given in Version 2.0.0. %\label{sec:interpol} %Both charge assignment and electric field interpolation are related to the interpolation %scheme used. A detailed discussion is given in~\ref{hockney}. %If $e_i$ is the charge of a particle, we can write the density at mesh point $\mathbf{k}_m$ as %$$\label{eq:discRho} %\rho(\mathbf{k}_m)^D = \sum_{i=1}^N e_i\cdot W(\mathbf{q}_i,\mathbf{k}_m), ~ m=1\dots M %$$ %where $W$ is a suitably chosen weighting function (with local support). %The simplest scheme is the nearest grid point (NGP) method, where the total particle charge is assigned to %the nearest grid point and the electric field is also evaluated at the nearest grid point. A more %elaborate scheme is called cloud in cell (CIC). It assigns the charge to the $2^d$ nearest grid points and %also interpolates the electric field from these grid points. The assigned density changes are continuous when %a particle moves across a cell boundary, although the first derivative is discontinuous. A schematic of the %CIC interpolation scheme is shown in Figure~\ref{CIC} for the two-dimensional case. %%total momentum consered. errors small at large particle separations. charge assignement schould vary smoothly as particle position changes. % \begin{figure} % \begin{center} %% \includegraphics[width=0.75\linewidth]{\figdir/cic.eps} % \caption[CIC interpolation scheme]{\label{fig:CIC} \it CIC interpolation scheme.} % \end{center} % \end{figure} \section{Iterative Space Charge Solver} \index{Field Solver!Iterative} 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 \ref{Adelmann:2009p543}. \section{Energy Binning} \index{Field Solver!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. \section{The \texttt{FIELDSOLVER} Command} \label{sec:fieldsolvercmd} \index{Field Solver!Command} \index{FIELDSOLVER} See Table~\ref{fieldsolvercmd} for a summary of the Fieldsolver command. \begin{table}[ht] \footnotesize \begin{center} \caption{Fieldsolver command summary} \label{tab:fieldsolvercmd} \begin{tabular}{|l|p{0.6\textwidth}|l|} \hline \tabhead Command & Purpose \\ \hline \texttt{FIELDSOLVER} & Specify a fieldsolver \\ \texttt{FSTYPE} & Specify the type of field solver: \texttt{FFT}, \texttt{FFTPERIODIC}, \texttt{MG}, \texttt{FMG} (AMR only) and \texttt{NONE} \index{FFT}\index{FFTPERIODIC}\index{MG}\index{AMR}\index{NONE} \\ \texttt{PARFFTX} & If \texttt{TRUE}, the dimension $x$ is distributed among the processors \\ \texttt{PARFFTY} & If \texttt{TRUE}, the dimension $y$ is distributed among the processors \\ \texttt{PARFFTZ} & If \texttt{TRUE}, the dimension $z$ is distributed among the processors \\ \texttt{MX} & Number of grid points in $x$ specifying rectangular grid \\ \texttt{MY} & Number of grid points in $y$ specifying rectangular grid \\ \texttt{MZ} & Number of grid points in $z$ specifying rectangular grid \\ \texttt{BCFFTX} & Boundary condition in $x$ [\texttt{OPEN}] \\ \texttt{BCFFTY} & Boundary condition in $y$ [\texttt{OPEN}] \\ \texttt{BCFFTZ} & Boundary condition in $z$ [\texttt{OPEN,PERIODIC}] \\ \texttt{GREENSF} & Defines the Greens function for the FFT Solver \\ \texttt{BBOXINCR} & Enlargement of the bounding box in \% \\ \texttt{GEOMETRY} & Geometry to be used as domain boundary \\ \texttt{ITSOLVER} & Type of iterative solver \\ \texttt{INTERPL} & Interpolation used for boundary points \\ \texttt{TOL} & Tolerance for iterative solver \\ \texttt{MAXITERS} & Maximum number of iterations of iterative solver \\ \texttt{PRECMODE} & Behavior of the preconditioner \\ \texttt{AMR\_MAXLEVEL} & Maximum adaptive mesh refinement level (single level if \texttt{AMR\_MAXLEVEL} is zero) \\ \texttt{AMR\_MAXGRID} & Maximum grid size.\ It has to be smaller than \texttt{MX}, \texttt{MY}, \texttt{MZ} when running in parallel \\ \texttt{AMR\_TAGGING} & Mesh-refinement strategy \\ \hline \end{tabular} \end{center} \end{table} \section{Define the Fieldsolver to be used} \label{sec:FSFSTYPE} \index{FSFSTYPE} 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. \section{Define Domain Decomposition} \label{sec:FSDomDEC} \index{FSDomDEC} The dimensions in $x$, $y$ and $z$ can be parallel (\texttt{TRUE}) or serial \texttt{FALSE}. The default settings are: parallel in $z$ and serial in $x$ and $y$. \section{Define Number of Grid Points} \label{sec:FSMX} \index{FSMX} Number of grid points in $x$, $y$ and $z$ for a rectangular grid. \section{Define Boundary Conditions} \label{sec:FSBC} \index{FSBC} Two boundary conditions can be selected independently among $x$, $y$ namely: \texttt{OPEN} and for $z$ \texttt{OPEN} \& \texttt{PERIODIC}\index{OPEN}\index{PERIODIC}. In the case you select for $z$ periodic you are about to model a DC-beam. \section{Define Greens Function} \label{sec:FSGREEN} \index{FSGREEN} Two Greens functions can be selected: \texttt{INTEGRATED}, \texttt{STANDARD}. The integrated Green's function is described in \ref{qiang2005, qiang2006-1, qiang2006-2}. Default setting is \texttt{INTEGRATED}. \section{Define Bounding Box Enlargement} \label{sec:FSBBOX} \index{FSBBOX} The bounding box defines a minimal rectangular domain including all particles. With \texttt{BBOXINCR} the bounding box can be enlarged by a factor given in percent of the minimal rectangular domain. \section{Define Geometry} \label{sec:GEOMETRY} \index{GEOMETRY} The list of geometries defining the beam line boundary. For further details see \chpref{geometry}. \section{Define Iterative Solver} \label{sec:ITSOLVER} \index{ITSOLVER} The iterative solver for solving the preconditioned system: \texttt{CG}\index{ITSOLVER!CG}, \texttt{BiCGSTAB}\index{ITSOLVER!BiCGSTAB} or \texttt{GMRES}\index{ITSOLVER!GMRES}. \section{Define Interpolation for Boundary Points} \label{sec:INTERPL} \index{INTERPL} The interpolation method for grid points near the boundary: \texttt{CONSTANT}\index{INTERPL!CONSTANT}, \texttt{LINEAR}\index{INTERPL!LINEAR} or \texttt{QUADRATIC}\index{INTERPL!QUADRATIC}. \section{Define Tolerance} \label{sec:TOL} \index{TOL} The tolerance for the iterative solver used by the \texttt{MG} solver. \section{Define Maximal Iterations} \label{sec:MAXITERS} \index{MAXITERS} The maximal number of iterations the iterative solver performs. \section{Define Preconditioner Behavior} \label{sec:PRECMODE} \index{PRECMODE} The behavior of the preconditioner can be: \texttt{STD}\index{PRECMODE!STD}, \texttt{HIERARCHY}\index{PRECMODE!HIERARCHY} or \texttt{REUSE}\index{PRECMODE!REUSE}. This argument is only relevant when using the \texttt{MG} solver and should \textbf{only be set if the consequences to simulation and solver are evident}. A short description is given in Table~\ref{preconditioner_behaviour}. \begin{table}[ht] \footnotesize \begin{center} \caption{Preconditioner behavior command summary} \label{tab:preconditioner_behaviour} \begin{tabular}{|l|p{0.6\textwidth}|} \hline \tabhead Value & Behavior \\ \hline \texttt{STD} & The preconditioner is rebuilt in every time step (enabled by default) \\ \texttt{HIERARCHY} & The hierarchy (tentative prolongator) is reused \\ \texttt{REUSE} & The preconditioner is reused \\ \hline \end{tabular} \end{center} \end{table} \section{Define the number of Energy Bins to use} \label{sec:FSENBINS} \index{FSENBINS} 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~\ref{cyclotron} the number of energy bins must be at minimum the number of neighboring bunches (\texttt{NNEIGHBB}) i.e. $\text{\texttt{ENBINS}} \le \text{\texttt{NNEIGHBB}}$. The variable \texttt{MINSTEPFORREBIN} defines the number of integration step that have to pass until all energy bins are merged into one. \section{Define AMR Solver} \label{sec:AMR} \index{AMR} The option \texttt{AMR=TRUE} enables further commands to be used in the Fieldsolver command.\ In order to run a proper AMR simulation one requires following ingredients: \begin{itemize} \item Upper bound of refinement, specified by \texttt{AMR\_MAXLEVEL} \item Maximum allowable grid size, specified by \texttt{AMR\_MAXGRID}.\ This size holds for all levels in all directions. \item Mesh refinement strategy, specified by \texttt{AMR\_TAGGING}.\ Depending on the tagging scheme, further keywords can be used.\ A summary is given in Table~\ref{tagging_strategies}. \end{itemize} \begin{table}[ht] \footnotesize \begin{center} \caption{Mesh-refinement strategies} \label{tab:tagging_strategies} \begin{tabular}{|l|p{0.6\textwidth}|} \hline \tabhead Value & Behavior \\ \hline \texttt{POTENTIAL} & Mark each cell if $\phi^{level}_{cell}\ge\alpha\max\phi^{level}$, where $\alpha\in[0, 1]$ is the scaling factor \texttt{AMR\_SCALING} \\ \texttt{EFIELD} & Mark each cell if the electric field component of any direction satisfies $E^{level}_{d, cell}\ge\alpha\max E_{d}^{level}$, where $d=x, y, z$ and $\alpha\in[0, 1]$ is the scaling factor \texttt{AMR\_SCALING} \\ \texttt{MOMENTA} & 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 $||\mathbf{p}|| \ge \alpha \max_{level} ||\mathbf{p}||$ is refined.\ The scalar $\alpha\in[0, 1]$ is a user-defined value \texttt{AMR\_SCALING}. \\ \texttt{CHARGE\_DENSITY} & If the charge density of a cell is greater or equal to the value specified with \texttt{AMD\_DENSITY} the cell is tagged for refinement\\ \texttt{MIN\_NUM\_PARTICLES} & Cells with equal or more particles are refined.\ The bound is specified with \texttt{AMR\_MIN\_NUM\_PART} \\ \texttt{MAX\_NUM\_PARTICLES} & Cells with equal or less particles are refined.\ The bound is specified with \texttt{AMR\_MAX\_NUM\_PART}\\ \hline \end{tabular} \end{center} \end{table} \index{Field Solver|)} \input{footer}