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

## Abstract

We present a new method for solving the three-dimensional gravitational potential of a density field on the Yin-Yang grid. Our algorithm is based on a multipole decomposition and completely symmetric with respect to the two Yin-Yang grid patches. It is particularly efficient on distributed-memory machines with a large number of compute tasks, because the amount of data being explicitly communicated is minimized. All operations are performed on the original grid without the need for interpolating data onto an auxiliary spherical mesh.

## 1 Introduction

Solving for Newtonian gravity in three dimensions (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, which 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 nowadays 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; Melson et al., 2015a; Lentz et al., 2015; Melson et al., 2015b; Roberts et al., 2016; Müller et al., 2017; Summa et al., 2018; Ott 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. Since 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 three-dimensional 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.

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 three-dimensional 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.

## 2 Solving Poisson’s equation on a spherical grid

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

Given the density distribution , the gravitational potential at the location is determined by solving Poisson’s equation, which can be expressed in its integral form as

(1) |

where is the gravitational constant and denotes the whole computational domain. In the following, we employ spherical coordinates and express spatial vectors as . A decomposition of into spherical harmonics yields (Müller & Steinmetz, 1995, Eq. 8)

(2) |

In the latter equation, and are defined as

(3) | ||||

(4) |

where 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 , , and a second for the potential outside, ,

(5) |

where are the associated Legendre polynomials. The integrals for these inner and outer contributions can be expressed as

(6) |

and

(7) |

The normalization factor in Eq. (5) is given by

(8) |

with

(9) |

To discretize these equations, we divide the spatial domain into grid cells, each spanning from to with , , and . 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.,

(10) |

for , , and . This assumption allows for simplifying Eqs. (6) and (7), which can then be written as

(11) | ||||

(12) |

with

(13) |

In our chosen coordinates, the gravitational potential is computed at the cell interfaces. The radial index introduced above is equal to the cell index if . The two integrals and can be evaluated analytically,

(14) |

and

(15) |

The radial integrals in Eqs. (11) and (12) can also be computed directly,

(16) | ||||

(17) |

The remaining integrals in Eqs. (11) and (12),

(18) |

can be evaluted efficiently using recurrence relations (see Appendix A).

## 3 Parallelization of the method

For efficiently computing the gravitational potential on distributed-memory systems, it is important to minimize the amout 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 taks. This, however, would require a large amount of data being communicated very inefficiently. More precisely, roughly floating-point numbers would be sent serially, where the factor of accounts for the collection of the data and their subsequent re-distribution. 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 Eqs. (11) and (12) as

(19) | ||||

(20) |

and define

(21) | ||||

(22) |

we immediately see that the summation over the angular domain is separated from the remaining calculation. After each compute task has solved Eqs. (21) and (22) locally, a global summation of and 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 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 communication steps to compute the global sum (Thakur et al., 2005).

With our method, we thus communicate floating-point numbers serially, where the factor of accounts for the two sums, and , and the third factor counts the number of different terms for and . 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 two degrees, our procedure reduces the amount of data being sent as long as . In the core-collapse supernova context, a realistic setting of would be considerably smaller than this threshold (e.g., ), because the gravitational potential is dominated by the proto-neutron star, whose density distribution is spherically symmetric to good approximation.

The global summation of and 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.

## 4 Algorithm and its parallelization for the Yin-Yang grid

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 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 Fig. 1.

Let the coordinates on Yin and Yang be denoted as (, , ) and (, , ), respectively, and defined in the intervals

where is the outer radius of the grid. and are small values that control the width of the buffer zones in the Yin-Yang overlap region. Both grid patches consist of cells and have equal coordinate axes. The angular coordinates of a certain point on Yin are transformed into the Yang system according to

(23) | ||||

(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. 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.

In order to efficiently compute the gravitational potential on the parallelized Yin-Yang grid, we developed a new method based on Eqs. (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 being for grid cells inside the overlap and 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

(25) |

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

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 and 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,

(26) | ||||

(27) |

Note that due to the symmetry between both grid patches and their identical local coordinates, the values of , , , and 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 , , , and . Four different potentials are calculated using Eqs. (19), (20), and (5), namely , , , and . Here, the subscript of denotes the grid part on which and are computed, while the subscript of determines in which coordinate system the angles and in Eqs. (5), (19), and (20) are expressed.

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

(28) | ||||

(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 Fig. 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 one degree. 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 tasks as implemented in the unparallelized code by Wongwathanarat et al. (2010) would result in a flat curve in Fig. 2. Our approach is thus suitable for large-scale simulations involving thousands of processors.

## 5 Test calculations

In order 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 one degree
angular resolution and its gravitational potential was compared to the result
gained from an axisymmetric gravity solver^{8}

We constructed an ellipsoid of constant density with semi-axes , , and with and . In Cartesian coordinates, a point is inside the ellipsoid if

(30) |

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

(31) |

where and

(32) | ||||

(33) | ||||

(34) | ||||

(35) |

If the point is located inside the ellipsoid, . Otherwise, is given as the positive root of the equation

(36) |

We evaluated the root and the integrals numerically up to a precision of .

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

(37) | ||||

(38) | ||||

(39) |

where we have introduced a rotation angle in order to break the symmetry between both Yin-Yang grid patches. To compute the Cartesian coordinates for Yang, the transformation given by Eqs. (23) and (24) has to be applied to express the Yang coordinates in the Yin system.

To estimate the correctness of our three-dimensional 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) to 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 . We varied both the grid resolution and the maximum multipole order . In the following, we denote the models with and an angular resolution of two degrees as “low resolution” (LR). The runs with twice as many grid points in all directions are referred to as “high resolution” (HR).

In Fig. 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 adaption of the gravity solver performs correctly on the Yin-Yang grid.

With increasing maximum multipole order , the relative error decreases in all setups, since 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 Fig. 4 for the HR case with in the plane . 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 Eq. (5) and the term there). Since the ellipsoidal surface is not sampled to machine precision, we do not expect this either for the gravitational potential in the central region.

## 6 Discussion and Conclusions

In this work, we have presented a new method for solving the three-dimensional 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. Since 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 of the spherical harmonics decomposition determines the accuracy of the gravitational potential. The computational cost increases roughly as . 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 of should be able to capture most of the asymmetries.

## Acknowledgements

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).

## Appendix A Recurrence relations for the integrals of the associated Legendre polynomials

### Footnotes

### References

- Buras, R., Rampp, M., Janka, H.-T., & Kifonidis, K. 2006, A&A, 447, 1049
- Chandrasekhar, S. 1969, Ellipsoidal figures of equilibrium (Yale Univ. Press)
- Falle, S. A. E. G., Vaidya, B., & Hartquist, T. W. 2017, MNRAS, 465, 260
- Fryxell, B., Müller, E., & Arnett, D. 1989, in Proceedings of the 5 Workshop on Nuclear Astrophysics, ed. W. Hillebrandt & E. Müller, 100
- He, Y., & Ding, C. H. Q. 2001, J. Supercomput., 18, 259
- Hunter, J. D. 2007, CSE, 9, 90
- Kageyama, A., & Sato, T. 2004, GGG, 5, 9005
- Komatitsch, D., & Tromp, J. 2002, GeoJI, 150, 303
- Lentz, E. J., Bruenn, S. W., Hix, W. R., et al. 2015, ApJL, 807, L31
- Liebendörfer, M., Rampp, M., Janka, H.-T., & Mezzacappa, A. 2005, ApJ, 620, 840
- Marek, A., Dimmelmeier, H., Janka, H.-T., Müller, E., & Buras, R. 2006, A&A, 445, 273
- Marek, A., & Janka, H.-T. 2009, ApJ, 694, 664
- Melson, T., Janka, H.-T., Bollig, R., et al. 2015a, ApJL, 808, L42
- Melson, T., Janka, H.-T., & Marek, A. 2015b, ApJL, 801, L24
- Mondal, P., & Korenaga, J. 2018, GeoJI, 212, 1859
- Müller, B., Janka, H.-T., & Dimmelmeier, H. 2010, ApJS, 189, 104
- Müller, B., Melson, T., Heger, A., & Janka, H.-T. 2017, MNRAS, 472, 491
- Müller, E., & Steinmetz, M. 1995, CoPhC, 89, 45
- Oliphant, T. E. 2007, CSE, 9, 10
- Ott, C. D., Roberts, L. F., da Silva Schneider, A., et al. 2018, ApJ, 855, L3
- Peng, X., Xiao, F., & Takahashi, K. 2006, QJRMS, 132, 979
- Perez, F., & Granger, B. E. 2007, CSE, 9, 21
- Rampp, M., & Janka, H.-T. 2002, A&A, 396, 361
- Roberts, L. F., Ott, C. D., Haas, R., et al. 2016, ApJ, 831, 98
- Simon, J. B., Armitage, P. J., Li, R., & Youdin, A. N. 2016, ApJ, 822, 55
- Summa, A., Hanke, F., Janka, H.-T., et al. 2016, ApJ, 825, 6
- Summa, A., Janka, H.-T., Melson, T., & Marek, A. 2018, ApJ, 852, 28
- Takiwaki, T., Kotake, K., & Suwa, Y. 2012, ApJ, 749, 98
- —. 2014, ApJ, 786, 83
- Thakur, R., Rabenseifner, R., & Gropp, W. 2005, IJHPCA, 19, 49
- Wongwathanarat, A., Hammer, N. J., & Müller, E. 2010, A&A, 514, A48
- Zwerger, T. 1995, PhD thesis, Technische Universität München