Issue 
A&A
Volume 514, May 2010



Article Number  A48  
Number of page(s)  14  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/200913435  
Published online  12 May 2010 
An axisfree overset grid in spherical polar coordinates for simulating 3D selfgravitating flows
A. Wongwathanarat  N. J. Hammer^{}  E. Müller
MaxPlanck Institut für Astrophysik, KarlSchwarzschildStraße 1, 85740 Garching, Germany
Received 9 October 2009 / Accepted 3 March 2010
Abstract
Aims. Three dimensional explicit hydrodynamic codes
based on spherical polar coordinates using a single spherical polar
grid suffer from a severe restriction of the time step size due to the
convergence of grid lines near the poles of the coordinate system. More
importantly, numerical artifacts are encountered at the symmetry axis
of the grid where boundary conditions have to be imposed that flaw the
flow near the axis. The first problem can be eased and the second one
avoided by applying an overlapping grid technique.
Methods. A type of overlapping grid in spherical
coordinates is adopted. This so called ``YinYang'' grid is a twopatch
overset grid proposed by Kageyama and Sato for geophysical simulations.
Its two grid patches contain only the lowlatitude regions of the usual
spherical polar grid and are combined together in a simple manner. This
property of the YinYang grid greatly simplifies its implementation
into a 3D code already employing spherical polar coordinates.
It further allows for a much larger time step in
3D simulations using high angular resolution ()
than that required in 3D simulations using a regular spherical
grid with the same angular resolution.
Results. The YinYang grid is successfully
implemented into a 3D version of the explicit Eulerian
gridbased code PROMETHEUS including selfgravity. The modified code
successfully passed several standard hydrodynamic tests producing
results which are in very good agreement with analytic solutions.
Moreover, the solutions obtained with the YinYang grid exhibit no
peculiar behaviour at the boundary between the two grid patches. The
code has also been successfully used to model astrophysically relevant
situations, namely equilibrium polytropes, a TaylorSedov explosion,
and RayleighTaylor instabilities. According to our results, the usage
of the YinYang grid greatly enhances the suitability and efficiency of
3D explicit Eulerian codes based on spherical polar
coordinates for astrophysical flows.
Key words: methods: numerical  hydrodynamics  gravitation  supernovae: general
1 Introduction
Three dimensional hydrodynamic simulations employing a single spherical polar grid are computationally expensive because of the convergence of grid lines towards the north and south pole. The converging grid lines imply a severe restriction of the time step size for any hydrodynamic code using explicit time discretization due to the CFL condition. This socalled ``pole problem'' bothers astrophysicists when simulating selfgravitating flow in three dimensions (e.g., convection in stars, or stellar explosions) where the spherical coordinate system is often preferable. In particular, simulations of corecollapse supernovae are a problem with which astrophysicists have been struggling. While observations show clear evidence of asymmetric (3D) complex structures in supernova ejecta, numerical simulations, in most cases, are carried out only in two spatial dimensions assuming axisymmetry (e.g., Blondin & Mezzacappa 2006; Scheck et al. 2006; Ohnishi et al. 2007). Three dimensional corecollapse supernova simulations are rare (e.g., Iwakami et al. 2008; Janka et al. 2005; Mezzacappa et al. 2006; Scheck 2006). In addition to the severe restriction of the time step size, boundary conditions that have to be imposed at the symmetry axis flaw the simulations near the axis by causing undesired numerical artifacts in 2D axisymmetric simulations, as e.g., jetxlike flow features (Kifonidis et al. 2003). In 3D simulations, the axis represents a coordinate singularity that almost unavoidably will leave its mark on the flow near or across the axis.There have been attempts to construct a new type of grid which is able to ease the pole problem. However, it is not possible to construct a single grid patch that can cover the entire surface of a sphere, is orthogonal, and at the same time does not contain any coordinate singularity except at the origin. Therefore, multipatch grid and overlapping (or overset) grid approaches are employed. They are widely used in the field of computational fluid dynamics where complex grid structures are common. For flows possessing an approximate global spherical symmetry, the ``cubed sphere'' grid (Ronchi et al. 1996) has been developed and is currently applied to several astrophysical problems (e.g. Koldoba et al. 2002; Fragile et al. 2009; Romanova et al. 2003; Zink et al. 2008). It is an overset grid consisting of six identical patches covering a solid angle of steradians. The ``YinYang'' grid has the latter property, too, but up to now it has not been used in astrophysical applications.
The YinYang grid was introduced by Kageyama & Sato (2004). It consists of two overlapping grid patches named ``Yin'' and ``Yang'' grid. In comparison with other types of overset grids in spherical geometry, the YinYang grid geometry is simple, as both the Yin and the Yang grid consist of a part of a usual spherical polar grid. The transformation of coordinates and vector components between the two patches is straightforward and symmetric, thus allowing for an easy and straightforward implementation of the grid into a 3D code already employing spherical polar coordinates. The YinYang grid is successfully used on massively parallel supercomputers in the field of geophysical science for simulations of mantle convection and the geodynamo. In these applications the thermal convection equation and the magnetohydrodynamic (MHD) equations are solved on the YinYang grid using a secondorder accurate finite difference method. Here, we also adopt the YinYang grid, and use it for astrophysically relevant (finitevolume) hydrodynamic simulations for the first time.
The paper is structured as follows. In Sect. 2, we describe the basics of the YinYang grid configuration including the transformations of coordinates and vectors between the Yin and Yang grid patches. In Sect. 3, we provide the details of the implementation of the YinYang grid into the PROMETHEUS hydrodynamic code, and also discuss the resulting necessary modifications of its 3D gravity solver that is based on spherical harmonics. In Sect. 4, we present the results of the test calculations we have performed including a test with selfgravity. In Sect. 5, we discuss the conservation problem arising when applying the YinYang grid. Then we report on the efficiency and performance gain obtained with the YinYang grid compared to a spherical polar grid in Sect. 6. Finally, we give the conclusions from our study in Sect. 7.
2 YinYang grid
The YinYang grid configuration is shown in Fig. 1. Both
the Yin and the Yang grid are simply a part of a usual spherical polar
grid and are identical in geometry. The angular domain of each grid
patch is given by
where and are the colatitude and azimuth, respectively. Note that it is necessary to add at least one extra buffer grid zone to both angular directions in order to ensure an appropriate overlap of the grids. The angular width of this buffer zone depends on the grid resolution, i.e., , where for simplicity we assumed equal angular spacing in  and direction. The angular domain is hereby extended by in both angular directions. The Yin and Yang grid are patched together in a specific manner forming a spherical shell with a small overlapping region covering approximately of a sphere's surface. Stacking up YinYang shells in radial direction results in a 3D grid that is identical to the usual spherical polar grid in radial direction. It is obvious that, unlike in the case of the spherical polar grid, the problematic high latitude sections of the sphere are avoided, and the angular zoning is almost equidistant.
Figure 1: An axisfree YinYang grid configuration plotted on a spherical surface. Both the Yin (red) and Yang (blue) grid are the low latitude part of the normal spherical polar grid and are identical in geometry. The Yang grid is obtained from the Yin grid by two rotations, and vice versa. 

Open with DEXTER 
The Cartesian coordinates
corresponding to the Yin grid, denoted by a superscript (n), and the Cartesian coordinates
corresponding to the Yang grid, denoted by a superscript (e), are related to each other through the transformation
where
This YinYang coordinate transformation can also be considered as two subsequent rotations. Accordingly, the transformation matrix M can be written as , where R_{x} and R_{z}are the transformation matrices of rotations by around the xaxis and by around the zaxis in counterclockwise direction, respectively. For the inverse transformation matrix M^{1} = M holds.
The relation between the spherical coordinates of the Yin and
Yang
grid patches can be derived directly from the transformation matrix
M. Because the YinYang coordinate transformation
involves only
rotations, it implies that the radial coordinate is identical on the
Yin and the Yang grid. The angular coordinates transform as
Note that the inverse transformation has the same form as (6) and (7) but exchanging the (grid) superscripts.
Vector components in spherical coordinates transform according
to
where
is the vector transformation matrix. When switching (grid) superscripts (e) and (n) in matrix P, the inverse vector transformation matrix is obtained. For a detailed derivation of the transformation matrix P, we refer to Sect. 3 of Kageyama & Sato (2004). Note that the vector transformation matrix P is singular at , but this singular point is rectifiable. In practice, one can always decompose vectors into their Cartesian components and perform the corresponding transformation.
3 Implementation
We have implemented the YinYang grid into our explicit finitevolume Eulerian hydrodynamics code, PROMETHEUS, which integrates the equations of multidimensional hydrodynamics using the piecewise parabolic method (PPM; Collela & Woodward 1984) and dimensional splitting. The code also includes a Poisson solver based on spherical harmonics to handle selfgravity.
3.1 Hydrodynamics solver
Firstly, the YinYang grid needs to be constructed. Since both the Yin
and the Yang grid are part of a spherical polar grid an analogous
spatial discretization in angular direction can be used. For example,
the
and
coordinates of the zone center of an angular
zone (j, k) of a YinYang grid,
having
zones in
direction
and
zones in direction,
are given
by
(10)  
(11) 
where
(12)  
(13) 
are the respective angular grid spacings.
The range of values for the colatitude and the azimuth angle are as given in (1), and for simplicity we set . In radial direction no modification is required. The geometric property of the YinYang grid allows us to make use of the coordinate arrays r_{i}, and twice by enforcing the same grid resolution for both grid patches. This approach avoids doubling the coordinate arrays.
Only simple modifications are needed concerning the data and program structure. Arrays with three spatial indices, e.g., i, j, and k, need an extra grid index, say, l. For example, the array for the density field will be instead of . As a consequence any triple loop running over indices i, j, and k in the program becomes a fourfold loop over i, j, k, and linstead. Otherwise, the YinYang grid allows one to exploit without any further modification any already implemented finitevolume scheme in spherical coordinates to solve the equations of hydrodynamics.
Different from the spherical polar grid, the YinYang grid requires no boundary conditions in angular directions. Each grid patch communicates with its neighboring patch using information from ghost zones that is obtained by interpolation of data between internal grid zones of the neighboring grid patch. Interpolation is only required in the two angular coordinates as the radial part of the YinYang grid is identical to that of a spherical polar grid. It is straightforward to determine the corresponding interpolation coefficients. The mapping of vector quantities between the Yin and Yang grid patches requires an additional step. After interpolating the vector components they must be transformed according to the transformation given in Eq. (8) from the Yin to the Yang angular coordinate system, and vice versa.
We tested two interpolation procedures. In the first one all primitive state variables (density, velocity, energy, pressure, temperature, abundances) are interpolated ignoring the resulting small thermodynamic inconsistencies. In the second procedure, we only interpolate the conserved quantities (density, momentum, total energy, and abundances), and compute the velocity and the remaining thermodynamic state variables consistently via the equation of state. Both procedures produce very similar results which differ at the level of the discretization errors. As the second procedure is more consistent we use it as the standard one in our code.
An example of overlapping situations which are encountered when using a YinYang grid is shown in Fig. 2. For simplicity, we use bilinear interpolation in order to prevent unwanted oscillation. Because the grid patches are fixed in both angular directions the interpolation coefficients for each ghost zone need to be calculated only once per simulation at the initialization step. After initialization, the coefficient map is stored in an array for later usage. Moreover, the symmetry property of the YinYang transformation allows one to make use of the interpolation coefficients twice for both grids.
Figure 2: A Mercator projection of an overlap region of the YinYang grid. In case of bilinear interpolation, four neighboring values of the underlying grid (red) will be used to determine the zonecentered value of a ghost zone in the grid on top (blue). The interpolation coefficients are determined by the relative distances, denoted by black lines, between the interpolation point (diamond) and the four neighboring points (crosses). 

Open with DEXTER 
Because the YinYang grid is an overlapping grid integral quantities such as the total mass or total energy on the computational domain cannot be obtained by just summing local quantities from every grid cells. Doing so will result in counting the contributions in the overlapping region twice. To circumvent this problem, weights are given to each grid zone during the summation. Suppose a grid zone has an overlapping volume fraction the cell will receive a weight . Zones in the nonoverlapping region receive the weight of 1.0, i.e., the entire zone contributes to the integral while, on the other hand, zones that are fully contained within the overlapping region have a weight 0.5. The volume fraction does not depend on the radial coordinate and can be thought of as an area fraction since the grid patches are not offset in radial direction. Prior to the area integration, one needs to determine for each zone interface of the underlying grid the points where the interface is intersected by the boundary lines of the other grid, e.g., points on the Yin grid intersected by the boundary lines of the Yang grid. The intersection points can be determined using the YinYang coordinate transformation in (6) and (7), respectively. The integration in the overlapping area is then carried out using the trapezoidal method. This procedure is also described in Peng et al. (2006). Once the area or volume fraction is calculated, the weights for each cell are obtained easily. Note that these weights need to be calculated only at the initialization step, and are stored for later usage in a coefficient map w(j,k), where j and k are the indices referring to the and coordinates, respectively. The coefficient map can be applied to both grids without any modification. Using the above described approach, the volume or surface area of the grids can be calculated with an accuracy up to machine precision.
3.2 Gravity solver
The 3D Newtonian gravitational potential is computed from Poisson's equation in its integral form using an expansion into spherical harmonics as described in Müller & Steinmetz (1995). Because the algorithm of these authors is based on a (single) spherical polar grid the density on the YinYang sphere has to be interpolated onto an auxiliary spherical polar grid. The interpolation used is firstorder accurate, and due to the simplicity of the YinYang grid configuration has to be performed only in the two angular dimensions. Concerning the resolution of the auxiliary grid, it is natural to employ the same grid resolution as that used for the YinYang grid in all three spatial dimensions. The orientation of the auxiliary grid can be chosen freely in principle. However, it is convenient to align it with one of the two grid patches (the Yingrid in our case). Once the density is interpolated onto the auxiliary spherical grid we compute the gravitational potential, as suggested by Müller & Steinmetz (1995), at zone interfaces instead of at zone centers on both the Yin and Yang grid. The gravitational acceleration at zone centers can then be obtained by central differencing the potential. Note that the interpolation coefficients for the density need to be calculated only once per simulation, because both the auxiliary grid and the YinYang grid are fixed in angular directions. In addition, all angular weights, Legendre polynomials, and their integrals required for the calculation of the gravitational potential are stored after the initialization step for later usage.
It is also possible to directly calculate the gravitational
acceleration at zone centers. The gravitational potential is given by
(see Eqs. (5)(7) in Müller
& Steinmetz 1995)
with
(15)  
(16) 
where Y^{lm} and Y^{lm*} are the spherical harmonics and their complex conjugates, is the density, and . The gravitational acceleration in radial direction is then
Writing the radial derivative in Eq. (17) as
and noticing that the first and third term on the right hand side of this expression cancel each other because of the identities
(19) 
and
(20) 
the gravitational acceleration in radial direction becomes
The corresponding expressions for the gravitational acceleration in the two angular directions are easy to obtain since the spherical harmonics Y^{lm} are the only angulardependent terms in Eq. (14). Therefore, we only need to consider the partial derivatives of the spherical harmonics with respect to the and coordinates. As the spherical harmonics are given by
(22) 
where N^{lm} is the normalization constant and P^{lm} the associated Legendre polynomial, one finds
(23) 
and
(24) 
The derivatives of the associated Legendre polynomials are easily obtained using the recurrence formula
(25) 
Thus, one finds for the gravitational acceleration in the angular directions
and
Obviously, the expressions for three components of the gravitational acceleration (see Eqs. (21), (26), and (27)) are similar to that for the gravitational potential itself (see Eq. (14)). Hence, besides computing derivatives of Legendre polynomials, our extended Poisson solver can provide without much additional effort both the gravitational potential and the corresponding acceleration.
Usage of the analytic expressions for the gravitational acceleration avoids the errors arising from the numerical differentiation of the gravitational potential. However, tests show that the results obtained using either the gravitational potential computed with the ``standard'' Poisson solver and subsequent numerical differentiation or directly the gravitational acceleration provided by the extended Poisson solver differ only very slightly (see next section). Thus, we decided to stick to the ``standard'' Poisson solver in our simulations and compute the gravitational acceleration by numerical differentiation, as it requires no modification of our code.
Figure 3: One dimensional profiles of density , pressure p, velocity in xdirection v_{x}, and specific internal energy e are shown along the xdirection at z^{(n)}=0.25 cm and y^{(n)}=0 cm (dasheddotted line in Fig. 4) for the shock tube simulation at every 0.1 s. Open and filled symbols represent data points on the Yin and Yang grid, respectively. Solid lines give the distributions calculated with an exact Riemann solver. 

Open with DEXTER 
4 Test suites
4.1 Sod shock tube
The first problem of our test suite is the planar Sod shock tube
problem, a classical hydrodynamic test problem (Sod
1978). We
simulated this (1D Cartesian) flow problem using spherical coordinates
and the YinYang grid. The initial state consists of two constant
states given by
where , p, and v_{x} are the density, pressure and the velocity in xdirection of the fluid, respectively. We assume the fluid to obey an ideal gas equation of state with an adiabatic index . The surface separating the two constant states is a plane orthogonal to the xaxis located at x^{(n)}=0.4 cm, the (positive) xaxis corresponding to a radial ray with angular coordinates and . Thus, this 1D planar Sod shock tube problem invokes all three spherical velocity components v_{r}, , and when simulating the flow in spherical polar coordinates. This allows us to test both the scalar and vector transformations as well as the interpolation between the Yin and Yang grid patches. The simulation was carried out on an equidistant YinYang grid of zones (i.e., with an angular resolution of one degree; see Eq. (1)). In radial direction the computational domain ranges from r=0.05 cm to r=1.0 cm. We impose a zerogradient boundary condition at both edges of the radial domain.
Figure 4: Snapshot of density contours in the meridional plane at t=0.15 s for the shock tube test problem. Dashed lines mark the YinYang grid boundary, while the dotted circular curves represent the inner and outer radial boundary of the computational domain, respectively. The one dimensional profiles shown in Fig. 3 are resampled along the dasheddotted line at z^{(n)}=0.25 cm. 

Open with DEXTER 
The solution of the shock tube problem is wellknown. We compare our results with the solution calculated using an exact Riemann solver (Toro 1997). For comparison, data are resampled along the xdirection with a spacing cm. Figure 3 shows one dimensional profiles of , p, v_{x}, and e (specific internal energy), respectively, along the xdirection at z^{(n)}=0.25 cm and y^{(n)}=0 cm (dasheddotted line in Fig. 4) at different times. Our results agree very well with the solution obtained with the exact Riemann solver. The grid resolution is sufficiently high to give a sharp shock front and contact discontinuity while the rarefaction wave is smooth. The shock position is correct at all time throughout the simulation. The resampled data yield an accuracy of approximately on average for shock positions. The shock wave and the contact discontinuity propagate smoothly across the YinYang boundary located at x^{(n)}=0.25 cm without any noticeable effect by the existence of the boundary. To illustrate this behavior, Fig. 4 shows lines of constant density in the meridional plane at time t=0.15 s. The isocontours are nearly perfectly straight lines perpendicular to the xaxis that are unaffected by the YinYang boundary (dashed line). The contour lines are slightly bent near the outer radial edge of the computational domain due to the zerogradient boundary condition we have imposed there.
In order to firmly demonstrate that the YinYang boundary does not cause numerical artifacts, we also computed this shock tube problem with a standard spherical polar grid using the same radial and angular resolution as for the YinYang grid described above, i.e., . We imposed reflecting boundary conditions in direction and periodic ones in direction. Figure 5 shows a comparison of the results obtained with both simulations. The two panels give the tangential velocity, defined as , in the meridional plane at time t=0.15 s for the YinYang grid (left), and the standard spherical polar grid (right), respectively. This velocity component should remain exactly zero because of the chosen initial conditions. Thus, it is a sensitive indicator whether the YinYang boundary works properly, which obviously is indeed the case as the left panel of Fig. 5 shows no hint of the location of that boundary. The modulus of the tangential velocity does nowhere exceed a value of 0.05 cm/s or approximately of the shock velocity (in xdirection) except near the outer radial edge of the grids, where the boundary condition causes larger numerical errors. Note that nonzero tangential velocities are encountered on both the YinYang grid and the standard spherical polar grid in the same grid regions at the same level. We thus conclude that they are the result of numerical errors that unavoidably occur when propagating a planar shock across a spherical polar grid, be it a standard one or a YinYang grid.
Figure 5: Color maps of the tangential velocity defined by in the meridional plane resulting from the (1D Cartesian) shock tube problem. The snapshots are computed using the YinYang grid ( left) and a standard spherical polar grid ( right) at a time t=0.15 s. On the left panel, red and blue lines mark the boundaries of the Yin and the Yang grid patches, respectively. On the right panel, the two red circles show the inner and outer boundary in the radial direction of the standard spherical polar grid. The labels at the color bars give the tangential velocity in units of cm/s. The color range is limited to 0.05 cm/s to emphasize the smallness of the tangential velocity far from the outer radial grid boundary. 

Open with DEXTER 
4.2 TaylorSedov explosion
As a second test for our code we consider the TaylorSedov explosion problem. We set up the initial state for the problem by mapping a spherically symmetric analytic solution (Landau & Lifshitz 1959) onto the computational grid. We choose the parameters of the problem to mimic a supernova explosion in an interstellar medium. Because the shock wave resulting from the explosion is spherically symmetric with respect to the center of the explosion, we assume the explosion center to be located at the point cm. Hence, this second test problem also involves a nonzero flux of mass, momentum, and energy across the YinYang boundary, and as the previous shock tube test, it probes whether that boundary causes any numerical artifacts.
The initial shock radius is cm orresponding to a time s past the onset of the explosion, and the explosion energy was set to E_{0} = 10^{51} erg. The ambient medium into which the shock wave is propagating is at rest. It has a constant density g/cm^{3}, and a constant pressure erg/cm^{3}. The fluid is described by an ideal gas equation of state with an adiabatic index , resulting in a density jump across the shock front of . We use a grid resolution of zones, a computational domain covering the radial interval cm, and employ a zerogradient boundary condition at both the inner and the outer radial boundary.
Our results are shown together with the analytic solution in Fig. 6. We have resampled our data and calculated radial profiles of the density , pressure p, and radial velocity v_{r} along a line in zdirection through the explosion center using a uniform radial spacing cm. As one can see the numerical results agree very well with the analytic solution. All flow quantities are smooth across the YinYang boundary, i.e., the shock wave passes that boundary without any noticeable numerical artifact. Due to the finite resolution the density jump across the shock front is slightly smaller in the simulation than the analytic value of four. However, the shock front is sharp throughout the whole simulation, and it propagates with the correct speed. One distinct feature of the TaylorSedov solution is its spherical symmetry. To illustrate that the YinYang grid does not destroy this symmetry of the solution, we show a set of lines of constant density in the meridional plane in Fig. 7. We also marked the line (dasheddotted) along which the data given in Fig. 6 are resampled. The contour lines, all of which are almost perfectly circular, are drawn at a simulation time s (i.e., time step number 1276) corresponding to an explosion time s.
Figure 6: Distributions of density ( top), pressure ( middle) and radial velocity ( bottom) versus radius from the explosion center (located at cm for the TaylorSedov explosion problem plotted at every 10^{11} s. Open symbols are data points from the Yin grid, while filled symbols represent sampled data from the Yang grid. The solid lines give the corresponding analytic solution. The data are resampled along the dasheddotted line shown in Fig. 7. 

Open with DEXTER 
Figure 7: Lines of constant density in the meridional plane obtained from our simulation of a TaylorSedov explosion. The snapshot is taken at a simulation time s which corresponds to an explosion time s. The dashed lines mark the YinYang boundary, while the two dotted circles represent the inner and outer radial boundary of the computational domain, respectively. The data presented in Fig. 6 are resampled along the dasheddotted line. 

Open with DEXTER 
Figure 8: Evolution of the mass within the overlap region for the TaylorSedov test case computed on the Yin grid, minus the mass computed on the Yang grid, , divided by the sum of these two masses. The dashed, dotted and solid lines give the solutions computed on a grid of zones, (i.e., angular resolution), zones (i.e., angular resolution), and zones (i.e., angular resolution), respectively. 

Open with DEXTER 
We further studied how the solution differs in the region where the Yin and Yang grid overlap. To this end we compare the total mass within the overlap region computed on the Yin and the Yang grid, respectively. Figure 8 shows the evolution of the relative mass difference, i.e., the mass within the overlap region computed on the Yin grid, minus the mass computed on the Yang grid, , divided by the sum of these two masses. We calculated this quantity for three different (angular) grids with zones (i.e., angular resolution), zones (i.e., angular resolution), and (i.e., angular resolution), respectively. For all three grid resolutions the relative mass difference has a value of about 10^{4}. Although its evolution with time is different in case of the simulation (because the coarse angular grid causes large errors when mapping the analytic initial data onto the grid which determine the further evolution), Fig. 8 shows that for an angular resolution better than the relative mass difference behaves similarly, its maximum value decreasing from at angular resolution to at angular resolution.
4.3 RayleighTaylor instability
We also simulated a single mode RayleighTaylor instability (RTI) on a
YinYang sphere. The initial configuration consists of a spherical
shell of a heavier fluid of density g/cm^{3}
that is
supported against a constant gravitational field g
= 1 cm/s^{2}pointing in negative radial
direction by a spherical shell of a
lighter fluid of density g/cm^{3}.
The boundary between
the two fluid shells is initially located at a radius r=0.5 cm.
To
balance the gravitational force, the initial (radial) pressure
distribution is set to
where P_{0}=1 erg/cm^{3}. A radial velocity varying in angular direction as the spherical harmonics with l=3and m=2 is used to perturb the initial configuration. The amplitude of the velocity perturbation is of the local sound speed c_{s}(r). Hence, the initial radial velocity is given by
The spherical harmonics are connected with the associated Legendre polynomials P_{l}^{m} via the expression
(31) 
The perturbation mode (l,m)=(3,2) yields a maximum radial velocity in the directions
(32) 
where . The remaining two velocity components of the perturbation mode are set equal to 0. The fluids are described by an ideal gas equation of state with an adiabatic index . The simulation is carried out on a YinYang grid of zones. To keep the fluid in hydrostatic equilibrium, a zerogradient boundary condition is used for both the inner and outer boundary in radial direction. The inner radial boundary is located at r=0.1 cm.
A snapshot of the resulting density distribution obtained with the YinYang grid is displayed in Fig. 9 at epoch t=2.85 s. The left panel shows color coded contour lines in 3D, and the right one a meridional cut at . The contour lines are drawn using different color tables for the Yin and Yang grid, respectively. Four distinct bubbles of rising low density fluid (Yin: blue; Yang: bright gray) are clearly visible that reflect the initial perturbation mode (l,m)=(3,2). High density fluid (Yin: yellow/red; Yang: dark gray/black) sinks down and settles at the inner part of the sphere. One can also notice KelvinHelmholtz instabilities developing at the surface of the bubbles. This is particularly obvious in the meridional cut (right panel). One of the RTI bubbles is within the Yang grid, while the three others reside on the Yin grid. It is obvious that the bubbles are distributed symmetrically following the perturbation pattern regardless of the grid patch. The 2D contour lines shown in the right panel of Fig. 9 emphasize this fact.
The RTI bubbles grow with nearly the same growth rate in all four (perturbation) directions, as can also be seen from Fig. 10 that displays the position of each bubble's head versus time. The four curves lie exactly on top of each other during the phase of linear growth. There are slight discrepancies between the four curves in the nonlinear regime, because the linear grid resolution in angular direction is slightly nonequidistant (due to its dependence). Two curves from the Yin grid coincide perfectly since they represent the two bubbles that lie symmetrically above and below the equator in the Yin grid. The results confirm that the YinYang grid does not favor any angular direction on the sphere. Since our aim was only to demonstrate this important fact, we do not further analyze the growth rate of the RTI.
4.4 Gravitational potential of homogeneous spheroids
We investigate the accuracy of our gravity solver by calculating the gravitational potential of homogeneous spheroids. We consider two homogeneous selfgravitating configurations: a prolate spheroid with an axis ratio of 0.7, and a sphere. The configurations have a constant density g/cm^{3}, and are embedded into a homogeneous background of much lower density g/cm^{3} in order to minimize the background's contribution to the gravitational potential. The semimajor axis of the spheroid aligns with the xaxis, while its center is placed at the origin of the YinYang grid. To provide a more difficult test for our multipole based gravity solver, we shift the center of the sphere off the origin of the computational grid by more than one sphere radius.
The analytical form of the gravitational potential for both type of configurations are known. The solution for the prolate spheroid can be found in Chap. 3 of Chandrasekhar (1969), and the sphere's potential can be easily calculated. Figure 11 shows contour lines of the gravitational potential for both cases in the meridional plane . The potential is calculated on a grid of zones with L=15, where L is the number of spherical harmonics taken into account (see Sect. 3.2). The contour lines are smooth across the YinYang boundary for both the prolate spheroid and the sphere.
Concerning the convergence behavior of the solver, we consider various grid resolutions and a number of spherical harmonics ranging up to L=25 for this convergence test. The grid resolutions used in the test are zones, zones, zones, and zones, respectively. The maximum and mean error of the gravitational potential are given as a function of L for both considered configurations in the middle and right panels of Fig. 11, respectively. Both errors show a convergence behavior with higher grid resolution, and tend to saturate at large values of L. This behavior is similar to what is described in Müller & Steinmetz (1995). In addition, for lower grid resolution the accuracy saturates at a lower number of spherical harmonics compared to calculations with a higher grid resolution. This is expected since higher order terms in the multipole expansion are not well represented on grids of lower angular resolution.
Figure 9: Surfaces of constant density in 3D ( left) and 2D ( right; meridional cut at ) resulting from the simulation of the RayleighTaylor instability described in the text at t=2.85 s. Contour lines on the Yin grid are shown using the blueyellow colors while contour lines on the Yang grid are displayed using the whiteblack colors. 

Open with DEXTER 
We also tested our extended Poisson solver discussed in Sect. 3.2. In Table 1 we compare the mean errors in the components of the gravitational acceleration for both the prolate spheroid and the sphere test case computed with the numerically differentiated gravitational potential given in Eq. (14) with those obtained from the analytic expression given in Eqs. (21), (26), and (27), respectively. We used a grid of zones and L=15 for this comparison.
Figure 10: Position of the heads of the RTI bubbles versus time. Red symbols (circles, triangles, and squares) show data from the Yin grid, while blue symbols (diamonds) represent data on the Yang grid. 

Open with DEXTER 
Figure 11: Contour lines of the gravitational potential ( left column) for two homogeneous selfgravitating configurations: a prolate spheroid ( top row) with an axis ratio of 0.7, and a sphere ( bottom row). The configurations are indicated by the darkgray shaded areas. Dashed lines show the YinYang boundary, while dotted lines indicate the outer radial boundary of the computational grid. The middle and right columns give the maximum and mean error of the numerically calculated gravitational potential for different grid resolutions as a function of the number of spherical harmonics used in our multipole gravity solver. The solid, dotted, dashed, and dasheddotted lines in both columns correspond to a grid resolution of zones, zones, zones, and zones, respectively. 

Open with DEXTER 
Table 1: Mean errors in the gravitational acceleration.
4.5 Selfgravitating polytropes
Using our YinYang grid based hydrocode we have also considered selfgravitating, nonrotating and rotating equilibrium polytropes. Both kinds of polytropes provide another test of the Poisson solver, and a test of how well our hydrodynamics code can keep a selfgravitating configuration in hydrostatic and stationary equilibrium, respectively. In addition, the rotating polytrope also serves to test the proper working of the YinYang boundary treatment, as it involves a considerable and systematic flow of mass, momentum and energy flux across that boundary due to the polytrope's rotation.
The polytropes have a polytropic index n=1, a polytropic constant , and a central density of g/cm^{3}. For our test runs we interpolated equilibrium polytropes calculated with the method of Eriguchi & Müller (1985) onto a YinYang sphere, and simulated their dynamic evolution (occurring as the interpolated configuration is not in perfect hydrostatic equilibrium). The central region (r < 1 km) of the polytrope is cut out and replaced by a corresponding point mass to allow for a larger time step.
We use an artificial atmosphere technique to handle those regions of the computational grid that lie outside the (rotating, i.e., nonspherical) polytrope. The density in the atmosphere is set equal to a value , where is the central density of the polytrope. Here, atmosphere denotes any grid zone whose density is less than the cutoff density . Furthermore, for all zones in the atmosphere the velocity is set to zero in order to keep the atmosphere quiet. This procedure is applied at the end of every time step throughout the simulation. A zerogradient boundary condition is imposed at the outer radial boundary, and a reflecting boundary condition at the inner one. The polytrope's evolution is followed for 10 ms corresponding to approximately 10 dynamic time scales in order to check how well the initial approximate equilibrium configuration is maintained by the YinYang code.
For the nonrotating polytrope, we employ a grid of zones. Note that we are able to use a relatively low angular resolution compared to the other tests, because the problem has spherical symmetry. Our results show that the polytrope stays perfectly spherically symmetric throughout the simulation, and that the nonradial velocities inside the polytrope remain zero. This demonstrates that the YinYang grid is able to preserve the initial spherical symmetry. Figure 12 shows the evolution of the central density (more precisely of the density of the innermost radial zone at r=1 km), which exhibits oscillations with an amplitude of the order of 10^{4} without any sign of a systematic trend. Comparing the initial radial distributions of the density (Fig. 13, upper panel) and the radial velocity (Fig. 13, lower panel) of the polytrope with those after 10 ms of evolution, we find no significant deviations. Relative changes in the density profile are of the order of 10^{4}, comparable to the size of the fluctuations of the central density. Only for zones near the edge of the polytrope the deviations can reach a level of up to , in particular in the zone next to the atmosphere. The figure also shows that data points from the Yin and the Yang grid lie on top of each other confirming that the code preserves the initial spherical symmetry of the polytrope very well. Except for the zones at the polytrope's surface, where the radial velocity is fluctuating at a level of approximately cm/s, the radial velocities are less than 10^{6} cm/s (i.e., less than of the local sound speed). Thus, we conclude that a nonrotating (n=1) equilibrium polytrope is correctly handled by our YinYang hydrocode.
Figure 12: Relative change of the central density of a nonrotating (nearly) equilibrium polytrope as a function of time. 

Open with DEXTER 
Figure 13: Density ( top) and radial velocity ( bottom) of a nonrotating n=1 equilibrium polytrope as a function of radius after t=10 ms of ``evolution''. In the top panel, the solid line shows the initial density profile. Red circles and blue triangles correspond to data from the Yin and the Yang grid, respectively. 

Open with DEXTER 
The rotating polytrope needs a higher grid resolution in direction, as it is no longer spherically symmetric. Thus, we used a grid resolution of zones for this simulation. The initial oblate equilibrium configuration has an axis ratio of 0.7. We, again, evolve the configuration for 10 ms to test the correct treatment of the situation by our YinYang hydrocode.
Figure 14 shows the relative variation of the central density as a function of time along an equatorial ray ( ; ) and along the pole ( ), respectively. One also recognizes a slight systematic trend in the behavior of the density fluctuation, which is steeper along the equator than at the pole. However, in both cases the relative increase of the central density is very small (10^{3}). The initial radial density profiles along the pole and the equator do not show any significant change during the 10 ms of evolution we have simulated with the YinYang code (Fig. 15, upper panel). The axis ratio has slightly increased to a value of 0.719. The radial velocities (Fig. 15, lower panel) are larger than in the nonrotating case by about an order of magnitude, because it is obviously more difficult to keep a rotating polytrope in equilibrium than a nonrotating (spherically symmetric) one. We again find the largest radial velocities (a few times 10^{8} cm/s) near the surface of the polytrope, especially along the equator. However, these velocities vary with time. When averaged over time (in the time interval t= [9, 10] ms) the profiles become flatter and the velocities smaller. This confirms that the polytrope is oscillating around its equilibrium configuration.
Figure 14: Same as Fig. 12 but for a rotating polytrope. The solid and dashed curves show the relative variation of the density along an equatorial ray ( ; ) and along the pole ( ), respectively. 

Open with DEXTER 
Figure 15: Density ( upper panel) and radial velocity ( lower panel) of a n=1 rotating polytrope in stationary equilibrium as a function of radius after t=10 ms of ``evolution''. In both panels the solid and dashed lines show the profiles along an equatorial ray ( , ) and along the pole ( ), respectively. Red circles and blue triangles in the upper panel correspond to data from the Yin and the Yang grid, respectively. In the lower panel, we show in addition time averaged (over the interval t= [9, 10] ms) velocity profiles along the equatorial ray (dotted) and the pole (dasheddotted). 

Open with DEXTER 
5 Conservation problem
The YinYang grid has a disadvantage common with other types of overlapping grids (see, e.g., Wang 1995; Wu et al. 2007; Chesshire & Henshaw 1994). The communication via interpolation between the two grid patches does not guarantee conservation of conserved quantities even though the finitevolume difference scheme employed on each grid patch is conservative. Nonconservation occurs when a flow across the YinYang boundary is present. This is the case in most of our tests except for the simulation of the nonrotating polytrope that involves only radial flow.
Nevertheless, we are still able to obtain sufficiently good results for all the test simulations discussed in the previous section. The degree of nonconservation is highly problem dependent. A simulation involving a considerable and systematic flow across the YinYang boundary, as e.g., in the case of the rotating polytrope, will result in a larger degree of nonconservation. We observe that the total mass increases by within 10 ms (or about ten dynamical timescales) in the case of the rotating polytrope. For the TaylorSedov test case, which is the cleanest test case in this respect (as it involves, e.g., no boundary effects like the shock tube, and e.g., no artificial atmosphere like the rotating polytrope), we find a mass loss of the order of 10^{5}, only. As Fig. 16 demonstrates this mass loss can be reduced by using a higher angular resolution.
Figure 16: Evolution of the relative mass loss, (MM_{0})/M_{0}, where M_{0} is the initial total mass, for the TaylorSedov test simulated on three different (angular) grids with zones (i.e., angular resolution; dashed line), zones (i.e., angular resolution; dotted line), and zones (i.e., angular resolution; solid line), respectively. 

Open with DEXTER 
Conservation of conserved scalar quantities can be obtained to machine precision by applying the algorithm described in detail in Peng et al. (2006), and summarized below. According to this algorithm scalar fluxes at the outer edges of boundary zones of both the Yin and the Yang grid are replaced by scalar fluxes computed using only ``interior'' fluxes from adjacent grid zones.
As an illustration, consider the YinYang grid overlap configuration in Fig. 17, where PQRS is a grid zone at the boundary of the Yang grid (blue) which overlaps with the underlying grid zone of the Yin grid (red). Fluxes referring to the Yin and the Yang grid are denoted by f and g, respectively.
The boundary flux
of the Yang grid is replaced by the flux
where and are the fluxes through the segments and , respectively.
The flux
in Eq. (33)
is calculated using
information from zone .
The evolution of a scalar quantity
of zone
is given by
Similarly, for the fraction of the zone defined by the polygon one has,
Assuming a piecewise constant state within the zone , Eqs. (34) and (35) lead to
where is the overlapping volume fraction (area) described in Sect. 3. Therefore,
Note that the flux is the only unknown in Eq. (37). Since the intersection points E and F are already known from the step to calculate the volume fraction , the lengths of all segments can be obtained. The flux is then given by
After obtaining the still missing flux in Eq. (33) by a similar procedure, the scalar quantity of the boundary zone PQRS is updated according to
(39) 
This procedure is then repeated to update all boundary grid zones.
After implementing the above algorithm we are able to conserve mass and total energy up to machine precision. However, the conservation of momentum is more complicated since the momentum equations in spherical coordinates involve not only flux (i.e., divergence) terms but also source terms (due to the presence of fictitious and pressure forces), and due to the ``mixing'' of momentum components as the Yin and Yang grid patches are rotated relative to each other (see Fig. 17).
As we have not yet devised and implemented a corresponding momentum conservation algorithm, momentum is not yet perfectly conserved in our code. For that reason we also refrain from using the scalar conservation algorithm described above, since in some simulations (e.g., in the TaylorSedov explosion simulation) we encountered a negative internal energy in some zones due to the inconsistency arising from the perfect conservation of mass and total energy on one hand and the imperfect conservation of momentum on the other hand. In our test runs the momentum violation is small, e.g., amounting to 0.24% (0.03%) angular momentum loss in the case of the rotating polytrope for a grid with three (one) degree angular resolution.
Figure 17: Illustration of the YinYang grid overlap configuration, where PQRS is a grid zone at the boundary of the Yang grid (blue) which overlaps with the underlying grid zone ABCD of the Yin grid (red). Fluxes referring to the Yin and the Yang grid are denoted by f and g, respectively. 

Open with DEXTER 
6 Performance and efficiency
One of the main purposes in implementing the YinYang grid is to ease the severe restriction imposed on the size of the time step for any explicit hydrodynamics scheme by the CFL condition in the polar regions of 3D simulations using a grid in spherical polar coordinates. In most applications the size of the time step is restricted most strongly by the size of the zones in direction, which is smaller than the size in direction by the factor assuming an equal angular resolution in both angular directions.
For a spherical polar grid the factor
implies (assuming
zone centered variables) a minimum zone size
(in radians) in direction for the first zone next to the pole. Typically, . On the other hand, applying the YinYang grid yields
for the size of the smallest zone in direction, which is typically about 0.7. Hence, for the YinYang grid the smallest zone size in azimuthal direction is larger by the ratio
compared to the spherical polar grid.
Table 2 gives the value of this ratio for grids of various angular resolution, and various computational domains. These numbers provide an estimate of the gain in computation time one can expect when using the YinYang grid instead of the spherical polar grid.
Table 2: Expected gain factor when using the YinYang grid.
However, the gain factor calculated from the relative grid
spacings
does not determine the gain in the size of the time step, as the latter
is given in a more complicated way by the CFL condition
where C, v_{r}, , , and c_{s} are the Courant factor, the flow velocities in radial, colatitude and azimuthal direction, and the local sound speed, respectively. The CFL condition shows that the increase in the size of the CFL time step is somewhat smaller than implied by the gain factor resulting from the ratio of the sizes of the smallest zones of the YinYang grid and the spherical polar grid. In addition, the increase of the time step is problem dependent.
Besides the performance gain due to the increased size of the
CFL time
step, the YinYang grid also requires less computational zones to
cover the full sphere, and thus less computational time. For an
angular resolution
the spherical polar grid needs
zones to cover the full sphere, while the YinYang grid requires only
zones. Hence, up to 25% fewer computational zones are required. The gain depends only weakly on angular resolution and is problem independent.
However, employing the YinYang grid also requires some extra amount of computation compared to the spherical polar grid (see Sect. 3). In the following we only consider the extra costs of calculations during the actual simulation, but not the extra costs arising during the initialization, since these are negligible. We emphasize again that there are two major extra sets of calculations necessary when applying the YinYang grid. The first set concerns the interpolation of the ghost zone values that are needed for the communication between the Yin and Yang grid patches. The second set arises from the interpolation of the density onto the auxiliary spherical polar grid grid and the interpolation of the gravitational potential back from the auxiliary grid onto the YinYang grid. Exploiting the algorithms described in Sect. 3, the computational cost for both parts is almost negligible compared to the total computing time. Interpolation of the ghost zone values requires only of the total computing time per cycle in simulations with selfgravity, while the interpolation of density and gravitational potential performed within the gravity solver accounts for of the computing time needed for the gravity solver. This corresponds to approximately of the computing time per cycle.
To obtain actual numbers for the gain, we performed several timing tests including simulations with and without selfgravity using four different grid resolutions. The tests were carried on an IBM Power6 using a single processor. According to these tests the computing time per cycle for the YinYang grid averaged over five cycles is approximately and smaller than for the spherical polar grid for simulations without selfgravity and with and angular resolution, respectively. For simulations including selfgravity, the gain factor decreases by approximately.
Concerning the gain from the less restrictive CFL condition, we consider the case of the rotating polytrope since the size of the time step does not vary much throughout the simulation. For an angular resolution of , we find a gain of approximately a factor of 63 when using the same Courant number both for the YinYang grid and the spherical polar grid.
7 Conclusion
A twopatch overset grid in spherical coordinates called the ``YinYang'' grid is successfully implemented into our 3D Eulerian explicit hydrodynamics code, PROMETHEUS, including in particular the treatment of selfgravitating flows. The YinYang grid eases the severe restriction of the time step size in the polar regions of the sphere, because each YinYang grid patch contains only the lowlatitude part of the usual spherical polar grid. From our experiences, the implementation steps are easy and straightforward for a hydrodynamics code using directional splitting and having 3D spherical polar coordinates already implemented due to the simplicity of the YinYang transformation and its symmetry property. Basically it involves doubling the state variable arrays, calling the 1D core hydrodynamics solver in angular directions for both the Yin and the Yang grid, and adding a subroutine that handles the YinYang transformation and the interpolation of variables between both grids.
We validated the code with several standard hydrodynamic tests. The test results show good agreement with analytic solutions if these are available. Furthermore, as demonstrated by three of our test problems  a planar shock crossing the YinYang grid boundary, an offcenter TaylorSedov explosion involving mass, momentum (all three components) and energy flux across the YinYang grid boundary, and a polytrope whose rotation leads to considerable and systematic mass, momentum (only angular components) and energy flux across the YinYang grid boundary  the YinYang grid does not introduce any numerical artifact at the internal YinYang boundary. The tests also confirm that the numerical solutions obtained with the YinYang grid do not show any evidence of a preferred radial direction, as it eliminates numerical axes artifacts which seriously flaw the flow near the coordinate symmetry axis when using a spherical polar grid. Besides successfully simulating a TaylorSedov explosion and selfgravitating (rotating and nonrotating) equilibrium polytropes the code has also passed another astrophysically relevant test involving the growth of RayleighTaylor instabilities.
Because the communication between the two grid patches involves interpolation, flows across the YinYang boundary cause some small amount of nonconservation of conserved quantities. However, even for the (in this respect) severe test case of the rotating polytrope involving large flows across the YinYang boundary, we observe only a small amount (less than one percent) of nonconservation.
The YinYang grid offers a large gain in computing time arising from two sources. Firstly, the number of computational zones needed is reduced by approximately depending on the angular resolution. This gain reduces the computing time per cycle and is problem independent. Secondly, the size of the CFL time step is considerably enhanced, because the polar regions with converging meridional coordinate lines are not present in case of the YinYang grid. The corresponding gain in time step size highly depends on the problem simulated. The extra costsfor interpolation between the two grid patches and the interpolation performed in the gravity solver are negligible compared to the gain in the time step size.
In conclusion, our implementation of the YinYang grid into the multidimensional hydrodynamics code PROMETHEUS brings about the possibility to simulate three dimensional selfgravitating hydrodynamic flows in spherical coordinates which, in most cases, have been computationally inaccessible up to now due to the prohibitively large computational costs. With the possibility to add more physics such as neutrino transport (work in progress), the new code version can be used to carry out, e.g., core collapse supernova simulations in 3D.
AcknowledgementsThis research was supported by the Deutsche Forschungsgemeinschaft through the Transregional Collaborative Research Centers SFB/TR 27 ``Neutrinos and Beyond'' and SFB/TR 7 ``Gravitational Wave Astronomy'', and the Cluster of Excellence EXC 153 ``Origin and Structure of the Universe''. The simulations were performed at the Rechenzentrum Garching (RZG) of the MaxPlanckSociety. We would like to thank Prof. A. Kageyama for sharing with us his YinYang interpolation routine for scalar/vector fields, and the anonymous referee for his/her useful and supportive comments.
References
 Blondin, J. M., & Mezzacappa, A. 2006, ApJ, 642, 401 [NASA ADS] [CrossRef] [Google Scholar]
 Chandrasekhar, S. 1969, Ellipsoidal figures of equilibrium (Yale Univ. Press) [Google Scholar]
 Chesshire, G., & Henshaw, W. 1994, SIAM J. Sci. Comput., 15, 819 [CrossRef] [Google Scholar]
 Collela, P., & Woodward, P. R. 1984, J. Comput. Phys., 54, 174 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Eriguchi, Y., & Müller, E. 1985, A&A, 147, 161 [NASA ADS] [Google Scholar]
 Fragile, P. C., Lindner, C. C., Anninos, P., & Salmonson, J. D. 2009, ApJ, 691, 482 [NASA ADS] [CrossRef] [Google Scholar]
 Iwakami, W., Kotake, K., Ohnishi, N., Yamada, S., & Sawada, K. 2008, ApJ, 678, 1207 [NASA ADS] [CrossRef] [Google Scholar]
 Janka, H.T., Scheck, L., Kifonidis, K., Müller, E., & Plewa, T. 2005, in The Fate of the Most Massive Stars, ed. R. Humphreys, & K. Stanek, ASP Conf. Ser., 332, 363 [Google Scholar]
 Kageyama, A., & Sato, T. 2004, Geochem. Geophys. Geosyst., 5 [Google Scholar]
 Kifonidis, K., Plewa, T., Janka, H.T., & Müller, E. 2003, A&A, 408, 621 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Koldoba, A. V., Romanova, M. M., Ustyugova, G. V., & Lovelace, R. V. E. 2002, ApJ, 576, L53 [NASA ADS] [CrossRef] [Google Scholar]
 Landau, L., & Lifshitz, E. 1959, Fluid mechanics, 6 (Pergamon) [Google Scholar]
 Mezzacappa, A., Bruenn, S., Blondin, J. M., Hix, W., & Messer, O. 2006, in AIP Conf. Proc., 924, 234 [Google Scholar]
 Müller, E., & Steinmetz, M. 1995, Comput. Phys. Commun., 89, 45 [NASA ADS] [CrossRef] [Google Scholar]
 Ohnishi, N., Kotake, K., & Yamada, S. 2007, ApJ, 667, 375 [NASA ADS] [CrossRef] [Google Scholar]
 Peng, X., Xiao, F., & Takahashi, K. 2006, Q.J.R. Meteorol. Soc., 132, 979 [Google Scholar]
 Romanova, M. M., Ustyugova, G. V., Koldoba, A. V., Wick, J. V., & Lovelace, R. V. E. 2003, ApJ, 595, 1009 [NASA ADS] [CrossRef] [Google Scholar]
 Ronchi, C., Iacono, R., & Paolucci, P. S. 1996, J. Comput. Phys., 124, 93 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Scheck, L. 2006, Ph.D. Thesis, Technical University Munich [Google Scholar]
 Scheck, L., Kifonidis, K., Janka, H.T., & Müller, E. 2006, A&A, 457, 963 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sod, G. A. 1978, J. Comput. Phys., 27, 1 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Toro, E. F. 1997, Riemann solvers and numerical methods for fluid dynamics  a practical introduction (Springer) [Google Scholar]
 Wang, Z. 1995, J. Comput. Phys., 122, 96 [NASA ADS] [CrossRef] [Google Scholar]
 Wu, Z.N., Xu, S.S., Gao, B., & Zhuang, L.S. 2007, Comput. Fluids, 36, 1657 [CrossRef] [Google Scholar]
 Zink, B., Schnetter, E., & Tiglio, M. 2008, Phys. Rev. D, 77 [Google Scholar]
Footnotes
 ... Hammer^{}
 Present address: MaxPlanckInstitut für Plasmaphysik, Boltzmannstraße 2, 85748 Garching, Germany.
All Tables
Table 1: Mean errors in the gravitational acceleration.
Table 2: Expected gain factor when using the YinYang grid.
All Figures
Figure 1: An axisfree YinYang grid configuration plotted on a spherical surface. Both the Yin (red) and Yang (blue) grid are the low latitude part of the normal spherical polar grid and are identical in geometry. The Yang grid is obtained from the Yin grid by two rotations, and vice versa. 

Open with DEXTER  
In the text 
Figure 2: A Mercator projection of an overlap region of the YinYang grid. In case of bilinear interpolation, four neighboring values of the underlying grid (red) will be used to determine the zonecentered value of a ghost zone in the grid on top (blue). The interpolation coefficients are determined by the relative distances, denoted by black lines, between the interpolation point (diamond) and the four neighboring points (crosses). 

Open with DEXTER  
In the text 
Figure 3: One dimensional profiles of density , pressure p, velocity in xdirection v_{x}, and specific internal energy e are shown along the xdirection at z^{(n)}=0.25 cm and y^{(n)}=0 cm (dasheddotted line in Fig. 4) for the shock tube simulation at every 0.1 s. Open and filled symbols represent data points on the Yin and Yang grid, respectively. Solid lines give the distributions calculated with an exact Riemann solver. 

Open with DEXTER  
In the text 
Figure 4: Snapshot of density contours in the meridional plane at t=0.15 s for the shock tube test problem. Dashed lines mark the YinYang grid boundary, while the dotted circular curves represent the inner and outer radial boundary of the computational domain, respectively. The one dimensional profiles shown in Fig. 3 are resampled along the dasheddotted line at z^{(n)}=0.25 cm. 

Open with DEXTER  
In the text 
Figure 5: Color maps of the tangential velocity defined by in the meridional plane resulting from the (1D Cartesian) shock tube problem. The snapshots are computed using the YinYang grid ( left) and a standard spherical polar grid ( right) at a time t=0.15 s. On the left panel, red and blue lines mark the boundaries of the Yin and the Yang grid patches, respectively. On the right panel, the two red circles show the inner and outer boundary in the radial direction of the standard spherical polar grid. The labels at the color bars give the tangential velocity in units of cm/s. The color range is limited to 0.05 cm/s to emphasize the smallness of the tangential velocity far from the outer radial grid boundary. 

Open with DEXTER  
In the text 
Figure 6: Distributions of density ( top), pressure ( middle) and radial velocity ( bottom) versus radius from the explosion center (located at cm for the TaylorSedov explosion problem plotted at every 10^{11} s. Open symbols are data points from the Yin grid, while filled symbols represent sampled data from the Yang grid. The solid lines give the corresponding analytic solution. The data are resampled along the dasheddotted line shown in Fig. 7. 

Open with DEXTER  
In the text 
Figure 7: Lines of constant density in the meridional plane obtained from our simulation of a TaylorSedov explosion. The snapshot is taken at a simulation time s which corresponds to an explosion time s. The dashed lines mark the YinYang boundary, while the two dotted circles represent the inner and outer radial boundary of the computational domain, respectively. The data presented in Fig. 6 are resampled along the dasheddotted line. 

Open with DEXTER  
In the text 
Figure 8: Evolution of the mass within the overlap region for the TaylorSedov test case computed on the Yin grid, minus the mass computed on the Yang grid, , divided by the sum of these two masses. The dashed, dotted and solid lines give the solutions computed on a grid of zones, (i.e., angular resolution), zones (i.e., angular resolution), and zones (i.e., angular resolution), respectively. 

Open with DEXTER  
In the text 
Figure 9: Surfaces of constant density in 3D ( left) and 2D ( right; meridional cut at ) resulting from the simulation of the RayleighTaylor instability described in the text at t=2.85 s. Contour lines on the Yin grid are shown using the blueyellow colors while contour lines on the Yang grid are displayed using the whiteblack colors. 

Open with DEXTER  
In the text 
Figure 10: Position of the heads of the RTI bubbles versus time. Red symbols (circles, triangles, and squares) show data from the Yin grid, while blue symbols (diamonds) represent data on the Yang grid. 

Open with DEXTER  
In the text 
Figure 11: Contour lines of the gravitational potential ( left column) for two homogeneous selfgravitating configurations: a prolate spheroid ( top row) with an axis ratio of 0.7, and a sphere ( bottom row). The configurations are indicated by the darkgray shaded areas. Dashed lines show the YinYang boundary, while dotted lines indicate the outer radial boundary of the computational grid. The middle and right columns give the maximum and mean error of the numerically calculated gravitational potential for different grid resolutions as a function of the number of spherical harmonics used in our multipole gravity solver. The solid, dotted, dashed, and dasheddotted lines in both columns correspond to a grid resolution of zones, zones, zones, and zones, respectively. 

Open with DEXTER  
In the text 
Figure 12: Relative change of the central density of a nonrotating (nearly) equilibrium polytrope as a function of time. 

Open with DEXTER  
In the text 
Figure 13: Density ( top) and radial velocity ( bottom) of a nonrotating n=1 equilibrium polytrope as a function of radius after t=10 ms of ``evolution''. In the top panel, the solid line shows the initial density profile. Red circles and blue triangles correspond to data from the Yin and the Yang grid, respectively. 

Open with DEXTER  
In the text 
Figure 14: Same as Fig. 12 but for a rotating polytrope. The solid and dashed curves show the relative variation of the density along an equatorial ray ( ; ) and along the pole ( ), respectively. 

Open with DEXTER  
In the text 
Figure 15: Density ( upper panel) and radial velocity ( lower panel) of a n=1 rotating polytrope in stationary equilibrium as a function of radius after t=10 ms of ``evolution''. In both panels the solid and dashed lines show the profiles along an equatorial ray ( , ) and along the pole ( ), respectively. Red circles and blue triangles in the upper panel correspond to data from the Yin and the Yang grid, respectively. In the lower panel, we show in addition time averaged (over the interval t= [9, 10] ms) velocity profiles along the equatorial ray (dotted) and the pole (dasheddotted). 

Open with DEXTER  
In the text 
Figure 16: Evolution of the relative mass loss, (MM_{0})/M_{0}, where M_{0} is the initial total mass, for the TaylorSedov test simulated on three different (angular) grids with zones (i.e., angular resolution; dashed line), zones (i.e., angular resolution; dotted line), and zones (i.e., angular resolution; solid line), respectively. 

Open with DEXTER  
In the text 
Figure 17: Illustration of the YinYang grid overlap configuration, where PQRS is a grid zone at the boundary of the Yang grid (blue) which overlaps with the underlying grid zone ABCD of the Yin grid (red). Fluxes referring to the Yin and the Yang grid are denoted by f and g, respectively. 

Open with DEXTER  
In the text 
Copyright ESO 2010
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.