3D Finite Difference

A 3D Finite Difference space charge solver is optionally available in synergia3.

Theory

Let us remind ourselves that the Poisson Equation in 3D is : $$\Delta\phi = -\rho / \epsilon_0$$ where $\phi$ refers to the potential, $\rho$ refers to the space charge density and $\epsilon_0$ is a normalization factor. The domain is: $-L_x/L_y/L_z < x,y,z < L_x/L_y/L_z$, with dirichlet boundary conditions (i.e. potential is 0 at all boundaries).

Note that the scaled equation is implemented as it reduces the condition number of the matrix with $$\textrm{scaling factor} = \textrm{area of unit cell} = 1/(L_x * L_y * L_z)$$

Options

Option	Type	Default	Description
shape	std::array<int, 3>	{32,32,64}	shape of the domain
scale_thresholds	std::array<double, 3>	{15,15,2.5}	threshold at which preconditioner is rebuilt
domain_fixed	bool	false	use a fixed domain?
n_sigma	double	8	size of domain is n_sigma times stddev along each dimension
kick_scale	double	1	scale kicks by this factor
comm_group_size	int	1	number of MPI ranks for each subcomminicator

To fix the domain, call the following function on the options object:

set_fixed_domain(std::array<double, 3> offset, std::array<double, 3> size)

with the indended offset and size.

Using this solver

To use the solver, one needs to have PETSc installed. The only dependencies required for PETSc are : an MPI installation and a set of BLAS/LAPACK libraries (Fortran compilers are not necessary!). PETSc can bootstrap these if needed as described here. For use on NVIDIA GPUs, please configure PETSc with CUDA support.

A second route is to use the spack to install it and other dependencies required by Synergia3. For convience, we provide docker containers with all dependencies installed for x86_64_v2 ISA here. These docker containers are generated from spack environments that can be changed as needed (available here).

Finally, note that the installations of synergia3 on Wilson Cluster alreaedy include an appropriate build of PETSc included with them.

Implementation details

A rough overview of the control flow for this solver is as follows:

graph LR;
    subgraph subcomm-operations
    B[Recompute matrix];
    end
    subgraph global-operations
    A[Update Domain]
    F[Gather charge density on all subcomms]
    end
    subgraph local-operations
    A[Update Domain]-->E(Deposit bunch)
    I(Electric fields calculation)
    J(Apply kick to particles in bunch)
    end
    A[Update Domain]-->B[Recompute matrix];
    subgraph subcomm-operations
    B[Recompute matrix]-->C[Reuse Preconditioner];
    B[Recompute matrix]-->D[Recompute Preconditioner];
    C[Reuse Preconditioner]-->G[Solve linear system]
    D[Recompute Preconditioner]-->G[Solve linear system];
    E(Deposit bunch)-->F[Gather charge density on all subcomms];
    F[Gather charge density on all subcomms]-->G[Solve linear system];
    G[Solve linear system]-->H[Scatter phi to constituent ranks];
    end
    subgraph local-operations
    H[Scatter phi to constituent ranks]-->I(Electric fields calculation);
    I(Electric fields calculation)-->J(Apply kick to particles in bunch);
    end

Loading

MPI communications within the solver

There are two data transferst that occur within each solve:

• Charge density accumunation on each subcommunicator
• Phi accumulation on each MPI rank within a subcommunicator.

Both of these transfers are implemented using PETSc Star Forest [1][2][3] via the VecScatter routines [4]. The first transfer involves performs only as much communication as is necessary when compared to the MPI_AllReduce based method used in the 3D Open Hockney solver.

Linear solver details

Algebraic multigrid preconditioning [5][6][7] is used to solve the linear system.

An attempt to use a direct solver (superlu-dist) was made, but owing to an inability of the solver to use matrix/vectors directly from GPU memory, it was not pursued:[8].