fieldsolvers.tex 28.3 KB
 snuverink_j committed Sep 06, 2017 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 \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  snuverink_j committed Sep 06, 2017 17 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.  snuverink_j committed Sep 06, 2017 18 19 20 21  \section{FFT Based Particle-Mesh (PM) Solver} \label{sec:fieldsolver:fftbased} \index{Field Solver!FFT based}  snuverink_j committed Sep 08, 2017 22 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}.  snuverink_j committed Sep 06, 2017 23 24 25 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.  snuverink_j committed Sep 08, 2017 26 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  snuverink_j committed Sep 06, 2017 27 28 29 30 (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$.  snuverink_j committed Sep 08, 2017 31 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}.  snuverink_j committed Sep 06, 2017 32 33 34 35 36 37 38 39 40 41 42  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,  snuverink_j committed Sep 08, 2017 43 \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'),  snuverink_j committed Sep 06, 2017 44 45 46 47 48 49  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  snuverink_j committed Sep 08, 2017 50 \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'),  snuverink_j committed Sep 06, 2017 51 52 53 54 55 56 57 58 \label{eq:convolutionsolution} where G(u,v,w)={\frac{1}{\sqrt{u^2+v^2+w^2}}}. \label{eq:isolatedgreenfunction}  snuverink_j committed Sep 06, 2017 59 A simple discretization of Equation~\ref{convolutionsolution}  snuverink_j committed Sep 06, 2017 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 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}  snuverink_j committed Sep 06, 2017 137 138 139 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}.  snuverink_j committed Sep 06, 2017 140 141 142 143 144 145 146  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.  snuverink_j committed Sep 06, 2017 147 %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,  snuverink_j committed Sep 06, 2017 148 149 150 151 %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,  snuverink_j committed Sep 06, 2017 152 %and which have a grid containing the Green function ({\it e.g.}, Equation~\ref{isolatedgreenfunction}, Equation~\ref{rgreenfunction}, or Equation~\ref{rintgreenfunction})  snuverink_j committed Sep 06, 2017 153 154 %for both positive and negative arguments.  snuverink_j committed Sep 06, 2017 155 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.  snuverink_j committed Sep 06, 2017 156 157 158 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$)  snuverink_j committed Sep 06, 2017 159 can be dropped, since it will be recovered implicitly through the symmetry of Equation~\ref{periodicgreenfunction}.  snuverink_j committed Sep 08, 2017 160 In that case the approach is identical to the Hockney method \ref{hockney, eastwoodandbrownrigg,hockneyandeastwood}.  snuverink_j committed Sep 06, 2017 161 162 163 %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.  snuverink_j committed Sep 08, 2017 164 %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  snuverink_j committed Sep 06, 2017 165 %center of the grid to the corner corresponding to the origin of grid points, produces a periodic Green function that is identical to the  snuverink_j committed Sep 12, 2017 166 %periodic Green function used in the Hockney method. Viewed this way, the only difference between an "ordinary" FFT-based convolution  snuverink_j committed Sep 06, 2017 167 168 169 170 171 172 173 174 175 176 177 178 %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.  snuverink_j committed Sep 06, 2017 179 180 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  snuverink_j committed Sep 06, 2017 181 182 183 184 185 186 187 188 189 190 191  \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.  snuverink_j committed Sep 07, 2017 192 \subsubsection{Algorithm used in \textit{OPAL}}  snuverink_j committed Sep 06, 2017 193   snuverink_j committed Sep 06, 2017 194 As a result, the solution of Equation~\ref{openbruteforceconvolution} is then given by  snuverink_j committed Sep 06, 2017 195 196 197 198 199 200 201 202 203  \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.  snuverink_j committed Sep 06, 2017 204 %In order to obtain ${\cal M}_2$ in Equation~\ref{splitOper1} we must solve Poisson's equation,  snuverink_j committed Sep 06, 2017 205 206 207 %where $\rho$ stands for the %charge density and $\phi$ for the scalar electrostatic potential: %\label{eq:Poisson}  snuverink_j committed Sep 08, 2017 208 %\nabla^2 \phi(\mathbf{q}) = - \frac{\rho(\mathbf{q})}{\epsilon_0}  snuverink_j committed Sep 06, 2017 209 %  snuverink_j committed Sep 08, 2017 210 %subject to open boundary conditions in all spatial directions: $\phi (\mathbf{q}) \rightarrow 0$ as $|\mathbf{q}| \rightarrow \infty$ or  snuverink_j committed Sep 12, 2017 211 %imposing periodic boundary conditions in longitudinal directions. The assumption of using an "isolated system" is physically  snuverink_j committed Sep 06, 2017 212 %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}.  snuverink_j committed Sep 06, 2017 213 214 215 216 %The computational domain $\Omega \subset \R^3$ is simply connected and has a cylindrical %or rectilinear shape. %The corresponding integral equation reads: %  snuverink_j committed Sep 08, 2017 217 %\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  snuverink_j committed Sep 06, 2017 218 219 220 % %where $G$ is the Green's function which gives the response to a unit source term. In 3D we have %  snuverink_j committed Sep 08, 2017 221 %G(\mathbf{q}-\mathbf{q}\primed) = \dis\frac{1}{4 \pi \,|\mathbf{q} - \mathbf{q}\primed|} \,.  snuverink_j committed Sep 06, 2017 222 223 224 % %The electric field then follows from the electrostatic potential %\label{eq:Electric}  snuverink_j committed Sep 08, 2017 225 %\mathbf{E} = - \nabla \phi.  snuverink_j committed Sep 06, 2017 226 227 228 229 230 231 232 233 234 % %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$.  snuverink_j committed Sep 08, 2017 235 %The solution of the discretized Poisson equation with $\mathbf{k}=(l,n,m,)$  snuverink_j committed Sep 06, 2017 236 %\label{eq:DiscretizedPoisson}  snuverink_j committed Sep 08, 2017 237 %\nabla^{2} \phi^D(\mathbf{k}) = - \frac{\rho^D(\mathbf{k})}{\epsilon_0}, \mathbf{k} \in \Omega^D.  snuverink_j committed Sep 06, 2017 238 239 240 241 242 243 244 % %$\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 %$ snuverink_j committed Sep 08, 2017 245 %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}  snuverink_j committed Sep 06, 2017 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 %$ %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$ \\  snuverink_j committed Sep 08, 2017 264 265 %\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$ \\  snuverink_j committed Sep 06, 2017 266 267 %\end{tabbing} %\subsubsection{Open and Periodic Boundary Conditions}  snuverink_j committed Sep 08, 2017 268 %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.  snuverink_j committed Sep 06, 2017 269 270 %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  snuverink_j committed Sep 08, 2017 271 %region of the calculation (see \ref{hockney} on p. 213).  snuverink_j committed Sep 06, 2017 272 \subsubsection{Interpolation Schemes}  snuverink_j committed Sep 06, 2017 273 \index{Field Solver!Interpolation Schemes} More details will be given in Version 2.0.0.  snuverink_j committed Sep 06, 2017 274 275 276  %\label{sec:interpol} %Both charge assignment and electric field interpolation are related to the interpolation  snuverink_j committed Sep 08, 2017 277 %scheme used. A detailed discussion is given in~\ref{hockney}.  snuverink_j committed Sep 08, 2017 278 %If $e_i$ is the charge of a particle, we can write the density at mesh point $\mathbf{k}_m$ as  snuverink_j committed Sep 06, 2017 279 %\label{eq:discRho}  snuverink_j committed Sep 08, 2017 280 %\rho(\mathbf{k}_m)^D = \sum_{i=1}^N e_i\cdot W(\mathbf{q}_i,\mathbf{k}_m), ~ m=1\dots M  snuverink_j committed Sep 06, 2017 281 282 283 284 285 286 287 % %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  snuverink_j committed Sep 06, 2017 288 %CIC interpolation scheme is shown in Figure~\ref{CIC} for the two-dimensional case.  snuverink_j committed Sep 06, 2017 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 %%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)  snuverink_j committed Sep 08, 2017 307  preconditioning. More details are given in \ref{Adelmann:2009p543}.  snuverink_j committed Sep 06, 2017 308 309 310 311 312 313 314  \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.  snuverink_j committed Sep 06, 2017 315  More details will be given in Version 2.0.0.  snuverink_j committed Sep 06, 2017 316   snuverink_j committed Sep 08, 2017 317 \section{The \texttt{FIELDSOLVER} Command}  snuverink_j committed Sep 06, 2017 318 319 320 \label{sec:fieldsolvercmd} \index{Field Solver!Command} \index{FIELDSOLVER}  snuverink_j committed Sep 06, 2017 321 See Table~\ref{fieldsolvercmd} for a summary of the Fieldsolver command.  snuverink_j committed Sep 06, 2017 322 323 324 325 326 327 \begin{table}[ht] \footnotesize \begin{center} \caption{Fieldsolver command summary} \label{tab:fieldsolvercmd} \begin{tabular}{|l|p{0.6\textwidth}|l|} \hline  snuverink_j committed Sep 11, 2017 328  \tabhead Command & Purpose \\  snuverink_j committed Sep 06, 2017 329  \hline  snuverink_j committed Sep 11, 2017 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353  \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 \\  snuverink_j committed Sep 06, 2017 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368  \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}  snuverink_j committed Sep 06, 2017 369 The dimensions in $x$, $y$ and $z$ can be parallel (\texttt{TRUE}) or serial \texttt{FALSE}. The  snuverink_j committed Sep 06, 2017 370 371 372 373 374 375 376 377 378 379 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}  snuverink_j committed Sep 06, 2017 380 Two boundary conditions can be selected independently among $x$, $y$ namely: \texttt{OPEN} and for $z$ \texttt{OPEN} \& \texttt{PERIODIC}\index{OPEN}\index{PERIODIC}.  snuverink_j committed Sep 06, 2017 381 382 383 384 385 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}  snuverink_j committed Sep 08, 2017 386 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}.  snuverink_j committed Sep 06, 2017 387 388 389 390  \section{Define Bounding Box Enlargement} \label{sec:FSBBOX} \index{FSBBOX}  snuverink_j committed Sep 06, 2017 391 The bounding box defines a minimal rectangular domain including all particles. With \texttt{BBOXINCR}  snuverink_j committed Sep 06, 2017 392 393 394 395 396 397 398 399 400 401 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}  snuverink_j committed Sep 06, 2017 402 The iterative solver for solving the preconditioned system: \texttt{CG}\index{ITSOLVER!CG}, \texttt{BiCGSTAB}\index{ITSOLVER!BiCGSTAB} or \texttt{GMRES}\index{ITSOLVER!GMRES}.  snuverink_j committed Sep 06, 2017 403 404 405 406  \section{Define Interpolation for Boundary Points} \label{sec:INTERPL} \index{INTERPL}  snuverink_j committed Sep 06, 2017 407 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}.  snuverink_j committed Sep 06, 2017 408 409 410 411  \section{Define Tolerance} \label{sec:TOL} \index{TOL}  snuverink_j committed Sep 06, 2017 412 The tolerance for the iterative solver used by the \texttt{MG} solver.  snuverink_j committed Sep 06, 2017 413 414 415 416 417 418 419 420 421  \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}  snuverink_j committed Sep 06, 2017 422 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}.  snuverink_j committed Sep 06, 2017 423 424 425 426 427 428 429  \begin{table}[ht] \footnotesize \begin{center} \caption{Preconditioner behavior command summary} \label{tab:preconditioner_behaviour} \begin{tabular}{|l|p{0.6\textwidth}|} \hline  snuverink_j committed Sep 11, 2017 430  \tabhead Value & Behavior \\  snuverink_j committed Sep 06, 2017 431  \hline  snuverink_j committed Sep 06, 2017 432 433 434  \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 \\  snuverink_j committed Sep 06, 2017 435 436 437 438 439 440 441 442  \hline \end{tabular} \end{center} \end{table} \section{Define the number of Energy Bins to use} \label{sec:FSENBINS} \index{FSENBINS}  snuverink_j committed Sep 11, 2017 443 Suppose $\mathrm{d} E$ the energy spread in the particle bunch is to large, the electrostatic approximation is no longer valid.  snuverink_j committed Sep 06, 2017 444 One solution to that problem is to introduce $k$ energy bins and perform $k$ separate field solves  snuverink_j committed Sep 11, 2017 445 in which $\mathrm{d} E$ is again small and hence the electrostatic approximation valid. In case of a cyclotron  snuverink_j committed Sep 06, 2017 446 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}}$.  snuverink_j committed Sep 06, 2017 447   snuverink_j committed Sep 06, 2017 448 The variable \texttt{MINSTEPFORREBIN} defines the number of integration step that have to pass until all energy bins are merged into one.  snuverink_j committed Sep 06, 2017 449 450 451 452 453  \section{Define AMR Solver} \label{sec:AMR} \index{AMR}  snuverink_j committed Sep 06, 2017 454 The option \texttt{AMR=TRUE} enables further commands to be used in the Fieldsolver command.\ In order to run a proper AMR  snuverink_j committed Sep 06, 2017 455 456 simulation one requires following ingredients: \begin{itemize}  snuverink_j committed Sep 06, 2017 457 458 459  \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  snuverink_j committed Sep 06, 2017 460  used.\ A summary is given in Table~\ref{tagging_strategies}.  snuverink_j committed Sep 06, 2017 461 462 463 464 465 466 467 468 469 \end{itemize} \begin{table}[ht] \footnotesize \begin{center} \caption{Mesh-refinement strategies} \label{tab:tagging_strategies} \begin{tabular}{|l|p{0.6\textwidth}|} \hline  snuverink_j committed Sep 11, 2017 470  \tabhead Value & Behavior \\  snuverink_j committed Sep 06, 2017 471  \hline  snuverink_j committed Sep 06, 2017 472 473 474  \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  snuverink_j committed Sep 06, 2017 475  $E^{level}_{d, cell}\ge\alpha\max E_{d}^{level}$, where $d=x, y, z$ and $\alpha\in[0, 1]$ is the scaling factor  snuverink_j committed Sep 06, 2017 476 477  \texttt{AMR\_SCALING} \\ \texttt{MOMENTA} & It performs a loop over all particles of a level and computes the dot product of the momenta.\ Every  snuverink_j committed Sep 08, 2017 478  cell that contains a particle with $||\mathbf{p}|| \ge \alpha \max_{level} ||\mathbf{p}||$ is refined.\ The scalar $\alpha\in[0, 1]$  snuverink_j committed Sep 06, 2017 479 480 481 482 483 484 485  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}\\  snuverink_j committed Sep 06, 2017 486 487 488 489 490 491 492  \hline \end{tabular} \end{center} \end{table} \index{Field Solver|)} \input{footer}`