A Generalized Solution Method for Parallelized Computation of the Three-dimensional Gravitational Potential on a Multi-patch grid in Spherical Geometry

TL;DR

This paper introduces a new algorithm for efficiently computing the 3D gravitational potential on multi-patch spherical grids, reducing computational cost and communication overhead while maintaining accuracy.

Contribution

The proposed method computes the gravitational potential directly across all patches using spherical harmonics, improving efficiency and parallel scalability over previous approaches.

Findings

Achieves the same accuracy as previous methods

Reduces computational cost significantly

Improves parallel scaling performance

Abstract

We present a generalized algorithm based on a spherical harmonics expansion method for efficient computation of the three-dimensional gravitational potential on a multi-patch grid in spherical geometry. Instead of solving for the gravitational potential by superposition of separate contributions from the mass density distribution on individual grid patch our new algorithm computes directly the gravitational potential due to contributions from all grid patches in one computation step, thereby reducing the computational cost of the gravity solver. This is possible by considering a set of angular weights which are derived from rotations of spherical harmonics functions defined in a global coordinate system that is common for all grid patches. Additionally, our algorithm minimizes data communication between parallel compute tasks by eliminating its proportionality to the number of subdomains in the grid configuration, making it suitable for parallelized computation on a multi-patch grid configuration with any number of subdomains. Test calculations of the gravitational potential of a tri-axial ellipsoidal body with constant mass density on the Yin-Yang two-patch overset grid demonstrate that our method delivers the same level of accuracy as a previous method developed for the Yin-Yang grid, while offering improved computation efficiency and parallel scaling behaviour.

Tables

Table 1. Table 1: Averaged wall-clock time for computation of the gravitational potential on the Yin-Yang grid with 800 radial zones and 1 ∘ superscript 1 1^{\circ} angular resolution for different number of MPI processes using our algorithm (second column) and the algorithm by AMJM18 (third column).

Number of MPI processes	Wall-clock time [s]
	This work	AMJM18
1	106.6	159.5
2	53.07	80.02
4	26.04	39.18
8	13.68	20.59
16	7.310	11.10
32	4.283	6.561
64	2.526	4.022
92	1.960	3.164
184	1.340	2.233
368	1.007	1.719
782	0.816	1.553

Equations

$$\Phi (\mathbf{r})=-G\int \text{displaylimits}{\mathrm{d}}^3{\mathbf{r}}^{\prime }\frac{\rho ({\mathbf{r}}^{\prime })}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}.$$

$$\Phi (\mathbf{r})=-G\sum _{\ell =0}^{\infty }\frac{4\pi }{2\ell +1}\sum _{m=-\ell }^{\ell }{Y}_{\ell }^m(\theta ,\varphi )\cdot \left[{\mathcal{A}}_{\ell m}(r)+{\mathcal{B}}_{\ell m}(r)\right]$$

$${\mathcal{A}}_{\ell m}(r)=\frac{1}{r^{\ell +1}}\iint {\text{displaylimits}}_{4\pi }\mathrm{d}{\Omega }^{\prime }{Y}_{\ell }^{m\ast }({\theta }^{\prime },{\varphi }^{\prime })\int {\text{displaylimits}}_0^r\mathrm{d}{r}^{\prime }(r^{\prime }{)}^{\ell +2}\rho ({\mathbf{r}}^{\prime })$$

$${\mathcal{B}}_{\ell m}(r)=r^{\ell }\iint {\text{displaylimits}}_{4\pi }\mathrm{d}{\Omega }^{\prime }{Y}_{\ell }^{m\ast }({\theta }^{\prime },{\varphi }^{\prime })\int {\text{displaylimits}}_r^{\infty }\mathrm{d}{r}^{\prime }(r^{\prime }{)}^{1-\ell }\rho ({\mathbf{r}}^{\prime })$$

$$Y_{\ell }^m(\theta ,\varphi )=\sqrt{\frac{2\ell +1}{4\pi }\frac{(\ell -m)!}{(\ell +m)!}}\cdot P_{\ell }^m(\cos \theta )\cdot e^{im\varphi }$$

$$P_{\ell }^{-m}=(-1{)}^m\frac{(\ell -m)!}{(\ell +m)!}{P}_{\ell }^m,$$

$$\begin{gathered} \Phi (\mathbf{r})=-G\sum _{\ell =0}^{\infty }\sum _{m=0}^{\ell }\frac{2}{{\delta }_m}\frac{(\ell -m)!}{(\ell +m)!}{P}_{\ell }^m(\cos \theta ) \\ \cdot \left[{\mathcal{G}}_{\ell m}(r,\varphi )+{\mathcal{H}}_{\ell m}(r,\varphi )\right] \end{gathered}$$

$$\begin{gathered} {\mathcal{G}}_{\ell m}(r,\varphi )=\frac{1}{r^{\ell +1}}\iint {\text{displaylimits}}_{4\pi }\mathrm{d}{\Omega }^{\prime }{P}_{\ell }^m(\cos {\theta }^{\prime }) \\ \times \cos (m(\varphi -{\varphi }^{\prime }))\int {\text{displaylimits}}_0^r\mathrm{d}{r}^{\prime }(r^{\prime }{)}^{\ell +2}\rho ({\mathbf{r}}^{\prime }) \end{gathered}$$

$$\begin{gathered} {\mathcal{H}}_{\ell m}(r,\varphi )=r^{\ell }\iint {\text{displaylimits}}_{4\pi }\mathrm{d}{\Omega }^{\prime }{P}_{\ell }^m(\cos {\theta }^{\prime }) \\ \times \cos (m(\varphi -{\varphi }^{\prime }))\int {\text{displaylimits}}_r^{\infty }\mathrm{d}{r}^{\prime }(r^{\prime }{)}^{1-\ell }\rho ({\mathbf{r}}^{\prime }). \end{gathered}$$

$${\delta }_m=\{\begin{array}{ll}2,& \mathrm{if}m=0,\\ 1,& \mathrm{otherwise}.\end{array}$$

$${\mathcal{G}}_{\ell m}(r_n^{+},\varphi )=\frac{1}{(r_n^{+}{)}^{\ell +1}}\sum _{i=1}^n\sum _{j=1}^{N_{\theta }}\sum _{k=1}^{N_{\varphi }}{\rho }_{ijk}{\mathcal{R}}_{\mathrm{in},i}^{(\ell )}{\mathcal{T}}_j^{(\ell m)}{\mathcal{F}}_k^{(m)}$$

$${\mathcal{H}}_{\ell m}(r_n^{+},\varphi )=(r_n^{+}{)}^{\ell }\sum _{i=n+1}^{N_r}\sum _{j=1}^{N_{\theta }}\sum _{k=1}^{N_{\varphi }}{\rho }_{ijk}{\mathcal{R}}_{\mathrm{out},i}^{(\ell )}{\mathcal{T}}_j^{(\ell m)}{\mathcal{F}}_k^{(m)}$$

$${\mathcal{F}}_k^{(m)}=\cos (m\varphi ){\mathcal{C}}_k^{(m)}+\sin (m\varphi ){\mathcal{S}}_k^{(m)}.$$

$${\mathcal{T}}_j^{(\ell m)}=\int {\text{displaylimits}}_{{\theta }_j^{-}}^{{\theta }_j^{+}}\mathrm{d}{\theta }^{\prime }\sin {\theta }^{\prime }{P}_{\ell }^m(\cos {\theta }^{\prime }),$$

$${\mathcal{C}}_k^{(m)}=\int {\text{displaylimits}}_{{\varphi }_k^{-}}^{{\varphi }_k^{+}}\mathrm{d}{\varphi }^{\prime }\cos (m{\varphi }^{\prime }),$$

$${\mathcal{S}}_k^{(m)}=\int {\text{displaylimits}}_{{\varphi }_k^{-}}^{{\varphi }_k^{+}}\mathrm{d}{\varphi }^{\prime }\sin (m{\varphi }^{\prime }).$$

$${\mathcal{R}}_{\text{in},i}^{(\ell )}=\int {\text{displaylimits}}_{r_i^{-}}^{r_i^{+}}\mathrm{d}{r}^{\prime }(r^{\prime }{)}^{\ell +2}$$

$${\mathcal{R}}_{\text{out},i}^{(\ell )}=\int {\text{displaylimits}}_{r_i^{-}}^{r_i^{+}}\mathrm{d}{r}^{\prime }(r^{\prime }{)}^{1-\ell }$$

$$A_{C,i}^{(\ell m)}=\sum _{j=1}^{N_{\theta }}\sum _{k=1}^{N_{\varphi }}{\rho }_{ijk}{\mathcal{T}}_j^{(\ell m)}{\mathcal{C}}_k^{(m)}$$

$$A_{S,i}^{(\ell m)}=\sum _{j=1}^{N_{\theta }}\sum _{k=1}^{N_{\varphi }}{\rho }_{ijk}{\mathcal{T}}_j^{(\ell m)}{\mathcal{S}}_k^{(m)}.$$

$$\sum _{i^{\prime }=1}^i{\mathcal{R}}_{\mathrm{in},i^{\prime }}^{(\ell )}{A}_{C/S,i^{\prime }}^{(\ell m)}={\mathcal{R}}_{\mathrm{in},i}^{(\ell )}{A}_{C/S,i}^{(\ell m)}+\sum _{i^{\prime }=1}^{i-1}{\mathcal{R}}_{\mathrm{in},i^{\prime }}^{(\ell )}{A}_{C/S,i^{\prime }}^{(\ell m)}$$

$$\sum _{i^{\prime }=i}^{N_r}{\mathcal{R}}_{\mathrm{out},i^{\prime }}^{(\ell )}{A}_{C/S,i^{\prime }}^{(\ell m)}={\mathcal{R}}_{\mathrm{out},i}^{(\ell )}{A}_{C/S,i}^{(\ell m)}+\sum _{i^{\prime }=i+1}^{N_r}{\mathcal{R}}_{\mathrm{out},i^{\prime }}^{(\ell )}{A}_{C/S,i^{\prime }}^{(\ell m)}.$$

$$\begin{gathered} {\mathcal{A}}_{\ell m}(r)=\frac{1}{r^{\ell +1}}\sum _{g=1}^{N_g}\iint {\text{displaylimits}}_{{\Omega }_g}\mathrm{d}{\Omega }^{\prime }w({\theta }^{\prime },{\varphi }^{\prime })Y_{\ell }^{m\ast }({\theta }^{\prime },{\varphi }^{\prime }) \\ \times \int {\text{displaylimits}}_0^r\mathrm{d}{r}^{\prime }(r^{\prime }{)}^{\ell +2}\rho ({\mathbf{r}}^{\prime }) \end{gathered}$$

$$\begin{gathered} {\mathcal{B}}_{\ell m}(r)=r^{\ell }\sum _{g=1}^{N_g}\iint {\text{displaylimits}}_{{\Omega }_g}\mathrm{d}{\Omega }^{\prime }w({\theta }^{\prime },{\varphi }^{\prime })Y_{\ell }^{m\ast }({\theta }^{\prime },{\varphi }^{\prime }) \\ \times \int {\text{displaylimits}}_r^{\infty }\mathrm{d}{r}^{\prime }(r^{\prime }{)}^{1-\ell }\rho ({\mathbf{r}}^{\prime }) \end{gathered}$$

$$\sum _{g=1}^{N_g}\iint {\text{displaylimits}}_{{\Omega }_g}\mathrm{d}{\Omega }^{\prime }w({\theta }^{\prime },{\varphi }^{\prime })=4\pi .$$

$${\xi }_g=r\sin {\vartheta }_g\cos {\phi }_g,$$

$${\eta }_g=r\sin {\vartheta }_g\sin {\phi }_g,$$

$${\zeta }_g=r\cos {\vartheta }_g.$$

$$r=\sqrt{{\xi }_g^2+{\eta }_g^2+{\zeta }_g^2},$$

$${\vartheta }_g=\arccos ({\zeta }_g/r),$$

$${\phi }_g=\arctan ({\eta }_g/{\xi }_g).$$

Full text

methods: numerical — gravitation

††software: VisIt (Childs et al., 2012)

1 Introduction

When modeling self-gravitating systems a gravity solver that calculates the gravitational potential by solving the Poisson’s equation is one of the central components in the simulations. For multi-dimensional hydrodynamical simulations of self-gravitating flows the Poisson’s equation is solved at every hydrodynamical time step, and thus can be responsible for a significant fraction of the computational cost. Hence, several techniques for an efficient gravity solver have been developed over the past decades. An algorithm of choice is often decided by the complexity of the mass density distribution on the computational domain and the numerical technique used to solve the coupled hydrodynamic equations. For instance, hierarchical tree-based algorithm (e.g., Appel, 1985; Jernigan, 1985; Porter, 1985; Barnes & Hut, 1986; Hernquist & Katz, 1989) is usually the preferred choice in particle-based codes due to its flexibility in considering arbitrary geometry. On the other hand, for grid-based codes a solver based on Fast Fourier Transform (FFT) can easily be applied when the grid spacing is uniform (Hockney, 1970; Boris & Roberts, 1969). For problems with large spatial dynamic ranges such as cosmological structure formation or star cluster formation an adaptive mesh refinement (AMR) technique is employed in order to achieve high effective spatial resolution. On such adaptive grids iterative multi-grid gravity solvers (e.g., Ricker, 2008) are usually used, but an application of the tree-based solver on AMR grids has also recently been investigated (Wünsch et al., 2018).

In stellar hydrodynamics, three-dimensional (3D) simulations have mostly been performed on a spherical polar grid. For this class of simulations a common choice of gravity solvers is the algorithm based on a spherical harmonics expansion of the Green’s function developed by Müller & Steinmetz (1995). This gravity solver has been employed, for example, in several 3D CCSNe simulations (e.g., Vartanyan et al., 2019; Glas et al., 2018; Wongwathanarat et al., 2017; Lentz et al., 2015) due to its high computation efficiency in a case where the gravitational potential is dominated by the monopole term contribution from a central quasi spherical body (e.g., the proto-neutron star in CCSNe). It has also been adapted for a calculation on a Cartesian mesh (Couch et al., 2013), and is now implemented as a modular component in recent versions of the Flash code (Fryxell et al., 2000; Dubey et al., 2009). On the other hand, an alternative algorithm which solves the discretized Poisson’s equation on a spherical polar grid using FFT has recently been developed by Müller & Chan (2018). This algorithm gives more accurate solutions of the gravitational potential than those obtained by the multipole expansion technique for extremely asymmetric density configurations such as an off-center point mass. Such an algorithm will particularly be more advantageous in cases where multiple components of mass concentration are present on the grid.

Because spatial discretization of the computational domain by a spherical polar grid introduces a severe time step constraint imposed by small grid zones in the polar regions, modern multi-patch grid techniques in spherical geometry such as the Yin-Yang grid (Kageyama & Sato, 2004) and the cubed-sphere grid (Ronchi et al., 1996) have recently been receiving more attention. These grid techniques avoid coordinate singularities at the poles, and therefore help to speed up simulations by increasing the allowed time step size. At the Garching supernova group, the Yin-Yang grid, which is a two-patch overset grid configuration, is implemented into the finite-volume neutrino-radiation hydrodynamic code Prometheus-Vertex (Fryxell et al., 1989; Rampp & Janka, 2002), and is now being used extensively by for performing state-of-the-art calculations of core-collapse supernovae (CCSNe) in 3D (Summa et al., 2018; Melson et al., 2015).

Computation of self-gravity on a multi-patch grid in spherical geometry is non-trivial. Extensions of the spherical harmonics solver by Müller & Steinmetz (1995) and the FFT-based solver by Müller & Chan (2018) for a multi-patch grid are not readily available. For the case of the Yin-Yang grid, Wongwathanarat et al. (2010) resorted to a simple approach. They interpolated the density field on the Yin-Yang grid onto an auxilliary spherical polar grid, and then applied the algorithm by Müller & Steinmetz (1995) without any modification. While this approach provides an easy workaround to the solution, it is not ideal since it introduces an additional source of errors through interpolation of the density field. Moreover, the method cannot be easily and efficiently paralellized for computation with a large number of processes because of complicated data communication pattern.

A more efficiently parallelized approach for solving the gravitational potential directly on the Yin-Yang grid based on the algorithm by Müller & Steinmetz (1995) has recently been proposed by Almanstötter et al. (2018, hereafter AMJM18). In their method, the gravitational potential is computed by adding contributions from the mass density distribution on the Yin and the Yang grid section, which are evaluated separately. Since each grid patch in the Yin-Yang grid configuration is simply the low-latitude part of the usual spherical polar grid**,** solving for the gravitational potential using this approach is straightforward. Although this method eliminates the need for interpolation of the density field onto a spherical polar grid and is easily parallelized on distributed-memory systems, it is still not optimal. On the one hand, the computational cost is increased because it calculates two sets of potential. On the other hand, the data communication volume among parallel compute tasks is also enlarged by a factor of two when compared to computation on a spherical polar grid with the same number of compute tasks. The latter imposes a limit to the parallel scalability of the algorithm. Furthermore and most importantly, an extension of this algorithm for other multi-patch configurations in spherical geometry with a larger number of grid sections (e.g., the cubed-sphere grid) would significantly decrease its computational efficiency because the additional computational cost and the size of data communication are proportional to the total number of grid patches.

In this paper, we derive a new algorithm for efficient computation of the 3D gravitational potential on a multi-patch grid in spherical geometry based on spherical harmonics expansion. Our method is a generalization of the method by Müller & Steinmetz (1995). It differs from the previous algorithm by AMJM18 in that our method calculates the sum of all contributions to the gravitational potential from all grid patches in one computation step. It takes full advantage of the symmetry property of the multi-patch grid configuration when calculating angular weights for the density distribution on each grid patch by utilizing rotational transformations of spherical harmonics. Data communication between compute tasks in parallel computation is minimized such that there is no dependency on the number of grid patches in the mesh configuration. Consequently, this gravity solver is suitable not only for the Yin-Yang grid with two grid sections but also for other multi-patch configurations in spherical geometry that consist of a larger number of grid patches.

Our paper is organized as follows: We begin by summarizing the basic algorithm by Müller & Steinmetz (1995) for solving the 3D gravitational potential on a spherical polar grid in Section 2. Then, we present our generalization of the basic algorithm for a multi-patch grid configuration in spherical geometry, and give explicit formulae of angular and radial weights needed for reconstruction of the gravitational potential in Section 3. In Section 4, we detail our implementation of the new algorithm for the case of computation on the Yin-Yang overset grid. Our algorithm is validated with a test computation on the Yin-Yang grid. The results are shown in Section 5 along with a performance analysis of our algorithm in comparison with the previous method by AMJM18. We conclude with discussions in Section 6.

2 Base Algorithm on a spherical polar grid

2.1 Basic equations

A brief summary of the efficient algorithm by Müller & Steinmetz (1995) for solving the integral form of the Poisson’s equation is as follows. The Poisson’s equation in its integral form reads

[TABLE]

Here, $G$ is the gravitational constant, $\boldsymbol{r}=(r,\theta,\phi)$ is a coordinate vector in spherical polar coordinates, and $\rho(\boldsymbol{r})$ is the density distribution function. To solve this equation the Green’s function $\lvert\boldsymbol{r}-\boldsymbol{r^{\prime}}\rvert^{-1}$ is expanded into spherical harmonics. Following this expansion, the gravitational potential at a point $\boldsymbol{r}$ can then be expressed as

[TABLE]

where the radius dependence functions $\mathcal{A}_{\ell m}$ and $\mathcal{B}_{\ell m}$ describing contributions to the gravitational potential from the mass distribution inside and outside of a radial coordinate $r$ are defined as

[TABLE]

and

[TABLE]

with $\mathrm{d}\Omega=\sin{\theta}\,\mathrm{d}\theta\,\mathrm{d}\phi$. Here we use the definition

[TABLE]

for spherical harmonics of degree $\ell$ and order $m$ where $P_{\ell}^{m}$ are the associated Legendre polynomials. Complex conjugates of spherical harmonics $Y_{\ell}^{m}$ are denoted as $Y_{\ell}^{m\,*}$.

Using the definition in Eq. (5) and an identity for $P_{\ell}^{m}$ with negative orders,

[TABLE]

Eq. (2) can be re-written as

[TABLE]

where

[TABLE]

and

[TABLE]

The coefficient $\delta_{m}$ is defined by

[TABLE]

2.2 Discretized formulae on a spherical polar grid

Let us construct a spherical polar grid of $N_{r}\times N_{\theta}\times N_{\phi}$ grid cells, and use indices $i$, $j$, and $k$ to label the $i^{\mathrm{th}}$, $j^{\mathrm{th}}$, and $k^{\mathrm{th}}$ grid cell in the $r$-, $\theta$-, and $\phi$-direction, respectively. For each grid cell $\mathcal{Z}_{ijk}$ coordinates of the (lower)higher side of the cell in each coordinate direction are denoted by $r_{i}^{(-)+}$, $\theta_{j}^{(-)+}$, and $\phi_{k}^{(-)+}$. To compute the gravitational potential on this grid numerically using Eq. (7) we make two assumptions. Firstly, we assume that the density distribution inside a grid cell $\mathcal{Z}_{ijk}$ is constant, and is approximated by the cell-averaged density $\rho_{ijk}$. Secondly, we truncate the summation over spherical harmonic degree $\ell$ at a chosen value $\ell_{\mathrm{max}}$.

Consider a case that the gravitational potential is to be calculated at grid zone interfaces in the radial direction, $\mathcal{G}_{\ell m}$ and $\mathcal{H}_{\ell m}$ can be numerically computed for the $n^{\mathrm{th}}$ radial grid interface as

[TABLE]

and

[TABLE]

where

[TABLE]

The angular weights $\mathcal{T}_{j}^{(\ell m)},\mathcal{C}_{k}^{(m)},$ and $\mathcal{S}_{k}^{(m)}$ in Eqs. (11)–(13) are defined by

[TABLE]

[TABLE]

and

[TABLE]

The integrals $\mathcal{T}_{j}^{(\ell m)}$ can be computed analytically with help of recurrence formulae for the associated Legendre polynomials, while the integrals $\mathcal{C}_{k}^{(m)}$ and $\mathcal{S}_{k}^{(m)}$ are elementary. Analytic solutions to these integrals are already explicitly given by Zwerger (1995) and recently also in the work by AMJM18. Finally, the radial integrals

[TABLE]

and

[TABLE]

are also easy to compute.

Following the implementation steps suggested by Müller & Steinmetz (1995) the angular parts of the summations in Eqs. (11) and (12) defined by

[TABLE]

and

[TABLE]

are computed first. Then, the radial summations can be evaluated efficiently by utilizing recurrence relations

[TABLE]

and

[TABLE]

3 Generalized Algorithm for a multi-patch grid in spherical geometry

When solving for the gravitational potential on a multi-patch grid configuration in spherical geometry, e.g., the Yin-Yang grid (Kageyama & Sato, 2004) and the cubed-sphere grid (Ronchi et al., 1996), the angular integrals over the $4\pi$ spherical surface area in Eqs. (3) and (4) are splitted into $N_{g}$ parts with $N_{g}$ being the number of grid patches in the considered grid configuration. Each integration part is performed over the solid angle $\Omega_{g}$ covered by the $g^{\mathrm{th}}$ grid patch. The functions $\mathcal{A}_{\ell m}$ and $\mathcal{B}_{\ell m}$ can then be re-written as

[TABLE]

and

[TABLE]

where the complex conjugates of spherical harmonics are multiplied by a weight function $w(\theta,\phi)$ to account for overlapping area in the case of an overlapping grid configuration. This surface weight function takes a value of $1/n_{g}(\theta,\phi)$ where $n_{g}(\theta,\phi)$ is the number of grid patches covering an angular position $(\theta,\phi)$. By definition, $w$ thus equals to 1 in the nonoverlapping region. An integration of this weight function over the entire surface of all grid patches yields $4\pi$ steradian, i.e.

[TABLE]

Typically, multi-patch grids in spherical geometry are designed such that each patch is geometrically identical to ease complications in the implementation of a numerical scheme on such grids. Meshes are often constructed using coordinate systems that are local to each grid patch. These local coordinate reference frames are related by rotational transformation about the common coordinate origin. Evaluation of the angular parts of the integrations in Eqs. (23) and (24) now become much more complicated than in the case of a spherical polar grid since each spherical harmonic mode takes a different functional form on each grid patch due to coordinate transformations. Nevertheless, because under rotational transformation a spherical harmonic of degree $\ell$ and order $m$ is simply a linear combination of spherical harmonics of the same degree defined in the rotated reference frame and, in addition, because each grid patch is geometrically identical the integrations in Eqs. (23) and (24) can be simplified considerably.

Let $(r,\vartheta_{g},\varphi_{g})$ be the spherical coordinates of a point $\boldsymbol{r}$ in a reference frame $[\xi\eta\zeta]^{(g)}$ defined for the construction of the $g^{\mathrm{th}}$ grid patch. Since coordinates in reference frames of all grid patches transform only by rotations about the coordinate origin the radial coordinate $r$ remains equal for all grid patches, and is thus denoted without the subscript $g$. Components of the Cartesian coordinates $(\xi_{g},\eta_{g},\zeta_{g})$ are related to the spherical coordinates in the same reference frame by the usual coordinate transformation,

[TABLE]

and

[TABLE]

The corresponding inverse transformation reads

[TABLE]

and

[TABLE]

The rotational transformation from the $[\xi\eta\zeta]^{(g)}$ coordinate system to the $[xyz]$ coordinate system can be decomposed into three elemental rotations with the amount of rotations given by the Euler angles $\alpha_{g}$, $\beta_{g}$, and $\gamma_{g}$. In the $z$-$y^{\prime}$-$z^{\prime\prime}$ convention, the sequence of rotation is a rotation by an angle $\alpha_{g}$ around the $\zeta_{g}$ axis, then by an angle $\beta_{g}$ around the rotated $\eta_{g}^{\prime}$ axis, and finally by an angle $\gamma_{g}$ around the rotated $\zeta_{g}^{\prime\prime}$ axis (see Fig.1). Correspondingly, the rotational transformation matrix relating the Cartesian coordinates $(\xi_{g},\eta_{g},\zeta_{g})$ to $(x,y,z)$ is given by

[TABLE]

where we abbreviate $\cos{\{\alpha_{g},\beta_{g},\gamma_{g}\}}$ and $\sin{\{\alpha_{g},\beta_{g},\gamma_{g}\}}$ as $c_{\{\alpha_{g},\beta_{g},\gamma_{g}\}}$ and $s_{\{\alpha_{g},\beta_{g},\gamma_{g}\}}$, respectively, in the above matrix equation for compactness.

Once the rotation operator is defined relations between spherical harmonics in the $[xyz]$ and $[\xi\eta\zeta]^{(g)}$ coordinate systems are given by

[TABLE]

where coefficients of the linear combinations are elements of the Wigner D-matrix, $D_{m^{\prime}m}^{\ell}$ (Wigner, 1931). Naturally, the Wigner D-matrix is a function of the three Euler angles characterizing the rotational transformation between the two coordinate systems. With the rotational transformation defined using the $z$-$y^{\prime}$-$z^{\prime\prime}$ convention, elements of the Wigner D-matrix are expressed as 111A number of different notations of the Wigner D-matrix are used in the literature. Here, we adopt the same notation as in Morrison & Parker (1987) where $D_{m^{\prime}m}^{\ell}(\alpha,\beta,\gamma)$ is equal to $D_{m^{\prime}m}^{\ell}(-\alpha,-\beta,-\gamma)$ with the notation that is employed in the original work by Wigner (1931).

[TABLE]

with the reduced Wigner d-matrix, $d_{m^{\prime}m}^{\ell}$ given by

[TABLE]

The summation index $s$ is an integer starting from 0 to $s_{\mathrm{max}}$ which is set by $s_{\mathrm{max}}=\min{(\ell-m,\ell-m^{\prime})}$ such that arguments of the factorials in the denominator are always positive.

However, care must be taken when evaluating $d_{m^{\prime}m}^{\ell}(\beta)$ numerically. Computation of $d_{m^{\prime}m}^{\ell}(\beta)$ directly using Eq. (35) is known to suffer from serious loss of precision at high degree $\ell$ due to cancellation of terms consisting of huge floating point numbers (Tajima, 2015). To circumvent this problem we calculate $d_{m^{\prime}m}^{\ell}(\beta)$ by Fourier decomposition, which results from factorization of the second elemental rotation of the transformation (Edmonds, 1964). The reduced Wigner-d matrix $d_{m^{\prime}m}^{\ell}(\beta)$ of any angle $\beta$ can be computed by (Trapani & Navaza, 2006)

[TABLE]

with the Fourier coefficients computed by utilizing recurrence formulae

[TABLE]

[TABLE]

and

[TABLE]

The starting condition for these recursions, i.e. the apex of the $d_{m^{\prime}m}^{\ell}(\frac{\pi}{2})$ matrix pyramid, is given by $d_{00}^{0}(\frac{\pi}{2})=1$.

We proceed in our derivation by taking the complex conjugate of Eq. (33), and then substituting the result in Eqs. (23) and (24). The functions $\mathcal{A}_{\ell m}$ and $\mathcal{B}_{\ell m}$ now take the forms

[TABLE]

and

[TABLE]

where $\mathrm{d}\omega_{g}=\sin{\vartheta_{g}}\,\mathrm{d}\vartheta_{g}\,\mathrm{d}\varphi_{g}$. It is important to note that, along with the transformation of spherical harmonics, the angular integrals on each grid section in Eqs. (40) and (41) have also been transformed into integrals in the $[\xi\eta\zeta]^{(g)}$ coordinate system that is local to each grid patch. This coordinate transformation will allow us to fully exploit the symmetry property of the grid configuration when we evaluate these angular integrals numerically.

Inserting these results into Eq. (2), and using identities for the reduced Wigner d-matrix (Edmonds, 1964)

[TABLE]

[TABLE]

and the identity in Eq. (6) for the associated Legendre polynomials to eliminate terms containing spherical harmonics with a negative order, we obtain, after algebraic rearrangement, an expression for the gravitational potential

[TABLE]

The functions $\mathcal{I}_{\ell m}$ and $\mathcal{J}_{\ell m}$ expand into

[TABLE]

and

[TABLE]

with the definitions

[TABLE]

and

[TABLE]

where each of the symbols $\mathcal{l}$ and $\mathcal{m}$ represent modes $C$ or $S$. The integrands $\mathcal{P}_{C,g}^{(\ell m^{\prime})}$ and $\mathcal{P}_{S,g}^{(\ell m^{\prime})}$ are defined by

[TABLE]

and

[TABLE]

And finally, the four modes normalization factors $\mathcal{N}_{CC,g}^{(\ell mm^{\prime})}$, $\mathcal{N}_{CS,g}^{(\ell mm^{\prime})}$ $\mathcal{N}_{SS,g}^{(\ell mm^{\prime})}$, and $\mathcal{N}_{SC,g}^{(\ell mm^{\prime})}$ are given by

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

with the coefficients

[TABLE]

and

[TABLE]

It is worth noting that the expression for the gravitational potential on a spherical polar grid derived in Section 2 (Eq. 7) can easily be recovered by setting $N_{g}=1$ with the Euler angles $\alpha_{g}=\beta_{g}=\gamma_{g}=0$.

4 Implementation for computation on the Yin-Yang grid

In this section, we demonstrate how the algorithm that we derived in the previous section is applied to compute the gravitational potential on the Yin-Yang overset grid configuration. The Yin-Yang grid in its most basic configuration consists of two geometrically identical overlapping grid patches. Each grid section, Yin or Yang, is simply the low-latitude part of the usual spherical polar grid, and therefore forms an orthogonal grid on the surface of a sphere. For this particular reason, the algorithm for computation of the gravitational potential on the Yin-Yang grid is an easy extension of the base algorithm derived for the case of a spherical polar grid in Section 2. A pseudo-algorithm providing guidance for implementing our new gravity solver is presented in Section 4.3.

4.1 Yin-Yang grid orientation and transformations

First of all, we construct the Yin and the Yang grid which spans the angular ranges

[TABLE]

and

[TABLE]

in the colatitude and azimuthal directions of their local coordinate reference frames, respectively. The angular resolution in both coordinate directions on both grid patches is denoted by $\Delta$. We choose a radial grid which is equidistant with a radial grid resolution $\Delta_{r}$. It spans the range

[TABLE]

where $R_{\mathrm{ib}}$ and $R_{\mathrm{ob}}$ are the inner and the outer radius of the computational domain. Transformation rules between Cartesian and spherical coordinates are given by Eqs. (26)–(31). Because the Yin-Yang grid configuration is symmetric Cartesian coordinates of a point in both the Yin and the Yang coordinate system are transformed to coordinates of the other grid patch by a matrix equation of the same form, i.e.

[TABLE]

The Yin-Yang coordinate transformation matrix translates to a rotation about an axis $\boldsymbol{\hat{S}}=(\xi_{\mathrm{Yin/Yang}},\eta_{\mathrm{Yin/Yang}},\zeta_{\mathrm{Yin/Yang}})=(0,\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}})$ by an angle $\pi$. Hence, instead of aligning the $[xyz]$ coordinate reference frame with either the reference frame of the Yin or the Yang grid we choose to align the polar axis $\boldsymbol{\hat{z}}$ with the axis $\boldsymbol{\hat{S}}$, and define the relative orientation of the Yin grid with respect to the $[xyz]$ coordinate reference frame such that

[TABLE]

Transformation of the Yang coordinates to the $[xyz]$ system can then be obtained by combining Eq. (60) and (61). This particular choice of orientation results in a grid that possesses rotational symmetry of order 2 about all coordinate axes of the $[xyz]$ coordinate system (see Fig. 2).

Once the choice of grid orientation for the Yin-Yang grid with respect to the $[xyz]$ coordinate system is defined, the three Euler angles describing the rotational transformation from the Yin and the Yang coordinate system to the $[xyz]$ system can be computed by solving trigonometric equations resulting from Eqs. (32),(60), and (61). This yields

[TABLE]

and

[TABLE]

One can see that only the rotation angle $\gamma$ of the third elemental rotation differs between the two sets of Euler angles. As a result, this allows us to simplify calculations of angular weights necessary for computation of the gravitational potential by evaluating only one set of weights, either for the Yin or the Yang grid, and obtain the weights for the other grid patch by multiplying with a factor 1 or -1. The multiplication factor depends on the spherical harmonic mode. This will be demonstrated in the following steps.

4.2 Discretized formulae on the Yin-Yang grid

As in the case of computation on a spherical polar grid we approximate the density distribution inside a grid cell $ijk$ on the Yin and the Yang patch by the cell-averaged density $\rho_{ijk,\mathrm{Yin/Yang}}$. In addition,we also truncate the summation series over spherical harmonics degree at a degree $\ell_{\mathrm{max}}$, and an assumption for the weight accounting for overlapping surface area $w$ is applied. When computing surface integrals we assume a constant weight within an angular grid zone $jk$, thereby replacing the weight function $w$ by a surface-averaged value $w_{jk}=1-0.5\alpha_{jk}$ where $\alpha_{jk}$ is the fraction of overlapping surface area (see e.g., Wongwathanarat et al., 2010, for details).

To calculate the gravitational potential at cell vertices of a grid zone $ijk$ on the Yin-Yang grid we rewrite the potential in a compact form as

[TABLE]

The gravitational potential expressed in this form reflects directly the actual implementation of the method in our numerical code. The prefactors $\mathcal{Q}_{\mathcal{l},jk,g}^{(\ell m)}$ are defined by

[TABLE]

and

[TABLE]

where $\theta_{jk,g}^{+}=\theta(\vartheta_{j,g}^{+},\varphi_{k,g}^{+})$ and $\phi_{jk,g}^{+}=\phi(\vartheta_{j,g}^{+},\varphi_{k,g}^{+})$, both of which can be computed easily by utilizing coordinate transformations. Furthermore, by coordinate transformation rules in Eqs. (60) and (61), one finds that

[TABLE]

and

[TABLE]

Thus $\mathcal{Q}_{\mathcal{l},jk,\mathrm{Yang}}^{(\ell m)}=(-1)^{m}\mathcal{Q}_{\mathcal{l},jk,\mathrm{Yin}}^{(\ell m)}$, thereby allowing us to simply store only one set of weights in an actual computation.

On the other hand, the radial weights for reconstruction of the gravitational potential at the $n^{\mathrm{th}}$ radial grid interface, $\mathcal{M}_{\mathcal{l},i}^{(\ell m)}$ are computed as

[TABLE]

with

[TABLE]

and $\mathcal{R}_{\mathrm{in},i^{\prime}}^{(\ell)}$ and $\mathcal{R}_{\mathrm{out},i^{\prime}}^{(\ell)}$ defined by Eqs. (17) and (18). The angular weights $\mathcal{U}_{\mathcal{l},jk,g}^{(\ell m)}$ are defined by

[TABLE]

and

[TABLE]

with $\mathcal{N}_{\mathcal{l}\mathcal{m},g}^{(\ell mm^{\prime})}$ defined by Eqs.(51–54). It is worth noting that the surface weights $w_{jk}$ and the integrals $\mathcal{T}^{(\ell m^{\prime})}_{j}$, $\mathcal{C}^{(m^{\prime})}_{k}$, and $\mathcal{S}^{(m^{\prime})}_{k}$ (Eqs.14–16) take the same values on both the Yin and the Yang grid because of the symmetry of the grid configuration and, in addition, because these angular integrals are performed using coordinates that are local on each grid patch.

Finally, it is also important to note that both the prefactors $\mathcal{Q}_{\mathcal{l},jk,g}^{(\ell m)}$ and the angular weights $\mathcal{U}_{\mathcal{l},jk,g}^{(\ell m)}$ need to be evaluated only once at an initialization step. These coefficients can be re-used to compute the gravitational potential of any mass distribution represented on the Yin-Yang grid. This is valid under the assumption that the angular grid remains fixed throughout the simulation.

4.3 Computation steps

Consider a case in which the Yin-Yang grid is decomposed only in the angular directions into smaller subdomains. Assume that the number of subdomains equals the number of compute tasks, $N_{\mathrm{tasks}}$, being used for computation of the gravitational potential on a distributed memory system. For this computational setup the parallelized algorithm to compute the gravitational potential using Eq. (64) can be summarized into the following steps:

Compute and store prefactors $\mathcal{Q}_{\mathcal{l},jk,g}^{(\ell m)}$ (Eqs. 65–66) and angular weights $\mathcal{U}_{\mathcal{l},jk,g}^{(\ell m)}$ (Eqs. 69–70). 2. 2.

Each compute task calculates angular summations $E_{\mathcal{l},i,g}^{(\ell m)}=\sum_{j=1}^{N_{\theta}}\sum_{k=1}^{N_{\phi}}\rho_{ijk,g}\mathcal{U}_{\mathcal{l},jk,g}^{(\ell m)}$ (Eq. 68). 3. 3.

Perform a summation of $E_{\mathcal{l},i,g}^{(\ell m)}$ across all compute tasks to obtain $B_{\mathcal{l},i}^{(\ell m)}$. 4. 4.

Each compute task calculates radial summations $\mathcal{M}_{\mathcal{l},i}^{(\ell m)}$ (Eq. 67) by using recurrence relations

[TABLE]

and

[TABLE] 5. 5.

Multiply $\mathcal{M}_{\mathcal{l},i}^{(\ell m)}$ with prefactors $\mathcal{Q}_{\mathcal{l},jk,g}^{(\ell m)}$ and add contributions from all $\ell$ and $m$ moments to reconstruct the gravational potential.

To ensure consistency of results when computing the summations in step 2 and 3 in parallel we evaluate these summations by using the two-sum algorithm (Møller, 1965; Knuth, 1981). The implementation of the two-sum algorithm for the summation between compute tasks (step 3) using the Message Passing Interface (MPI) library follows that of He & Ding (2001).

The computation steps summarized above are very similar to the paralellized algorithm proposed by AMJM18. The differences between the two algorithms are as follows: First of all and most importantly, definitions of the angular weights $\mathcal{U}_{\mathcal{l},jk,g}^{(\ell m)}$ are different. In AMJM18, the angular weights are simply the integrations of spherical harmonics defined in the Yin or the Yang coordinate system, which are rotated with respect to each other. These weights take similar forms as those computed for a spherical polar grid, but are multiplied by the surface weight factor $w$ to account for grid overlaps. Because of this the radial weights $E_{\mathcal{l},i,g}^{(\ell m)}$ for each multipole moment of the expansion that results from these angular integrations cannot be directly added into one set of radial weights $B_{\mathcal{l},i}^{(\ell m)}$. The algorithm of AMJM18 thus computes two sets of potential from these radial weights in the subsequent steps, and then these are added together in the final computation step. That is, the gravitational potential is calculated by considering two sources corresponding to the mass density distribution on each grid patch separately. On the other hand, the angular weights $\mathcal{U}_{\mathcal{l},jk,g}^{(\ell m)}$ derived in our algorithm consider integrations of spherical harmonics defined in a global coordinate system that is common for all grid patches. These spherical harmonics functions are transformed into linear combinations of spherical harmonics defined in the local coordinate system of each grid patch, and are integrated. This transformation is directly reflected by the appearance of summations over all spherical harmonics order $m^{\prime}$ in Eqs. (69) and (70). As a result, our algorithm computes only one set of radial weights $B_{\mathcal{l},i}^{(\ell m)}$ which is used by all grid patches in subsequent steps to reconstruct the gravitational potential.

These fundamental differences lead to an improvement of the computational efficiency of the new gravity solver. It is easy to see that the operation counts in steps 4 and 5 of the algorithm are reduced by a factor of two compared with AMJM18. In addition, the size of data communication between MPI processes at step 3 decreases by half in comparison to AMJM18. The algorithm by AMJM18 exchanges $2\times N_{\textrm{grid}}\times N_{r}\times\frac{1}{2}[(\ell_{\textrm{max}}+1)^{2}+\ell_{\textrm{max}}+1]\times\log_{2}{N_{\textrm{tasks}}}$ floating-point numbers in total, assuming that the recursive doubling algorithm (Thakur et al., 2005) is used for the global reduction operation. In contrast, our new method eliminates the factor $N_{\textrm{grid}}$ from the expression.

5 Numerical test and analysis

5.1 Gravitational potential of a homogeneous ellipsoid

As a test case for our algorithm we calculate the gravitational potential of a homogeneous ellipsoidal body on a Yin-Yang grid setup with an angular resolution $\Delta$ of $1^{\circ}$ and an equidistant radial grid of 800 zones. The radius of the inner and outer grid boundary, $R_{\mathrm{ib}}$ and $R_{\mathrm{ob}}$, are set to 0 and 2, respectively. For the purpose of comparison with results computed by AMJM18 we employ the same parameters as listed in their work for this test setup. The surface of the tri-axial ellipsoid is defined by the equation

[TABLE]

with parameters for the semi-axes set to $a=1,b=1.5$ and $c=2$. Coordinates $X,Y$ and $Z$ relate to Cartesian coordinates in the Yin and the Yang grid system by

[TABLE]

[TABLE]

and

[TABLE]

i.e. the shortest principle axis of the ellipsoid is tilted with respect to the $\xi$-axis of the Yin grid by an angle $\frac{\pi}{8}$. The density $\rho$ for any given point inside the ellipsoidal surface is set to $\rho_{0}=1$, while $\rho=0$ outside of the ellipsoid.

The analytical solution of the gravitaional potential of this homogeneous ellipsoidal body is given by (Chandrasekhar, 1969). The solution at a point $\boldsymbol{R}=(X,Y,Z)$ reads

[TABLE]

with the functions $\mathfrak{A}(\boldsymbol{R}),\mathfrak{B}(\boldsymbol{R}),\mathfrak{C}(\boldsymbol{R})$ and $\mathfrak{D}(\boldsymbol{R})$ defined by

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

The value $u_{0}$ determining the lower limit of the integrations is $u_{0}=0$ for a point $\boldsymbol{R}$ that lies outside of the ellipsoidal surface. On the other hand, in the case of a point $\boldsymbol{R}$ inside of the ellipsoidal surface, $u_{0}$ is given by the positive root of the equation

[TABLE]

To compute this solution semi-analytically we solve Eq. (82) for $u_{0}$ by using the bisection method with the tolerance error of $10^{-14}$. Then, the functions $\mathfrak{A}(\boldsymbol{R}),\mathfrak{B}(\boldsymbol{R}),\mathfrak{C}(\boldsymbol{R})$ and $\mathfrak{D}(\boldsymbol{R})$ are integrated numerically using the Simpson’s rule as implemented in the Fortran subroutine qsimp (Press et al., 1986). The upper limit for the integrations is set to $10^{16}$, and the integrations are evaluated up to a fractional accuracy of $5\times 10^{-14}$.

Our result for this test is shown in Figure 3, which displays the color-coded distribution of relative errors in the equatorial slice through the Yin-Yang grid. The figure shows errors for the case computed with $\ell_{\mathrm{max}}=80$. The error distribution shows a maximum value at the surface of the ellipsoid, which is marked in Figure 3 by the black solid line to guide the eyes. This is because of poor representation of the ellipsoidal surface on the Yin-Yang grid with limited spatial resolution. We also compute the maximum relative error as a function of $\ell_{\mathrm{max}}$, and show the results in Figure 4. The maximum relative error rapidly decreases until it approaches an asymptotic value at $\ell_{\mathrm{max}}\sim 70$ since the spatial discretization error play a more dominant role at high values of $\ell_{\mathrm{max}}$. In addition, we implemented the algorithm by AMJM18 into our gravity solver, and perform the same test with their algorithm. As one can see from the curve of the maximum error versus $\ell_{\mathrm{max}}$ in Figure 4 our solver shows excellent agreement with the method by AMJM18 at all values of $\ell_{\mathrm{max}}$. We observe the maximum relative difference between results computed with our and their algorithm at a level of $\sim 10^{-9}$ only, even at high values of $\ell_{\mathrm{max}}$.

5.2 Performance and scaling efficiency

Table 1 lists the wall-clock time to solution averaged over 20 calculations of the gravitational potential on the Yin-Yang grid with $1^{\circ}$ angular resolution and 800 radial zones using different numbers of MPI processes. For comparison we show both timing data for our method and for our implementation of the method by AMJM18. The data is also plotted in Fig. 5 (upper panel) along with the strong scaling efficiency (bottom panel). These numbers are measured using the Intel Xeon Gold 6148 Processors equipped on the Cobra high-performance computing system at the Max Planck Computing and Data Facility.

Our data shows that by applying our new algorithm for computation on the Yin-Yang grid using a single CPU core the computational efficiency is increased by about 30% when compared with the method by AMJM18. A more detailed analysis of both methods reveals that about 70% of the total computing time is, in fact, spent to compute the radial weights $E_{\mathcal{l},i,g}^{(\ell m)}$ in step 2 due to additional costs associated with the usage of the two-sum algorithm. Although the operation count at this step is equal for both algorithms, we already observe a 30% gain here. This gain factor results from the fact that our algorithm computes only one set of radial weights instead of two sets as required by AMJM18 method. Consequently, the number of load/store instructions in the angular summation loops, which is the computational bottleneck in our implementations, is reduced. This demonstrates that it can be misleading to compare operation counts when gauging the relative computational efficiency between two algorithms.

As we increase the number of MPI compute tasks we begin to observe benefits from smaller data communication volume required by our algorithm when compared with the method by AMJM18. While the scaling efficiency of our method is improved only slightly relative to the algorithm of AMJM18, the wall-clock time to solution is reduced by almost a factor of two when using 782 MPI processes with respect to the method by AMJM18. We also point out that although the strong scaling efficiency we report in this work is very different from that which is shown by AMJM18, cross-comparison of the scaling efficiency should be taken with cautions since it depends strongly on details of how an algorithm is implemented and also on details of the system running the algorithm.

6 Discussions and conclusions

In this work, we have presented a generalization of the multipole expansion based gravity solver by Müller & Steinmetz (1995) for efficient computation of the 3D gravitational potential on a multi-patch grid configuration in spherical geometry. We derive explicit formulae of angular and radial weights for reconstruction of the gravitational potential by considering integrals of spherical harmonics defined in a global coordinate system that is common to all subdomains in the multi-patch grid configuration. These spherical harmonics functions are transformed into linear combinations of spherical harmonics defined in the local coordinate reference frame of each individual grid patch. This transformation eases complications of having to integrate different functions on different grid patches when evaluating these angular and radial weights numerically. Linear coefficients for the rotational transformation of spherical harmonics are given by elements of the well-known Wigner D-matrix (Wigner, 1931) that can be evaluated efficiently by recursion relations for any set of Euler angles characterizing the transformation between the local and global coordinate system.

We have applied our new algorithm for calculations of the 3D gravitational potential on the Yin-Yang overset grid. Validation of our algorithm is done by comparison of the numerical solution to a semi-analytical solution of the gravitational potential of a tri-axial ellipsoidal body with homogeneous mass density. For this test we computed using the maximum degree of the multipole expansion $\ell_{\textrm{max}}$ of up to 80. At this value of $\ell_{\textrm{max}}$ the numerical error of the gravitational potential is dominated by the spatial discretization error associated with the chosen grid resolution. It is important to note that a suitable choice of $\ell_{\textrm{max}}$ is problem- and resolution dependent. Judging from our experiences, in 3D simulations of CCSNe, which is one of the application areas of our new method, an $\ell_{\textrm{max}}$ of $\sim$20 should already be adequate for typical angular grid resolutions of 1–2 degree.

Our results demonstrate that our algorithm yields a solution that is as accurate as that obtained by the recent algorithm of AMJM18 proposed for the Yin-Yang grid. Performance wise, our algorithm benefits from reduced computational cost and smaller data communiation volume between parallel compute tasks, thus yielding a faster gravity solver with better parallel scaling efficiency in comparison with the previous method.

Our new algorithm is easy to implement into an exisiting solver that is based on the multipole expansion method because it involves only minor modification to the calculations of angular weights at an initialization step of the gravity solver. We present detailed implementation steps of the algorithm for the case of the Yin-Yang grid configuration in Section 4.3. These implementation steps can be applied also for computations on other multi-patch grids in spherical geometry. In the case of a non-orthogonal angular grid, computation of angular weights, which involves integrations of spherical harmonics, are more complicated than the computation on the Yin-Yang grid configuration that considers orthogonal angular meshes. Nevertheless, these integrations can either be approximated or evaluated by numerical integrations. Once the angular weights are computed, the remaining steps of the algorithm remain unchanged.

In a future work, we plan to implement this algorithm into our newly developed high-order finite-volume hydrodynamic code, Apsara (Wongwathanarat et al., 2016), which is capable of dealing with general multi-block structured grids in curvilinear coordinates. We also plan to investigate how our algorithm can be re-formulated such that it yields higher-order of accuracy of the solution.

The author is grateful to Ewald Müller for a careful reading of the manuscript, and to Ninoy Rahman and Tobias Melson for fruitful discussions. The author thanks also the anonymous referee for his/her constructive comments. Computations are carried out on the Cobra high-performance computing system at the Max Planck Computing and Data Facility.