Parallelized Solution Method of the Three-dimensional Gravitational Potential on the Yin–Yang Grid

Solving for Newtonian gravity in 3D is relevant in a large number of astrophysical and geophysical problems. For instance, the propagation of seismic waves on Earth (Komatitsch & Tromp 2002), the formation of planetesimals in protoplanetary disks (Simon et al. 2016), shock propagation in protostellar clouds (Falle et al. 2017), and the formation of the Earth's core (Mondal & Korenaga 2018) are 3D situations that require self-gravity to be taken into account. One example, which is in the focus of our interest, are core-collapse supernovae.

The explosion mechanism of core-collapse supernovae is one of the long-standing riddles in stellar astrophysics. Thanks to growing supercomputing power, 3D neutrino-hydrodynamics simulations can now be performed to study the physical processes responsible for the onset of the explosion with highly optimized codes on distributed-memory architectures (Takiwaki et al. 2012, 2014; Lentz et al. 2015; Melson et al. 2015a, 2015b; Roberts et al. 2016; Müller et al. 2017; Ott et al. 2018; Summa et al. 2018). The Garching group uses the Prometheus-Vertex package (Rampp & Janka 2002), which extends the finite-volume hydrodynamics module Prometheus (Fryxell et al. 1989) with a state-of-the-art neutrino transport and interaction treatment. Its latest code version applies the Yin–Yang grid (Kageyama & Sato 2004)—a composite spherical mesh—to discretize the spatial domain.

Until now, self-gravity of the stellar plasma in 3D simulations has been treated in spherical symmetry on an averaged density profile in Prometheus, because the gravitational potential is dominated by the spherical proto-neutron star in the center. As stellar collapse to neutron stars is only a mildly relativistic problem, the Garching group applies Newtonian hydrodynamics but uses a correction of the monopole of the gravitational potential (Marek et al. 2006), which has been turned out to yield results that are well compatible with fully relativistic calculations (Liebendörfer et al. 2005; Müller et al. 2010). However, as 3D core-collapse supernova simulations become more elaborate by taking more and more physical aspects into account, it is highly desirable to include a realistic 3D gravitational potential to treat large-scale asymmetries correctly. This holds true, in particular, in cases where the collapsing star develops a global deformation, e.g., due to centrifugal effects in the case of rapid rotation.

Müller & Chan (2018) recently presented a method for solving Poisson's equation on spherical polar grids using 3D Fast Fourier transforms. Although their algorithm yields an accurate solution of the gravitational potential even for highly aspherical density configurations, it is currently not available for the Yin–Yang grid, and its parallel efficiency needs to be improved to serve our purposes.

In this work, we present a method for efficiently computing the gravitational potential on the Yin–Yang grid based on the gravity solver of Müller & Steinmetz (1995). Also, Wongwathanarat et al. (2010) applied a 3D gravity solver on Yin–Yang data; however, their code was not parallelized for distributed-memory systems. It relied on mapping the data to an auxiliary spherical polar grid. With our approach presented here, we compute the gravitational potential on the Yin–Yang grid directly.

In Section 2, we will briefly summarize the algorithm developed by Müller & Steinmetz (1995) for 3D spherical grids. A method for its efficient parallelization will be discussed in Section 3, followed by a detailed explanation of the modifications for the Yin–Yang grid in Section 4. Test calculations will be presented in Section 5.

In this section, we briefly summarize the procedure for computing the 3D gravitational potential on a spherical polar grid as shown by Müller & Steinmetz (1995).

Given the density distribution, $\varrho ({\boldsymbol{r}})$, the gravitational potential at the location ${\boldsymbol{r}}$ is determined by solving Poisson's equation, which can be expressed in its integral form as

$$\begin{eqnarray}&&{\rm{\Phi }}({\boldsymbol{r}})=-G{\int }_{V}\displaystyle \frac{\varrho ({\boldsymbol{r}}^{\prime} )}{| {\boldsymbol{r}}-{\boldsymbol{r}}^{\prime} | }{d}^{3}{\boldsymbol{r}}^{\prime} ,\end{eqnarray} \tag{ 1 }$$

where G is the gravitational constant, and V denotes the whole computational domain. In the following, we employ spherical coordinates and express spatial vectors as ${\boldsymbol{r}}$ = (r,ϑ,φ). A decomposition of $| {\boldsymbol{r}}-{\boldsymbol{r}}^{\prime} {| }^{-1}$ into spherical harmonics Y_ℓm yields (Müller & Steinmetz 1995, Equation (8))

$$\begin{eqnarray}{\rm{\Phi }}({\boldsymbol{r}}) & = & -G\displaystyle \sum _{{\ell }=0}^{\infty }\displaystyle \sum _{m=-{\ell }}^{{\ell }}\displaystyle \frac{4\pi }{2{\ell }+1}{Y}_{{\ell }m}(\vartheta ,\varphi )\,\\ & & \times \,\left(\displaystyle \frac{1}{{r}^{{\ell }+1}}{C}_{{\ell }m}(r)+{r}^{{\ell }}{D}_{{\ell }m}(r)\right).\end{eqnarray} \tag{ 2 }$$

In the latter equation, C_ℓm and ${D}_{{\ell }m}$ are defined as

$$\begin{eqnarray}{C}_{{\ell }m}(r) & := & {\int }_{0}^{2\pi }{\int }_{0}^{\pi }d\varphi ^{\prime} d\vartheta ^{\prime} \sin \vartheta ^{\prime} \,\overline{{Y}_{{\ell }m}}(\vartheta ^{\prime} ,\varphi ^{\prime} )\\ & & \times {\int }_{0}^{r}{dr}^{\prime} {(r^{\prime} )}^{{\ell }+2}\varrho ({\boldsymbol{r}}^{\prime} ),\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}{D}_{{\ell }m}(r) & := & {\int }_{0}^{2\pi }{\int }_{0}^{\pi }d\varphi ^{\prime} d\vartheta ^{\prime} \sin \vartheta ^{\prime} \,\overline{{Y}_{{\ell }m}}(\vartheta ^{\prime} ,\varphi ^{\prime} )\\ & & \times {\int }_{r}^{\infty }{dr}^{\prime} {(r^{\prime} )}^{1-{\ell }}\varrho ({\boldsymbol{r}}^{\prime} ),\end{eqnarray} \tag{ 4 }$$

where $\overline{{Y}_{{\ell }m}}$ are the complex conjugates of the spherical harmonics. After inserting their definition and rearranging the terms, Müller & Steinmetz (1995) wrote the gravitational potential as a sum of two contributions: one for the gravitational potential inside a sphere of radius r, ${{\rm{\Phi }}}_{\mathrm{in}}^{({\ell }m)}({\boldsymbol{r}})$, and a second for the potential outside, ${{\rm{\Phi }}}_{\mathrm{out}}^{({\ell }m)}({\boldsymbol{r}})$,

$$\begin{eqnarray}&&{\rm{\Phi }}({\boldsymbol{r}})=-G\displaystyle \sum _{{\ell }=0}^{\infty }\displaystyle \sum _{m=0}^{{\ell }}{N}^{({\ell }m)}{P}_{{\ell }}^{m}(\cos \vartheta )[{{\rm{\Phi }}}_{\mathrm{in}}^{({\ell }m)}({\boldsymbol{r}})+{{\rm{\Phi }}}_{\mathrm{out}}^{({\ell }m)}({\boldsymbol{r}})],\end{eqnarray} \tag{ 5 }$$

where ${P}_{{\ell }}^{m}$ are the associated Legendre polynomials. The integrals for these inner and outer contributions can be expressed as

$$\begin{eqnarray}{{\rm{\Phi }}}_{\mathrm{in}}^{({\ell }m)}({\boldsymbol{r}}) & = & \displaystyle \frac{1}{{r}^{{\ell }+1}}{\int }_{0}^{2\pi }{\int }_{0}^{\pi }d\varphi ^{\prime} d\vartheta ^{\prime} \sin \vartheta ^{\prime} {P}_{{\ell }}^{m}(\cos \vartheta ^{\prime} )\\ & & \times \,\cos (m(\varphi -\varphi ^{\prime} )){\int }_{0}^{r}{dr}^{\prime} {(r^{\prime} )}^{{\ell }+2}\varrho ({\boldsymbol{r}}^{\prime} )\end{eqnarray} \tag{ 6 }$$

and

$$\begin{eqnarray}{{\rm{\Phi }}}_{\mathrm{out}}^{({\ell }m)}({\boldsymbol{r}}) & = & {r}^{{\ell }}{\int }_{0}^{2\pi }{\int }_{0}^{\pi }d\varphi ^{\prime} d\vartheta ^{\prime} \sin \vartheta ^{\prime} {P}_{{\ell }}^{m}(\cos \vartheta ^{\prime} )\\ & & \times \,\cos (m(\varphi -\varphi ^{\prime} )){\int }_{r}^{\infty }{dr}^{\prime} {(r^{\prime} )}^{1-{\ell }}\varrho ({\boldsymbol{r}}^{\prime} ).\end{eqnarray} \tag{ 7 }$$

The normalization factor in Equation (5) is given by

$$\begin{eqnarray}&&{N}^{({\ell }m)}:=\displaystyle \frac{({\ell }-m)!}{({\ell }+m)!}\displaystyle \frac{2}{{\delta }^{(m)}},\end{eqnarray} \tag{ 8 }$$

with

$$\begin{eqnarray}{\delta }^{(m)}:=\,\left\{\begin{array}{ll}2, & \mathrm{if}\quad m=0,\\ 1, & \mathrm{if}\quad m\gt 0.\end{array}\right.\end{eqnarray} \tag{ 9 }$$

To discretize these equations, we divide the spatial domain into ${n}_{R}\times {n}_{\vartheta }\times {n}_{\varphi }$ grid cells, each spanning from $({r}_{i}^{-},{\vartheta }_{j}^{-},{\varphi }_{k}^{-})$ to $({r}_{i}^{+},{\vartheta }_{j}^{+},{\varphi }_{k}^{+})$ with $i=1,\,\ldots ,\,{n}_{R}$, $j=1,\,\ldots ,\,{n}_{\vartheta }$, and $k=1,\,\ldots ,\,{n}_{\varphi }$. The finite-volume method as being used in the Prometheus code assumes that in each cell, the density is given as a cell average, i.e.,

$$\begin{eqnarray}&&\varrho ({\boldsymbol{r}})={\varrho }_{{ijk}},\end{eqnarray} \tag{ 10 }$$

for ${r}_{i}^{-}\leqslant r\leqslant {r}_{i}^{+}$, ${\vartheta }_{j}^{-}\leqslant \vartheta \leqslant {\vartheta }_{j}^{+}$, and ${\varphi }_{k}^{-}\leqslant \varphi \leqslant {\varphi }_{k}^{+}$. This assumption allows for simplifying Equations (6) and (7), which can then be written as

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{in}}^{({\ell }m)}({\boldsymbol{r}})=\displaystyle \frac{1}{{r}^{{\ell }+1}}\displaystyle \sum _{i=1}^{{n}_{r}}\displaystyle \sum _{j=1}^{{n}_{\vartheta }}\displaystyle \sum _{k=1}^{{n}_{\varphi }}{\varrho }_{{ijk}}\,{{ \mathcal R }}_{\mathrm{in},{\rm{i}}}^{({\ell })}\,{{ \mathcal T }}_{j}^{({\ell }m)}\,{{ \mathcal F }}_{k}^{(m)}(\varphi ),\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{out}}^{({\ell }m)}({\boldsymbol{r}})={r}^{{\ell }}\displaystyle \sum _{i={n}_{r}+1}^{{n}_{R}}\displaystyle \sum _{j=1}^{{n}_{\vartheta }}\displaystyle \sum _{k=1}^{{n}_{\varphi }}{\varrho }_{{ijk}}\,{{ \mathcal R }}_{\mathrm{out},{\rm{i}}}^{({\ell })}\,{{ \mathcal T }}_{j}^{({\ell }m)}\,{{ \mathcal F }}_{k}^{(m)}(\varphi ),\end{eqnarray} \tag{ 12 }$$

with

$$\begin{eqnarray}&&{{ \mathcal F }}_{k}^{(m)}(\varphi ):=\cos (m\varphi ){{ \mathcal C }}_{k}^{(m)}+\sin (m\varphi ){{ \mathcal S }}_{k}^{(m)}.\end{eqnarray} \tag{ 13 }$$

In our chosen coordinates, the gravitational potential is computed at the cell interfaces. The radial index, n_r, introduced above is equal to the cell index, i, if $r={r}_{i}^{+}$. The two integrals ${{ \mathcal C }}_{k}^{(m)}$ and ${{ \mathcal S }}_{k}^{(m)}$ can be evaluated analytically:

$$\begin{eqnarray}{{ \mathcal C }}_{k}^{(m)} & := & {\int }_{{\varphi }_{k}^{-}}^{{\varphi }_{k}^{+}}d\varphi ^{\prime} \cos (m\varphi ^{\prime} )\\ & = & \left\{\begin{array}{ll}{\varphi }_{k}^{+}-{\varphi }_{k}^{-}={\rm{\Delta }}{\varphi }_{k}, & \mathrm{if}\,m=0,\\ \displaystyle \frac{1}{m}[\sin (m{\varphi }_{k}^{+})-\sin (m{\varphi }_{k}^{-})], & \mathrm{if}\,m\gt 0,\end{array}\right.\end{eqnarray} \tag{ 14 }$$

and

$$\begin{eqnarray}{{ \mathcal S }}_{k}^{(m)} & := & {\int }_{{\varphi }_{k}^{-}}^{{\varphi }_{k}^{+}}d\varphi ^{\prime} \sin (m\varphi ^{\prime} )\\ & = & \left\{\begin{array}{ll}0, & \mathrm{if}\,m=0,\\ \tfrac{1}{m}[\cos (m{\varphi }_{k}^{-})-\cos (m{\varphi }_{k}^{+})], & \mathrm{if}\,m\gt 0.\end{array}\right.\end{eqnarray} \tag{ 15 }$$

The radial integrals in Equations (11) and (12) can also be computed directly:

$$\begin{eqnarray}&&{{ \mathcal R }}_{\mathrm{in},i}^{({\ell })}:={\int }_{{r}_{i}^{-}}^{{r}_{i}^{+}}{dr}^{\prime} {(r^{\prime} )}^{{\ell }+2}=\displaystyle \frac{1}{{\ell }+3}({({r}_{i}^{+})}^{{\ell }+3}-{({r}_{i}^{-})}^{{\ell }+3}),\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}{{ \mathcal R }}_{\mathrm{out},i}^{({\ell })} & := & {\int }_{{r}_{i}^{-}}^{{r}_{i}^{+}}{dr}^{\prime} {(r^{\prime} )}^{1-{\ell }}\\ & = & \left\{\begin{array}{ll}\mathrm{ln}({r}_{i}^{+})-\mathrm{ln}({r}_{i}^{-}), & \mathrm{if}\,{\ell }=2,\\ \tfrac{1}{2-{\ell }}[{({r}_{i}^{+})}^{2-{\ell }}-{({r}_{i}^{-})}^{2-{\ell }}], & \mathrm{if}\,{\ell }\ne 2.\end{array}\right.\end{eqnarray} \tag{ 17 }$$

The remaining integrals in Equations (11) and (12),

$$\begin{eqnarray}&&{{ \mathcal T }}_{j}^{({\ell }m)}:={\int }_{{\vartheta }_{j}^{-}}^{{\vartheta }_{j}^{+}}d\vartheta ^{\prime} \sin \vartheta ^{\prime} {P}_{{\ell }}^{m}(\cos \vartheta ^{\prime} ),\end{eqnarray} \tag{ 18 }$$

can be evaluated efficiently and analytically, i.e., without numerical integration errors, using recurrence relations (see the Appendix).

Solving Equation (5) numerically requires setting an upper bound for the summation over ℓ, which we denote as ${{\ell }}_{\max }$. We will discuss the choice of ${{\ell }}_{\max }$ below in Section 6.

For efficiently computing the gravitational potential on distributed-memory systems, it is important to minimize the amount of data being explicitly exchanged between compute tasks. The easiest parallel solution of the aforementioned equations would be to collect the entire density field from all computing units, calculate the gravitational potential in serial, and send the result back to all tasks. This, however, would require a large amount of data being communicated very inefficiently. More precisely, roughly ${n}_{R}\times {n}_{\vartheta }\times {n}_{\varphi }\times 2$ floating-point numbers would be sent serially, where the factor of two accounts for the collection of the data and their subsequent redistribution. The serial calculation of the gravitational potential would cause an additional load imbalance between the tasks. Especially if the number of grid cells is very large, this strategy is unrewarding.

Here, we present an approach for solving the gravitational potential accurately on a spherical grid, whose angular domain is decomposed among a large number of tasks, while the full radial slice is retained. This is precisely the setup used in our neutrino-hydrodynamics code, Prometheus-Vertex. If we rewrite Equations (11) and (12) as

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{in}}^{({\ell }m)}({\boldsymbol{r}})=\displaystyle \frac{1}{{r}^{{\ell }+1}}\displaystyle \sum _{i=1}^{{n}_{r}}{{ \mathcal R }}_{\mathrm{in},{\rm{i}}}^{({\ell })}(\cos (m\varphi ){A}_{{\rm{C}},i}^{({\ell }m)}+\sin (m\varphi ){A}_{{\rm{S}},i}^{({\ell }m)}),\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{out}}^{({\ell }m)}({\boldsymbol{r}})={r}^{{\ell }}\displaystyle \sum _{i={n}_{r}+1}^{{n}_{R}}{{ \mathcal R }}_{\mathrm{out},{\text{}}{\rm{i}}}^{({\ell })}(\cos (m\varphi ){A}_{{\rm{C}},i}^{({\ell }m)}+\sin (m\varphi ){A}_{{\rm{S}},i}^{({\ell }m)}),\end{eqnarray} \tag{ 20 }$$

and define

$$\begin{eqnarray}&&{A}_{{\rm{C}},i}^{({\ell }m)}:=\displaystyle \sum _{j=1}^{{n}_{\vartheta }}\displaystyle \sum _{k=1}^{{n}_{\varphi }}{\varrho }_{{ijk}}\,{{ \mathcal T }}_{j}^{({\ell }m)}\,{{ \mathcal C }}_{k}^{(m)},\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&{A}_{{\rm{S}},i}^{({\ell }m)}:=\displaystyle \sum _{j=1}^{{n}_{\vartheta }}\displaystyle \sum _{k=1}^{{n}_{\varphi }}{\varrho }_{{ijk}}\,{{ \mathcal T }}_{j}^{({\ell }m)}\,{{ \mathcal S }}_{k}^{(m)},\end{eqnarray} \tag{ 22 }$$

we immediately see that the summation over the angular domain is separated from the remaining calculation. After each compute task has solved Equations (21) and (22) locally, a global summation of ${A}_{{\rm{C}},i}^{({\ell }m)}$ and ${A}_{{\rm{S}},i}^{({\ell }m)}$ is performed. This is the only step, where explicit data communication is necessary. Parallel computing libraries provide efficient built-in methods to compute the sum across n_tasks compute tasks. In case of the commonly used Message Passing Interface (MPI), such a function is MPI_Allreduce, which is usually an implementation of the "recursive doubling" algorithm using ${\mathrm{log}}_{2}{n}_{\mathrm{tasks}}$ communication steps to compute the global sum (Thakur et al. 2005).

With our method, we thus communicate $2\times {n}_{R}\,\times \tfrac{1}{2}[{({{\ell }}_{\max }+1)}^{2}+({{\ell }}_{\max }+1)]\times {\mathrm{log}}_{2}{n}_{\mathrm{tasks}}$ floating-point numbers serially, where the factor of two accounts for the two sums, ${A}_{{\rm{C}},i}^{({\ell }m)}$ and ${A}_{{\rm{S}},i}^{({\ell }m)}$, and the third factor counts the number of different terms for $0\leqslant {\ell }\leqslant {{\ell }}_{\max }$ and $0\leqslant m\leqslant {\ell }$. The fourth factor is based on the assumption that the recursive doubling algorithm is used implicitly by the parallel computing library, as explained above. Considering a computational mesh with an angular resolution of 2°, our procedure reduces the amount of data being sent as long as ${{\ell }}_{\max }\lesssim 55$. In the core-collapse supernova context, a realistic setting of ${{\ell }}_{\max }$ would be considerably smaller than this threshold (e.g., $\sim 25$), because the gravitational potential is dominated by the proto-neutron star, whose density distribution is spherically symmetric to good approximation.

The global summation of ${A}_{{\rm{C}},i}^{({\ell }m)}$ and ${A}_{{\rm{S}},i}^{({\ell }m)}$ can be dominated by round-off errors, even in double floating-point precision. The reproducibility of the results is not guaranteed, because the summation order depends on the number of compute tasks and on implementation details of the message-passing library used for data communication. To tackle this issue, the "double-double summation" algorithm described by He & Ding (2001) is applied, which treats each number as a pair of the value itself and its error being accumulated during floating-point operations. With this method, we achieve reproducibility of the results up to machine precision.

Generally, only the terms that explicitly depend on the density have to be calculated in each step of a time-dependent simulation. All other terms remain constant and therefore need to be computed only once, as long as the grid coordinates do not change with time.

In its latest version, the Prometheus-Vertex code solves the hydrodynamics on the Yin–Yang grid (Kageyama & Sato 2004; Wongwathanarat et al. 2010), which is a combined mesh consisting of two low-latitude parts of a spherical polar grid. Its main advantage over a spherical polar grid covering the full 4π sphere is the exclusion of the grid singularities at the poles with their numerical difficulties. We show the geometry of the Yin–Yang grid in Figure 1.

Zoom In Zoom Out Reset image size

Let the coordinates on Yin and Yang be denoted as (r_Yin, ${\vartheta }_{\mathrm{Yin}}$, ${\varphi }_{\mathrm{Yin}}$) and (r_Yang, ${\vartheta }_{\mathrm{Yang}}$, ${\varphi }_{\mathrm{Yang}}$), respectively, and defined in the intervals

$$\begin{eqnarray*}0\leqslant r & \leqslant & {r}_{\max },\\ \displaystyle \frac{\pi }{4}-{\delta }_{\vartheta } & \leqslant & \vartheta \leqslant \displaystyle \frac{3\pi }{4}+{\delta }_{\vartheta },\\ -\,\displaystyle \frac{3\pi }{4}-{\delta }_{\varphi } & \leqslant & \varphi \leqslant \displaystyle \frac{3\pi }{4}+{\delta }_{\varphi },\end{eqnarray*}$$

where r_max is the outer radius of the grid. ${\delta }_{\vartheta }$ and ${\delta }_{\varphi }$ are small values that control the width of the buffer zones in the Yin–Yang overlap region. Both grid patches consist of ${n}_{R}\times {n}_{\vartheta }\times {n}_{\varphi }$ cells, and have equal coordinate axes. The angular coordinates of a certain point on Yin are transformed into the Yang system according to

$$\begin{eqnarray}&&{\vartheta }_{\mathrm{Yang}}=\arccos (\sin {\vartheta }_{\mathrm{Yin}}\sin {\varphi }_{\mathrm{Yin}}),\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&{\varphi }_{\mathrm{Yang}}=\arctan \left(-\displaystyle \frac{\cos {\vartheta }_{\mathrm{Yin}}}{\sin {\vartheta }_{\mathrm{Yin}}\cos {\varphi }_{\mathrm{Yin}}}\right),\end{eqnarray} \tag{ 24 }$$

where the inverse transformation is obtained by exchanging the subscripts.

The computation of the gravitational potential on the Yin–Yang grid is not trivial, because both grid patches are rotated against each other and are partially overlapping. Wongwathanarat et al. (2010) also used a version of Prometheus with the Yin–Yang grid. They computed the gravitational potential by interpolating the density from the individual Yin–Yang grid parts onto a new "auxiliary" spherical polar grid that covers the entire domain. The algorithm by Müller & Steinmetz (1995) could then be applied without any modifications. After that, the gravitational potential was mapped back to the Yin–Yang grid patches. Although the orientation and grid resolution of the auxiliary grid is arbitrary, Wongwathanarat et al. (2010) aligned it with the Yin part for simplicity.

The method by Wongwathanarat et al. (2010), however, has two major drawbacks for us: their code runs on shared-memory systems only, while our implementation of the gravity solver should be optimized for distributed-memory architectures; and collecting the density data from all compute tasks to a single worker would be very inefficient, due to the amount of explicit data communication and the serialization of the calculation (see Section 3). Furthermore, interpolation errors of the density and the calculated gravitational potential occur on the grid patch that is not aligned with the auxiliary spherical polar grid (Yang), while the data on the Yin part can simply be copied. This causes an artificial asymmetry between the Yin–Yang grid patches.

To efficiently compute the gravitational potential on the parallelized Yin–Yang grid, we developed a new method based on Equations (21) and (22). It benefits from the linearity of integrals, i.e., integrals over the entire computational domain can be written as a sum of two integrals computed separately on both Yin and Yang. Only the overlap regions between the grid patches have to be treated differently by introducing a surface weight, w_jk, being 0.5 for grid cells inside the overlap and one for cells outside. The weight factor for cells that are intersected by a boundary line, and therefore only partially overlapped, is calculated numerically (Peng et al. 2006). The sum of all surface elements multiplied with the surface weight yields

$$\begin{eqnarray}&&\displaystyle \sum _{j=1}^{{n}_{\vartheta }}\displaystyle \sum _{k=1}^{{n}_{\varphi }}(\cos {\vartheta }_{j}^{-}-\cos {\vartheta }_{j}^{+})({\varphi }_{k}^{+}-{\varphi }_{k}^{-})\,{w}_{{jk}}=2\pi .\end{eqnarray} \tag{ 25 }$$

Note that this summation is done on one grid patch (either Yin or Yang), which therefore yields 2π.

Similar to the parallelization procedure described in Section 3 being valid for spherical polar grids, our new approach for the Yin–Yang grid is based on computing ${A}_{{\rm{C}},i}^{({\ell }m)}$ and ${A}_{{\rm{S}},i}^{({\ell }m)}$ locally, followed by a global summation over all compute tasks. In the Yin–Yang setup, each task operates on data either from Yin or Yang, i.e., it calculates the following sums locally,

$$\begin{eqnarray}&&{A}_{{\rm{C}},i,\mathrm{Yin}/\mathrm{Yang}}^{({\ell }m)}=\displaystyle \sum _{j=1}^{{n}_{\vartheta }}\displaystyle \sum _{k=1}^{{n}_{\varphi }}{\varrho }_{{ijk},\mathrm{Yin}/\mathrm{Yang}}\,{w}_{{jk}}\,{{ \mathcal T }}_{j}^{({\ell }m)}\,{{ \mathcal C }}_{k}^{(m)},\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&{A}_{{\rm{S}},i,\mathrm{Yin}/\mathrm{Yang}}^{({\ell }m)}=\displaystyle \sum _{j=1}^{{n}_{\vartheta }}\displaystyle \sum _{k=1}^{{n}_{\varphi }}{\varrho }_{{ijk},\mathrm{Yin}/\mathrm{Yang}}\,{w}_{{jk}}\,{{ \mathcal T }}_{j}^{({\ell }m)}\,{{ \mathcal S }}_{k}^{(m)}.\end{eqnarray} \tag{ 27 }$$

Note that due to the symmetry between both grid patches and their identical local coordinates, the values of w_jk, ${{ \mathcal T }}_{j}^{({\ell }m)}$, ${{ \mathcal C }}_{k}^{(m)}$, and ${{ \mathcal S }}_{k}^{(m)}$ are equal on Yin and Yang.

After evaluating the latter equations, the global sums are computed separately for both grid parts resulting in different values for ${A}_{{\rm{C}},i,\mathrm{Yin}}^{({\ell }m)}$, ${A}_{{\rm{S}},i,\mathrm{Yin}}^{({\ell }m)}$, ${A}_{{\rm{C}},i,\mathrm{Yang}}^{({\ell }m)}$, and ${A}_{{\rm{S}},i,\mathrm{Yang}}^{({\ell }m)}$. Four different potentials are calculated using Equations (19), (20), and (5), namely ${{\rm{\Phi }}}_{\mathrm{Yin}}({{\boldsymbol{r}}}_{\mathrm{Yin}})$, ${{\rm{\Phi }}}_{\mathrm{Yin}}({{\boldsymbol{r}}}_{\mathrm{Yang}})$, ${{\rm{\Phi }}}_{\mathrm{Yang}}({{\boldsymbol{r}}}_{\mathrm{Yin}})$, and ${{\rm{\Phi }}}_{\mathrm{Yang}}({{\boldsymbol{r}}}_{\mathrm{Yang}})$. Here, the subscript of Φ denotes the grid part on which ${A}_{{\rm{C}},i}^{({\ell }m)}$ and ${A}_{{\rm{S}},i}^{({\ell }m)}$ are computed, while the subscript of ${\boldsymbol{r}}$ determines in which coordinate system the angles ϑ and φ in Equations (5), (19), and (20) are expressed.

The final gravitational potential for all compute tasks on Yin and Yang is, respectively,

$$\begin{eqnarray}&&{\rm{\Phi }}({{\boldsymbol{r}}}_{\mathrm{Yin}})={{\rm{\Phi }}}_{\mathrm{Yin}}({{\boldsymbol{r}}}_{\mathrm{Yin}})+{{\rm{\Phi }}}_{\mathrm{Yang}}({{\boldsymbol{r}}}_{\mathrm{Yin}}),\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&{\rm{\Phi }}({{\boldsymbol{r}}}_{\mathrm{Yang}})={{\rm{\Phi }}}_{\mathrm{Yin}}({{\boldsymbol{r}}}_{\mathrm{Yang}})+{{\rm{\Phi }}}_{\mathrm{Yang}}({{\boldsymbol{r}}}_{\mathrm{Yang}}).\end{eqnarray} \tag{ 29 }$$

Hence, the contributions from both grid parts are transformed into the local coordinates of the task and added. This procedure is completely symmetric with respect to both Yin–Yang grid patches and does not introduce a preferred direction. An interpolation onto an auxiliary grid is not required thus avoiding further numerical errors.

In Figure 2, we present the strong scaling behavior of our implementation of the Yin–Yang gravity solver. The gravitational potential was computed for a density distribution on the Yin–Yang grid with an angular resolution of 1°. While the grid configuration was kept constant, we altered the number of processors and measured the speedup given as the current wall-clock time normalized to the wall-clock time of the slowest run. The strong scaling deviates from the theoretically expected perfect behavior with increasing number of processors, because communication operations start to dominate. However, we emphasize that the computation of the gravitational potential by one single compute task as implemented in the unparallelized code by Wongwathanarat et al. (2010) would result in a flat curve in Figure 2. Our approach is thus suitable for large-scale simulations involving thousands of processors.

Zoom In Zoom Out Reset image size

To validate our approach for solving the gravitational potential on the Yin–Yang grid, we performed various tests. First, the gravitational potential of a spherically symmetric density distribution was verified to be equal to the corresponding analytical solution up to machine precision. Second, an axisymmetric density distribution was mapped to the Yin–Yang grid with 1° angular resolution, and its gravitational potential was compared with the result gained from an axisymmetric gravity solver³ yielding accurate agreement of 10⁻⁴. Third, the potential of a homogeneous ellipsoid was compared to the semi-analytical solution. We will elaborate on the latter test here in detail.

We constructed an ellipsoid of constant density, ϱ₀, with semi-axes a, b, and c with b = 1.5a and c = 2a. In Cartesian coordinates, a point is inside the ellipsoid if

$$\begin{eqnarray}&&{\left(\displaystyle \frac{x}{a}\right)}^{2}+{\left(\displaystyle \frac{y}{b}\right)}^{2}+{\left(\displaystyle \frac{z}{c}\right)}^{2}\leqslant 1.\end{eqnarray} \tag{ 30 }$$

The gravitational potential of this configuration is given by (Chandrasekhar 1969)

$$\begin{eqnarray}&&{\rm{\Phi }}({\boldsymbol{r}})=\pi G\,{\varrho }_{0}\,{abc}(A({\boldsymbol{r}})\,{x}^{2}+B({\boldsymbol{r}}){y}^{2}+C({\boldsymbol{r}}){z}^{2}-D({\boldsymbol{r}})),\end{eqnarray} \tag{ 31 }$$

where ${\boldsymbol{r}}=(x,y,z)$ and

$$\begin{eqnarray}&&A({\boldsymbol{r}})\,:={\int }_{{u}_{0}({\boldsymbol{r}})}^{\infty }{du}{[{({a}^{2}+u)}^{3}({b}^{2}+u)({c}^{2}+u)]}^{-1/2},\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&B({\boldsymbol{r}})\,:={\int }_{{u}_{0}({\boldsymbol{r}})}^{\infty }{du}{[({a}^{2}+u){({b}^{2}+u)}^{3}({c}^{2}+u)]}^{-1/2},\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}&&C({\boldsymbol{r}})\,:={\int }_{{u}_{0}({\boldsymbol{r}})}^{\infty }{du}{[({a}^{2}+u)({b}^{2}+u){({c}^{2}+u)}^{3}]}^{-1/2},\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&D({\boldsymbol{r}})\,:={\int }_{{u}_{0}({\boldsymbol{r}})}^{\infty }{du}{[({a}^{2}+u)({b}^{2}+u)({c}^{2}+u)]}^{-1/2}.\end{eqnarray} \tag{ 35 }$$

If the point ${\boldsymbol{r}}$ is located inside the ellipsoid, ${u}_{0}({\boldsymbol{r}})=0$. Otherwise, ${u}_{0}({\boldsymbol{r}})$ is given as the positive root of the equation

$$\begin{eqnarray}&&\displaystyle \frac{{x}^{2}}{{a}^{2}+{u}_{0}}+\displaystyle \frac{{y}^{2}}{{b}^{2}+{u}_{0}}+\displaystyle \frac{{z}^{2}}{{c}^{2}+{u}_{0}}=1.\end{eqnarray} \tag{ 36 }$$

We evaluated the root and the integrals numerically up to a precision of 10⁻¹³.

The transformation from the Yin spherical polar coordinates to Cartesian coordinates is given by

$$\begin{eqnarray}&&x=r\sin {\vartheta }_{\mathrm{Yin}}\cos ({\varphi }_{\mathrm{Yin}}+{\rm{\Delta }}\varphi ),\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}&&y=r\sin {\vartheta }_{\mathrm{Yin}}\sin ({\varphi }_{\mathrm{Yin}}+{\rm{\Delta }}\varphi ),\end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray}&&z=r\cos {\vartheta }_{\mathrm{Yin}},\end{eqnarray} \tag{ 39 }$$

where we have introduced a rotation angle ${\rm{\Delta }}\varphi =\pi /8$ to break the symmetry between both Yin–Yang grid patches. To compute the Cartesian coordinates for Yang, the transformation given by Equations (23) and (24) has to be applied to express the Yang coordinates in the Yin system.

To estimate the correctness of our 3D gravity solver for the case of the homogeneous ellipsoid, we compared the calculated gravitational potential in the entire computational domain (i.e., inside and outside the ellipsoid) with the semi-analytical solution. This was done on the Yin–Yang grid, but also on the standard spherical polar grid covering the whole sphere. Both mesh configurations had an equidistant radial spacing and an outer radius r_max = c. We varied both the grid resolution and the maximum multipole order ${{\ell }}_{\max }$. In the following, we denote the models with n_R =  400 and an angular resolution of 2° as low resolution (LR). The runs with twice as many grid points in all directions are referred to as high resolution (HR).

In Figure 3, we show the maximum relative error of the calculated gravitational potential with respect to the semi-analytical solution for the different mesh configurations and grid resolutions. As expected, the relative errors of the HR models are always smaller than in the LR cases, because the gravitational potential is sampled more accurately in the former runs and the surface of the ellipsoid is better resolved. While for the LR case the results computed on the Yin–Yang grid and on the spherical polar grid deviate from each other, both setups yield basically identical solutions in the HR case. This demonstrates that our adaptation of the gravity solver performs correctly on the Yin–Yang grid.

Zoom In Zoom Out Reset image size

With increasing maximum multipole order ${{\ell }}_{\max }$, the relative error decreases in all setups, as the shape of the ellipsoid is approximated more accurately. Ideally, an infinite number of multipole terms would be required to precisely reproduce an ellipsoidal configuration. Including more higher-order terms thus reduces the deviation from the semi-analytical solution.

The spatial distribution of the relative error on the Yin–Yang grid is shown in Figure 4 for the HR case with ${{\ell }}_{\max }=80$ in the plane z = 0. It can clearly be seen that the error is dominated by the cells close to the surface of the ellipsoid. Their deviation from the semi-analytical solution is two to three orders of magnitude larger than for cells near the grid origin. Machine precision is not reached in the center, because the gravitational potential there is also dependent on the density distribution outside (see Equation (5) and the term ${{\rm{\Phi }}}_{\mathrm{out}}^{({\ell }m)}$ there). As the ellipsoidal surface is not sampled to machine precision, we do not expect this either for the gravitational potential in the central region.

Zoom In Zoom Out Reset image size

In this work, we have presented a new method for solving the 3D gravitational potential of a density field on the Yin–Yang grid based on the algorithm developed by Müller & Steinmetz (1995). Our approach is particularly efficient on distributed-memory machines with a large number of compute tasks, because the amount of data being explicitly communicated is minimized and the work load is evenly distributed across all tasks resulting in a good scaling behavior.

In contrast to the strategy by Wongwathanarat et al. (2010), who interpolated the density from the Yin–Yang grid onto an auxiliary spherical mesh, we perform all operations on the original Yin–Yang grid. Our algorithm is completely symmetric with respect to the Yin–Yang grid patches and does not introduce a preferred direction. The validity of the results have been verified by performing test calculations that have an analytical or a semi-analytical solution.

Our ellipsoidal test configuration mimics the situation near the surface of the newly formed neutron star during the later post-bounce evolution of a core-collapse supernova. The increasingly steeper density gradient between the neutron star, which may be centrifugally deformed due to rapid rotation, and its surrounding medium is represented by the boundary of the ellipsoid. Because our solver can accurately determine the gravitational potential in this considered test setup with its extremely "sharp" surface, proper results can also be expected for the more moderate situation met during the birth of neutron stars during stellar core collapse and explosion.

The maximum multipole order ${{\ell }}_{\max }$ of the spherical harmonics decomposition determines the accuracy of the gravitational potential. The computational cost increases roughly as ${{\ell }}_{\max }^{2}$. A suitable balance between accuracy and computing time has to be found depending on the characteristics of the investigated physical problem. In core-collapse supernova models, the gravitational potential is dominated by the approximately spherically symmetric proto-neutron star. Therefore, a moderate ${{\ell }}_{\max }$ of $\sim 25$ should be able to capture most of the asymmetries.

We thank Alexander Summa and Annop Wongwathanarat for fruitful discussions. This project was supported by the European Research Council through grant ERC-AdG No. 341157-COCO2CASA and by the Deutsche Forschungsgemeinschaft through the Excellence Cluster "Universe" (EXC-153). The test calculations were performed on the Hydra system of the Max Planck Computing and Data Facility (MPCDF). The code scaling was tested on SuperMUC at the Leibniz Supercomputing Center with resources granted by the Gauss Centre for Supercomputing (LRZ project ID: pr53yi).

Software: Prometheus-Vertex (Fryxell et al. 1989; Rampp & Janka 2002; Buras et al. 2006), NumPy and Scipy (Oliphant 2007), IPython (Perez & Granger 2007), Matplotlib (Hunter 2007).

The integrals ${{ \mathcal T }}_{j}^{({\ell }m)}$ defined in Equation (18) can be calculated using the following recurrence relations (Zwerger 1995), which are evaluated in the order given below.

$$\begin{eqnarray}&&{{ \mathcal T }}_{j}^{(00)}={P}_{1}^{0}(\cos {\vartheta }_{j}^{-})-{P}_{1}^{0}(\cos {\vartheta }_{j}^{+}),\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}{{ \mathcal T }}_{j}^{({\ell }0)} & = & \displaystyle \frac{1}{2{\ell }+1}[{P}_{{\ell }-1}^{0}(\cos {\vartheta }_{j}^{+})-{P}_{{\ell }-1}^{0}(\cos {\vartheta }_{j}^{-})\\ & & +\,{P}_{{\ell }+1}^{0}(\cos {\vartheta }_{j}^{-})-{P}_{{\ell }+1}^{0}(\cos {\vartheta }_{j}^{+})]\\ & & \quad \,(\mathrm{only}\ \mathrm{for}\ {\ell }\gt 0),\end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray}&&{{ \mathcal T }}_{j}^{(11)}=\displaystyle \frac{1}{2}\left[{\vartheta }_{j}^{-}-{\vartheta }_{j}^{+}+\displaystyle \frac{1}{2}(\sin (2{\vartheta }_{j}^{+})-\sin (2{\vartheta }_{j}^{-}))\right],\end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray}&&{{ \mathcal T }}_{j}^{(22)}=\displaystyle \frac{1}{4}(\cos (3{\vartheta }_{j}^{+})-\cos (3{\vartheta }_{j}^{-}))-9(\cos {\vartheta }_{j}^{+}-\cos {\vartheta }_{j}^{-}),\end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray}{{ \mathcal T }}_{j}^{({\ell }{\ell })} & = & \displaystyle \frac{1}{{\ell }+1}[\cos {\vartheta }_{j}^{-}\,{P}_{{\ell }}^{{\ell }}(\cos {\vartheta }_{j}^{-})-\cos {\vartheta }_{j}^{+}\,{P}_{{\ell }}^{{\ell }}(\cos {\vartheta }_{j}^{+})\\ & & +\,{\ell }(2{\ell }-3)(2{\ell }-1){{ \mathcal T }}_{j}^{({\ell }-2,{\ell }-2)}],\end{eqnarray} \tag{ 44 }$$

$$\begin{eqnarray}&&{{ \mathcal T }}_{j}^{({\ell },{\ell }-1)}=-\displaystyle \frac{1}{{\ell }+1}[\sin {\vartheta }_{j}^{+}\,{P}_{{\ell }}^{{\ell }}(\cos {\vartheta }_{j}^{+})-\sin {\vartheta }_{j}^{-}\,{P}_{{\ell }}^{{\ell }}(\cos {\vartheta }_{j}^{-})],\end{eqnarray} \tag{ 45 }$$

$$\begin{eqnarray}{{ \mathcal T }}_{j}^{({\ell }m)} & = & \displaystyle \frac{1}{({\ell }+1)({\ell }-m)}[(2{\ell }-1)({\sin }^{2}{\vartheta }_{j}^{+}\,{P}_{{\ell }-1}^{m}(\cos {\vartheta }_{j}^{+})\\ & & -\,{\sin }^{2}{\vartheta }_{j}^{-}\,{P}_{{\ell }-1}^{m}(\cos {\vartheta }_{j}^{-}))\\ & & +\,({\ell }-2)({\ell }+m-1){{ \mathcal T }}_{j}^{({\ell }-2,m)}].\end{eqnarray} \tag{ 46 }$$