fieldsolvers.asciidoc

Back to Main Page

Table of Contents

• 1. Field Solver
  • 1.1. FFT Based Particle-Mesh (PM) Solver
    • 1.1.1. FFT-based Convolutions and Zero Padding
    • 1.1.2. Algorithm used in OPAL
    • 1.1.3. Interpolation Schemes
  • 1.2. Iterative Space Charge Solver
  • 1.3. Energy Binning
  • 1.4. The FIELDSOLVER Command
  • 1.5. Define the Fieldsolver to be used
  • 1.6. Define Domain Decomposition
  • 1.7. Define Number of Grid Points
  • 1.8. Define Boundary Conditions
  • 1.9. Define Greens Function
  • 1.10. Define Bounding Box Enlargement
  • 1.11. Define Geometry
  • 1.12. Define Iterative Solver
  • 1.13. Define Interpolation for Boundary Points
  • 1.14. Define Tolerance
  • 1.15. Define Maximal Iterations
  • 1.16. Define Preconditioner Behavior
  • 1.17. Define the number of Energy Bins to use
  • 1.18. Define AMR Solver
    • 1.18.1. Hardware-Architecture Independent AMR Poisson Solver
    • 1.18.2. Use of AMR
    • 1.18.3. AMR Example
  • 1.19. References

1. Field Solver

Space charge effects are included in the simulation by specifying a field solver described in this chapter and attaching it to the track command as described in Chapter Tracking. 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.

1.1. FFT Based Particle-Mesh (PM) Solver

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 [1]. 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 [1], [2] and [3]. 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 [4] , [5].

The solution of the Poisson equation,

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

for the scalar potential,

$\phi$

, due to a charge density,

$\rho$

, and appropriate boundary conditions, can be expressed as,

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

where

$G(x,x',y,y',z,z')$

is the Green function, subject to the appropriate boundary conditions, describing the contribution of a source charge at location

$(x',y',z')$

to the potential at an observation location

$(x,y,z)$

.

For an isolated distribution of charge this reduces to

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

where

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

A simple discretization of Convolution solution on a Cartesian grid with cell size

$(h_x,h_y,h_z)$

leads to,

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

where

$\rho_{i,j,k}$

and

$G_{i-i',j-j',k-k'}$

denote the values of the charge density and the Green function, respectively, defined on the grid

$M$

.

1.1.1. 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}$$

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.

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,

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~\mathrm{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.

In the physical (unpadded) region,

j\in \left[0,J-1\right]

, so the quantity

j-n

in Final equation 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, FFT convolution, matches the convolution in Brute force convolution.

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.

The above FFT-based approach – zero-padding the charge density array, and circular-shifting the Green function in accordance with Periodic Green function – 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

) can be dropped, since it will be recovered implicitly through the symmetry of Periodic Green function. In that case the approach is identical to the Hockney method [1], [2] and [3]. Lastly, note that the above proof that the convolution, FFT convolution, is identical to Brute force convolution 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 FFT convolution. 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}

simply by changing the direction of the Fourier transform of the Green function and changing the direction of the final Fourier transform.

1.1.2. Algorithm used in OPAL

As a result, the solution of Open brute force convolution is then given by

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

where the notation has been introduced that

\text{FFT}\{ . \}

denotes a forward FFT in all 3 dimensions, and

\text{FFT}^{-1}\{ . \}

denotes a backward FFT in all 3 dimensions.

1.1.3. Interpolation Schemes

More details will be given in Version 2.0.0.

1.2. Iterative Space Charge Solver

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 [6].

1.3. Energy Binning

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

1.4. The FIELDSOLVER Command

See Table 1 for a summary of the Fieldsolver command.

Table 1. Fieldsolver command summary

Command	Purpose
FIELDSOLVER	Specify a fieldsolver
FSTYPE	Specify the type of field solver: FFT, FFTPERIODIC, SAAMG and NONE. Further arguments are enabled with the AMR solver (cf. Define AMR Solver).
PARFFTX	If TRUE, the dimension x is distributed among the processors
PARFFTY	If TRUE, the dimension y is distributed among the processors
PARFFTT	If TRUE, the dimension z is distributed among the processors
MX	Number of grid points in x specifying rectangular grid
MY	Number of grid points in y specifying rectangular grid
MT	Number of grid points in z specifying rectangular grid
BCFFTX	Boundary condition in x [OPEN]. (FFT + AMR_MG only)
BCFFTY	Boundary condition in y [OPEN]. (FFT + AMR_MG only)
BCFFTZ	Boundary condition in z [OPEN,PERIODIC]. (FFT + AMR_MG only)
GREENSF	Defines the Greens function for the FFT Solver. (FFT only)
BBOXINCR	Enlargement of the bounding box in %
GEOMETRY	Geometry to be used as domain boundary. (SAAMG only)
ITSOLVER	Type of iterative solver. (SAAMG + AMR_MG only)
INTERPL	Interpolation used for boundary points. (SAAMG only)
TOL	Tolerance for iterative solver. (SAAMG only)
MAXITERS	Maximum number of iterations of iterative solver. (SAAMG only)
PRECMODE	Behavior of the preconditioner. (SAAMG only)

1.5. Define the Fieldsolver to be used

At present FFT based solver is the most stable solver. Other solvers include Finite Element solvers and a Multigrid solver with Shortley-Weller boundary conditions for irregular domains.

1.6. Define Domain Decomposition

The dimensions in

x

,

y

and

z

can be parallel (TRUE) or serial FALSE. The default settings are: parallel in

z

and serial in

x

and

y

.

1.7. Define Number of Grid Points

Number of grid points in

x

,

y

and

z

for a rectangular grid.

1.8. Define Boundary Conditions

Two boundary conditions can be selected independently among

x

,

y

namely: OPEN and for

z

OPEN & PERIODIC. In the case you select for

z

periodic you are about to model a DC-beam.

1.9. Define Greens Function

Two Greens functions can be selected: INTEGRATED, STANDARD. The integrated Green’s function is described in [7], [8], [9]. Default setting is INTEGRATED.

1.10. Define Bounding Box Enlargement

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

1.11. Define Geometry

The list of geometries defining the beam line boundary. For further details see Chapter Geometry.

1.12. Define Iterative Solver

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

1.13. Define Interpolation for Boundary Points

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

1.14. Define Tolerance

The tolerance for the iterative solver used by the SAAMG solver. Default value is 1e-8.

1.15. Define Maximal Iterations

The maximal number of iterations the iterative solver performs. Default value is 100.

1.16. Define Preconditioner Behavior

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

Table 2. Preconditioner behavior command summary

Value	Behavior
STD	The preconditioner is rebuilt in every time step
HIERARCHY	The hierarchy (tentative prolongator) is reused (enabled by default)
REUSE	The preconditioner is reused

1.17. Define the number of Energy Bins to use

Suppose

\mathrm{d} E

the energy spread in the particle bunch is to large, the electrostatic approximation is no longer valid. One solution to that problem is to introduce

k

energy bins and perform

k

separate field solves in which

\mathrm{d} E

is again small and hence the electrostatic approximation valid. In case of a cyclotron see Cyclotron the number of energy bins must be at minimum the number of neighboring bunches (NNEIGHBB) i.e.

\mathrm{ENBINS} \le \mathrm{NNEIGHBB}

.

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

1.18. Define AMR Solver

After compiling OPAL with ENABLE_AMR=ON, the option AMR=TRUE enables further commands to be used in the Fieldsolver command (cf. Table 3).

The current AMR interface in OPAL is based on AMReX [10] which provides the multigrid solvers FMG (till release 18.07) and ML. We also provide a Trilinos based solver which is described in the section Hardware-Architecture Independent AMR Poisson Solver.

In order to setup an AMR simulation, the user has to specify the number of AMR levels. The hierarchy is built based on the values of the refinement ratios, blocking factors and maximum grid sizes of the base level. The maximum grid size of the next higher level is given by the division with the refinement ratio. The same is true for the blocking factors. The mesh refinement follows user-dependent rules that are provided by AMR_TAGGING. The various mesh refinement methods are provided in Table 4.

Table 3. AMR extensions to the FIELDSOLVER command

Command	Purpose
FSTYPE	AMR Poisson solvers: FMG (BoxLib), ML (AMReX) and AMR_MG
AMR_MAXLEVEL	Maximum adaptive mesh refinement level (single level if AMR_MAXLEVEL is zero)
AMR_REFX	Grid cell refinement ratio in x, only AMR_REFX=2 is tested
AMR_REFY	Grid cell refinement ratio in y, only AMR_REFY=2 is tested
AMR_REFZ	Grid cell refinement ratio in z, only AMR_REFZ=2 is tested
AMR_MAXGRIDX	Maximum grid size of base level in x. It has to be smaller than MX when running in parallel
AMR_MAXGRIDY	Maximum grid size of base level in y. It has to be smaller than MY when running in parallel
AMR_MAXGRIDZ	Maximum grid size of base level in z. It has to be smaller than MZ when running in parallel
AMR_BFX	AMReX blocking factor in x. AMR_MAXGRIDX has to be a multiple.
AMR_BFY	AMReX blocking factor in y. AMR_MAXGRIDY has to be a multiple.
AMR_BFZ	AMReX blocking factor in z. AMR_MAXGRIDZ has to be a multiple.
AMR_TAGGING	Mesh-refinement strategy [CHARGE_DENSITY, POTENTIAL, EFIELD, MOMENTA, MAX_NUM_PARTICLES, MIN_NUM_PARTICLES]
AMR_DENSITY	Charge density refinement threshold when AMR_TAGGING=CHARGE_DENSITY. See also Table 4.
AMR_MAX_NUM_PART	Refinement threshold for AMR_TAGGING=MAX_NUM_PARTICLES. See also Table 4.
AMR_MIN_NUM_PART	Refinement threshold for AMR_TAGGING=MIN_NUM_PARTICLES. See also Table 4.
AMR_SCALING	Refinement threshold for AMR_TAGGING=POTENTIAL and AMR_TAGGING=EFIELD. See also Table 4.
AMR_DOMAIN_RATIO	Box ratio of the AMR computation domain. It is an array of size 3. The entries define the ratio in (x, y, z) direction. For example, [1, 0.75, 0.75] yields a computation domain of [-1, 1]\times [-0.75, 0.75]\times [-0.75, 0.75] .

Table 4. Mesh refinement strategies

Value	Behavior
POTENTIAL	Mark each cell if \phi^{level}_{cell}\ge\alpha\max\phi^{level} , where \alpha\in[0, 1] is the scaling factor AMR_SCALING
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 AMR_SCALING
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 AMR_SCALING.
CHARGE_DENSITY	If the charge density of a cell is greater or equal to the value specified with AMD_DENSITY the cell is tagged for refinement
MIN_NUM_PARTICLES	Cells with equal or more particles are refined. The bound is specified with AMR_MIN_NUM_PART
MAX_NUM_PARTICLES	Cells with equal or less particles are refined. The bound is specified with AMR_MAX_NUM_PART

1.18.1. Hardware-Architecture Independent AMR Poisson Solver

This solvers is enabled with ENABLE_AMR_MG_SOLVER=ON and requires a working Trilinos installation with at least the following packages:

• Tpetra [11]

• Ifpack2 [12]

• Amesos2 and Belos [13]

• MueLu [14]

Some base level linear solvers may require additional third party libraries. These are SUPERLU [15], UMFPACK [16], PARDISO_MKL [17] [18], MUMPS [19] and LAPACK [20]. A full list of the arguments and linear solvers is given in Table 5. The solver is tested with Trilinos version 12.14.1. The corresponding paper is [21].

Table 5. Extension of the FIELDSOLVER command for the AMR Poisson Solver

Command	Purpose
AMR_MG_SMOOTHER	The pre- and post-smoother which is either [GS (Gauss-Seidel), JACOBI or SGS (symmetric Gauss-Seidel)]. Default: GS
AMR_MG_NSWEEPS	The number of smoothing steps. Default: 8
AMR_MG_PREC	The preconditioner of the bottom linear solver. Please see [12] for further information. [NONE, ILUT, CHEBYSHEV, RILUK, SA, JACOBI, BLOCK_JACOBI, GS, BLOCK_GS]. Default: NONE
AMR_MG_INTERP	The prolongator of the level solution to next higher level [PC (piecewise constant), TRILINEAR, LAGRANGE (coarse-fine interface only)]. Default: PC
AMR_MG_NORM	The norm of the convergence criterion [L1, L2, LINF]. Default: LINF
AMR_MG_VERBOSE	Boolean to enable/disable solver output. It writes an SDDS file. Default: FALSE
AMR_MG_REBALANCE	Boolean to rebalance the smoothed aggregation preconditioner or linear solver (if TRUE, it uses subcommunicators to reduce the communication overhead). Default: FALSE
AMR_MG_REBALANCE	Reuse type of the smoothed aggregation multigrid solver of MueLu. Please see [14] for further information. [NONE, RP, RAP, S, FULL]. Default: RAP
AMR_MG_TOL	The tolerance of the solver. Default: 10^{-10}
ITSOLVER	The solver of the linear system of equations on the base level. [BICGSTAB, MINRES, PCPG, CG, GMRES, STOCHASTIC_CG, RECYCLING_CG, RECYCLING_GMRES, KLU2, SUPERLU, UMFPACK, PARDISO_MKL, MUMPS, LAPACK, SA]. Please see [13] for further information. Default: CG

1.18.2. Use of AMR

AMR is only available in the multi-bunch mode (cf. [sec.opalcycl.MultiBunch]) of OPAL-cycl. As soon as more than one bunch is in the simulation, the AMR hierarchy is built. Note that AMReX handles the MPI parallelism, hence, other Poisson solvers of OPAL (cf. FSTYPE in Table 1) are not available.

1.18.3. AMR Example

In the example below we show an AMR simulation setup. It has a

[128]^3

base grid and creates two AMR levels. The maximum grid size of the base level is 32. Hence, we can have 4 MPI-processes per spatial dimension on level zero.

// global options
OPTION, AMR              = TRUE;
OPTION, AMR_REGRID_FREQ  = 10;
OPTION, AMR_YT_DUMP_FREQ = 100000;
...

// initialize kinetic energy etc.
REAL Edes      = .072;
REAL turns     = 8;
REAL nstep     = 360;
REAL frequency = 50.650;
...

// define the cyclotron
ring: Cyclotron, ...
rf0: RFCavity, ...
rf1: RFCavity, ...
rf2: RFCavity, ...
rf3: RFCavity, ...
rf4: RFCavity, ...
l1:   Line = (ring, rf0, rf1, rf2, rf3, rf4);

// define the distribution
Dist1: DISTRIBUTION, ...

// fieldsolver command
Fs1:FIELDSOLVER, FSTYPE=AMR_MG,
                 MX=128, MY=128, MT=128,
                 PARFFTX=true, PARFFTY=true, PARFFTT=true,
                 BCFFTX=open, BCFFTY=open, BCFFTT=open,
                 BBOXINCR=20, AMR_MAXLEVEL=2,
                 AMR_MAXGRIDX=32, AMR_MAXGRIDY=32, AMR_MAXGRIDZ=32,
                 AMR_BFX=16, AMR_BFY=16, AMR_BFZ=16,
                 AMR_REFX=2, AMR_REFY=2, AMR_REFZ=2,
                 AMR_DOMAIN_RATIO={1.0, 0.75, 0.75},
                 AMR_TAGGING=CHARGE_DENSITY, AMR_DENSITY=1.0e-9,
                 AMR_MG_VERBOSE=true, AMR_MG_REBALANCE=true, AMR_MG_REUSE=FULL,
                 ITSOLVER=SA, AMR_MG_NORM=LINF, AMR_MG_NSWEEPS=12;

beam1: BEAM, PARTICLE=PROTON, pc=P0, NPART=1e5, BCURRENT=2.0E-3, CHARGE=1.0, BFREQ= frequency;

Select, Line=l1;

TRACK,LINE=l1, BEAM=beam1, MAXSTEPS=nstep*turns, STEPSPERTURN=nstep,TIMEINTEGRATOR="RK-4";
    run, method = "CYCLOTRON-T", beam=beam1, fieldsolver=Fs1, distribution=Dist1,
    MBMODE=FORCE, TURNS=11, MB_BINNING=GAMMA, MB_ETA=0.25;

endtrack;

Stop;

1.19. References

[1] R. W. Hockney, The potential calculation and some applications, Methods Comput. Phys. 9, 136–210 (1970).

[2] J. W. Eastwood and D. R. K. Brownrigg, Remarks on the solution of poisson’s equation for isolated systems, J. Comput. Phys. 32, 24–38 (1979).

[3] R. W. Hockney and J. W. Eastwood, Computer simulation using particles (Taylor & Francis Group, 1988).

[4] J. Qiang and R. Ryne, Parallel 3d poisson solver for a charged beam in a conducting pipe, Comput. Phys. Commun. 138, 18–28 (2001).

[5] J. Qiang and R. L. Gluckstern, Three-dimensional poisson solver for a charged beam with large aspect ratio in a conducting pipe, Comput. Phys. Commun. 160, 120–128 (2004).

[6] A. Adelmann, P. Arbenz, and Y. Ineichen, A fast parallel poisson solver on irregular domains applied to beam dynamic simulations, arXiv (2009).

[7] J. Qiang et al., A three-dimensional quasi-static model for high brightness beam dynamics simulation, tech. rep. LBNL-59098 (LBNL).

[8] J. Qiang et al., Three-dimensional quasi-static model for high brightness beam dynamics simulation, Phys. Rev. ST Accel. Beams 9 (2006).

[9] J. Qiang et al., Erratum: three-dimensional quasi-static model for high brightness beam dynamics simulation, Phys. Rev. ST Accel. Beams 10, 12990 (2007).

[10] W. Zhang, A. Almgren, et al., AMReX: a framework for block-structured adaptive mesh refinement, Journal of Open Source Software, 4(37):1370 (2019).

[11] C. G. Baker and M. A. Heroux, Tpetra, and the use of generic programming in scientific computing, Scientific Programming, 20(2):115–128 (2012).

[12] A. Prokopenko, C. M. Siefert, J. J. Hu, M. Hoemmen, and A. Klinvex, Ifpack2 User’s Guide 1.0, Technical Report SAND2016-5338, Sandia National Labs (2016).

[13] E. Bavier, M. Hoemmen, S. Rajamanickam, and H. Thornquist, Amesos2 and Belos: Direct and iterative solvers for large sparse linear systems, Scientific Programming, 20(3):241–255 (2012).

[14] L. Berger-Vergiat, C. A. Glusa, J. J. Hu, M. Mayr, A. Prokopenko, C. M. Siefert, R. S. Tuminaro, and T. A. Wiesner, MueLu User’s Guide, Technical Report SAND2019-0537, Sandia National Laboratories (2019).

[15] X. Li, J. Demmel, J. Gilbert, L. Grigori, M. Shao, and I. Yamazaki, SuperLU Users' Guide, Technical Report LBNL-44289, Lawrence Berkeley National Laboratory (1999).

[16] T. A. Davis, Algorithm 832: UMFPACK V4.3—an Unsymmetric-Pattern Multifrontal Method, ACM Trans. Math. Softw., 30(2):196–199 (2004).

[17] O. Schenk and K. Gärtner, Solving unsymmetric sparse systems of linear equations with PARDISO, Future Generation Computer Systems, 20(3):475–487 (2004).

[18] Developer Guide for Intel® oneAPI Math Kernel Library for Linux, last updated: February 7, 2020.

[19] P. R. Amestoy, I. S. Duff, J.-Y. L’Excellent, and J. Koster, A Fully Asynchronous Multifrontal Solver Using Distributed Dynamic Scheduling, SIAM Journal on Matrix Analysis and Applications, 23(1):15–41 (2001).

[20] E. Anderson, Z. Bai, C. Bischof, L. S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen, LAPACK Users' Guide, SIAM, third edition, 1999, ISBN 978-0-89871-447-0.

[21] M. Frey, A. Adelmann, and U. Locans, On architecture and performance of adaptive mesh refinement in an electrostatics Particle-In-Cell code, Computer Physics Communications, 247:106912 (2020).