fieldsolvers.tex 28.3 KB
Newer Older
snuverink_j's avatar
snuverink_j committed
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's avatar
snuverink_j committed
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's avatar
snuverink_j committed
18 19 20 21

\section{FFT Based Particle-Mesh (PM) Solver}
\label{sec:fieldsolver:fftbased}
\index{Field Solver!FFT based}
snuverink_j's avatar
snuverink_j committed
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's avatar
snuverink_j committed
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's avatar
snuverink_j committed
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's avatar
snuverink_j committed
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's avatar
snuverink_j committed
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's avatar
snuverink_j committed
32 33 34 35 36 37 38 39 40 41 42

The solution of the Poisson equation,


\begin{equation}
%\label{ }
\nabla^2\phi=-\rho/\epsilon_0,
\end{equation}
for the scalar potential, $\phi$, due to a charge density, $\rho$, and appropriate boundary conditions, can be expressed as,

\begin{equation}
snuverink_j's avatar
snuverink_j committed
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's avatar
snuverink_j committed
44 45 46 47 48 49
\end{equation}
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

\begin{equation}
snuverink_j's avatar
snuverink_j committed
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's avatar
snuverink_j committed
51 52 53 54 55 56 57 58
\label{eq:convolutionsolution}
\end{equation}
where

\begin{equation}
G(u,v,w)={\frac{1}{\sqrt{u^2+v^2+w^2}}}.
\label{eq:isolatedgreenfunction}
\end{equation}
snuverink_j's avatar
snuverink_j committed
59
A simple discretization of Equation~\ref{convolutionsolution}
snuverink_j's avatar
snuverink_j committed
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,

\begin{equation}
\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}
\end{equation}
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,

\begin{equation}
\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}
\end{equation}
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$,
\begin{equation}
\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.
\end{equation}
Define a periodic Green function, $G_m$, as follows,
\begin{equation}
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
\end{equation}
Now consider the sum
\begin{equation}
{\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}
\end{equation}
where $W=e^{-2\pi i/N}$. This is just the FFT-based convolution of $\{\rho_k\}$ with $\{G_m\}$.
Then,
\begin{equation}
{\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.
\end{equation}
Now use the relation
\begin{equation}
\sum_{k=0}^{N-1} W^{(m+n-j)k}= N \delta_{m+n-j,iN}~~~~~(i~\rm an~integer).
\end{equation}

It follows that
\begin{equation}
{\phi}_j=\sum_{n=0}^{K-1}~\bar{\rho}_n G_{j-n+iN}
~~~~~~0 \le j \le N-1.
\end{equation}
But $G$ is periodic with period $N$. Hence,
\begin{equation}
{\phi}_j=\sum_{n=0}^{K-1}~\bar{\rho}_n G_{j-n}
~~~~~~0 \le j \le N-1.
\label{eq:finaleqn}
\end{equation}
snuverink_j's avatar
snuverink_j committed
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's avatar
snuverink_j committed
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's avatar
snuverink_j committed
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's avatar
snuverink_j committed
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's avatar
snuverink_j committed
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's avatar
snuverink_j committed
153 154
%for both positive and negative arguments.

snuverink_j's avatar
snuverink_j committed
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's avatar
snuverink_j committed
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's avatar
snuverink_j committed
159
can be dropped, since it will be recovered implicitly through the symmetry of Equation~\ref{periodicgreenfunction}.
snuverink_j's avatar
snuverink_j committed
160
In that case the approach is identical to the Hockney method \ref{hockney, eastwoodandbrownrigg,hockneyandeastwood}.
snuverink_j's avatar
snuverink_j committed
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's avatar
snuverink_j committed
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's avatar
snuverink_j committed
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
166
%periodic Green function used in the Hockney method. Viewed this way, the only difference between an "ordinary" FFT-based convolution
snuverink_j's avatar
snuverink_j committed
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's avatar
snuverink_j committed
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's avatar
snuverink_j committed
181 182 183 184 185 186 187 188 189 190 191
\begin{equation}
\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}
\end{equation}
simply by changing the direction of the Fourier transform of the Green function and changing the direction of the final Fourier transform.

192
\subsubsection{Algorithm used in \textit{OPAL}}
snuverink_j's avatar
snuverink_j committed
193

snuverink_j's avatar
snuverink_j committed
194
As a result, the solution of Equation~\ref{openbruteforceconvolution} is then given by
snuverink_j's avatar
snuverink_j committed
195 196 197 198 199 200 201 202 203

\begin{equation}
\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}
\end{equation}
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's avatar
snuverink_j committed
204
%In order to obtain  ${\cal M}_2$ in Equation~\ref{splitOper1} we must solve Poisson's equation,
snuverink_j's avatar
snuverink_j committed
205 206 207
%where $\rho$ stands for the
%charge density and $\phi$ for the scalar electrostatic potential:
%\begin{equation}\label{eq:Poisson}
snuverink_j's avatar
snuverink_j committed
208
%\nabla^2 \phi(\mathbf{q}) = - \frac{\rho(\mathbf{q})}{\epsilon_0}
snuverink_j's avatar
snuverink_j committed
209
%\end{equation}
snuverink_j's avatar
snuverink_j committed
210
%subject to open boundary conditions in all spatial directions: $\phi (\mathbf{q}) \rightarrow 0$ as $|\mathbf{q}| \rightarrow \infty$ or
211
%imposing periodic boundary conditions in longitudinal directions. The assumption of using an "isolated system" is physically
snuverink_j's avatar
snuverink_j committed
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's avatar
snuverink_j committed
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:
%\begin{equation}
snuverink_j's avatar
snuverink_j committed
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's avatar
snuverink_j committed
218 219 220
%\end{equation}
%where $G$ is the Green's function which gives the response to a unit source term. In 3D we have
%\begin{equation}
snuverink_j's avatar
snuverink_j committed
221
%G(\mathbf{q}-\mathbf{q}\primed) = \dis\frac{1}{4 \pi \,|\mathbf{q} - \mathbf{q}\primed|} \,.
snuverink_j's avatar
snuverink_j committed
222 223 224
%\end{equation}
%The electric field then follows from the electrostatic potential
%\begin{equation}\label{eq:Electric}
snuverink_j's avatar
snuverink_j committed
225
%\mathbf{E} = - \nabla \phi.
snuverink_j's avatar
snuverink_j committed
226 227 228 229 230 231 232 233 234
%\end{equation}

%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's avatar
snuverink_j committed
235
%The solution of the discretized Poisson equation with $\mathbf{k}=(l,n,m,)$
snuverink_j's avatar
snuverink_j committed
236
%\begin{equation}\label{eq:DiscretizedPoisson}
snuverink_j's avatar
snuverink_j committed
237
%\nabla^{2} \phi^D(\mathbf{k}) = - \frac{\rho^D(\mathbf{k})}{\epsilon_0}, \mathbf{k} \in \Omega^D.
snuverink_j's avatar
snuverink_j committed
238 239 240 241 242 243 244
%\end{equation}
%$\phi^D$ then is given by convolution with the appropriate discretized Green's function $G_D$:
%\begin{equation}
%\phi^D = \rho^D * G^D.
%\end{equation}
%In Fourier space (hats) the convolution becomes a simple multiplication, with
%\[
snuverink_j's avatar
snuverink_j committed
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's avatar
snuverink_j committed
246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263
%\]
%and we get:
%\begin{equation}\label{eq:FourierPoisson}
%\widehat{\phi}^D = \widehat{\rho}^D \cdot \widehat{G}^D.
%\end{equation}
%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's avatar
snuverink_j committed
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's avatar
snuverink_j committed
266 267
%\end{tabbing}
%\subsubsection{Open and Periodic Boundary Conditions}
snuverink_j's avatar
snuverink_j committed
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's avatar
snuverink_j committed
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's avatar
snuverink_j committed
271
%region of the calculation (see \ref{hockney} on p. 213).
snuverink_j's avatar
snuverink_j committed
272
\subsubsection{Interpolation Schemes}
snuverink_j's avatar
snuverink_j committed
273
\index{Field Solver!Interpolation Schemes} More details will be given in Version 2.0.0.
snuverink_j's avatar
snuverink_j committed
274 275 276

%\label{sec:interpol}
%Both charge assignment and electric field interpolation are related to the interpolation
snuverink_j's avatar
snuverink_j committed
277
%scheme used. A detailed discussion is given in~\ref{hockney}.
snuverink_j's avatar
snuverink_j committed
278
%If $e_i$ is the charge of a particle, we can write the density at mesh point $\mathbf{k}_m$ as
snuverink_j's avatar
snuverink_j committed
279
%\begin{equation}\label{eq:discRho}
snuverink_j's avatar
snuverink_j committed
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's avatar
snuverink_j committed
281 282 283 284 285 286 287
%\end{equation}
%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's avatar
snuverink_j committed
288
%CIC interpolation scheme is shown in Figure~\ref{CIC} for the two-dimensional case.
snuverink_j's avatar
snuverink_j committed
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's avatar
snuverink_j committed
307
  preconditioning.  More details are given in \ref{Adelmann:2009p543}.
snuverink_j's avatar
snuverink_j committed
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's avatar
snuverink_j committed
315
 More details will be given in Version 2.0.0.
snuverink_j's avatar
snuverink_j committed
316

snuverink_j's avatar
snuverink_j committed
317
\section{The \texttt{FIELDSOLVER} Command}
snuverink_j's avatar
snuverink_j committed
318 319 320
\label{sec:fieldsolvercmd}
\index{Field Solver!Command}
\index{FIELDSOLVER}
snuverink_j's avatar
snuverink_j committed
321
See Table~\ref{fieldsolvercmd} for a summary of the Fieldsolver command.
snuverink_j's avatar
snuverink_j committed
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's avatar
snuverink_j committed
328
      \tabhead Command & Purpose \\
snuverink_j's avatar
snuverink_j committed
329
      \hline
snuverink_j's avatar
snuverink_j committed
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's avatar
snuverink_j committed
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's avatar
snuverink_j committed
369
The dimensions in  $x$, $y$ and $z$ can be parallel (\texttt{TRUE})  or serial \texttt{FALSE}. The
snuverink_j's avatar
snuverink_j committed
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's avatar
snuverink_j committed
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's avatar
snuverink_j committed
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's avatar
snuverink_j committed
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's avatar
snuverink_j committed
387 388 389 390

\section{Define Bounding Box Enlargement}
\label{sec:FSBBOX}
\index{FSBBOX}
snuverink_j's avatar
snuverink_j committed
391
The bounding box defines a minimal rectangular domain including all particles. With \texttt{BBOXINCR}
snuverink_j's avatar
snuverink_j committed
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's avatar
snuverink_j committed
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's avatar
snuverink_j committed
403 404 405 406

\section{Define Interpolation for Boundary Points}
\label{sec:INTERPL}
\index{INTERPL}
snuverink_j's avatar
snuverink_j committed
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's avatar
snuverink_j committed
408 409 410 411

\section{Define Tolerance}
\label{sec:TOL}
\index{TOL}
snuverink_j's avatar
snuverink_j committed
412
The tolerance for the iterative solver used by the \texttt{MG} solver.
snuverink_j's avatar
snuverink_j committed
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's avatar
snuverink_j committed
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's avatar
snuverink_j committed
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's avatar
snuverink_j committed
430
      \tabhead Value & Behavior \\
snuverink_j's avatar
snuverink_j committed
431
      \hline
snuverink_j's avatar
snuverink_j committed
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's avatar
snuverink_j committed
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's avatar
snuverink_j committed
443
Suppose $\mathrm{d} E$ the energy spread in the particle bunch is to large, the electrostatic approximation is no longer valid.
snuverink_j's avatar
snuverink_j committed
444
One solution to that problem is to introduce  $k$ energy bins  and perform $k$ separate field solves
snuverink_j's avatar
snuverink_j committed
445
in which $\mathrm{d} E$ is again small and hence the electrostatic approximation valid. In case of a cyclotron
snuverink_j's avatar
snuverink_j committed
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's avatar
snuverink_j committed
447

snuverink_j's avatar
snuverink_j committed
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's avatar
snuverink_j committed
449 450 451 452 453


\section{Define AMR Solver}
\label{sec:AMR}
\index{AMR}
snuverink_j's avatar
snuverink_j committed
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's avatar
snuverink_j committed
455 456
simulation one requires following ingredients:
\begin{itemize}
snuverink_j's avatar
snuverink_j committed
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's avatar
snuverink_j committed
460
    used.\ A summary is given in Table~\ref{tagging_strategies}.
snuverink_j's avatar
snuverink_j committed
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's avatar
snuverink_j committed
470
      \tabhead Value & Behavior \\
snuverink_j's avatar
snuverink_j committed
471
      \hline
snuverink_j's avatar
snuverink_j committed
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's avatar
snuverink_j committed
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's avatar
snuverink_j committed
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's avatar
snuverink_j committed
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's avatar
snuverink_j committed
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's avatar
snuverink_j committed
486 487 488 489 490 491 492
      \hline
    \end{tabular}
  \end{center}
\end{table}

\index{Field Solver|)}
\input{footer}