A Fast, Second-Order Accurate Poisson Solver in Spherical Polar Coordinates

We present an efficient and accurate algorithm for solving the Poisson equation in spherical polar coordinates with a logarithmic radial grid and open boundary conditions. The method employs a divide-and-conquer strategy, decomposing the computational domain into hierarchical units with varying cell sizes. James's algorithm is used to compute the zero-boundary potentials of lower-level units, which are systematically integrated to reconstruct the zero-boundary potential over the entire domain. These calculations are performed efficiently via matrix-vector operations using various precomputed kernel matrices. The open-boundary potential is then obtained by applying a discrete Green's function to the effective screening density induced at the domain boundaries. The overall algorithm achieves a computational complexity of $\mathcal{O}(N^3 \log N)$, where $N$ denotes the number of cells in one dimension. We implement the solver in the FARGO3D magnetohydrodynamics code and demonstrate its performance and second-order accuracy through a series of test problems.