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

TL;DR

This paper introduces a fast, second-order accurate Poisson solver for calculating gravitational potentials in 3D Cartesian and cylindrical coordinates, suitable for isolated systems with vacuum boundary conditions.

Contribution

It develops a novel eigenfunction expansion and boundary screening method, including a discrete Green's function in cylindrical coordinates, integrated into Athena++ for improved accuracy and efficiency.

Findings

The solver achieves second-order accuracy in tests.

It demonstrates good parallel performance.

It effectively computes gravitational potential for isolated systems.

Abstract

We present an accurate and efficient method to calculate the gravitational potential of an isolated system in three-dimensional 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 {\tt Athena++} magnetohydrodynamics code, and perform various tests to check that our solver is second-order accurate and exhibits good parallel performance.

Tables

Table 1. Table 1: Parameters for the uniform sphere test

Coordinate System	Domain Size	Resolution	$x_0$ or $R_0$	$y_0$ or ${\varphi }_0$	$z_0$
Cartesian	$[-0.5,0.5]\times [-0.5,0.5]\times [-0.5,0.5]$	$64\times 64\times 64$	$0.25$	$0.1$	$-0.04$
uniform cylindrical	$[0.5,1]\times [0,2\pi ]\times [-0.25,0.25]$	$64\times 256\times 64$	$0.72$	$0.63$	$-0.04$
logarithmic cylindrical	$[{10}^{-2},1]\times [0,2\pi ]\times [-0.25,0.25]$	$128\times 64\times 64$	$0.27$	$0.38$	$-0.04$

Table 2. Table 2: Parameters for the convergence test

Coordinate system	$x_1$ or $R_1$	$x_2$ or $R_2$	$y_1$ or ${\varphi }_1$	$y_2$ or ${\varphi }_2$	$z_1$	$z_2$
Cartesian	0.0625	0.4375	-0.375	0	-0.125	0.25
uniform cylindrical	0.625	0.90625	0	1.570796327	-0.0625	0.1875
logarithmic cylindrical	0.1	0.7498942093	0	1.570796327	-0.0625	0.1875

Equations

$${\nabla }^2\Phi =4\pi G\rho ,$$

$$\Phi (\mathbf{x})=\iiint {\mathcal{G}}_{\infty }(\mathbf{x},{\mathbf{x}}^{\prime })\rho ({\mathbf{x}}^{\prime })d^3{x}^{\prime },$$

$$\left({\Delta }_x^2+{\Delta }_y^2+{\Delta }_z^2\right){\Phi }_{i,j,k}=4\pi G{\rho }_{i,j,k},$$

$${\Delta }_x^2{\Phi }_{i,j,k}$$

$${\Delta }_y^2{\Phi }_{i,j,k}$$

$${\Delta }_z^2{\Phi }_{i,j,k}$$

$$\left({\Delta }_R^2+{\Delta }_{\varphi }^2+{\Delta }_z^2\right){\Phi }_{i,j,k}=4\pi G{\rho }_{i,j,k},$$

$${\Delta }_R^2{\Phi }_{i,j,k}=$$

$${\Delta }_{\varphi }^2{\Phi }_{i,j,k}=$$

$$R_i\equiv \frac{{\int }_{R_{i-1/2}}^{R_{i+1/2}}{R}^{2}dR}{{\int }_{R_{i-1/2}}^{R_{i+1/2}}RdR}=\frac{2(f^2+f+1)}{3(f+1)}{f}^{i-1}{R}_{\min },$$

$${\Delta }_R^2{\Phi }_{i,j,k}=\frac{{\Phi }_{i-1,j,k}-2{\Phi }_{i,j,k}+{\Phi }_{i+1,j,k}}{(R_i\ln f{)}^2}.$$

$${\mathcal{X}}_i^l=\sin \left(\frac{\pi li}{N_x+1}\right),$$

$${\mathcal{Y}}_j^m=\sin \left(\frac{\pi mj}{N_y+1}\right),$$

$${\mathcal{Z}}_k^n=\sin \left(\frac{\pi nk}{N_z+1}\right),$$

$$\lambda_{x}^{l}=-k_{l}^{2}\left[\sin\left(\frac{\pi l}{2(N_{x}+1)}\right)\bigg{/}\left(\frac{\pi l}{2N_{x}}\right)\right]^{2},$$

$$\lambda_{y}^{m}=-k_{m}^{2}\left[\sin\left(\frac{\pi m}{2(N_{y}+1)}\right)\bigg{/}\left(\frac{\pi m}{2N_{y}}\right)\right]^{2},$$

$$\lambda_{z}^{n}=-k_{n}^{2}\left[\sin\left(\frac{\pi n}{2(N_{z}+1)}\right)\bigg{/}\left(\frac{\pi n}{2N_{z}}\right)\right]^{2},$$

$$\frac{2}{N_x+1}\sum _{i=1}^{N_x}{\mathcal{X}}_i^l{\mathcal{X}}_i^{l^{\prime }}={\delta }_{l{l}^{\prime }}\text{and}\frac{2}{N_x+1}\sum _{l=1}^{N_x}{\mathcal{X}}_i^l{\mathcal{X}}_{i^{\prime }}^l={\delta }_{i{i}^{\prime }},$$

$${\tilde{\Phi }}_{i,j,k}=\frac{8}{(N_x+1)(N_y+1)(N_z+1)}\sum _{l=1}^{N_x}\sum _{m=1}^{N_y}\sum _{n=1}^{N_z}{\tilde{\Phi }}^{lmn}{\mathcal{X}}_i^l{\mathcal{Y}}_j^m{\mathcal{Z}}_k^n,$$

$${\rho }_{i,j,k}=\frac{8}{(N_x+1)(N_y+1)(N_z+1)}\sum _{l=1}^{N_x}\sum _{m=1}^{N_y}\sum _{n=1}^{N_z}{\rho }^{lmn}{\mathcal{X}}_i^l{\mathcal{Y}}_j^m{\mathcal{Z}}_k^n,$$

$${\tilde{\Phi }}^{lmn}$$

$${\rho }^{lmn}$$

$${\tilde{\Phi }}^{lmn}=\frac{4\pi G{\rho }^{lmn}}{{\lambda }_x^l+{\lambda }_y^m+{\lambda }_z^n}.$$

$${\mathcal{P}}_j^m=\exp \left[\frac{2\pi \sqrt{-1}mj}{N_{\varphi }}\right],$$

$$\lambda^{m}_{\phi}=-\frac{m^{2}}{R_{i}^{2}}\left[\sin\left(\frac{\pi m}{N_{\phi}}\right)\bigg{/}\left(\frac{\pi m}{N_{\phi}}\right)\right]^{2}.$$

$$\frac{1}{N_{\varphi }}\sum _{j=1}^{N_{\varphi }}({\mathcal{P}}_j^m{)}^{\ast }{\mathcal{P}}_j^{m^{\prime }}={\delta }_{m{m}^{\prime }}\text{and}\frac{1}{N_{\varphi }}\sum _{m=1}^{N_{\varphi }}({\mathcal{P}}_j^m{)}^{\ast }{\mathcal{P}}_j^{m^{\prime }}={\delta }_{j{j}^{\prime }}.$$

$${\tilde{\Phi }}_{i,j,k}=\frac{2}{N_{\varphi }(N_z+1)}\sum _{m=1}^{N_{\varphi }}\sum _{n=1}^{N_z}{\tilde{\Phi }}_i^{mn}{\mathcal{P}}_j^m{\mathcal{Z}}_k^n,$$

$${\rho }_{i,j,k}=\frac{2}{N_{\varphi }(N_z+1)}\sum _{m=1}^{N_{\varphi }}\sum _{n=1}^{N_z}{\rho }_i^{mn}{\mathcal{P}}_j^m{\mathcal{Z}}_k^n,$$

$${\tilde{\Phi }}_i^{mn}$$

$${\rho }_i^{mn}$$

$$\left({\Delta }_R^2+{\lambda }_{\varphi }^m+{\lambda }_z^n\right){\tilde{\Phi }}_i^{mn}=4\pi G{\rho }_i^{mn},$$

$$\left[\frac{1}{(\delta R{)}^2}-\frac{1}{2R_i\delta R}\right]{\tilde{\Phi }}_{i-1}^{mn}+\left[{\lambda }_{\varphi }^m+{\lambda }_z^n-\frac{2}{(\delta R{)}^2}\right]{\tilde{\Phi }}_i^{mn}+\left[\frac{1}{(\delta R{)}^2}+\frac{1}{2R_i\delta R}\right]{\tilde{\Phi }}_{i+1}^{mn}=4\pi G{\rho }_i^{mn}$$

$$\frac{1}{(R_i\ln f{)}^2}{\tilde{\Phi }}_{i-1}^{mn}+\left[{\lambda }_{\varphi }^m+{\lambda }_z^n-\frac{2}{(R_i\ln f{)}^2}\right]{\tilde{\Phi }}_i^{mn}+\frac{1}{(R_i\ln f{)}^2}{\tilde{\Phi }}_{i+1}^{mn}=4\pi G{\rho }_i^{mn}$$

$${\sigma }_{i,j}(\mathrm{bot})$$

Full text

hydrodynamics — magnetohydrodynamics (MHD) — methods: numerical

††software: Athena++ (J. M. Stone et al. 2019, in preparation, http://princetonuniversity.github.io/athena/), IPython (Perez & Granger, 2007), FFTW (www.fftw.org), Plimpton’s parallel transpose (https://www.sandia.gov/~sjplimp/docs/fft/README.html)

1 INTRODUCTION

There are a number of astronomical systems, such as galactic and protostellar disks, where self-gravity and rotation play an essential role in dynamical evolution. For instance, starburst activity occurring in massive circumnuclear disks can not only inflate the natal disks to form thick tori surrounding active galactic nuclei (AGN) (Wada & Norman, 2002; Wada et al., 2009) but also drive large-scale galactic winds and outflows (Strickland & Stevens, 2000; Strickland et al., 2004a, b; Schneider & Robertson, 2018; Schneider et al., 2018). Accretion disks around AGN may be gravitationally unstable at some radii to form stars (Goodman, 2003; Goodman & Tan, 2004; Levin, 2007; Nayakshin et al., 2007; Jiang & Goodman, 2011). Self-gravity is also important in formation of large-scale spiral structure (Goldreich & Lynden-Bell, 1965; Baba et al., 2013; D’Onghia et al., 2013) and giant molecular clouds (Kim et al., 2003; Dobbs, 2008; Tasker & Tan, 2009) on larger scales in galactic disks. In addition, recent observations of young stellar objects indicate that at least in the early stage of evolution, protostellar disks are massive enough to be self-gravitating (Kratter & Lodato, 2016; Tobin et al., 2016). Gravitational instability of such disks may form trailing spirals that can redistribute the mass and angular momentum and induce heat via shocks (Mejía et al., 2005; Evans et al., 2015) and may be responsible for the formation of giant planets (Boss, 2007; Zhu et al., 2012).

To follow evolution of self-gravitating disks, one needs to solve the Poisson equation

[TABLE]

in cylindrical coordinates $(R,\phi,z)$ subject to a proper boundary condition. In Equation (1), $\Phi$, $\rho$, and $G$ refer to the gravitational potential, mass density, and gravitational constant, respectively. For an isolated system, $\Phi$ has to satisfy vacuum (or “open”) boundary conditions (i.e., $\Phi$ vanishes at infinite distances), for which the formal solution of Equation (1) is given by

[TABLE]

where ${\cal G}_{\infty}({\bf x,x^{\prime}})\equiv-G/|{\bf x-x^{\prime}}|$ is the gravitational potential per unit mass due to a point source situated at ${\bf x^{\prime}}$. Hereafter, we call ${\cal G}_{\infty}({\bf x,x^{\prime}})$ the continuous Green’s function (CGF) to distinguish it from the discrete Green’s function (DGF) based on the discrete Laplace operator (e.g., Burkhart 1997) discussed in Section 4.3.

In simulating dynamics of geometrically thin disks, it has been customary to assume that the disk density along the vertical direction follows a simple function such as Dirac’s delta function for a razor-thin disk and a Gaussian function for a slightly extended disk (e.g., Kalnajs 1971; Miller 1976; Li et al. 2009; Wang et al. 2015). In this case, the integral along the $z$-direction in Equation (2) can be performed analytically, and finding $\Phi(R,\phi)$ at $z=0$ reduces to numerical evaluation of the remaining two-dimensional (2D) integral in the $R$–$\phi$ plane. For example, Miller (1976) solved the gravitational potential of an infinitesimally-thin disk by using a fast Fourier transform (FFT) technique along the azimuthal direction, while directly summing the individual contributions from concentric rings. He introduced a constant softening factor in order to avoid singularity at $\bf x=x^{\prime}$ of the CGF. Li et al. (2009) applied this method to develop an efficient gravity solver for disks with finite thickness on a 2D uniform polar grid. They improved the parallel efficiency of their method by cutting off the high azimuthal Fourier modes based on the energy criterion.

When the grid spacing is logarithmic in the radial direction, a suitable change of variables recasts the integral in Equation (2) to a 2D convolution (Kalnajs, 1971; Binney & Tremaine, 2008), for which the standard FFT convolution method works efficiently (Hockney & Eastwood, 1988). For example, Baruteau & Masset (2008) applied this technique to a razor-thin disk by taking a softening factor proportional to $R$ to avoid divergence of the CGF. Shetty & Ostriker (2008) extended this method to a slightly vertically-extended disk, in which finite disk thickness naturally provides the required softening. Noting that softening reduces the accuracy of a gravity solver, Wang et al. (2015) avoided singularity by using the force kernels integrated over cells, and achieved a second-order accuracy for self-gravity of a razor-thin disk.

Although the methods described above are useful and efficient, they are all limited to 2D polar geometry in the $R$–$\phi$ plane. To our knowledge, there is no efficient method available for fully three-dimensional (3D) cylindrical systems with the vacuum (open) boundary conditions. This is presumably because the Green’s function integral takes a convolution form only along the azimuthal and vertical directions: there is no variable transformation that can cast the integral to a full 3D convolution. One may still attempt to perform the radial integral in Equation (2) by direct summation, while applying the FFT convolution along the azimuthal and vertical directions. But, the associated computational cost is of order ${\cal O}(N^{4}+N^{3}\log N)$, with $N$ being the typical number of cells in one spatial dimension (Pfenniger & Friedli, 1993; Sellwood & Valluri, 1997), making the method computationally prohibitive.

In many cases, it is computationally more efficient to solve Equation (1) directly, rather than evaluating the integral in Equation (2). For example, Gupta et al. (1997) discretized Equation (1) using a fourth-order scheme in 2D Cartesian coordinates and employed a V-cycle multigrid method to solve the resulting linear system. Lai & Tseng (2007) adopted another fourth-order formula to discretize Equation (1) in cylindrical coordinates and solved the resulting linear system using FFT combined with a varient of the Bi-Conjugate Gradient iterative method. Perhaps, the most efficient and robust method to solve a discretized Poisson equation may be a full multigrid algorithm (e.g., Matsumoto & Hanawa, 2003), which can in principle be implemented in either Cartesian or cylindrical coordinates.

However, all the methods mentioned above in turn require provision of appropriate potentials at the domain boundaries in advance. Since Equation (2) naturally satisfies vacuum boundary conditions, it is reasonable to use it to find the desired boundary potentials for Equation (1). Still, the computational cost of ${\cal O}(N^{4}+N^{3}\log N)$ would be inevitable if radial summation is employed for the boundary potentials. One way to reduce the computational cost is to expand the Green’s function in eigenfunction series and truncate it at some point. For instance, the so-called “multipole expansion method” (Black & Bodenheimer, 1975; Stone & Norman, 1992; Boley & Durisen, 2008; Katz et al., 2016) in spherical polar coordinates costs ${\cal O}(l_{\rm max}m_{\rm max}N^{3})$ operations for the boundary potential calculations, where $l_{\rm max}$ and $m_{\rm max}$ refer to the maximum meridional and azimuthal mode numbers, respectively. Although this method appears feasible for small $l_{\rm max}$ and $m_{\rm max}$, the computational cost would increase to ${\cal O}(l_{\rm max}m_{\rm max}N^{4})$ for a flattened mass distribution with surfaces lying close to domain boundaries. An additional $N$ factor in the computational cost arises from the fact that the interior and exterior multipole moments for such a flattened mass distribution are different at most boundary points (see, e.g., Cohl & Tohline 1999).

Cohl & Tohline (1999) derived an alternative expansion of the Green’s function in cylindrical coordinates, which they termed the compact cylindrical Green’s function (CCGF). Their CCGF method can perfectly resolve a highly flattened mass distribution, effectively using $l_{\rm max}=\infty$. Coupled with FFT, the CCGF method requires ${\cal O}(m_{\rm max}N^{3}+N^{3}\log N)$ operations for the boundary potential calculations, and has been applied to 3D simulations of self-gravitating disks (e.g., Fromang, 2005; Mellon & Li, 2008; Marcello & Tohline, 2012; Motl et al., 2017). When the mass distribution is highly non-axisymmetric and/or dominated by small-scale modes, however, the simulation outcomes depend rather sensitively on the choice of $m_{\rm max}$ (Fromang et al., 2004; Fromang, 2005). In such situations, an accurate force evaluation requires $m_{\rm max}$ comparable to $N$ (Fromang, 2005), and the associated computational cost of ${\cal O}(m_{\rm max}N^{3})$ as well as the memory requirement to store four-dimensional arrays would be overwhelming for high-resolution simulations.

A very powerful method to deal with vacuum boundary conditions is the four-step algorithm developed by James (1977) (see also Magorrian, 2007, for more compact description). In the first step, a preliminary solution $\Psi$ is computed based on the interior density distribution and subject to the $\Psi=0$ boundary condition. In the second step, $\Psi$ is used to compute a screening surface charge, $\sigma=\Delta^{2}\Psi/(4\pi G)$, and in the third step a boundary potential $\Theta$ generated by the surface charge and consistent with vacuum boundary conditions is computed via convolution with the DGF. In the fourth step, the final interior solution $\Phi$ is computed, making use of $\Theta$ to enforce a new (nonzero) Dirichlet boundary condition.

The key physical principle underlying James’s method is best understood via an electrostatic analogy. Consider the potential $\Phi$ produced by an isolated, charged box with a metal surface. If the box is then grounded, charges would flow in and be distributed over the surface to enforce zero potential. The potential $\Psi$ of the grounded box is then the sum $\Psi=\Phi+\Theta$ of the potential produced by the original internal charge distribution plus the potential generated by the surface charges. On the surface of the box, the desired potential produced by the original interior charge distribution is then given by $\Phi^{\rm B}=-\Theta$, since by definition $\Psi=0$ on the surface. While direct calculation of $\Phi$ in the whole interior via a volume integration with a Green’s function that enforces vacuum BCs would be computationally very expensive (and would also require storing large arrays), calculation of $\Theta$ on the surface requires only integration over the surface. Once the surface potential $\Phi^{\rm B}$ is known, it can be used to compute the interior $\Phi$ via an efficient numerical method.

The DGF and its convolution presented in James (1977) are valid for Cartesian coordinates, but it is straightforward to extend James’s method to cylindrical coordinates. Snytnikov (2011) was the first to adapt James’s method for 3D cylindrical problems, but he used the CGF in place of DGF. Since the difference between CGF and DGF is quite large at small $\bf|x-x^{\prime}|$ (see Appendix B), this forced him to make a computational domain larger than the original volume in order to improve accuracy. We note, however, that the domain cannot be extended arbitrarily across the inner radial boundary $R_{\rm min}$ in cylindrical coordinates, especially when $R_{\rm min}$ is small. The James algorithm with CGF loses accuracy when the ratio of the outer to inner radial boundary is large, which is frequently encountered in various astronomical applications. In such situations, it is desirable to use the DGF rather than the CGF for accurate potential computations.

In this paper, we develop an efficient, accurate, and scalable algorithm to calculate the gravitational potential of an isolated system in cylindrical coordinates. Our algorithm is fully 3D, and considers both uniform and logarithmic cylindrical grids. For completeness, and to connect with James’s original method, we also present the method for Cartesian coordinates. Our method utilizes the James algorithm with the DGF for the boundary potential from screening charges. For the interior solver, for cylindrical grids we employ a hybrid method that incorporates the eigenfunction expansion in azimuthal and vertical directions and a tridiagonal solver in the radial direction. For Cartesian grids our interior solver employs standard Fourier methods. We implement our algorithm in the Athena++ code framework (e.g., White et al. 2016), and parallelize it on a distributed memory platform using the message passing interface (MPI). Using various test problems, we confirm that our method is efficient and second-order accurate.

The remainder of the paper is organized as follows. In Section 2, we introduce our 3D computational domain in Cartesian and cylindrical coordinates and discretize the Poisson equation. In Section 3, we describe the hybrid method we adopt to solve the discrete Poisson equation in the interior of the computational domain for a given Dirichlet boundary condition. In Section 4, we introduce the James algorithm in Cartesian coordinates and its extension to cylindrical coordinates. In Section 5, we present the results of our Poisson solver on various test problems to demonstrate its accuracy and efficiency. In Section 6, we summarize and discuss the present work.

2 DISCRETE POISSON EQUATION

A standard way to solve Equation (1) is a finite-difference method in which the differential operator $\nabla^{2}$ is replaced by the difference operator $\Delta^{2}$, yielding the discrete Poisson equation $\Delta^{2}\Phi=4\pi G\rho$. The definition of $\Delta^{2}$ depends on the coordinates and the desired level of approximation. This section defines our computational domains and the finite-difference representations of $\Delta^{2}$, which are second-order accurate, in uniform Cartesian, uniform cylindrical, and logarithmic cylindrical coordinate systems.

2.1 Uniform Cartesian Grid

In Cartesian coordinates, we discretize the computational domain $[x_{\rm min},x_{\rm max}]\times[y_{\rm min},y_{\rm max}]\times[z_{\rm min},z_{\rm max}]$ with size $L_{x}\times L_{y}\times L_{z}$ uniformly into $N_{x}\times N_{y}\times N_{z}$ cells. We define the face-centered coordinates as $x_{i+1/2}=x_{\rm min}+i\delta x$ where $\delta x\equiv L_{x}/N_{x}$ and similarly for $y_{j+1/2}$ and $z_{k+1/2}$. We also define the cell-centered coordinates as $x_{i}=(x_{i-1/2}+x_{i+1/2})/2$ with index $i$ running from $1$ to $N_{x}$, and similarly for $y_{j}$ with $j=1,2,\cdots,N_{y}$ and $z_{k}$ with $k=1,2,\cdots,N_{z}$.

We denote the cells inside the nominal index range given above as “active cells”, because these are the places where other fluid variables are updated by a hydrodynamics solver. To the active cells, we add one extra layer of “ghost cells” to the boundaries of the computational domain, with their cell-centered coordinates are denoted, for example, by $i=0$ and $i=N_{x}+1$ in the $x$-direction. The boundary conditions for the Poisson equation and other equations of hydrodynamics are provided using these ghost cells. We similarly define the ghost cells in the other coordinate systems described below.

The second-order accurate, finite-difference approximation to Equation (1) can be written as

[TABLE]

where $\Phi_{i,j,k}$ and $\rho_{i,j,k}$ are the cell-centered potential-density pair and the difference operators $\Delta_{x}^{2}$, $\Delta_{y}^{2}$, and $\Delta_{z}^{2}$ are defined by

[TABLE]

2.2 Uniform Cylindrical Grid

In uniform cylindrical coordinates, we discretize the computational domain $[R_{\rm min},R_{\rm max}]\times[\phi_{\rm min},\phi_{\rm max}]\times[z_{\rm min},z_{\rm max}]$ with size $L_{R}\times L_{\phi}\times L_{z}$ uniformly into $N_{R}\times N_{\phi}\times N_{z}$ cells. We require that $L_{\phi}=\phi_{\rm max}-\phi_{\rm min}$ should be an integer fraction of 2$\pi$ to impose periodic boundary condition along the azimuthal direction. We define the face-centered radial and azimuthal coordinates as $R_{i+1/2}=R_{\rm min}+i\delta R$ and $\phi_{j+1/2}=\phi_{\rm min}+j\delta\phi$, where $\delta R\equiv L_{R}/N_{R}$ and $\delta\phi\equiv L_{\phi}/N_{\phi}$. Unlike in Cartesian coordinates, the definition of the radial cell-center is ambiguous in cylindrical coordinates because the geometric center does not coincide with the volumetric center. When the radial grid is uniform, finite difference of quantities defined at the geometric centers can retain second-order accuracy. We thus define the cell-centered coordinates as $R_{i}=(R_{i-1/2}+R_{i+1/2})/2$ with $i=1,2,\cdots,N_{R}$ and $\phi_{j}=(\phi_{j-1/2}+\phi_{j+1/2})/2$ with $j=1,2,\cdots,N_{\phi}$. Discretization in the vertical direction is the same as in the uniform Cartesian coordinates.

The second-order finite-difference approximation to Equation (1) can be written as

[TABLE]

where the difference operators $\Delta_{R}^{2}$ and $\Delta_{\phi}^{2}$ are defined by

[TABLE]

while $\Delta_{z}^{2}$ is defined through Equation (6).

2.3 Logarithmic Cylindrical Grid

In logarithmic cylindrical coordinates, we discretize a cylindrical computational domain in the same way as in the uniform cylindrical grid, but with logarithmic radial spacing. We define the face-centered radial coordinates as $R_{i+1/2}=f^{i}R_{\rm min}$, with a common multiplication factor $f\equiv(R_{\rm max}/R_{\rm min})^{1/N_{R}}>1$. Since the radial zone width, given by $R_{i+1/2}-R_{i-1/2}=(f-1)f^{i-1}R_{\rm min}$ shrinks toward small radii, a logarithmic cylindrical grid is advantageous in resolving the central regions of a disk with high accuracy. We also define the cell-centered radial coordinates using the volumetric centers as

[TABLE]

for $i=1,2,\cdots,N_{R}$. Note that the radial cell spacing $\delta R_{i}\equiv R_{i+1}-R_{i}=(f-1)R_{i}$ increases with $R_{i}$.

The second-order finite-difference approximation to Equation (1) takes the same form as Equation (7), but with the radial difference operator defined as

[TABLE]

Appendix A shows that Equation (10) makes the finite difference approximation second-order accurate.

3 Calculation of Interior Potential for Dirichlet Boundary Conditions

In this section, we provide the general method that we use to obtain the interior potential within the original domain, given Dirichlet boundaries for the potential on the surface. We describe our methods for both Cartesian grids (Section 3.1) and cylindrical grids (Section 3.2). We note that in principle, alternative fast and efficient solvers (such as multigrid) could be employed for computing the interior potential given Dirichlet boundary conditions for the potential.

The interior solver is employed in three different instances in our method. The first instance is our use of the interior solver to obtain the numerical DGF, as described in Appendix B; this is done once at the beginning of any simulation. The other two instances are in the first and fourth step of the James’s method; each time the Poisson solution is required, two calls to the interior solver are made.

In the first step of James’s method, the interior potential from the original density distribution is computed subject to the zero boundary condition. In the fourth step of James’s method, the interior potential must be computed subject to a known boundary potential $\Phi^{\rm B}$; a method to obtain $\Phi^{\rm B}$ will be presented in Section 4. Formally, we define $\Phi^{\rm B}$ as having nonzero value only in a single layer of ghost zones immediately outside the active domain. Then we can write the desired potential as $\Phi=\widetilde{\Phi}+\Phi^{\rm B}$ where the required boundary condition is $\widetilde{\Phi}=0$. From the definition of $\Phi^{\rm B}$, $\Delta^{2}\Phi^{\rm B}$ will be nonzero only in the single layer of active zones adjoining the domain boundaries. We can thus define a modified density distribution $\rho\rightarrow\rho-\Delta^{2}\Phi^{\rm B}/(4\pi G)$ which is the same as the original density distribution everywhere except in the layer just inside the domain boundaries. We then employ this modified density distribution following the procedure of Section 3.1 or Section 3.2 to compute $\widetilde{\Phi}$. Within the interior, where $\Phi^{\rm B}=0$, this solution is then the desired final solution $\Phi$.

Note that in Sections 3.1 and 3.2 below, $\rho_{i,j,k}$ is any arbitrary density distribution on the grid, and in fact represents a different quantity for each of the three instances where we solve for the interior potential.

3.1 Cartesian Grid Solution with Zero Boundary Value

It is conventional to utilize the eigenfunctions of a differential operator in solving an elliptic partial differential equation. The same technique can be applied to the discretized Poisson equation, if the eigenfunctions of the corresponding discrete Laplace operator can be found.

Let ${\cal X}^{l}_{i}$, ${\cal Y}^{m}_{j}$, and ${\cal Z}^{n}_{k}$ be the eigenfunctions of the discrete Laplace operators $\Delta_{x}^{2}$, $\Delta_{y}^{2}$, and $\Delta_{z}^{2}$ satisfying $\Delta_{x}^{2}{\cal X}^{l}_{i}=\lambda_{x}^{l}{\cal X}^{l}_{i}$, $\Delta_{y}^{2}{\cal Y}^{m}_{j}=\lambda_{y}^{m}{\cal Y}^{m}_{j}$, and $\Delta_{z}^{2}{\cal Z}^{n}_{k}=\lambda_{z}^{n}{\cal Z}^{n}_{k}$, with respective eigenvalues $\lambda^{l}_{x}$, $\lambda_{y}^{m}$, and $\lambda_{z}^{n}$. It is straightforward to show that

[TABLE]

[TABLE]

[TABLE]

are the desired eigenfunctions satisfying the zero boundary condition at the ghost cells. The corresponding eigenvalues are

[TABLE]

[TABLE]

[TABLE]

where $k_{l}\equiv\pi l/L_{x}$, $k_{m}\equiv\pi m/L_{y}$, and $k_{n}\equiv\pi n/L_{z}$. In the limit of $l/N_{x},m/N_{y},n/N_{z}\ll 1$, the discrete eigenvalues reduce to the counterpart of the continuous Laplace operator ($-k^{2}$).

The discrete analog of the Sturm-Liouville theory (e.g., Hildebrand 1968, chap. 1.10–1.16; see also, Atkinson 1964) guarantees that the eigenfunctions given in Equation (11)–(13) satisfy discrete orthogonality relations, for example,

[TABLE]

where the symbol $\delta$ denotes the Kronecker delta. These orthogonality relations allow us to expand $\widetilde{\Phi}_{i,j,k}$ and $\rho_{i,j,k}$ as

[TABLE]

where the expansion coefficients $\widetilde{\Phi}^{lmn}$ and $\rho^{lmn}$ are given by

[TABLE]

Plugging Equations (18)–(19) in Equation (3), we obtain a simple algebraic relation

[TABLE]

Therefore, the Poisson equation in Cartesian coordinates can be solved by the following three steps:

Perform a forward transform $\rho_{i,j,k}\to\rho^{lmn}$ using Equation (21) : ${\cal O}(N_{x}N_{y}N_{z}\log_{2}[N_{x}N_{y}N_{z}])$. 2. 2.

Convert $\rho^{lmn}\to\widetilde{\Phi}^{lmn}$ using the kernel in Equation (22) : ${\cal O}(N_{x}N_{y}N_{z})$. 3. 3.

Perform a backward transform $\widetilde{\Phi}^{lmn}\to\widetilde{\Phi}_{i,j,k}$ using Equation (18) : ${\cal O}(N_{x}N_{y}N_{z}\log_{2}[N_{x}N_{y}N_{z}])$.

In practice, the discrete transforms in Equations (18)–(21) can be performed efficiently with an FFT algorithm. The public FFTW library111http://www.fftw.org/ performs sine transforms of various kinds, among which we use FFTW_R0DFT00 consistent with the zero boundary condition. To perform 3D FFTs in parallel, we decompose the computational domain into 2D pencils along, for example, the $y$-axis and execute 1D sine transforms locally in each pencils (e.g., Li & Laizet 2010). We then transpose the pencils to the $z$- and $x$-axes sequentially, each time by executing corresponding sine transforms, which completes a 3D FFT. The parallel transpose among different pencil decompositions are done with the remap_3d function in Steve Plimpton’s parallel FFT package222https://www.sandia.gov/~sjplimp/docs/fft/README.html. We note that, since mass density and the gravitational potential in general are distributed as blocks rather than pencils in real applications, we have to transpose between block and pencil decompositions at the input and output stage of the Poisson solver. Plimpton’s remap routine provides this functionality as well.

3.2 Cylindrical Grid Solution with Zero Boundary Value

A natural boundary condition in the azimuthal direction is that both mass density and gravitational potential are periodic, with period $L_{\phi}$. This holds true even when the problem under study has $P$-fold symmetry in the $\phi$-direction, with a domain size $L_{\phi}=2\pi/P$. The algorithm presented below is applicable for such systems as long as the $P$-fold symmetry is considered in the boundary condition for the DGF (Equation (B3)).

The eigenfunction ${\cal P}^{m}_{j}$ for the discrete Laplace operator $\Delta_{\phi}^{2}$ and the corresponding eigenvalue $\lambda_{\phi}^{m}$ are given by

[TABLE]

[TABLE]

Note that $\lambda^{m}_{\phi}\rightarrow-m^{2}/R^{2}$ for $m/N_{\phi}\ll 1$. The eigenfunction ${\cal P}^{m}_{j}$ satisfies the discrete orthogonality relation

[TABLE]

We expand $\widetilde{\Phi}_{i,j,k}$ and $\rho_{i,j,k}$ only along the azimuthal and vertical directions as

[TABLE]

where the expansion coefficients $\widetilde{\Phi}^{mn}_{i}$ and $\rho^{mn}_{i}$ satisfy the inverse transforms

[TABLE]

One cannot analytically expand $\widetilde{\Phi}_{i,j,k}$ and $\rho_{i,j,k}$ along the radial direction because radial eigenfunction ${\cal R}^{l}_{i}$, defined through $\Delta_{R}^{2}{\cal R}^{l}_{i}=\lambda^{l}_{R}{\cal R}^{l}_{i}$, has no closed-form expression and is not compatible with FFT.333The radial eigenfunction ${\cal R}^{l}_{i}$ can instead be obtained numerically by solving the eigenvalue problem, and the resulting eigenfunction may be called the discrete Bessel function. Since it satisfies the exact discrete orthogonality relation, it may also serve as discrete kernel for the discrete Hankel transform (Johnson, 1987; Baddour & Chouinard, 2015).

Plugging Equation (26)–(27) into Equation (7) yields

[TABLE]

which, using Equations (8) and (10), can be written as

[TABLE]

in uniform cylindrical coordinates, and

[TABLE]

in logarithmic cylindrical coordinates. Note that Equations (31) and (32) are tridiagonal matrix equations subject to the zero boundary conditions of $\widetilde{\Phi}_{0}^{mn}=\widetilde{\Phi}_{N_{R}+1}^{mn}=0$, which can easily be solved via the Thomas algorithm involving back substitutions (e.g., Press et al. 2007).

Therefore, the Poisson equation in cylindrical coordinates can be solved by the following three steps:

Perform a forward transform $\rho_{i,j,k}\to\rho^{mn}_{i}$ using Equation (29) : ${\cal O}(N_{R}N_{\phi}N_{z}\log_{2}[N_{\phi}N_{z}])$. 2. 2.

Solve Equation (31) or (32) for $\widetilde{\Phi}^{mn}_{i}$ : ${\cal O}(N_{R}N_{\phi}N_{z})$. 3. 3.

Perform a backward transform $\widetilde{\Phi}^{mn}_{i}\to\widetilde{\Phi}_{i,j,k}$ using Equation (26) : ${\cal O}(N_{R}N_{\phi}N_{z}\log_{2}[N_{\phi}N_{z}])$.

In actual computation, the discrete transforms in Equations (26)–(29) can be carried out efficiently with an FFT algorithm. For transforms involving ${\cal P}^{m}_{j}$, we use the real-to-complex transform in FFTW, which halves the size of the output by utilizing the Hermitian symmetry. For transforms involving ${\cal Z}^{n}_{k}$, we use the sine transform as in Cartesian coordinates. For parallel computations, we employ the pencil decomposition technique along with the Steve Plimpton’s parallel transpose routines, similarly to the Cartesian solver.

4 CALCULATION OF THE BOUNDARY POTENTIAL

4.1 Overview of the James Algorithm

We adopt the James algorithm to calculate the boundary potential $\Phi^{\rm B}$ which is second-order accurate. As explained in Introduction, James’s method first solves for the preliminary potential $\Psi=\Phi+\Theta$ with zero boundary condition, where $\Phi({\bf x})=-\int G\rho({\bf x^{\prime}})/|{\bf x-x^{\prime}}|\,d^{3}x^{\prime}$ is the gravitational potential generated from the original density distribution (i.e., the desired solution) and $\Theta({\bf x})=-\oint G\sigma({\bf x^{\prime}})/|{\bf x-x^{\prime}}|\,d^{2}x^{\prime}$ is that from the screening charges. Since $\Psi^{\rm B}=0$ by definition, one only needs to compute $\Theta^{\rm B}=-\oint G\sigma({\bf x^{\prime}})/|{\bf x}^{\rm B}-{\bf x^{\prime}}|\,d^{2}x^{\prime}$ to obtain $\Phi^{\rm B}=-\Theta^{\rm B}$.

Once the preliminary gravitational potential $\Psi$ with the zero boundary condition is obtained (using the method of Section 3 and the original density distribution), the screening charges are found by applying the discrete Laplace operators at the ghost cells. It should be noted that at the ghost cells, the discrete Laplace operator calls for the value of $\Psi$ outside the ghost cells, which is set to zero in James (1977). Figure 1 depicts this situation, where ${\mathscr{R}}$ and $\partial{\mathscr{R}}$ represent the active and the ghost cells defined in Section 2. Encompassing this, one can imagine the infinite domain ${\mathscr{R}}_{\infty}$ where $\Psi=\rho=0$ everywhere exterior to $\partial{\mathscr{R}}$. It is evident that the $\Psi$ satisfies the discrete Poisson equation at every cell in ${\mathscr{R}}_{\infty}$ if one adds $\sigma=\Delta^{2}\Psi/(4\pi G)$ at $\partial{\mathscr{R}}$. Since $\Delta^{2}\Psi=4\pi G(\rho+\sigma)$ in $\mathscr{R}_{\infty}$, one can solve $\Delta^{2}\Theta=4\pi G\sigma$ to obtain $\Phi=\Psi-\Theta$ that satisfies the original Poisson equation $\Delta^{2}\Phi=4\pi G\rho$ subject to the vacuum boundary condition (i.e., $\Phi$ is the desired solution for an isolated mass distribution). Note that one needs $\Theta$ only at the domain boundary, which can be efficiently calculated using the Green’s function for the discrete Laplace operator subject to vacuum boundary condition (i.e., the DGF). Then, $\Phi^{\rm B}=-\Theta^{\rm B}$ gives a new boundary condition for the final step of the James’s method.

To summarize, the James algorithm consists of the following four steps: (1) Solve the Poisson equation with the zero boundary condition to obtain $\Psi$; (2) Evaluate the screening charge $\sigma$ by applying the discrete Laplace operator to the ghost cells ($\sigma=\Delta^{2}\Psi/4\pi G$); (3) Use the DGF to calculate the gravitational potential $\Theta^{\rm B}$ at the domain boundary due to $\sigma$; (4) Solve the Poisson equation with the Dirichlet boundary condition $\Phi^{\rm B}=-\Theta^{\rm B}$. In section 3 we have presented the method we adopt for the interior Poisson solver, which is employed in Steps (1) and (4) and also is used to pre-compute the DGF (see Appendix B), which enters in Step 3. In what follows, we describe Steps (2) and (3) in more detail.

4.2 Computation of the Screening Charges

Once the preliminary gravitational potential $\Psi$ with the zero boundary condition is obtained (using the method of Section 3 and the original density distribution), one can readily apply the discrete Laplace operators at the ghost cells in each boundary to calculate the screening charges.

4.2.1 Cartesian Grid

A Cartesian grid has six boundary surfaces consisting of the loci of ghost zones immediately outside the problem domain: bottom (bot; $k=0$), top (top; $k=N_{z}+1$), south (sth; $j=0$), north (nth; $j=N_{y}+1$), west (wst; $i=0$), and east (est; $i=N_{x}+1$). With $\Psi=0$ in both the first and second layer of ghost zones outside the domain, the screening charges on these boundary surfaces are given by

[TABLE]

where $\sigma_{i,j}({\rm bot})$ denotes the screening charge on the bottom boundary, etc. Note that the screening charges have units of mass density rather than surface density, because the charge is assumed to fill a volume $\delta x\delta y\delta z$.

4.2.2 Cylindrical Grid

In cylindrical coordinates, one needs to deal with only four boundary surfaces: bottom (bot; $k=0$), top (top; $k=N_{z}+1$), inner (inn; $i=0$), and outer (out; $i=N_{R}+1$): the azimuthal direction is assumed periodic. Using the discrete Laplace operators given in Section 2, one can show that the screening charges are calculated as

[TABLE]

in a uniform cylindrical grid, and

[TABLE]

in a logarithmic cylindrical grid.

4.3 Discrete Green’s Function and the Potential Generated by Screening Charges

The gravitational potential $\Theta$ that results from the screening charges can be obtained by convolving $\sigma$ with the Green’s function of the operator that determines $\sigma$. Since $\sigma$ is obtained through the application of the discrete Laplace operator, the corresponding Green’s function should be the DGF rather than the CGF. The proper operation of the James’s method thus relies on the accurate calculation of the DGF, yet finding its analytic expressions in cylindrical coordinates is a daunting task. To our knowledge, the analytic DGF is available only in 2D Cartesian coordinates (Buneman, 1971). In 3D Cartesian coordinates, Burkhart (1997) addressed the definition, existence, and uniqueness of the DGF and derived asymptotic expansion formulae, applicable at distances far from the source. In this section, we provide a working definition of the DGF and a numerical method to calculate $\Theta$ in Cartesian and cylindrical coordinates. We refer the reader to Appendix B for our method for the DGF.

4.3.1 Cartesian Grid

The DGF, ${\cal G}_{i-i^{\prime},j-j^{\prime},k-k^{\prime}}$, in Cartesian coordinates is the gravitational potential per unit mass due to a discrete point mass at $(i^{\prime},j^{\prime},k^{\prime})$ and ought to satisfy

[TABLE]

where the symbol $\delta_{ii^{\prime}}$ is the Kronecker delta and ${\cal V}=\int_{z_{k^{\prime}-1/2}}^{z_{k^{\prime}+1/2}}\int_{y_{j^{\prime}-1/2}}^{y_{j^{\prime}+1/2}}\int_{x_{i^{\prime}-1/2}}^{x_{i^{\prime}+1/2}}dxdydz=\delta x\delta y\delta z$ is the volume of the $(i^{\prime},j^{\prime},k^{\prime})$-th cell. Note that in writing the indices of ${\cal G}_{i-i^{\prime},j-j^{\prime},k-k^{\prime}}$, we implicitly allow for the translational symmetry on a Cartesian grid. In Appendix B.1, we follow James (1977) to calculate the Cartesian DGF numerically.

The gravitational potential $\Theta_{i,j,k}$ generated by the screening charges $\sigma_{i,j,k}$ is given by

[TABLE]

Substituting Equation (37) into Equation (3), one can easily check that $\Theta_{i,j,k}$ and $\sigma_{i,j,k}$ is a valid potential-density pair.

Because Equation (37) involves a discrete convolution, it is efficient to use FFTs for computations. Making use of certain symmetries of the problem for a hollow charge distribution, James (1977) devised a formulation that uses sine and cosine transforms to expresses the potential on the each surface as the sum of seven terms. We refer the reader to Equations (4.7)–(4.20) of James (1977) for the description of this formulation, which costs ${\cal O}(N^{3}+N^{2}\log N)$ operations.

The gravitational potential $\Phi^{\rm B}$ at the domain boundary is obtained by

[TABLE]

which provides the required Dirichlet boundary condition for the interior solver (Section 3). Note that $\Psi^{\rm B}_{i,j,k}=0$ by definition.

4.3.2 Cylindrical Grid

The cylindrical DGF, ${\cal G}_{i,i^{\prime},j-j^{\prime},k-k^{\prime}}$, satisfies

[TABLE]

where ${\cal V}_{i^{\prime}}=\int_{z_{k^{\prime}-1/2}}^{z_{k^{\prime}+1/2}}\int_{\phi_{j^{\prime}-1/2}}^{\phi_{j^{\prime}+1/2}}\int_{R_{i^{\prime}-1/2}}^{R_{i^{\prime}+1/2}}R\,dRd\phi dz=\tfrac{1}{2}(R_{i^{\prime}+1/2}^{2}-R_{i^{\prime}-1/2}^{2})\delta\phi\delta z$ is the volume of the $(i^{\prime},j^{\prime},k^{\prime})$-th cell. Note that the cylindrical DGF has four indices due to lack of the translational symmetry along the radial direction. In Appendix B.2, we present the method to calculate the cylindrical DGF and compare it with the continuous counterpart.

The gravitational potential $\Theta_{i,j,k}$ generated by the screening charges $\sigma_{i,j,k}$ takes a form of

[TABLE]

One can readily verify that $\Theta_{i,j,k}$ satisfies the discrete Poisson equation in cylindrical coordinates. Since all functions are periodic in the azimuthal direction, it is natural to apply a discrete Fourier transform such that

[TABLE]

and similarly for ${\cal G}$ and $\sigma$. Then, Equation (40) can be cast into a more compact form

[TABLE]

Since $\sigma$ is nonzero only at the ghost cells and we need to evaluate $\Theta$ also only at the ghost cells, we do not have to perform full double summations for all $m,i,k$ indices in Equation (42). By collecting the individual contributions of the surface charges, Equation (42) yields the Fourier-transformed potentials at the four boundaries

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where we use symbolic notations such that ${\cal G}^{m}_{i,k^{\prime}}({\rm inn\to top})={\cal G}^{m}_{i,0,N_{z}+1-k^{\prime}}$, $\sigma^{m}_{i^{\prime}}({\rm top})=\sigma^{m}_{i^{\prime},N_{z}+1}$, etc.

Finally, we apply an inverse Fourier transform to $\Theta^{m}_{i}({\rm top}),\Theta^{m}_{i}({\rm bot}),\Theta^{m}_{k}({\rm inn})$, and $\Theta^{m}_{k}({\rm out})$ to obtain the boundary potential $\Theta_{i,j,k}^{\rm B}$ due to the surface charges. Then, the desired boundary potential $\Phi^{\rm B}$ due to the original charge $\rho$ is given by

[TABLE]

which gives the required Dirichlet boundary condition for the interior solver. Note that the boundary potential calculations explained above involve FFTs on 2D arrays (e.g., $\sigma^{m}_{i^{\prime}}({\rm top})$) together with the summations of the Green’s function amounting to ${\cal O}(N^{3})$ operations. Therefore, the overall computational cost of the boundary potential calculation is of order ${\cal O}(N^{3}+N^{2}\log N)$, similarly to the case with a Cartesian grid.

We note that the above formulation for the boundary potential is valid for a mass distribution under $P$-fold symmetry in $\phi$. In Appendix C, we directly demonstrate that this is really the case as long as ${\cal G}_{i,i^{\prime},j-j^{\prime},k-k^{\prime}}$ in Equation (40) properly accounts for the contributions to the boundary potential from all periodic images of the mass density.

5 TEST RESULTS

We implement our Poisson solver in Athena++ which is a state-of-art astrophysical magnetohydrodynamics (MHD) code with very flexible coordinate and grid options. Using Cartesian and uniform/logarithmic cylindrical grids, we test our solver on a few test problems to check its accuracy, convergence, and parallel performance. We also run time-dependent simulations of a gravitationally-unstable isothermal ring to check if the gravity module combines well with the MHD solver of Athena++ to produce the expected results of ring fragmentation. For all tests presented below, we set the gravitational constant to $G=1$.

5.1 Uniform Sphere Test

To test the accuracy of our Poisson solver, we consider a uniform sphere with radius $r_{0}$ and density $\rho_{0}$. The analytic gravitational potential of such a sphere is given by

[TABLE]

where $r$ denotes the distance from the center of the sphere. We take $r_{0}=0.2$ and $\rho_{0}=1$, and place it at an off-centered position $(x_{0},y_{0},z_{0})$ in Cartesian coordinates and $(R_{0},\phi_{0},z_{0})$ in cylindrical coordinates, and calculate the gravitational potential $\Phi$ numerically. Table 1 lists the grid dimension, resolution, and sphere position in each coordinate system adopted.

As a measure of accuracy, we define the relative error between the numerical solution and the analytic solution as

[TABLE]

evaluated at the cell centers. Figures 2–4 plot the test results on the Cartesian, uniform cylindrical, and logarithmic cylindrical grid, respectively. Panels (a) and (b) plot the one-dimensional (1D) cut profiles of $\Phi$ and $\epsilon$ along the $x$- or $R$-direction, respectively, while panel (c) gives the 2D distribution of $\epsilon$ in the $z=0$ plane. Overall, the numerical results are exceedingly close to the analytic potential, with the mean relative error less than 0.1%. The errors are largest near the sphere boundary whose exact shape is not well resolved by any of the adopted grids.

5.2 Convergence Test

The discrete Poisson equation used for the interior solver in Section 3 and for the boundary condition in Section 4 are second-order accurate by construction. If our implementation of the Poisson solver is correct, therefore, the relative errors should be inversely proportional to the square of the grid spacing. To check if this is indeed the case, we repeat the uniform sphere tests by varying the number of cells from $16^{3}$ to $512^{3}$. Figure 5 plots as circles the mean relative errors $\left\langle{\epsilon}\right\rangle$ from the sphere tests as functions of $N_{z}$. Overall, the errors decrease roughly at a second-order rate with increasing $N_{z}$, but exhibit some fluctuations. Katz et al. (2016) noted that these fluctuations of the errors are caused not by the truncation errors of the finite-difference scheme but by inability of an adopted grid to perfectly resolve a spherical mass distribution. This is true even for a spherical grid when the sphere center offsets from the origin.

To delineate the truncation errors alone, it is thus necessary to design a solid figure whose shape is identical to the cell shape of an adopted grid. In addition, the size and mass of the solid figure should be unchanged with varying resolution. For this purpose, we consider a uniform cube with density $\rho=1$, located at $x_{1}\leq x\leq x_{2}$, $y_{1}\leq y\leq y_{2}$, and $z_{1}\leq z\leq z_{2}$ in a Cartesian grid. In a cylindrical grid, we consider a rectangular torus with density $\rho=1$, occupying the regions with $R_{1}\leq R\leq R_{2}$, $\phi_{1}\leq\phi\leq\phi_{2}$, and $z_{1}\leq z\leq z_{2}$. Table 2 lists the parameters of the solid figures that we adopt: these values ensure that the mass distribution does not change with resolution.

We use our Poisson solver to calculate the gravitational potentials of the solid figures by varying resolution from $N_{z}=16$ to $512$, while keeping the domain sizes the same as in Section 5.1. As the reference potential, we take Equation (20) of Katz et al. (2016) for the gravitational potential of a uniform cube. There is no algebraic expression for the potential of a rectangular torus, but Huré et al. (2014) provided a closed-form expression, in terms of line integrals with smooth integrands, in their Equation (29). We use the Romberg’s method with a relative tolerance $10^{-10}$ to ensure that the numerical integrations are accurate enough to serve as a reference solution. The squares in Figure 5 plot the resulting mean relative errors $\left<{\epsilon}\right>$ for the cube and rectangular torus. Note that $\left<{\epsilon}\right>$ against $N_{z}$ follows almost a straight line with slope of $-1.993$, $-2.003$, and $-2.003$ in the Cartesian, uniform cylindrical, and logarithmic cylindrical grid, respectively, confirming that our implementation of the Poisson solver retains a second-order accuracy.

5.3 Performance Test

To check the parallel performance of our implementation, we conduct a weak scaling test of our Poisson solver on the TigerCPU linux cluster at Princeton University444https://researchcomputing.princeton.edu/systems-and-services/available-systems/tiger. The TigerCPU cluster consists of 408 nodes, with each node comprised of 40 $2.4\,{\rm GHz}$ intel Skylake processors. We measure the wall clock time using the MPI_Wtime function, after a call to MPI_Barrier to synchronize all processors. We run the job with SLURM --exclusive option to make sure that other jobs do not interfere with ours.

For the weak scaling test, we divide the whole computational domain into $N_{\rm core}$ subdomains consisting of $64^{3}$ cells each. In the Athena++ terminology, the whole domain and subdomain are referred to as Mesh and MeshBlock, respectively. We assign each MeshBlock to a single processor, while varying $N_{\rm core}$ from $1$ to $4096$. The corresponding size of Mesh varies from $64^{3}$ to $1024^{3}$. In each run, we call our Poisson solver, together with the MHD solver for comparison, 100 times, and then measure the wall clock time per cycle $t_{\rm wall}$ for various steps. We repeat the calculations five times in order to avoid unusual runs due to stale nodes.

Figure 6 plots the mean values $\left\langle t_{\rm wall}\right\rangle$ of the wall clock times per cycle as functions of $N_{\rm core}$ for the Cartesian (left) and cylindrical (right) grids. As noted earlier, our Poisson solver requires to run the interior solver twice and the boundary solver once. The total time taken by the Poisson solver (triangles) is dominated by the interior solver (squares) rather than the boundary solver (circles). While the time taken by the MHD solver is comparable between Cartesian and cylindrical grids, the Poisson solver is more efficient in the cylindrical grid. This is because Cartesian coordinates has no inherent periodic direction and thus requires more operations to implement the open boundary conditions. In addition, the Cartesian grid has two more boundaries than the cylindrical grid and thus needs more boundary-to-boundary interactions in the boundary solver. Notwithstanding these differences, the Poisson solver in both Cartesian and cylindrical grids takes less time than the MHD solver (stars) at least up to 4096 processors, leading us to conclude that our Poisson solver is very efficient and does not contribute much to the total computational cost of self-gravitating MHD simulations.555The weak scaling test shown in Figure 6 hints some performance degradation from $N_{\rm core}=1$ to $64$ relative to $\langle t_{\rm wall}\rangle\propto\ln N_{\rm core}$ expected for the theoretical FFT. Our parallel FFT utilizes a “transpose algorithm” known to be efficient when a data size for communication is larger than the critical size that depends on the latency/bandwidth of the interconnecting network device and the network topology (e.g., Foster & Worley, 1997). An alternative “binary exchange algorithm” may work efficiently for a small data size (e.g., Müller & Chan, 2019).

5.4 Ring Fragmentation

As a final test of our Poisson solver, we run time-dependent simulations for fragmentation of a self-gravitating isothermal ring using Athena++. As an initial condition, we consider a rigidly-rotating ring at angular velocity $\Omega_{0}$ and sound speed $c_{s}$, surrounded by a tenuous hot external medium. The ring is initially in hydrosteady equilibrium under both self-gravity and external gravity. Assuming that the external gravity ${\bf g}_{\rm ext}=-\Omega_{e}^{2}{\bf R}$ alone makes the ring rotate at angular frequency $\Omega_{e}$, the equation for such an equilibrium reads

[TABLE]

together with Equation (1), where $\Omega_{s}\equiv(\Omega_{0}^{2}-\Omega_{e}^{2})^{1/2}$ is the angular velocity due to self-gravity alone (Kim & Moon, 2016). One can show that the equilibrium configurations are completely specified by two dimensionless parameters: $\alpha\equiv c_{s}^{2}/(GR_{A}^{2}\rho_{c})$ and $\widehat{\Omega}_{s}\equiv\Omega_{s}/(G\rho_{c})^{1/2}$, where $\rho_{c}$ and $R_{A}$ denote the maximum density and the maximum radial extent of an equilibrium object, respectively.

Using the self-consistent field method of Hachisu (1986), we solve Equations (1) and (50) alternatively and iteratively to find the equilibrium configuration for $\alpha=0.015$ and $\widehat{\Omega}_{s}=0.22$. We then boost the angular velocity of the ring to $\widehat{\Omega}_{0}=0.3$ to account for ${\bf g}_{\rm ext}$. The external medium is set to be hotter than the ring by two orders of magnitude. In order to check whether the configuration we set is really in equilibrium, we evolve it over time on a logarithmic cylindrical grid with size $\widehat{R}\equiv R/R_{A}\in[0.75,1.05]$, $\phi\in[0,2\pi]$, and $z/R_{A}\in[-0.15,0.15]$, without imposing any perturbations. The number of cells used is $N_{R}=64$, $N_{\phi}=1024$, and $N_{z}=64$. Figure 7 compares the radial profiles of the dimensionless density $\widehat{\rho}\equiv\rho/\rho_{c}$ and pressure $\widehat{p}=\alpha\widehat{\rho}$ in the equatorial plane as functions of $\widehat{R}$, at $t/T_{\rm orb}=0$, 1, and 5, where $T_{\rm orb}=2\pi/\Omega_{0}$. The density and pressure profiles are almost unchanged over time except near the contact discontinuity between the ring and external medium, demonstrating that the initial steady configuration is well maintained over many orbital periods.

By performing a linear-stability analysis, Kim & Moon (2016) found that the above equilibrium is gravitationally unstable to non-axisymmetric perturbations for a range of the azimuthal mode number $m$. The most unstable modes were found to have $m=9$ and $10$, with an almost equal growth rate ${\rm Im}(\omega)\approx 0.8(G\rho_{c})^{1/2}\approx 17\,T_{\rm orb}^{-1}$ and a phase speed ${\rm Re}(\omega)/m\approx 0.3(G\rho_{c})^{1/2}$, corresponding to overstability. To check if our Poisson solver can pick up these unstable modes and capture their growth, we generate random density perturbations with amplitude $10^{-5}$ and apply them to the equilibrium density $\rho$. We then evolve the system and monitor how various modes grow. Figure 8 plots the azimuthal profiles of the dimensionless density at $\widehat{R}=0.91$ and $\widehat{z}=0$ at a few selected epochs as well as the surface density $\Sigma=\int\rho\,dz$ at $t(G\rho_{c})^{1/2}=20.2$. It is apparent that the perturbations grow as they propagate along the $\phi$-direction. In the highly nonlinear stage, the density distribution is dominated equally by the $m=9$ and $m=10$ modes, indicating that these are two fastest growing modes of the instability.

To measure the growth rate and the phase speed of each mode in the numerical simulation, we calculate the Fourier transform ${\cal L}_{m}e^{-i\vartheta_{m}}\equiv\int\bar{\rho}(\phi)e^{-im\phi}\,d\phi$ of the radially- and vertically-integrated density $\bar{\rho}(\phi)=\iint\rho\,dRdz$. Figure 9 plots the temporal variations of the Fourier amplitude ${\cal L}_{m}$ and the phase angle $\vartheta_{m}$ for $m=9$ and $m=10$. The two modes grow exponentially as they propagate at a constant phase speed. The dimensionless growth rate $d\ln{\cal L}_{m}/dt/(G\rho_{c})^{1/2}$ and the phase speed $d{\vartheta_{m}}/dt/[m(G\rho_{c})^{1/2}]$ are measured to be $0.80$ and $0.30$, respectively, consistent with the results of the linear stability analysis.

Taken together, all the test results presented in this section demonstrate that our 3D cylindrical Poisson solver is reliable, accurate to second order, and very efficient.

6 SUMMARY AND DISCUSSION

To study dynamical evolution of self-gravitating, rotating disks, it is desirable to use a fully 3D Poisson solver in a cylindrical geometry subject to the open boundary condition. In this paper, we have presented an accurate and efficient algorithm for such a Poisson solver that works in Cartesian, uniform cylindrical, and logarithmic cylindrical coordinates. Our algorithm adopts “surface screening charge” method introduced by James, employing the DGF to calculate the boundary potential consistent with the open boundary condition (Section 4), and utilizes the eigenfunction expansion method to solve the interior potential (Section 3). The computational cost of our algorithm is of order ${\cal O}(N^{3}+N^{2}\log N)$, with $N$ being the number of cells in one dimension. The results of the various test problems presented in Section 5 confirm that our Poisson solver is second-order accurate and takes less computational cost than the MHD solver in Athena++ up to 4096 cores.

The second-order accuracy of our Poisson solver is enabled by the usage of the cylindrical DGF. Although Serafini et al. (2005) and Snytnikov (2011) implemented the James algorithm in Cartesian and cylindrical coordinates, respectively, they simply used the CGF instead of the DGF after enlarging the computational domain. Since the difference between CGF and DGF is quite small at large distances from the source (see Appendix B), the potential based on the CGF would be similar to that with the DGF if the domain is expanded sufficiently. As mentioned in Introduction, however, domain expansion in cylindrical coordinates is very limited toward the inner radial boundary, so that the potential with the CGF would become less accurate for larger $R_{\rm max}/R_{\rm min}$, where $R_{\rm max}$ and $R_{\rm min}$ denote the outer and inner radial boundaries of the domain, respectively. To demonstrate this, we recalculate the gravitational potential of a rectangular torus on a logarithmic cylindrical grid presented in Section 5.2, but this time by using the CGF. For fare comparison, we apply domain expansion technique similar to Snytnikov (2011) when using CGF. Specifically, we use the same extended domain for DGF calculation defined in Appendix B.2 for the enlarged domain. Figures 10 plots the resulting mean relative errors as red triangles against $R_{\rm max}/R_{\rm min}$, in comparison with the cases based on the DGF plotted as blue circles. The CGF works as well as the DGF for $R_{\rm max}/R_{\rm min}\lesssim 10^{2}$. But, it fails to give second-order convergence for systems with $R_{\rm max}/R_{\rm min}\gtrsim 10^{2}$, which are common in astronomical applications (e.g., Kuiper et al., 2010; Seo & Kim, 2013; Zhu et al., 2012; Bae et al., 2014; Ju et al., 2016; Kim & Elmegreen, 2017). In these cases, it is necessary to use the cylindrical DGF for accurate potential calculations.

SM wishes to thank Chang-Goo Kim, Jeong-Gyu Kim and Kengo Tomida for their helpful discussions and advice. The work of SM was supported by NRF (National Research Foundation of Korea) Grant funded by the Korean Government (NRF-2017-Fostering Core Leaders of the Future Basic Science Program/Global Ph.D. Fellowship Program). The work of WTK was supported by the grant (2017R1A4A1015178) of National Research Foundation of Korea. The work of ECO on this project is supported by grant 510940 from the Simons Foundation. The computation of this work was supported by the Supercomputing Center/Korea Institute of Science and Technology Information with supercomputing resources including technical support (KSC-2018-C3-0015) and the PICSciE TIGRESS High Performance Computing Center at Princeton University.

Appendix A SECOND-ORDER FINITE-DIFFERENCE IN LOGARITHMIC CYLINDRICAL COORDINATES

Here we derive a second-order finite-difference approximation to the radial part of the Laplace operator in logarithmic cylindrical coordinates. With the change of the variables $u\equiv\ln R$, the radial part of the Laplacian becomes

[TABLE]

Since the logarithmic grid in $R$ ($R_{i}=R_{0}f^{i}$) corresponds to a uniform grid in $u$ ($u_{i}=u_{0}+i\ln f$), we can apply a centered difference scheme in the $u$-space to obtain

[TABLE]

Noting that $\delta u=\ln f=N_{R}^{-1}\ln(R_{\rm max}/R_{\rm min})$, the above expression can be expressed in the $R$-space as

[TABLE]

It is evident that the remainder decreases at a second-order rate with increasing $N_{R}$. We adopt Equation (A3) as our discrete Laplace operator in the logarithmic cylindrical grid (see Equation (10)).

Appendix B COMPUTATION OF THE DISCRETE GREEN’S FUNCTION

The DGF is needed in order to calculate the surface potential associated with the surface screening charges. The DGF is pre-computed once at the beginning of any simulation. In this Appendix we outline the numerical method used to evaluate the DGF.

B.1 Cartesian Grid

We calculate the DGF ${\cal G}_{i-i^{\prime},j-j^{\prime},k-k^{\prime}}$ in Cartesian coordinates by solving Equation (36) for a point source with unit mass located in a cell at $(i^{\prime},j^{\prime},k^{\prime})=(1,1,1)$. Once we obtain ${\cal G}_{i-1,j-1,k-1}$ for $1\leq i\leq N_{x}$, $1\leq j\leq N_{y}$, and $1\leq k\leq N_{z}$, the values for other indices can be computed from the symmetry requirement

[TABLE]

for $|i-i^{\prime}|\leq N_{x}-1$, $|j-j^{\prime}|\leq N_{y}-1$, and $|k-k^{\prime}|\leq N_{z}-1$.

As we use the method presented in Section 3 to obtain the numerical solution of Equation (36), we must supply an appropriate boundary condition a priori.

It is reasonable to assume that far from the source, the DGF asymptotes to the CGF:

[TABLE]

Since high-order derivatives of the CGF are non-negligible close to the source, the near-field DGF that results from a second-order finite-difference approximation to the Poisson equation would deviate greatly from the CGF. How far should the boundary be away from the source to safely apply Equation (B2) as a proper boundary condition? The answer of James (1977) to this question was 16 cells. In fact, one can compare the first two terms in the asymptotic expansion of the DGF in Burkhart (1997) to verify that the relative deviation between the DGF and CGF is $0.1\%$ at a distance of $16$ cells, regardless of the total number of cells or the grid spacing.

To ensure that all boundaries are sufficiently away from the point source at $(i^{\prime},j^{\prime},k^{\prime})=(1,1,1)$, for the calculation of the DGF only, we extend the computational domain by adding 16 additional cells to one side in each direction. The newly added cells (except for the ghost cells) are numbered as $i=-15,-14,\cdots,-1$ in the $x$-direction (and similarly in the $y$- and $z$-directions). We then apply the boundary condition (Equation (B2)) to $i=-15$ and $N_{x}+1$ (and similarly for $j$ and $k$) and solve Equation (36) in the extended domain to obtain ${\cal G}_{i-1,j-1,k-1}$ for $i=-15,-14,\cdots,N_{x}+1$ (and similarly for $j$ and $k$). In practice, following the method of Section 3 we use the boundary condition to define a modified interior charge, and then employ sine transforms.

After obtaining the solution, we discard the portion pertaining to the extended part of the domain, and use Equation (B1) to calculate the Cartesian DGF for whole indices. Since the DGF is calculated once and for all in the initialization step, its contribution to the computational cost of an entire simulation is almost negligible. Figure 11 plots the resulting DGF as a function of the grid distance from the source, in comparison with the CGF. The CGF diverges at the source position, whereas the DGF remains finite everywhere. Although the DGF deviates from the CGF close to the source, their difference becomes smaller as the grid distance increases, consistent with the asymptotic expansion of Burkhart (1997).

B.2 Cylindrical Grid

In a cylindrical grid, we need to store only 16 kinds of the DGFs, such as ${\cal G}_{i,k^{\prime}}^{m}(\rm inn\rightarrow top)$, appearing in Equation (4.3.2)–(4.3.2), each of which represents the Fourier-transformed potential at one boundary due to point masses located at the same or different boundaries. These are calculated as follows.

To calculate the potential generated by point sources at the top vertical boundary, we place a point source with unit mass at $(i^{\prime},j^{\prime},k^{\prime})=(i^{\prime},1,N_{z}+1)$ for any $i^{\prime}\in[1,N_{R}]$ and solve Equation (39) numerically using the method of Section 3. This yields ${\cal G}_{i,i^{\prime},j-1,k-N_{z}-1}$ and its Fourier transform ${\cal G}_{i,i^{\prime},k-N_{z}-1}^{m}$ after applying Equation (41), but we only keep four 3D arrays ${\cal G}^{m}_{i,i^{\prime}}({\rm top\to top})={\cal G}^{m}_{i,i^{\prime},0}$, ${\cal G}^{m}_{i,i^{\prime}}({\rm top\to bot})={\cal G}^{m}_{i,i^{\prime},-N_{z}-1}$, ${\cal G}^{m}_{k,i^{\prime}}({\rm top\to inn})={\cal G}^{m}_{0,i^{\prime},k-N_{z}-1}$, and ${\cal G}^{m}_{k,i^{\prime}}({\rm top\to out})={\cal G}^{m}_{N_{R}+1,i^{\prime},k-N_{z}-1}$, corresponding to four DGFs due to the point sources at the top boundary. We repeat the above calculations by varying $i^{\prime}\in[1,N_{R}]$ to fill all the components of the four DGFs.

For the potential generated by point sources at the inner radial boundary, we place a point mass at $(i^{\prime},j^{\prime},k^{\prime})=(0,1,k^{\prime})$ for any $k^{\prime}\in[1,N_{z}]$. We then solve Equation (39) and apply Fourier transform to obtain ${\cal G}_{i,0,k-k^{\prime}}^{m}$. We keep only four 3D arrays ${\cal G}^{m}_{i,k^{\prime}}({\rm inn\to top})={\cal G}^{m}_{i,0,N_{z}+1-k^{\prime}}$, ${\cal G}^{m}_{i,k^{\prime}}({\rm inn\to bot})={\cal G}^{m}_{i,0,-k^{\prime}}$, ${\cal G}^{m}_{k-k^{\prime}}({\rm inn\to inn})={\cal G}^{m}_{0,0,k-k^{\prime}}$, and ${\cal G}^{m}_{k-k^{\prime}}({\rm inn\to out})={\cal G}^{m}_{N_{R}+1,0,k-k^{\prime}}$666Although ${\cal G}^{m}_{k-k^{\prime}}({\rm inn\to inn})$ and ${\cal G}^{m}_{k-k^{\prime}}({\rm inn\to out})$ can be stored in 2D arrays, we store them as 3D arrays for simple coding., corresponding to four DGFs due to the point sources at the inner radial boundary. We repeat the above calculations for all $k^{\prime}\in[1,N_{z}]$ to completely fill the elements of the arrays. We follow a similar procedure to obtain the remaining eight DGFs due to points sources at the bottom and outer radial boundaries.

The boundary condition in solving Equation (39) can be obtained by requiring that the DGF at a large distance from the source is approximately equal to

[TABLE]

The summation over $p$ in Equation (B3) is to add all the contributions from the periodic images of the point mass when the mass distribution holds $P$-fold symmetry in the $\phi$-direction.

As discussed in Section B.1, Equation (B3) remains valid as long as the distance between the cells $(i,j,k)$ and $(i^{\prime},j^{\prime},k^{\prime})$ is sufficiently large. Similarly to the Cartesian case, we extend the computational domain, only for the calculation of the DGF, by adding extra cells to each of the four boundaries. In the vertical and outer radial directions, 16 cells are wide enough to ensure that Equation (B3) is a good approximation to the DGF. In the inner radial direction, however, the 16-cell criterion based on a uniform grid spacing does not guarantee that Equation (B3) is a valid approximation especially when $f=(R_{\rm max}/R_{\rm min})^{1/N_{R}}$ is large. We numerically confirmed that using Equation (B3) as a Dirichlet boundary condition at the extended inner radial boundary causes large errors in the computation of the cylindrical DGF and ultimately gravitational potential returned from our Poisson solver when $R_{\rm max}/R_{\rm min}>10^{2}$.

It turns out that using the radial gradient of Equation (B3) as a Neumann boundary condition is a cure to this inner boundary problem. The Neumann condition on the extended inner boundary is not very accurate, either, but we empirically find that it enables the DGF to converge rapidly to the desired values and results in the accurate DGF at least in the original computational domain, if the domain is sufficiently expanded. We found that adding $N_{R}-16$ cells in the inner radial boundary ($16$ cells are kept for the outer radial boundary) produces enough accuracy for the DGF in the original domain. To illustrate this directly, we set up a logarithmic cylindrical grid with $64^{3}$ cells spanning $R\in[10^{-4},1]$, $\phi\in[0,2\pi]$, $z\in[-0.25,0.25]$. We then distribute additional $N_{R}=64$ cells, by adding 16 cells with $i=65,66,\cdots,80$ to outside of the outer radial boundary and $N_{R}-16=48$ cells with $i=-47,\cdots,-1,0$ to inside of the inner radial boundary, and calculate the DGF for a point source located at $(i^{\prime},j^{\prime},k^{\prime})=(0,1,1)$ for both Neumann and Dirichlet conditions at the extended inner boundary. Figure 12(a) compares the resulting ${\cal G}_{i,0,0,0}$ as functions of $R_{i}$ in logarithmic scale, with crosses and squares corresponding to the cases with Neumann and Dirichlet conditions, respectively. Figure 12(b) zooms in the section with $|R_{i}-R_{-1}|\leq 10^{-4}$ into linear scale, with the inset plotting the DGF from the Neumann condition for $-9.4\times 10^{-5}\leq R_{i}-R_{-1}\leq-9.0\times 10^{-5}$. It is clear that the DGF resulting from the Dirichlet condition is not symmetric with respect to the source, and cannot thus be considered correct. On the other hand, the DGF from the Neumann condition retains symmetry with respect to the source and converges to the CGF at large distances toward the outer boundary. The test problems presented in Section 5 verify that the DGF under the Neumann condition works very well for our purposes, yielding accurate gravitational potentials in the original domain even when $R_{\rm max}/R_{\rm min}=10^{5}$ (see Figure 10).

In the case of a uniform cylindrical grid, it is sufficient to add only $16$ cells to inside of the inner radial boundary either using Neumann or Dirichlet boundary condition777We use Neumann condition for consistency.. This is because $16$ cells are wide enough at both inner and outer boundary with uniform spacing. Since the radial coordinate of the extended inner boundary should be positive, this requires $R_{\rm min}>16\delta R$, or equivalently, $R_{\rm max}/R_{\rm min}<1+N_{R}/16$. This limitation is not be severe, given that the logarithmic cylindrical grid is more appropriate for large $R_{\rm max}/R_{\rm min}$.

Appendix C For Mass Distribution Under $P$-fold Azimuthal Symmetry

We often meet a problem that is periodic in the $\phi$-direction, with period $2\pi/P$ for a positive integer $P$. For example, gas flows in barred galaxies may have $P=2$ symmetry (e.g., Seo & Kim 2013), while $P=4$ for four-armed spiral galaxies (e.g., Dobbs et al., 2006; Shetty & Ostriker, 2008). In such situations, one can save computational time by restricting the domain to $\phi\in[0,2\pi/P]$, with a periodic boundary condition.

The interior solver presented in Section 3.2 works without modification even when $L_{\phi}=2\pi/P$ does not cover a full $2\pi$ domain, because the $P$-fold symmetry is automatically taken into account in Equation (23). This holds true also for the boundary solver (Section 4) as long as the periodic boundary condition with period $L_{\phi}$ is imposed in solving Equation (39). The cylindrical DGF defined as such then represents the gravitational potential from $P$ identical point masses lying along the azimuth with a uniform angular separation of $L_{\phi}$. When applying the boundary condition to the cylindrical DGF, therefore, one should add all the contributions from the image masses located at $\phi\in[L_{\phi},2\pi]$ to the gravitational potential of a real point mass located in the original domain with $\phi\in[0,L_{\phi}]$, as in Equation (B3). We note that the inclusion of the image masses is implicit in the DGF, so that all calculations are done in the original computational domain with size $L_{\phi}$, enabling a factor of $P$ reduction in the computational time as well as memory compared to the cases with a full $2\pi$-periodic domain.

To verify that our method handles $P$-fold azimuthal symmetry, we consider four uniform spheres located at $\phi_{0}$, $\phi_{0}+\pi/2$, $\phi_{0}+\pi$, and $\phi_{0}+3\pi/2$ in a uniform cylindrical grid presented in Table 1. We calculate the gravitational potential from the resulting mass distribution that clearly has $P=4$ symmetry over $\phi\in[0,2\pi]$. Figure 13 plots as black solid lines the mass density and the gravitational potential along the $\phi$-direction at $R=0.75$ and $z=0.0039$. We then recompute the gravitational potential of a single sphere at $\phi_{0}$ by reducing the computational domain to $\phi\in[0,\pi/2]$, and finally the gravitational potential of two spheres at $\phi_{0}+\pi$ and $\phi_{0}+3\pi/2$ in the domain covering $\phi\in[\pi,2\pi]$. The resulting potentials, plotted as red and blue dashed lines, are identical within machine precision to the potential from the full $2\pi$ domain. This confirms that our Poisson solver in a restricted $\phi\in[0,2\pi/P]$ domain correctly deals with mass distributions under $P$-fold azimuthal symmetry.