Numerical simulation of core convection by a multi-layer semi-implicit spherical spectral method

Journal of Computational Physics 230 (2011) 8698–8712

Contents lists available at SciVerse ScienceDirect

Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp

Numerical simulation of core convection by a multi-layer semi-implicit spherical spectral method Tao Cai a,⇑, Kwing L. Chan a, Licai Deng b a b

Department of Mathematics, The Hong Kong University of Science and Technology, Clear Water Bay, Hong Kong, China National Astronomical Observatories, Chinese Academy of Science, Beijing, China

a r t i c l e

i n f o

Article history: Received 23 March 2011 Received in revised form 17 August 2011 Accepted 18 August 2011 Available online 5 September 2011 Keywords: Spherical spectral method Multi-layer scheme Turbulence Core convection

a b s t r a c t A semi-implicit multi-layer spherical spectral method for simulating stellar core convection is described. The fully compressible three-dimensional hydrodynamic equations with rotation and energy generation are solved. Prognostic variables are expressed as finite sums of spherical harmonics in the horizontal directions and handled by the finite difference method in the radial direction. The stratified approximation is used to simplify the nonlinearity to quadratic. A multi-layer scheme is employed to overcome the time step problem arising from shrinking grid sizes in the physical space near the center of the star. Despite of the different spectral truncations in different layers, round-off conservation of the total mass and total angular momentum of the whole domain can be maintained, and were confirmed numerically. The code is parallelized; with 12 processors the speedup factor is about 9. The solutions of model core convection with and without rotation are discussed. Ó 2011 Elsevier Inc. All rights reserved.

1. Introduction Core convection is important for the structure and evolution of stars. The energy generated by the nuclear burning core is transported efficiently by convection, which makes the temperature gradient close to adiabatic. Theoretical description of convection faces difficulties due to the nonlinearity and complexity of the turbulent motion. In stellar evolution calculations, the turbulent effects are usually treated by the one dimensional mixing length assumption. The basic idea of mixing length theory is that the fluid elements travel over a distance on the order of a pressure scale height Hp, and then blend into the surroundings. The mixing length approach cannot be used at the center of a star as the pressure scale height becomes infinite when the radius r approaches zero. Another drawback is that it has difficulty in describing the overshooting process. Convective fluid elements generally do not stop abruptly on the boundary between convectively unstable and stable regions. Instead, they penetrate into the convectively stable region and then bounce back. Fresh nuclear fuel can be supplied to the burning core through the mixing process, thereby changing the lifetimes of stars in different evolutionary phases. By taking account of the nonlocal effects of convection in stellar cores, Maeder [1] found that as stellar mass increases the effect of overshooting on stellar evolution becomes more significant. In his work, the overshooting distance was assumed to be a fraction of the pressure scale height lov = aovHp. However, determination of the overshooting parameter aov is still a challenging problem. An empirical value of aov for massive stars is 0.25. For low mass stars the overshooting distance needs to be tuned smaller because of the small convective core. Using the same value of aov would predict an unrealistically large overshooting distance. For intermediate and low mass stars aov is usually chosen to be 0.125 and 0, respectively [2]. More fundamentally, ⇑ Corresponding author. Tel.: +852 23587453; fax: +852 23581643. E-mail addresses: [email protected] (T. Cai), [email protected] (K.L. Chan), [email protected] (L. Deng). 0021-9991/$ - see front matter Ó 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.jcp.2011.08.014

T. Cai et al. / Journal of Computational Physics 230 (2011) 8698–8712

8699

whether or not overshooting is non-negligible is still under debate. Although quite a number of studies agree that core overshooting is needed to explain observation [3,2,4,5], Stothers and Chin [6] argued that when the new opacity tables of Rogers and Iglesias [7] are used the extent of overshooting is small. Beyond the standard mixing length theory, Xiong [8] pointed out that convection should be treated by a global non-local approach. He developed a time-dependent Reynolds stress method based on taking the horizontal averages of the hydrodynamic equations. A similar approach was adopted by Canuto and Dubovikov [9], and Kupka [10]. The Reynolds stress method is more physical than the mixing length theory, but there are still uncertain parameters left to be calibrated. In addition, one dimensional convection theory cannot effectively estimate the effects introduced by differential rotation. Numerical simulations can offer more consistent information about convection. The finite-difference method and the spectral method are two popular schemes used to solve the nonlinear hydrodynamic equations. By focusing on a local region, the finite-difference method can afford higher resolution, however, it has difficulty with the poles of a spherical shell as the small grid size there makes the time steps extremely small. An alternative method is the spectral approach, but the operation count of the interaction coefficients grows very fast with resolution. The development of the spectral transform method [11] made the spectral approach more affordable and became competitive to the finite difference approach. The spectral approach avoids the pole problem and in principle has a very high order of accuracy. For moderate horizontal resolutions, it can be more efficient than the finite-difference approach. The spectral approach is widely used for carrying out climate simulations [12–15] in which the hydrostatic approximation can be used and aliasing can be avoided in transforms. For stellar convection problems, the situation, however, is different. Length scales of motions are comparable to the pressure scale height so that the quasi-hydrostatic approximation is inappropriate. Navier–Stokes equations should be solved directly, but the time step is strongly restricted by the CFL condition associated with waves. By setting the time derivative of density to zero (i.e. @ tq = 0), the anelastic approximation introduced by Ogura and Charney [16] filters out sound waves and thus avoids the associated CFL condition. Large time steps can be used without the hydrostatic approximation. When there is no background density stratification, the anelastic approximation degenerates to the Boussinesq approximation. Scale analysis of the anelastic equations by Ogura and Phillips [17] showed that the intrinsic error is related to the entropy gradient. To reduce the error, modification of the momentum or energy equation is usually used. Ogura and Phillips [17] suggested modifying the momentum equation by taking the isentropic state as the reference state. Lipps and Hemler [18] concluded that the anelastic approximation is consistent if the potential temperature is a slowly varying function. In general, however, radiative transport processes are important and the potential temperature may vary substantially. Gough [19] derived the anelastic approximation by taking into account the radiative processes. By his suggestion, the reference state is set as the horizontal average which is a function of the radius. Latour et al. [20] realized that the reference state may vary when the layer thermally evolves and thus should be treated as a function of time and radius. Combining the previous results, Gilman and Glatzmaier [21] adopted a reference state which is nearly adiabatic and a function of radius only. Pioneered by the work of Glatzmaier [22–24], the anelastic approach was successfully applied to study stellar convective dynamos with the inclusion of the magnetic fields. The anelastic spherical harmonic (ASH) code was then parallelized by Clune et al. [25] and a series of works about stellar convection were published [26–30]. As mentioned above, the anelastic approximation eliminates sound waves, but @ tq = 0 also affects the description of gravity waves. Such elimination and inaccuracy may not be desirable in certain applications. As these waves are linear, CFL conditions arising from them can be overcome by an implicit treatment of the linear terms of the hydrodynamic equations in spectral space. The nonlinear terms can be handled explicitly in real space with the transform method. If the Mach number (Ma) is sufficiently small, the CFL condition associated with the fluid motions would not be so restrictive. Based on such observation, Chan et al. [31] developed the semi-implicit stratified numerical spectral model (NSM) for deep atmospheres. The stratified approximation, which replaces the density in advection terms by its horizontally averaged value, was introduced to limit the nonlinearity to quadric. This facilitates non-alias transformation between the spectral and real spaces (for suppression of nonlinear instability). This approximation is valid for Ma  1 (accurate up to Ma2), and one important property is that the resulting equations preserve mass and angular momentum conservations (confirmed both analytically and numerically). For the numerical simulation of convection, the center of the sphere is usually excluded from the computational domain for numerical reasons, namely, the pseudo singularity at r = 0 and the very small time steps associated with the decreasing mesh sizes [28]. The effects of an artificially introduced solid core on the results are not yet clear. Using a two-dimensional finite volume approach, Evonuk and Glatzmaier [32] investigated the effects of solid cores in giant planets. They found that it is difficult to distinguish between a giant planet without a solid core and one with a 10% radius solid core, but the conclusion is different when the core size increases to 35% radius. As two- and three-dimensional turbulence behave differently [33], such result cannot be generalized to a three dimensional situation directly. To simulate core convection, an approach that can reach the center of the sphere is needed. In this paper, we describe a spherical spectral approach that can achieve this objective. On the basis of the work of Chan et al. [31], a multi-layer semi-implicit stratified spectral method is developed. Starting with very low order harmonic expansions at the center, successively higher order harmonics are used in the outer layers. The time-step restrictions are thus kept quite uniform across the layers and the singularity problem at the center is well handled.

8700

T. Cai et al. / Journal of Computational Physics 230 (2011) 8698–8712

2. The numerical model Our objective is to simulate stellar core convection, thus there are a number of requirements to consider. First, it is necessary to solve the fully compressible three-dimensional nonlinear hydrodynamic equations with rotation and energy generation by nuclear reactions. Second, the computational domain should include the center of the star. Third, the time steps cannot be too small as the time period required for thermal and dynamical relaxation is generally very long. Fourth, even when approximations are applied, the resulting numerical scheme should be consistent with the conservation of certain model characterizing quantities (in particular total mass and total angular momentum). This is to ensure that these parameters do not drift substantially when the computation becomes long. 2.1. The equations Following the work of [31], the linear wave terms of the hydrodynamic equations stay intact. Temporal marching of the linear terms can be treated implicitly in the spectral space to avoid the restrictive CFL conditions associated with acoustic and gravity waves. The stratified approximation is applied to the nonlinear terms; it facilitates the transformation of these terms between the spectral and physical spaces. The contribution of these terms to the temporal evolution of the prognostic variables are handled by an explicit method in the physical space. By the fact that the distributions of thermodynamic variables are dominated by the stratification, a thermodynamic variable can be expressed as the sum of a mean value (horizontal average) and a small term representing its horizontal variation. Both terms are time-dependent. For example, the density can be expressed as

qðr; h; /; tÞ ¼ q0 ðr; tÞ þ q1 ðr; h; /; tÞ;

ð1Þ

where r is the radial distance from the center; t is the time; h and / stand for the colatitude and the azimuthal angle, respectively; the subscripts 0 and 1 stand for the horizontal average and the deviation from it, respectively. The following assumptions are made: (i) In the momentum equation, the horizontal variation of density is ignored in the nonlinear advection terms and viscous terms. (ii) In the energy equation, terms containing products of two or more horizontal variations of the thermodynamic variables are ignored. In terms of the Mach number, the formal errors are of order higher than O(Ma2) [19,34,35]. The resulting equations have only quadratic nonlinearity, thereby alias-free transformations can be easily implemented. These assumptions are also used by the anelastic approximation. The stratified approximation, however, does not need to suppress the temporal variation of density in the continuity equation which introduces a relative error @ t q=ðq$  VÞ  OðMaÞ when sound waves are significant. The equations of the semi-implicit stratified spherical spectral method for a spherical shell are given in [31]. Here we recast these equations in a form that facilitates implementation of boundary conditions at the center of the sphere. The reduced fully compressible hydrodynamic equations under the above assumptions can be written as

~ ¼ Wq @t q fr ¼ W r þ Dr þ C r þ Nr @t M @ t ~d ¼ W d þ Dd þ C d þ Nd @ t ~f ¼ þDn þ C f þ Nf ~ ¼ W p þ Dp þ Np @t p where @ t is the partial derivative with respect to time t, and W, D, C, and N represent the wave, diffusion, Coriolis, and nonlinear terms, respectively. The wave terms are given by:

fr Þ  r 2 ~d W q ¼ @ k ðw1 M 2

2

ð2Þ

1

~  q ~g W r ¼ r @ k ½ðr wÞ p 2

ð3Þ

2 ~ Hp

W d ¼ r r

ð4Þ

fr w @ k ðlog p  C0 log q Þ þ C0 @ k ðw M fr Þ þ C0 r d; W p ¼ p0 =q0 ½ M 0 0 1

1

2 ~

ð5Þ

where @ k is the partial derivative with respect to the radial grid coordinate k. The prognostic variables are fr ¼ M r r 2 w; ~ ~ ¼ pr 2 w  w ¼ dr=dk is the transformation factor between the grid and q~ ¼ qr2 w; M d ¼ dr 4 w; ~f ¼ fr 4 w, and p the radial coordinates; q is the density; p is the pressure; M = qV is the momentum density; V is the velocity; Mr is the radial component of M; MH is the horizontal component of M; s = sinh; d = rH  MH = (rs)1(@ hsMh + @ /M/) is the horizontal divergence of MH; f = rH  MH = (rs)1(@ h(sM/)  @ /Mh) is the curl of MH along the radial direction; r2H ¼ ðrsÞ2 ðs@ h ðs@ h Þ þ @ 2/ Þ is the horizontal Laplacian; g is the magnitude of gravitational acceleration; C = (dlnp/dlnq)ad is the adiabatic exponent; the tilde variable f M is defined as M r2w where r2w is the volume factor. The diffusion terms are:

T. Cai et al. / Journal of Computational Physics 230 (2011) 8698–8712

fr Þ þ l r 2 @ k ðr 4 w1 dfV Þ Dr ¼ 2@ k ½lr r 2 w1 @ k ðr 2 w1 V r þ lH ð r

2f HVr



2r 3 df V



fr Þ 4r 2 V

2

2

1

þ r @ k ½kr w

8701

ð6Þ fr Þ @ k ðw1 V

þ kr

4

w1 dfV 

ð7Þ

2f 2f 2 1 2 f 2 f fr Þ þ r2 dfV  dV Þ þ k½@ k ðr 2 w1 r2H V Dd ¼ @ k ½lr r4 w1 @ k ðr 4 w1 df V Þ þ lH r w rH V r  þ 2lH ðrH dV þ r rH V r þ r H

ð8Þ

fV Þ þ lH ðr2H f fV þ 2r2 f fV Þ Df ¼ @ k ½lr r4 w1 @ k ðr 4 w1 f

ð9Þ

~Þ þ jqr @ k ðw1 r 2 q ~ Þg þ rad CðjpH r2H p ~ þ jqH r2H q ~ Þ; Dp ¼ rad C@ k fr 2 w1 ½jpr @ k ðw1 r2 p

ð10Þ

f f f ~ f ~ where d V ¼ d=q0 ; fV ¼ f=q0 ; and V r ¼ M r =q0 ; lr and lH are the radial and horizontal viscosities respectively; k is the second coefficient of viscosity; rad = (dlnT/dlnp)ad is the adiabatic temperature gradient; j is the coefficient of conduction. The Coriolis terms are:

C r ¼ 2Xs g M/

ð11Þ

fr Þ C d ¼ 2X½v~f  rðs g M/ þ @ / M

ð12Þ

fh  s@ h M fr  2v M fr Þ; C f ¼ 2X½v~d  rðs M

ð13Þ

where X is angular velocity and v = cosh. The nonlinear terms are:

fr V r Þ  $H  ð M g fr VH Þ þ r 1 M Nr ¼ @ k ðw1 M H VH

ð14Þ

fr VH Þ þ $H  ð~fVH Þ  $H  ð~dVH Þ  r 2 r2 ð M g Nd ¼ @ k ½r 2 w1 $H  ð M H V H Þ=2 H

ð15Þ

fr VH Þ  $H  ð~fVH Þ  $H  ð~dVH Þ Nf ¼ @ k ½r 2 w1 $H  ð M

ð16Þ

1 2 2 fr Þ þ r2 df f1 V=q0  þ ðC0  1Þp0 $  ð q f1 V=q0 Þ p1 ½@ k ðw1 V Np ¼  r 2 w$  ½ðr 2 wÞ1 f p1 V  ðC0  1Þ f V  þ r w$  ½ðr wÞ p0 q

þ r 2 wrad Ced þ Q;

ð17Þ

where Q is the energy source term. ed is the dissipation function

ed ¼ 2flr ð@ r V r Þ2 þ lH r2 ½ð@ h V h þ V r Þ2 þ ðs1 @ / V / þ V r þ V h cot hÞ2 g þ lH ½ðrsÞ1 @ / V h þ r1 s@ h ðV / =sÞ2 þ ½lH ðrsÞ1 @ / V r þ lr r@ r ðV / =rÞ  ½ðrsÞ1 @ / V r þ r@ r ðV / =rÞ þ ½lr r@ r ðV h =rÞ þ lH r 1 @ h V r   ½r@ r ðV h =rÞ þ r 1 @ h V r  þ kð$  VÞ2 :

ð18Þ

fh and g ~ 0 . The two diagnostic variables M In the above expressions (Eqs. (14)–(18)), V can be approximated by f M=q M / are evaluated by the formulas:

fh Þm ¼ r 1 fcðn; mÞ~dm þ cðn þ 1; mÞ~dm þ im½nðn þ 1Þ1~fm g ðs M n1 nþ1 n n

ð19Þ

1 ~m 1 ~m ~m ðs g M / Þm n ¼ r fcðn; mÞfn1 þ cðn þ 1; mÞfnþ1 þ im½nðn þ 1Þ dn g

ð20Þ

where the function c(n, m) is defined as

cðn; mÞ ¼ n1 ½ðn2  m2 Þ=ð4n2  1Þ1=2 :

ð21Þ

The Coriolis terms can be put in the form 1 ~m 1 ~m ~m Cm rn ¼ 2Xr fcðn;mÞfn1 þ cðn þ 1;mÞfnþ1  im½ðnðn þ 1ÞÞ dn g e m rg C m ¼ 2Xfðn þ 1Þcðn;mÞ~fm þ ncðn þ 1;mÞ~fm þ im½nðn þ 1Þ1 ~dm  im M dn

Cm fn

n1

nþ1

n

rn

1 ~m em em ~m ¼ 2Xfðn þ 1Þcðn;mÞ~dm n1  ncðn þ 1;mÞdnþ1 þ im½nðn þ 1Þ fn þ nðn þ 1Þ½cðn;mÞ M rn1 þ cðn þ 1;mÞ M rnþ1 rg:

ð22Þ ð23Þ ð24Þ

For simplicity of notation, we use y to represent the column array of prognostic variables. The hydrodynamic equations can be written in the form:

@ t y ¼ ðW þ D þ CÞy þ N½y:

ð25Þ

Note that the first three terms of this equation are linear. They can be most conveniently handled in spectral space. It is only necessary to apply the transform method to the last term.

8702

T. Cai et al. / Journal of Computational Physics 230 (2011) 8698–8712

2.2. Numerical method Since all the prognostic variables are scalars, they can be approximated by truncated series of spherical harmonics:

y¼

X

m im/ ym n ðrÞP n ðvÞe

ð26Þ

Pm n

where are the normalized associated Legendre functions. We adopt a triangular truncation of the n, m indices with 4(nmax) = {(n, m)j0 6 m 6 n, 0 6 n 6 nmax} where nmax is the maximum degree used in the truncated expansion. The operator im/ im/ r2H that appears in the linear terms can be easily handled since r2H pm ¼ nðn þ 1Þr2 pm . In the radial direction, n ðvÞe n ðvÞe ~ are located on the k ~; ~ the second order finite difference method is used for the partial derivatives with a staggered grid. q d; ~f; p fr is located on the k + 1/2 levels. In the horizontal direction, the grid points on each latitude circle are equally levels while M spaced, while Gaussian grids are used on the longitude circles. There are two kinds of singular points associated with the spherical coordinate system. One is on the poles, the other is at the center. Gaussian grids do not contain any point on the poles, therefore the spherical harmonics approach is not affected by the polar singularity. The singular problem caused by the center can be avoided as the tilde variables (as well as r2lr. . . etc) all go to zero at the center. The steps for computing the nonlinear terms consist of the following: (i) Transform the variables from the spectral domain to the physical domain by synthesizing the truncated spherical harmonics series; (ii) Evaluate the nonlinear terms at the grid points of the physical space; (iii) Backward transform the nonlinear terms from the physical domain to the spectral domain. To ensure non-alias transformation of the quadric nonlinear terms. The number of grid points in the latitudinal (Imax) and azimuthal (Jmax) directions need to satisfy

Imax P ð3nmax þ 1Þ=2

ð27Þ

J max P 3nmax þ 1

ð28Þ

To facilitate the fast Fourier transformation, Jmax is chosen to be a power of two. In order to keep the grid spacing more or less uniform (to avoid severe restriction of time steps by the CFL condition), the domain of computation is decomposed into a number of layers over different radial ranges. The maximum degree of spherical harmonics nmax varies from layer to layer and decreases from outside to inside. As the radius shrinks, the number of horizontal grid points also decreases. As a result the grid spacings do not vary by a big factor, and the time steps need not become very small. A semi-implicit method is adopted for time marching. The contributions from the Coriolis and nonlinear terms 4yC+N are evaluated by an explicit second-order method, and then an implicit step solving the following equation is performed.

ðyqþ1  yq Þ ¼ 4tbðW þ DÞðyqþ1  yq Þ þ 4tðW þ DÞyq þ 4yCþN þ ðb  1=2ÞOð4t 2 Þ þ Oð4t 3 Þ

ð29Þ

where b is a parameter that describes the degree of implicitness. Second order accuracy in time can be obtained by choosing b = 0.5, but in real calculations it is chosen to be a little larger than 0.5 to ensure numerical stability. The linear terms W and D can be treated implicitly in a simple way as they couple vertical levels but do not couple different (n, m) modes. The C term couples the (n, m) mode to the (n ± 1, m) modes. In a multi-layer situation where nmax varies across layers, the solution of the implicit equation requires some special arrangements. For example, let us consider a three-layer case. Suppose that the interfaces of the layers are at the vertical levels k1 and k2. In the region 1 6 k 6 k1, the maximum spectral degree is n1max; in the region k1 + 1 6 k 6 k2, the maximum spectral degree is n2max; and in the region k2 + 1 6 k 6 kmax, the maximum spectral degree is n3max. The spectral degrees satisfy n1max < n2max < n3max. The implicit equation above is split into three matrix equations by decomposing the harmonics-radius domain into three regions:

B1 : 4ðn1max Þ  ½1; kmax 

ð30Þ

B2 : ½4ðn2max Þ  4ðn1max Þ  ½k1 þ 1; kmax  B3 : ½4ðn3max Þ  4ðn2max Þ  ½k2 þ 1; kmax :

ð31Þ ð32Þ

Only the smallest spectral triangle 4(n1max) penetrates all the layers (B1). Each successive region only contains the additional spherical harmonics introduced across a layer interface. An illustration of the structure of the populated region in the (n, m, r) space is shown in Fig. 1. Note that the view is from a location below the nm-plane. For the calculations presented here, the impenetrable and stress-free boundary conditions are applied at the top. At the center, regularity implies that all the prognostic variables and the r2 flux terms vanish. The boundary conditions of the B regions are implemented as follows. At k = kmax, the impenetrable and stress-free boundary conditions are used. At the lower 4 1 f fr ; @ k ðr 4 w1 d f ~; M ~; p interface of a layer (k1 or k2), the spectral components of q w fV Þ of the layer above are set to V Þ, and @ k ðr zero. This is compatible with the increase of truncation from upper to lower layers. Higher resolution components of an upper layer do not interact with their counter parts that do not exist in the lower layer. With these boundary conditions, the total relative angular momentum can be shown conserved as discussed by Chan et al. [31]. The difference is that the first R term on the right hand side of Eq. (18) of that paper becomes @ k ½lr r4 w1 @ k ðr 4 w1 f fV Þdkd- where d- = sdh d/; the k integral is split into three parts by the three layers. The boundary conditions adopted here ensure that this integral vanishes. Total mass conservation can be similarly proven.

T. Cai et al. / Journal of Computational Physics 230 (2011) 8698–8712

8703

Fig. 1. Illustration of the multi-layer mesh in the spectral-radius domain.

2.3. Parallelized computation The code is parallelized with the OpenMP language. OpenMP achieves parallelization by splitting do loops into the multithreading processors without major modification. Thus a serial code can be easily converted to take advantage of a multiprocessor environment. Another advantage of OpenMP is that it uses shared memory, hence the cumbersome task of programming distributed memories can be avoided. To maximize the performance of all the processors, it is necessary to distribute the tasks into each processor evenly. Details of the scheme are described as follows. In solving the block tridiagonal matrices associated with the W + D terms, the spectral domain is split among processors. As the W and D terms do not couple harmonic components with different n or m values, the separate components can be solved independently in the slave processors. As described in Eqs. (30)–(32), the spectral domain is decomposed into several subdomains. In the program, each subdomain is divided into Npro (number of processors) smaller blocks. Computation of the non-linear terms in physical space uses a different workload distribution method. Each layer is divided radially into Npro blocks and distributed among the processors. The information needed across blocks or layers are passed to the neighbors through the buffer or extension layers. The C terms couple harmonic components with different n values. The computation is split the same way as for the N terms, but without transformation to the physical space. The workload of each processor is balanced by the following method. Consider for example the computation of N terms for the layer region [k1 + 1, k2]. Let Nsize be the maximum integer satisfying the condition NsizeNpro 6 k2  k1. If the two sides are exactly equal, Nsize of the k-levels are allocated into each processor. Otherwise Nsize + 1 of them are allocated to each of (Npro  1) processors, and the remaining levels are then put into the remaining processor. Therefore each processor compute the N terms for up to Nsize + 1 levels. The strategy for the allocation of W and D calculations is similar, except that the division is made in the (m, n) domain. The total clock time of the parallelized computation can be approximated by the formula ttot = ts + tp/Npro. Here ts and tp are the CPU times of the serial part and the parallelized part of the code, respectively. The layer configurations and numerical resolutions of four test cases are summarized in Table 1 (the physical setup is discussed in Section 3). The total number of vertical grids is kmax = 101 for all the cases. The test is done on a 12-processor Intel Xeon X5670 computer. The clock times versus the inverse number of processors are plotted in Fig. 2. The computational speed increases almost linearly with the number of processors. The data for Cases A and B can be fitted by the functions ttot = 0.0673 + 1.957/Npro and ttot = 0.364 + 12.012/Npro, respectively. The parallelized part of the code is about 97% of the total. This high ratio means that parallelization can reduce the clock time effectively. The speed up factor is about 9 when 12 processors are used.

8704

T. Cai et al. / Journal of Computational Physics 230 (2011) 8698–8712

Table 1 Configurations of the test cases. Case Case Case Case Case

A B C D

Radial intervals

nmax values

Radial grid numbers

Angular velocity

(0, 0.1)/(0.1, 0.2)/(0.2, 1) (0, 0.1)/(0.1, 0.2)/(0.2, 1) (0, 0.1)/(0.1, 0.2)/(0.2, 1) (0, 0.1)/(0.1, 0.2)/(0.2, 1)

4/10/42 4/10/84 4/10/170 4/10/84

10/10/81 10/10/81 10/10/81 10/10/81

0 0 0 0.1

14

12

time in seconds

10

8

6

4

2

0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

1/N pro Fig. 2. Total CPU times of Case A (square) and Case B (cross) versus the reciprocal number of processors.

Compared with Case A, the time used by Case B is 6.1 times longer. Case B has doubled nmax value at the outermost layer. The complexity of the fast Fourier transform (FFT) is n2max log nmax , while for the Legendre transform (LT) it is n3max . When nmax is doubled, the complexities of FFT and the LT are increased by factors of 4.74 and 8, respectively. To better understand the efficiency of the code, we have measured the CPU times spent in LT and FFT (both forward and backward) in a single step with one processor in Cases A and B. The times of LT are 0.98 s and 7.55 s, while those of FFT are 0.41 s and 1.83 s, respectively. Compared with Case A, LT of Case B takes 7.7 times longer, while FFT takes 4.4 times longer. The results are close to the theoretical analysis. In terms of operation counts, Cases A and B take 1.1  109 and 6.7  109 operations per time step, respectively. Despite of the semi-implicit advantage (suppression of the wave CFL condition for large time steps), the spherical spectral method is competitive only when nmax is not too large. To achieve much higher resolutions, it is anticipated that distributed memory computers [25,36] as well as other algorithms such as the finite difference method [37,38] or the spectral element method [39,40] need to be used. 3. Numerical models of 3D core convection 3.1. The model setup We use the one-dimensional hydrostatic structure of a self-gravitating idealized star as the initial background state of the test cases. The initial distribution of the gas is polytropic and in a hydrostatic state with

T=T t ¼ g

ð33Þ

q=qt ¼ gnp p=pt ¼ gnp þ1

ð34Þ ð35Þ

T. Cai et al. / Journal of Computational Physics 230 (2011) 8698–8712

8705

where np is the polytropic index. The temperature distribution function g can be derived form the Lane–Emden equation for the self-gravitating spherically symmetric polytropic fluid:

1 d 2 dH ðn Þ þ Hnp ¼ 0 n dn dn H ¼ 1 at n ¼ 0

ð36Þ ð37Þ

dH ¼ 0 at n ¼ 0 dn H ¼ 1=ðZ þ 1Þ at n ¼ nt

ð38Þ ð39Þ 1=2

where n ¼ r½4pGq2b =ððn þ 1Þpb Þ and H = T/Tb are the scaled radius and temperature variables, respectively; G is the gravitational constant; Z = Tb/Tt  1 describes the depth of stratification. The subscripts b and t denote values at the bottom and the top respectively. nt is the scaled radius at the top and is dependent on Z. By solving the Lane–Emden equation, the temperature distribution can be found as

g ¼ ZH

ð40Þ

n=nt ¼ r=r t

ð41Þ

For simplicity of notation and computation, all the quantities are scaled by values at the top so that Tt = qt = pt = rt = 1. The computational domain is a sphere with the range of radius 0 6 r 6 1. The initial gravitational acceleration is g = (np + 1)dg/ R dr. Although the gravity can be updated by the formula 4pG qr2 dr=r2 , here it is set fixed in time for simplicity. This approximation is valid if the perturbation in the gravitational potential is small. The condition is satisfied if the adjustment of mean density distribution and the density fluctuations are small. The later condition is generally satisfied as the relative density fluctuation is on the order of Mach number square [31]. The first condition depends on how close the initial density distribution approximates the final relaxed state. It can be managed by choosing the appropriate polytropic index. The fluid is taken to be a perfect ideal gas with specific heat C = 5/3. The viscosity is isotropic. The value of m is 1  104, and the second coefficient of viscosity k is 2/3l. Energy is supplied by the heat generation term

e ¼ ke T 4

ð42Þ 7

where ke is taken to be a constant with the value 10

1 d 2 ðr f Þ ¼ e: r2 dr

. In an equilibrium situation, the energy flux f is related to e by

ð43Þ

The luminosity L is given by 4pr2f. The domain of computation is split into two zones through adjusting the radiation conductivity. In the radiation zone (RZ, upper 10% radius) radiation can carry all the energy and is convectively stable. In the convection zone (CZ, lower 80% radius), the radiation conductivity is suppressed so that most of the generated energy needs to be delivered by convection. The radiation conductivity is set as

jT ¼ cT ðf =dT=drÞinitial

ð44Þ

where f is the initial energy flux (obtain by Eq. (43)) and cT = 1 and 0.001 in the RZ and the CZ, respectively. In the CZ, an additional diffusion term is added to ensure stability of the energy equation, it is chosen to be proportional to the entropy gradient:

f S ¼ jS qT rS

ð45Þ

where S is the specific entropy (see [37] for a discussion of this term). Both the radiation and entropy diffusion terms can be ~ and p ~ used in the fluid equations (see Eq. (10)). jS describes expressed in terms of the gradients of the prognostic variables q energy transport by the sub-grid-scale turbulence and is taken to satisfy m/jS = Pr = 1/3 in the convection zone. Pr represents the effective Prandtl number of the convective turbulence. The polytropic indices for the CZ and RZ are 1.5 and 4, respectively. For efficient convection, the thermally relaxed state of the convection zone is close to the adiabatic state adopted here as initial condition (at least as far as gravity is concerned). The change in the radiative zone is also small [41]. 3.2. The numerical runs Here we present the results of four idealized cases of stellar core convection. The numerical layouts are given in Table 1. Cases A, B, and C compare different resolutions. Case D tests the ability of the code in handling rotation. Otherwise, all the models have the same physical configuration. The RZ (upper 10% radius) is subadiabatic so that it is convectively stable. It acts as a buffer between the CZ (lower 90% radius) and the rigid upper boundary. The depth of the stratification is given by a Z value of 20. It contains a total of 9 pressure scale heights (psh) with 3 psh in the subadiabatic zone and 6 psh in the convection zone. The model values of the mean flux at the top are about f = 2.7  104.

8706

T. Cai et al. / Journal of Computational Physics 230 (2011) 8698–8712

root mean square velocities

0.025

0.02

0.015

0.01

2 1/2 θ 2 1/2 φ 1/2



0.005

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

radius Fig. 3. The distribution of rms vertical and horizontal velocities of case B.

pffiffiffiffiffiffiffiffiffiffiffiffi The CFL number based on the sound speed is N CFL ¼ Cp=q 4 t= minð4r; r 4 h; rs 4 /Þ. The minimum grid spacing occurs at the innermost layer. For nmax P 2, min(4r, r4h, rs4/) = 2p2rmin/(ImaxJmax) where rmin is the radius of the first grid level from the center. Thus, according to Eqs. 27 and 28, the limitation on time stepping with an explicit method is inversely proportional to n2max of the innermost layer. For our semi-implicit calculations, the sound speed CFL number can reach 117. Even with such efficiency, the runs took very long time (over 2  106 steps) to reach rough thermal equilibrium in which the enR ergy leaving the top (total luminosity) is within 10% of the heat generation integral ( er2 dr d-). 3.2.1. Statistical results The results presented in this subsection are mainly from Cases B and D as their statistics have been accumulated over the longest periods of time. The quantities are averaged both temporally and horizontally (over a spherical surface). Fig. 3 shows the radial distributions of root-mean-square (rms) velocities of Case B. In the CZ, the vertical velocity is larger than the horizontal velocities except near the upper and lower boundaries. In contrast, the rms horizontal velocities are nearly equal everywhere, which is in agreement with the expectation of horizontal isotropy. The horizontal velocities maximize at the top of the CZ and near the center of the sphere. At the upper boundary of the CZ, vertical flows are sharply deflected. At the center, vertical and horizontal flows are mixed up and the rms velocity components become the same. From the graphs, one can see that the velocities are quite small (on the order of 0.02). As the unit of velocity is C1/2 times the 1.4 1.2

L

energy luminosities

e

1

Lk

0.8

Ld L

t

0.6 0.4 0.2 0 −0.2 −0.4 −0.6

0

0.2

0.4

0.6

0.8

1

radius Fig. 4. The distributions of enthalpy, kinetic, diffusive and total luminosities of the case B. All of the variables are normalized by the total luminosity at r = 1.

8707

T. Cai et al. / Journal of Computational Physics 230 (2011) 8698–8712

sound speed at the top of the domain, the Mach number of the flow is even smaller. Therefore, the ‘stratified approximation’ should hold very well. The Reynolds number is on the order of 400. Fig. 4 shows the profiles of different forms of energy flows, expressed here in terms of the luminosities defined as following:

Le ¼ Lk ¼ Ld ¼

Z Z Z

ðe þ pÞV r r 2 d- ðenthalpyÞ

ð46Þ

1 qV 2 V r r2 d- ðkinetic energyÞ 2

ð47Þ

ðjpr @ r p þ jqr @ r qÞr 2 d- ðdiffusionÞ

ð48Þ

One can see that convection is very effective; in most of the CZ it is the dominant energy carrier. The diffusion luminosity in the convection zone reflects sub-grid-scale turbulence mixing (Eq. (45)) and is less than 10% of the total luminosity in most of the CZ (except very near the convection–radiation boundary where radiation picks up). The kinetic energy luminosity is negative and its magnitude can reach up to 48% of the total luminosity; this amount is comparable to the high values ob-

0.045 1/2 θ 2 1/2

root mean square velocities

0.04



0.035

2 1/2



0.03 0.025 0.02 0.015 0.01 0.005 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

radius Fig. 5. The distribution of rms vertical and horizontal velocities of case D.

1.4 1.2

L

e

L

energy luminosities

1

k

L

d

0.8

L

t

0.6 0.4 0.2 0 −0.2 −0.4

0

0.2

0.4

0.6

0.8

1

radius Fig. 6. The distributions of enthalpy, kinetic, diffusive and total luminosities of the case D. All of the variables are normalized by the total luminosity at r = 1.

8708

T. Cai et al. / Journal of Computational Physics 230 (2011) 8698–8712

Fig. 7. Global maps of instantaneous radial velocity. Rows represent cases A, B, C and D. Columns represent top and middle.

tained by simulations in rectangular geometry [42]. The enthalpy luminosity reaches a maximum of 137% of the total luminosity. In the lower part of the RZ, a small region with negative enthalpy luminosity exists because of upward overshooting. In this region, the temperature fluctuation is anti-correlated with the vertical velocity, resulting in the negative enthalpy luminosity. At the same time, the diffusion luminosity increases above the total luminosity to compensate for this negative luminosity. As the thermal relaxation time is very long, the current computation can achieve a balance between energy generation and total luminosity to only within 10%. A small wiggle in the total luminosity appears near the convection–radiation boundary. This is partly caused by the sharp jump in radiation conductivity there.

T. Cai et al. / Journal of Computational Physics 230 (2011) 8698–8712

8709

Fig. 8. Global maps of instantaneous temperature fluctuation. Rows represent cases A, B, C and D. Columns represent top and middle.

Fig. 5 shows the rms velocities for Case D. The profiles are markedly different from the non-rotating case. The rms vertical velocity decreases very substantially from the upper CZ to the center. Compared with the non-rotating case, the magnitude of the vertical velocity is smaller. This is compatible with the expectation that rotation tends to hinder the process of convection. The impact is more significant in the lower part of the CZ where the turbulence is weaker. The horizontal velocities have similar behavior towards the center. They decrease monotonically rather than forming a peak there. The zonal velocity (V/) has higher magnitudes than the meridional velocity (Vh). This is caused by the presence of differential rotation (see next subsection). The differential rotation is a mean zonal velocity that superimposes on the turbulence. Fig. 6 shows the profiles of different forms of energy flows for the rotating case. In the lower region of the CZ (r < 0.7), diffusion (Eq. (11)) transports

8710

T. Cai et al. / Journal of Computational Physics 230 (2011) 8698–8712

more energy than the non-rotating case as the super-adiabaticity needs to be raised to enhance the driving of convection [43]. The magnitudes of the enthalpy and kinetic energy luminosities are maximized to only 123% and 33% of the total luminosity respectively. This again illustrates that rotation hinders convection. 3.2.2. Flow patterns Figs. 7 and 8 show the horizontal distributions (in upper CZ, r = 0.8, and in the middle, r = 0.50) of the vertical velocity and the temperature fluctuation, respectively, at single instances. Cases A–D are arranged in the order top to bottom. Common to all the cases, the horizontal length scales of the flow structures near the top are smaller than those lower down. As both levels are in the outermost layer and have the same horizontal resolution (with largest nmax value), this reflect a real physical phenomenon. It is compatible with earlier findings concerning the behavior of deep compressible convection [44,38]) and justifies the lowering of nmax in deeper layers. These figures also show that adequate resolution is needed for resolving the finer flow structures in the upper region of the convection zone. Compared to the higher resolution Cases B (T84) and C (T170), Case A (T42) produces smaller cells with broader dark lanes (down flows). The top left panel of Fig. 7 also indicates the presence of some high wavenumber oscillatory features; they are spurious oscillations caused by the exceedingly high grid Reynolds number. The lanes of Case C look sharper than those of Case B, but the general patterns are similar. The connected lanes are present only in a shallow region beneath the interface of the CZ and RZ. In the middle region of the CZ, the down flows become spotty. The up and down flows are highly asymmetric. At the center of the sphere, the turbulence becomes isotropic. In the upper and middle part of the CZ, there is a strong positive correlation between the vertical velocity and the temperature fluctuation. Fluid parcels with higher enthalpy contents are moving upward, resulting in the positive enthalpy flux. Besides affecting the statistical quantities, rotation also modifies the flow pattern. The mean zonal velocity is negative (towards west) near the equator but positive in higher latitudes. Effectively, the equatorial region rotates slower than the polar region. This shear exerts a westward pull on the cells near the equator; evidence of this can be seen in the bottom rows of Figs. 7 and 8. Fig. 9 shows an instantaneous distribution of zonal velocity on the equatorial plane of Case D. The outer part of the region is mainly occupied by retrograde flows, but the distribution is not symmetric with respect to the center (the rotation axis). The inner part contains prograde flows. The shear flow is further illustrated on the left panel of Fig. 10 which shows the time averaged zonal velocity distribution on the meridional plane. The contours approximately align with the rotation axis as can be explained by the Taylor–Proudman theorem. The value of Coriolis number (Co = Xd/V) is 5.1, where d is the length scale (radius of the sphere) and V is the rms total velocity. According to the results of [45], this is a relatively small Coriolis number, and as a result the differential rotation is ‘negative’. The pseudo streamlines of the time averaged meridional circulation is shown on the right panel of the same figure. Two cells are formed symmetrically in the northern and southern hemispheres. The direction is counterclockwise in the northern hemisphere, and clockwise in the southern hemisphere. The flows rise at the equator and sink at the poles. It is compatible with the conclusion of [46]. 4. Conclusion In this paper, a multi-layer semi-implicit spherical spectral method for stellar core simulation is developed. Unlike the anelastic approximation, the wave terms of the fully compressible hydrodynamic equations are unaltered. The wave and diffusion terms are treated implicitly, hence the time steps for integration are not restricted by the sonic and gravity waves. The

Fig. 9. The zonal velocity on the equatorial plane at one instant.

T. Cai et al. / Journal of Computational Physics 230 (2011) 8698–8712

8711

Fig. 10. The left panel shows the time averaged zonal velocity. The right panel shows the pseudo streamline of meridional circulation.

CFL restriction only involves the fluid velocity. The computational domain is extended all the way to the center so that the fluid motions would not be influenced by the presence of a solid core as in most other spectral methods. The multi-layer scheme can avoid the problem of shrinking grid size towards the center. Thus it is possible to use reasonably large time steps for integration all over the sphere. As the tests demonstrate, the time steps are indeed only restricted by the CFL condition associated with the subsonic fluid motions, in accordance with our anticipation. The stratified approximation is employed in reducing the degree of nonlinearity. Only quadratics appear in the nonlinear terms so that de-aliasing can be easily attained by the transforms between the spectral and physical spaces. Without de-aliasing, the alias errors can accumulate through the integration steps and produce secular numerical instabilities. When the boundary conditions are without external exchange of mass and momentum, the multi-layer semi-implicit spherical spectral scheme preserves the conservation of total mass and total angular momentum both analytically and numerically. The change of spectral truncation across layers does not influence the alias-free and conservative properties. The present code is applied to model the process of energy generation and transportation in a stellar core with convection. In the convection zone energy is mainly delivered by the enthalpy flux. The negative kinetic energy flux has a substantial impact on the flux balance and cannot be neglected. A shallow overshoot layer can be seen above the convection zone. The effect of rotation on turbulent convection is also tested. The results are compatible with the expectation that rotation tends to hinder convection. Rotation strongly reduced velocities near the center. Convection cells near the top of the convection zone can be seen bent by the shear of the differential rotation. The preliminary study of core convection here mainly serves as a test of the numerical code. The study of real stellar cores will be a future undertaking. Acknowledgments We thank Prof. Da Run Xiong and Dr. Friedrich Kupka for helpful discussions. This research is supported by the Hong Kong Research Grants Council (Project numbers 600306, 600309). L. Deng would like to thank National Science Foundation of China for support through Grants 10973015 & 11061120454. References [1] A. Maeder, Stellar evolution. V – Evolutionary models of population I stars with or without overshooting from convective cores, Astronomy and Astrophysics 47 (1976) 389–400. [2] A. Maeder, G. Meynet, Grids of evolutionary models from 0.85 to 120 solar masses – Observational tests and the mass limits, Astronomy and Astrophysics 210 (1989) 155–173. [3] A. Bressan, C. Chiosi, G. Bertelli, Mass loss and overshooting in massive stars, Astronomy and Astrophysics 102 (1981) 25–30. [4] M. Alongi, G. Bertelli, A. Bressan, C. Chiosi, Effects of envelope overshooting on stellar models, Astronomy and Astrophysics 244 (1) (1991) 95–106. [5] C. Chiosi, G. Bertelli, A. Bressan, New developments on understanding the HR diagram, Annual Review of Astronomy and Astrophysics 30 (1992) 235– 285. [6] R.B. Stothers, C. Chin, Metal opacities and convective core overshooting in Population I stars, The Astrophysical Journal 381 (1991) L67–L70. [7] F. Rogers, C. Iglesias, Radiative atomic Rossland mean opacity tables, The Astrophysical Journal Supplement Series 79 (1992) 507–568. [8] D.R. Xiong, Convective overshooting in stellar internal models, Astronomy and Astrophysics 150 (1985) 133–138. [9] V. Canuto, M. Dubovikov, Stellar turbulent convection. I. Theory, The Astrophysical Journal 493 (1998) 834–847. [10] F. Kupka, Convection in stars, in: Proceedings IAU Symposium 224, 2004. [11] S.A. Orszag, Transform method for the calculation of vector-coupled sums: Application to the spectral form of the vorticity equation, Journal of Atmospheric Sciences 27 (1970) 890–895.

8712

T. Cai et al. / Journal of Computational Physics 230 (2011) 8698–8712

[12] W. Bourke, An efficient, one-level, primitive-equation spectral model, Monthly Weather Review 100 (9) (1972) 683–691. [13] W. Bourke, A multi-level spectral model. i. formulation and hemispheric integrations, Monthly Weather Review 102 (1974) 687–701. [14] R. Daley, J. Girard, J. Henderson, I. Simmonds, Short-term forecasting with a multi-level spectral primitive equations model, Atmosphere 14 (1976) 98– 116. [15] C.T. Gordon, W.F. Stern, A description of the GFDL global spectral model, Monthly Weather Review 110 (1982) 625–644. [16] Y. Ogura, J. Charney, A Numerical Model of Thermal Convection in the Atmosphere, in: Proceedings of the International Symposium on Numerical Weather Prediction, Tokyo, 1962. [17] Y. Ogura, N. Phillips, Scale analysis of deep and shallow convection in the atmosphere, Journal of Atmospheric Sciences 19 (1962) 173–179. [18] F.B. Lipps, R.S. Hemler, A scale analysis of deep moist convection and some related numerical calculations, Journal of Atmospheric Sciences 39 (1982) 2192–2210. [19] D.O. Gough, The anelastic approximation for thermal convection, Journal of Atmospheric Sciences 26 (1969) 448–456. [20] J. Latour, E. Spiegel, J. Toomre, J. Zahn, Stellar convection theory. I – The anelastic modal equations, The Astrophysical Journal 207 (1976) 223–243. [21] P. Gilman, G. Glatzmaier, Compressible convection in a rotating spherical shell. I. Anelastic equations, The Astrophysical Journal Supplement Series 45 (1981) 335–388. [22] G.A. Glatzmaier, Numerical simulations of stellar convective dynamos. I – The model and method, Journal of Computational Physics 55 (1984) 461– 484. [23] G.A. Glatzmaier, Numerical simulations of stellar convective dynamos. II – Field propagation in the convection zone, Astrophysical Journal 291 (1985) 300–307. [24] G.A. Glatzmaier, Numerical simulations of stellar convective dynamos. III – At the base of the convection zone, Geophysical and Astrophysical Fluid Dynamics 31 (1985) 137–150. [25] T.C. Clune, J.R. Elliott, M.S. Miesch, J. Toomre, G.A. Glatzmaier, Computational aspects of a code to study rotating turbulent convection in spherical shells, Parallel Computing 25 (4) (1999) 361–380. [26] J. Elliott, M. Miesch, J. Toomre, Turbulent solar convection and its coupling with rotation: The effect of Prandtl number and thermal boundary conditions on the resulting differential rotation, The Astrophysical Journal 533 (2000) 546–556. [27] M. Miesch, J. Elliott, J. Toomre, T. Clune, G. Glatzmaier, P. Gilman, Three-dimensional spherical simulations of solar convection. I. Differential rotation and pattern evolution achieved with laminar and turbulent stars, The Astrophysical Journal 532 (2000) 593–615. [28] M. Browning, S. Brun, J. Toomre, Simulation of core convection in rotating a-type stars: differential rotation and overshooting, The Astrophysical Journal 601 (2004) 512–529. [29] A. Brun, M. Browning, J. Toomre, Simulation of core convection in rotating A-type stars: magnetic dynamo action, The Astrophysical Journal 629 (2005) 461–481. [30] M. Miesch, A. Brun, J. Toomre, Solar differential rotation influenced by latitudinal entropy variations in the tachocline, The Astrophysical Journal 641 (2006) 618–625. [31] K.L. Chan, H.G. Mayr, J.G. Mengel, I. Harris, A ‘‘Stratified’’ spectral model for stable and convective atmospheres, Journal of Computational Physics 113 (2) (1994) 165–176. [32] M. Evonuk, G.A. Glatzmaier, A 2D study of the effects of the size of a solid core on the equatorial flow in giant planets, Icarus 181 (2) (2006) 458–464. [33] A. Nordlund, A. Brandenburg, R.L. Jennings, M. Rieutord, J. Ruokolainen, R.F. Stein, I. Tuominen, Dynamo action in stratified convection with overshoot, The Astrophysical Journal 392 (1992) 647–652. [34] K.L. Chan, S. Sofia, Turbulent compressible convection in a deep atmosphere. IV – Results of three-dimensional computations, The Astrophysical Journal 336 (1989) 1022–1040. [35] K.L. Chan, The ‘Stratified’ approximation for computing geophysical flows, Southeast Asian Bulletin of Mathematics 20 (1996) 65–70. [36] L. Rivier, R. Loft, L.M. Polvani, An efficient spectral dynamical core for distributed memory computers, Monthly Weather Review 130 (2002) 1384– 1396. [37] K.L. Chan, S. Sofia, Turbulent compressible convection in a deep atmosphere. III – Tests on the validity and limitation of the numerical approach, The Astrophysical Journal 307 (1986) 222–241. [38] N. Brummell, N. Hurlburt, J. Toomre, Turbulent compressible convection with rotation. i. flow structure and evolution, The Astrophysical Journal 473 (1996) 494–513. [39] M. Iskandarani, D.B. Haidvogel, J.P. Boyd, A staggered spectral element model with application to the oceanic shallow water equations, International Journal for Numerical Methods in Fluids 20 (1995) 393–414. [40] E.N. Curchitser, M. Iskandarani, D.B. Haidvogel, A spectral element solution of the shallow-water equations on multiprocessor computers, Journal of Atmospheric and Oceanic Technology 15 (1998) 510–521. [41] J.P. Cox, R.T. Giuli, Principles of Stellar Structure, 1968. [42] K.L. Chan, S. Sofia, Turbulent compressible convection in a deep atmosphere. V. Higher order statistical moments for a deeper case, The Astrophysical Journal 466 (1996) 372–383. [43] K.L. Chan, Rotating convection in f-boxes: faster rotation, Astronomische Nachrichten 328 (2007) 1059–1061. [44] F. Cattaneo, N.H. Brummell, J. Toomre, A. Malagoli, N.E. Hurlburt, Turbulent compressible convection, The Astrophysical Journal 370 (1991) 282–294. [45] K.L. Chan, ‘Negative’ surface differential rotation in stars having low Coriolis numbers (slow rotation or high turbulence), in: A.G. Kosovichev, A.H. Andrei, J.-P. Roelot (Eds.), IAU Symposium, IAU Symposium, vol. 264, 2010, pp. 219–221. [46] B.R. Durney, Meridional motions and the angular momentum balance in the solar convection zone, The Astrophysical Journal 528 (2000) 486–492.