A Fast Poisson Solver of Second-order Accuracy for Isolated Systems in Three-dimensional Cartesian and Cylindrical Coordinates

We present an accurate and efficient method to calculate the gravitational potential of an isolated system in 3D Cartesian and cylindrical coordinates subject to vacuum (open) boundary conditions. Our method consists of two parts: an interior solver and a boundary solver. The interior solver adopts an eigenfunction expansion method together with a tridiagonal matrix solver to solve the Poisson equation subject to the zero boundary condition. The boundary solver employs James's method to calculate the boundary potential due to the screening charges required to keep the zero boundary condition for the interior solver. A full computation of gravitational potential requires running the interior solver twice and the boundary solver once. We develop a method to compute the discrete Green's function in cylindrical coordinates, which is an integral part of the James algorithm to maintain second-order accuracy. We implement our method in the Athena++ magnetohydrodynamics code and perform various tests to check that our solver is second-order accurate and exhibits good parallel performance.